基因集的转录因子富集
为什么是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
