查看原文
其他

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

生信技能树 生信技能树 2022-06-06

一般来说,大家拿到了感兴趣的基因集后,通常是做超几何分布检验看看富集到了什么生物学功能数据库,比如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的。


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存