识别转录因子及互作基因的关联
<h1>背景</h1><p>Enrichment for Molecular Concept Map, to identify significant association within lists of TF genes。</p>
<p>通过转录因子靶基因的集合，用Fisher's exact test分析转录因子间的coupling（耦合、共调控），在形式上模仿paper里的这种图来展示结果：</p>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/30/milnmymnvw8a9od6n9l.png" alt="" /></p>
<p>出自PMID: 27288520文章</p>
<h1>应用场景</h1>
<p>分析任意集合之间的关系，例如转录因子之间（示例数据）、差异表达基因跟通路之间、染色质开放跟基因转录之间。</p>
<p>场景一：找转录因子之间的耦合关系（共同调控靶基因），例如示例数据(Interactions.RData)是201个转录因子，以及每个转录因子对应的多个有相互作用的基因。</p>
<p>场景二：差异表达基因在通路里的富集，就把差异表达基因跟通路里的基因都放到Interactions.RData里。例如：</p>
<p>$up_regulated_genes<br />
"ZBTB38" "SRA1" "IPO13""ALX4" "EP300""CREBBP" "IPO13"</p>
<p>$MAPK_pathway<br />
"ESR2" "ZNF688" "BANP" "ISYNA1" "TBC1D7"</p>
<p>场景三：分析染色质开放程度跟基因转录调控之间的关系，就把ATAC-seq peak附近的基因跟差异表达基因都放到Interactions.RData里。例如：</p>
<p>$up_regulated_genes<br />
"ZBTB38" "SRA1" "IPO13""ALX4" "EP300""CREBBP" "IPO13"</p>
<p>$ATAC-seq_peak<br />
"ESR2" "ZNF688" "BANP" "ISYNA1" "TBC1D7"</p>
<h1>环境设置</h1>
<p>使用国内镜像安装包</p>
<pre><code class="language-r">options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
install.packages("rentrez")
</code></pre>
<p>加载包</p>
<pre><code class="language-{r}">library(rentrez)
library(RColorBrewer)
library(corrplot)
library(ggplot2)
library(ggthemes)
source("my_function.R")
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
</code></pre>
<h1>输入文件</h1>
<p>如果你有自己的相互作用的基因list，保存在Interactions.RData里，就可以跳过这步，直接进入“Fisher's Exact Test”。</p>
<p>easy_input_TF.csv，此处用201个转录因子的gene symbol，可以替换成你感兴趣的基因的gene symbol，例如差异表达基因、差异表达lncRNA等。</p>
<p>easy_input_bg.csv，背景gene sybol，可以替换成与你感兴趣的基因相对应的背景基因，例如基因组上的全部基因、全部lncRNA等。</p>
<pre><code class="language-{r}">TFlist <- read.csv("easy_input_TF.csv", header = T)
head(TFlist)
dim(TFlist)
bglist <- read.table("easy_input_bg.csv", header = T)
head(bglist)
dim(bglist)
</code></pre>
<h2>从NCBI获取基因interaction</h2>
<p>用rentrez包获取NCBI gene数据库里的Interactions，即基因与它有相互作用的基因的list，保存到Interactions.RData。The general interactions in this section are provided, without review by Gene staff, by the external sources listed in https://ftp.ncbi.nlm.nih.gov/gene/GeneRIF/interaction_sources。详情看这里：https://www.ncbi.nlm.nih.gov/books/NBK3841/#EntrezGene.Interactions</p>
<pre><code class="language-r">#Download TF related interactions genes from NCBI Gene
#来源：https://github.com/ropensci/rentrez/wiki/Find-genes-known-to-interact-with-a-given-gene
res <- c()
for (i in 1:nrow(TFlist)) {
output <- print(as.character(TFlist$GeneSym),quote=FALSE)
gene_search <- entrez_search(db="gene",term=paste0("(",output,") AND (Homo sapiens)"))
# if you just want the interacting genes you can use this function, it's huge Xpath query
interactions_from_gene <- function(gene_id){
xmlrec <- entrez_fetch(db="gene", id=gene_id, rettype="xml", parsed=TRUE)
XML::xpathSApply(xmlrec,
"//Gene-commentary]//Other-source]//Other-source_anchor",
XML::xmlValue)
}
res1=interactions_from_gene(gene_search$ids)
res1=toupper(res1)
res1=list(res1)
names(res1)=output
res<- c(res,res1)
}
save(res, file = "Interactions.RData")
</code></pre>
<h1>Fisher's Exact Test</h1>
<p>Produce a geneset (lists) for hypergeometric test p-value</p>
<pre><code class="language-{r}">(load("Interactions.RData"))
head(res)
# Match and filter geneset by background genelist list size range (number of genes) from 150 - 2000
# 也可以跳过这行，不做筛选
res1 <- GeneSetsFilterByBackground(res,bglist$GeneSym,100,2000)
names(res1)
length(names(res1))
##########################################################
################# hypergeometric test ####################
##########################################################
# hyper and hyperTest function to caculate enrichment score and its pvalue by using hypergeometric test for each given genelist
# usage ## input 1. the list of GeneSet 2. length of background genes
# usage HyperGeoTest <- hyperTest(GeneSet, length(background.genes))
hyper <- function(X,Y,N,alpha=1){
K<-length(X)
M<-length(Y)
if(K==0 || M ==0){
return(1)
}
both<-length(intersect(X,Y));
XOnly<-length(setdiff(X,Y));
YOnly<-length(setdiff(Y,X));
Fisher.table<-N-both-XOnly-YOnly;
tab<-matrix(c(Fisher.table,YOnly, XOnly, both),2,2);
dimnames(tab)<-list(0:1,0:1);
f<-fisher.test(tab,alternative="greater");
pVal<-f$p.value;
estimate<-f$estimate;
return(c(pVal,estimate));
}
## s A list of gene vectors.
## N as a number to provide by using length (total background genes).
hyperTest<-function(s, N, alpha=1){
n<-length(s);
pVal<-matrix(0,n,n);
estimate<-matrix(0,n,n);
for(i in 1:n){
print(paste("Processing", i));
for(j in 1:n){
h<-hyper(s[],s[],N,alpha);
pVal<-h;
estimate<-h;
}
}
dimnames(pVal)<-list(names(s),names(s));
dimnames(estimate)<-list(names(s),names(s));
return(list(pVal=pVal,estimate=estimate));
}
# hypergeometric test p-value
HyperGeoTest <- hyperTest(res1, length(bglist$GeneSym))
rm(res1)
# Extract estimate value
est <- HyperGeoTest$estimate
# Remove self comparing results
est = NA
head(est)
# Extract pvalue
p <- HyperGeoTest$pVal
</code></pre>
<p>借用corrplot做clustring，为后面画图做准备</p>
<pre><code class="language-{r}"># set color
cols=rev(colorRampPalette(brewer.pal(9,"RdBu"))(199))
#cols=brewer.pal(n=9, name="RdBu")
p1 <- corrplot(as.matrix(log2(est+1)),
type = "upper", #只显示上三角
method="color",
order="hclust",
hclust.method="ward.D2",
col=cols,
tl.col ="black", #文字颜色
tl.cex = 0.5, #文字大小
tl.srt =45, #Text label color and rotation
is.corr = FALSE, #非相关系数矩阵。对于一般的矩阵，必须使用is.corr = FALSE
diag = F,#不展示相关系数
p.mat = as.matrix(p), #p值的矩阵
sig.level = c(1e-25,1e-35,1e-50), #显著水平
insig = c("label_sig"), #显著水平以***、**、*表示
pch.cex = 0.5, #标注的星号*的大小
font = 3) #斜体
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/30/milny5vdrlght25t7ca.png" alt="" /></p>
<h1>开始画图</h1>
<p>作者提供了三种画法展示以上结果，这里展示画法一，画法二三见压缩包里的2_3_ggplot2.R文件</p>
<h2>准备工作</h2>
<p>前面用corrplot做了聚类，提取顺序；只需要画一半，也就是三角形，需要把另一半变成NA。</p>
<pre><code class="language-{r}"># Remove self comparing p vaule
# p =NA
p=NA #or p = NA
p <- NA
# cluster ordering info by corrplot function
o <- rownames(p1$corr)
# p value reordered by ward.D2 hclust
p <- reorder(p, o)
# get -log10 fdr value
q <- -log10(p.adjust(p, method="fdr"))
# get -log10 p value
p <- -log10(p)
head(p)
# estimate value reordered by ward.D2 hclust
o <- rownames(p1$corr)
r <- reorder(est,o)
r <- NA
</code></pre>
<h2>画图</h2>
<pre><code class="language-{r}">pdf("Enrichment_baseplot.pdf", 11, 10)
par(bty="n",
mar=c(4,4,4,8)+.1, #四周留空
las=2,# the style of axis labels
tcl=-.33) #The length of tick marks as a fraction of the height of a line of text
m <- nrow(est)
n <- ncol(est)
# we should check range(est[!est==Inf]) for col and breaks setting
# we should check data distribution hist(est[!est==Inf]) for col and breaks setting
max_value<- range(r[!is.na(r)])
brks<- c(0,seq(1,20,l=10),seq(21,max_value,l=4))
cols2<- rev(colorRampPalette(brewer.pal(9,"RdBu"))(length(brks)-1))
# baseplot for heatmap
image(x=1:n, y=1:m, r, col=cols2, breaks=brks, xaxt="n", yaxt="n", xlab="",ylab="", xlim=c(0, n+4), ylim=c(0, n+1))
# add gene name
mtext(side=2, at=1:n, o, font=3, col="black") #left
mtext(side=3, at=1:n, o, font=3, col="black") #top
# add white border
abline(h=0:n+.5, col="white", lwd=.5)
abline(v=0:n+.5, col="white", lwd=.5)
# add title
#text(x=n/2, y=m+1, "Enrichment based on Fisher's Exact Test (greater)", pos=3)
q_range<- range(q[!is.na(q)])
# significant labels q>50
w = arrayInd(which(q > q_range/2), rep(m,2))
points(w, pch="*", col="black", cex=1)
# significant labels q>35
w = arrayInd(setdiff(which(q > q_range/3),which(q > q_range/2)),rep(m,2))
points(w, pch=3, col="black", cex=1) # "+""
# significant labels q>25
w = arrayInd(setdiff(which(q > q_range/4),which(q > q_range/3)), rep(m,2))
points(w, pch="-", col="black", cex=1)
# set up legend color bar
image(y = 1:16 +6, x=rep(n,2)+c(2.5,3)+1, z=matrix(c(1:16), nrow=1), col=cols2, add=TRUE)
# add legend color bar scale value
brks2 <- round(c(0,seq(1,20,l=10), seq(21,max_value,l=4)))
axis(side = 4, at = seq(1,15) + 6.5,tcl=-.15, label=brks2, las=1, lwd=.5)
points(x=rep(n,3)+3.5, y=1:3, pch=c("-","+","*"))
text(x=n+2, y=15, "enrichment estimate value",pos=1,srt=90)
mtext(side=4, at=c(1,2,3,4), c("-log10(FDR) > 25","-log10(FDR) > 35","-log10(FDR) > 50","ns"), line=0.2)
dev.off()
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/30/milo6480a5m13tewz1e.png" alt="" /></p>
<p>文中所用数据可以关注公众号“生信喵实验柴”<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/09/07/mf9nu72s90ofyseu9hl.png" alt="" /><br />
发送关键词“20250510”获取</p>
页:
[1]