前面我們通過RcisTarget包的 cisTarget()函數，一句代碼就完成了我們的hypoxiaGeneSet.txt文本文件的171個基因的轉錄因子注釋。見：基因集的轉錄因子富集分析
通過學習，我們知道這個RcisTarget包內置的motifAnnotations_hgnc是16萬行，可以看到每個基因有多個motif。而且下載好的 hg19-tss-centered-10kb-7species.mc9nr.feather 文件，也是 24453個motifs的基因排序信息。但是我們留下來了一個懸念，如何從幾萬個注釋結果裡面挑選到最後100個富集成功的motif呢？
首先批量計算AUC值如果是單細胞轉錄組數據裡面，每個單細胞都是有一個geneLists，那麼就是成千上萬個這樣的calcAUC分析，非常耗費計算資源和時間，就需要考慮並行處理，我們這裡暫時不需要，所以直接 nCores=1 即可。
其中geneLists和motifRankings來自於前面的教程，見：基因集的轉錄因子富集分析
motifs_AUC <- calcAUC(geneLists, motifRankings, nCores=1)
motifs_AUC可以看到是 24453個motifs的AUC值都被計算了：
> motifs_AUC
AUC for 1 gene-sets and 24453 motifs.
AUC matrix preview:
jaspar__MA1023.1
geneListName 0.03819963
taipale_cyt_meth__IRX3_NACGYRNNNNNNYGCGTN_eDBD_meth
geneListName 0.05625049
taipale__DBP_DBD_NRTTACGTAAYN
geneListName 0.0640565
cisbp__M4240
geneListName 0.02816458
scertf__macisaac.ACE2
geneListName 0.03124153
>
挑選統計學顯著的motifauc <- getAUC(motifs_AUC)[1,]
hist(auc, main="hypoxia", xlab="AUC histogram",
breaks=100, col="#ff000050", border="darkred")
nes3 <- (3*sd(auc)) + mean(auc)
abline(v=nes3, col="red")可以看到 24453個motifs的AUC值看起來滿足正態分布，一般來說，對正態分布，我們會挑選 mean+2sd範圍外的認為是統計學顯著，但是作者卡的比較嚴格，是 mean+3sd ，示意圖如下：
看看Area Under the Curve (AUC)如何計算這個時候就需要一個取捨了，我們是否需要知道每個細節，比如GSEA分析，我也多次講解：
但實際上，絕大部分讀者並沒有去細看這個統計學原理，也不需要知道gsea分析的nes值如何計算，或者說這個Area Under the Curve (AUC)如何計算。
我這裡也不想耗費時間去深究，去講解了。不理解原理並不影響大家使用，知道這個概念，知道如何根據AUC值去判斷結果就好。其實這個包的核心在於motifRankings變量，資料庫文件來自於前面的教程，見：基因集的轉錄因子富集分析，也是很容易製作的，選取人類的不到2000個TF的全部chip-seq數據的peaks文件的bed，把人類的2萬個基因的啟動子區域的該TF的信號強度排序即可。
然後看看motif的詳情這個RcisTarget包內置的motifAnnotations_hgnc是16萬行，可以看到每個基因有多個motif，我們挑選出來了105個moif，去這個表格裡面篩選一下，就只剩下82個了。
data(motifAnnotations_hgnc)
motifAnnotations_hgnc
cg=auc[auc>nes3]
names(cg)
cgmotif=motifAnnotations_hgnc[match(names(cg),motifAnnotations_hgnc$motif),]
cgmotif=na.omit(cgmotif)
高級分析之可視化motif前面的教程，一句代碼就完成了motif的富集 ，見：基因集的轉錄因子富集分析
這個時候，我們可以根據 addLogo 函數，對我們富集好的轉錄因子的motif加上一些可視化圖表，代碼如下：
motifEnrichmentTable_wGenes
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))這些motif都是資料庫已知的，其可視化如下：
高級分析之網絡圖這裡面的R代碼技巧還是蠻值得細細品讀的：
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:3]
incidenceMatrix <- getSignificantGenes(geneLists$hypoxia,
motifRankings,
signifRankingNames=signifMotifNames,
plotCurve=TRUE, maxRank=5000,
genesFormat="incidMatrix",
method="aprox")$incidMatrix
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)一個簡陋的網絡圖就出來了：
PPI調控網絡圖確實有點老套了我有預感，這個轉錄因子調控網絡圖應該是在未來5年內會逐步替代PPI調控網絡圖，直到轉錄因子調控網絡圖也變得俗氣為止。