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

生信技能树

发布于 2020-12-03 14:56:53

2.8K0

发布于 2020-12-03 14:56:53

文章被收录于专栏：生信技能树

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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-11-23，如有侵权请联系 cloudcommunity@tencent.com 删除

数据库

sql

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据库

sql

登录后参与评论

0 条评论

热度