基因集的转录因子富集分析

基因集的转录因子富集
为什么是AUC值而不是GSEA来挑选转录因子呢

library("RcisTarget")
# 官网的genelist
geneList1 <- read.table(file.path(system.file('examples', package='RcisTarget'), "hypoxiaGeneSet.txt"), 
                        stringsAsFactors=FALSE)[,1]
geneLists <- list(hypoxia=geneList1 ) 
geneLists


#motifAnnotations_hgnc里面包含了16万行，可以看到每个基因有多少个motif
data(motifAnnotations_hgnc)
motifAnnotations_hgnc


#前面我们下载好的 hg19-tss-centered-10kb-7species.mc9nr.feather 文件，需要存储在当前工作目录文件夹里面：
#这个文件下载了，在生信技能树的SCENIC文件夹里面
motifRankings <- importRankings("C:/Users/18700/Desktop/123/hg19-tss-centered-10kb-7species.mc9nr.feather")
# Motif enrichment analysis:geneSets必须是list
#最后结果可以看到100个motif被富集到，同时还计算了AUC和nes值。
motifEnrichmentTable_wGenes <- cisTarget(geneLists, motifRankings,
                                         motifAnnot=motifAnnotations_hgnc)



#因为通过RcisTarget包的 cisTarget()函数，一句代码完成的转录因子富集分析，等价于下面的3个步骤。
# 1. Calculate AUC
motifs_AUC <- calcAUC(geneLists, motifRankings)

# 2. Select significant motifs, add TF annotation & format as table
motifEnrichmentTable <- addMotifAnnotation(motifs_AUC, 
                                           motifAnnot=motifAnnotations_hgnc)
# 3. Identify significant genes for each motif
motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, 
                                                   geneSets=geneLists,
                                                   rankings=motifRankings, 
                                                   nCores=1,
                                                   method="aprox")

可以看到富集到了100个motif，可以看到每一个motif富集到的 highly ranked 基因

image.png
image.png

#首先批量计算AUC
motifs_AUC <- calcAUC(geneLists, motifRankings, nCores=1)
motifs_AUC#查看这个基因集有多少个motif参与了计算


#挑选统计学显著的motif
auc <- getAUC(motifs_AUC)[1,]
#画图判断auc是否符合正态分布
hist(auc, main="hypoxia", xlab="AUC histogram",
     breaks=100, col="#ff000050", border="darkred")
#一般来说，对正态分布，我们会挑选 mean+2sd范围外的认为是统计学显著，
#但这里选择的是 mean+3sd，更加严格 
nes3 <- (3*sd(auc)) + mean(auc)
abline(v=nes3, col="red")

data(motifAnnotations_hgnc) #读入所有的motif注释文件
motifAnnotations_hgnc
cg=auc[auc>nes3]#显示符合要求的有多少个motif
names(cg)
cgmotif=motifAnnotations_hgnc[match(names(cg),motifAnnotations_hgnc$motif),]
cgmotif=na.omit(cgmotif)#去除NA，剩余的就是挑选得到的motif

image.png 最后得到82个motif

#高级分析之可视化motif
motifEnrichmentTable_wGenes
#addLogo函数添加可视化图表
motifEnrichmentTable_wGenes_wLogo <- addLogo(motifEnrichmentTable_wGenes)
library(DT)
datatable(motifEnrichmentTable_wGenes_wLogo[,-c("enrichedGenes", "TF_lowConf"), with=FALSE], 
          escape = FALSE, # To show the logo
          filter="top", options=list(pageLength=5))



##高级分析之网络图
#数据库直接注释和同源基因推断的TF为高可信靶基因TF_highConf，
#使用motif序列相似性注释的TF为低可信度结果。
#将高可信靶基因拉出来，共51个
anotatedTfs <- lapply(split(motifEnrichmentTable_wGenes$TF_highConf,
                            motifEnrichmentTable$geneSet),
                      function(x) {
                        genes <- gsub(" \\(.*\\). ", "; ", x, fixed=FALSE)
                        genesSplit <- unique(unlist(strsplit(genes, "; ")))
                        return(genesSplit)
                      })

anotatedTfs$hypoxia
signifMotifNames <- motifEnrichmentTable$motif[1:4]#取前4个motif绘制网络图

#Identify which genes (of the gene-set) are highly ranked for each motif.
#incidMatrix里会有enrStatsh和incidMatrix两个，选择incidMatrix
#incidMatrix是靶基因和哪个TF相关的信息
incidenceMatrix <- getSignificantGenes(geneLists$hypoxia, 
                                       motifRankings,
                                       signifRankingNames=signifMotifNames,
                                       plotCurve=TRUE, maxRank=5000, 
                                       genesFormat="incidMatrix",
                                       method="aprox")$incidMatrix

#将绘图需要的edge和node文件进行创建
library(reshape2)
edges <- melt(incidenceMatrix)
edges <- edges[which(edges[,3]==1),1:2]
colnames(edges) <- c("from","to")
library(visNetwork)
motifs <- unique(as.character(edges[,1]))
genes <- unique(as.character(edges[,2]))
nodes <- data.frame(id=c(motifs, genes),   
                    label=c(motifs, genes),    
                    title=c(motifs, genes), # tooltip 
                    shape=c(rep("diamond", length(motifs)), rep("elypse", length(genes))),
                    color=c(rep("purple", length(motifs)), rep("skyblue", length(genes))))
visNetwork(nodes, edges) %>% visOptions(highlightNearest = TRUE, 
                                        nodesIdSelection = TRUE)

image.png image.png