基因集的转录因子富集分析 / 开普饭

一般来说，大家拿到了感兴趣的基因集后，通常是做超几何分布检验看看富集到了什么生物学功能数据库，比如KEGG或者GO数据库，或者走gsea/gsva这样的富集分析，也是注释生物学功能数据库。大家读我的表达芯片的公共数据库挖掘系列推文应该是够多了：

最近看到，单细胞领域有一个分析，可能是可以迁移到大家的常规表达量矩阵分析里面，就是最好拿到了基因集后：

使用RcisTarget包进行转录因子富集分析

下载及安装就不多说了，加载好RcisTarget包后，可以查看其详细的文档：

# Explore tutorials in the web browser: browseVignettes(package="RcisTarget") vignette("RcisTarget") # open

核心就是 cisTarget()函数的使用，需要3个输入数据：

首先是配套数据库文件

不同物种不一样，在 https://resources.aertslab.org/cistarget/ 查看自己的物种，按需下载，比如我这里就下载了人类和小鼠的数据：

# https://resources.aertslab.org/cistarget/ dbFiles <- c("https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather", "https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather") # dir.create("cisTarget_databases"); setwd("cisTarget_databases") # if needed dbFiles <- c("https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-7species.mc9nr.feather", "https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-7species.mc9nr.feather") # mc9nr: Motif collection version 9: 24k motifs dbFiles for(featherURL in dbFiles) { download.file(featherURL, destfile=basename(featherURL)) # saved in current dir # (1041.7 MB) # }

每个文件都是1G，下载好了存放在共享文件夹，可以多台电脑传输使用，没有必要每次都下载。

然后是基因集

一般来说呢，基因集可以是一次差异分析，挑选的符合统计学检验的，或者呢，直接挑选表达量或者离散度比较大的前500或者1000个基因，这个取决于大家的生物学研究目标。这里直接使用作者的提供的：

library(RcisTarget)

# Load gene sets to analyze. e.g.: geneList1 <- read.table(file.path(system.file('examples', package='RcisTarget'), "hypoxiaGeneSet.txt"), stringsAsFactors=FALSE)[,1] geneLists <- list(hypoxia=geneList1 ) geneLists

可以看到是 171个基因来自于hypoxiaGeneSet.txt文本文件，以symbol为名字：

> geneLists $geneListName [1] "ADM" "ADORA2B" "AHNAK2" "AK4" "AKAP12" "ALDOC" [7] "ANG" "ANGPTL4" "ANKZF1" "ARTN" "ASPH" "ATF3" [13] "ATG14" "ATXN1" "B3GNT4" "BBX" "BCOR" "BHLHE40" [19] "BNIP3" "BNIP3L" "C12orf24" "C5orf13" "C7orf68" "CA9" [25] "CADM1" "CAV1" "CCNG2" "CD59" "CITED2" "CSGALNACT1" [31] "CSRP2" "CXCR4" "CYB5A" "CYP1B1" "CYR61" "DAAM1"

接着需要一个ID转换列表

也是直接使用RcisTarget包提供的即可

# Select motif database to use (i.e. organism and distance around TSS) data(motifAnnotations_hgnc) motifAnnotations_hgnc

其中这个RcisTarget包内置的motifAnnotations_hgnc是16万行，可以看到每个基因有多个motif：

如果是其它物种，就需要自行去看文档，如何搞定这个ID转换列表啦。

最后就直接运行 `cisTarget()`函数

记住，前面我们下载好的 hg19-tss-centered-10kb-7species.mc9nr.feather 文件，需要存储在当前工作目录文件夹cisTarget_databases里面哦：

motifRankings <- importRankings("cisTarget_databases/hg19-tss-centered-10kb-7species.mc9nr.feather") # Motif enrichment analysis: motifEnrichmentTable_wGenes <- cisTarget(geneLists, motifRankings, motifAnnot=motifAnnotations_hgnc)

这个motifRankings是读取了数据库文件后的S4对象，最重要的信息是这个S4对象里面的rankings表格，节选如下:

解析结果

最后得到的结果是motifEnrichmentTable_wGenes，一个简单的数据框，总共是100个motif被富集到了，节选如下：

可以看到，NES概念大家都不陌生了，就是GSEA分析的结果而已，AUC的概念比较新颖。

其中EPAS1 (directAnnotation)这个基因命中了两个转录因子：

hocomoco__EPAS1_HUMAN.H11MO.0.B homer__GCACGTACCC_HIF2a

可以去前面的数据框查询看看，首先在RcisTarget包内置的motifAnnotations_hgnc，16万行数据信息里面可以找寻到这两个转录因子。

然后是在motifRankings查看：

motifs=unlist(as.data.frame(motifRankings@rankings[,1])) match('homer__GCACGTACCC_HIF2a',motifs) match('hocomoco__EPAS1_HUMAN.H11MO.0.B',motifs)

因为通过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 # (i.e. genes from the gene set in the top of the ranking) # Note: Method 'iCisTarget' instead of 'aprox' is more accurate, but slower motifEnrichmentTable_wGenes <- addSignificantGenes(motifEnrichmentTable, geneSets=geneLists, rankings=motifRankings, nCores=1, method="aprox")

后面我们再具体讲解，这3个步骤是如何挑选到这100个motif的。

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