单细胞转录组高级分析二:转录调控网络分析
SCENIC简介
组织内细胞异质性的基础是细胞转录状态的差异,转录状态的特异性又是由转录因子主导的基因调控网络(GRNs)决定并维持稳定的。因此分析单细胞的GRNs有助于深入挖掘细胞异质性背后的生物学意义,并为疾病的诊断、治疗以及发育分化的研究提供有价值的线索。然而单细胞转录组数据具有背景噪音高、基因检出率低和表达矩阵稀疏性的特点,给传统统计学和生物信息学方法推断高质量的GRNs带来了挑战。
Single-cell regulatory network inference and clustering (SCENIC)是一种专为单细胞数据开发的GRNs算法,它的创新之处在于引入了转录因子motif序列验证统计学方法推断的基因共表达网络,从而识别高可靠性的由转录因子主导的GRNs。SCENIC相关的文章2017年首先发表于nature methods,2020年又将流程整理后发表于nature protocls。需要深入了解分析原理和流程的朋友可以参考这两篇文章:
SCENIC : single-cell regulatory network inference and clustering
A scalable SCENIC workflow for single-cell gene regulatory network analysis
SCENIC分析流程
官方介绍的主要分析有四步:
- GENIE3/GRNBoost:基于共表达情况鉴定每个TF的潜在靶点;
- RcisTarget:基于DNA-motif 分析选择潜在的直接结合靶点;
- AUCell:分析每个细胞的regulons活性;
- 细胞聚类:基于regulons的活性鉴定稳定的细胞状态并对结果进行探索
分析流程可以细化为以下步骤:
图片SCENIC安装
详细的安装教程参阅官网:
安装代码如下
在安装SCENIC之前,请按照以下代码安装一些依赖包:
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::version()
# If your bioconductor version is previous to 3.9, see the section bellow
## Required
BiocManager::install(c("AUCell", "RcisTarget"))
BiocManager::install(c("GENIE3")) # Optional. Can be replaced by GRNBoost
## Optional (but highly recommended):
# To score the network on cells (i.e. run AUCell):
BiocManager::install(c("zoo", "mixtools", "rbokeh"))
# For various visualizations and perform t-SNEs:
BiocManager::install(c("DT", "NMF", "pheatmap", "R2HTML", "Rtsne"))
# To support paralell execution (not available in Windows):
BiocManager::install(c("doMC", "doRNG"))
# To export/visualize in http://scope.aertslab.org
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)
检查一下核心依赖包的版本,并确保版本符合以下要求:
AUCell >=1.4.1 (minimum 1.2.4);
RcisTarget>=1.2.0 (minimum 1.0.2);
GENIE3>=1.4.0 (minimum 1.2.1).
packageVersion("AUCell")
packageVersion("RcisTarget")
packageVersion("GENIE3")
安装SCENIC包
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCENIC")
packageVersion("SCENIC")
还需要下载一些数据库
可以在此网址用浏览器下载:https://resources.aertslab.org/cistarget/
也可以使用以下代码在R环境中直接下载:
##1, For human:
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")
# mc9nr: Motif collection version 9: 24k motifs
##2, For mouse:
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
##3, For fly:
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/drosophila_melanogaster/dm6/flybase_r6.02/mc8nr/gene_based/dm6-5kb-upstream-full-tx-11species.mc8nr.feather")
# mc8nr: Motif collection version 8: 20k motifs
##4, download
dir.create("cisTarget_databases"); #创建一个文件夹保存数据库
setwd("cisTarget_databases")
#如果3个参考数据库都想下载,每次设置变量dbFiles后,都要运行以下代码
for(featherURL in dbFiles)
{
download.file(featherURL, destfile=basename(featherURL)) # saved in current dir
}
SCENIC分析准备
SCENIC分析需要输入单细胞表达矩阵和细胞meta信息,正式分析之前还要设置运行线程数,参考数据库目录、版本,细胞meta信息存放目录等配置信息。
library(Seurat)
library(tidyverse)
library(patchwork)
library(SCENIC)
rm(list=ls())
##==分析准备==##
dir.create("SCENIC")
dir.create("SCENIC/int")
scRNA <- readRDS("scRNA.rds")
setwd("~/project/10xDemo2/SCENIC")
##准备细胞meta信息
cellInfo <- data.frame(scRNA@meta.data)
colnames(cellInfo)[which(colnames(cellInfo)=="orig.ident")] <- "sample"
colnames(cellInfo)[which(colnames(cellInfo)=="seurat_clusters")] <- "cluster"
colnames(cellInfo)[which(colnames(cellInfo)=="celltype_Monaco")] <- "celltype"
cellInfo <- cellInfo[,c("sample","cluster","celltype")]
saveRDS(cellInfo, file="int/cellInfo.Rds")
##准备表达矩阵
#为了节省计算资源,随机抽取1000个细胞的数据子集
subcell <- sample(colnames(scRNA),1000)
scRNAsub <- scRNA[,subcell]
saveRDS(scRNAsub, "scRNAsub.rds")
exprMat <- as.matrix(scRNAsub@assays$RNA@counts)
##设置分析环境
mydbDIR <- "./cisTarget"
mydbs <- c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather",
"hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather")
names(mydbs) <- c("500bp", "10kb")
scenicOptions <- initializeScenic(org="hgnc",
nCores=8,
dbDir=mydbDIR,
dbs = mydbs,
datasetTitle = "HNSCC")
saveRDS(scenicOptions, "int/scenicOptions.rds")
scenicOptions <- initializeScenic()这一步要注意自己的情况设置参数:
-
nCores=8,代表开启8个线程计算,你的参数要根据自己电脑CPU的情况设置;
-
mydbDIR <- "./cisTarget",设置的是我存放数据库的目录,你要填写自己存放数据库的目录;
-
dbs = mydbs,我在前面的变量设置中,把hg38的数据库文件赋值给了mydbs,你用hg19数据库或小鼠的数据库需要相应调整。
共表达网络计算
相关性回归分析
SCENIC正式分析的第一步是计算转录因子与每个基因的相关性,此步骤消耗的计算资源非常大,作者建议两个策略处理大数据集:
-
采用GENIE3推断共表达模块时随机抽取少量细胞计算,计算regulons的活性时,所有细胞都代入运算;
-
使用python版本的SCENIC,作者强烈推荐。
为了减少初学者的学习负担,本文采用第一种策略。
##==转录调控网络推断==##
##基因过滤
#过滤标准是基因表达量之和>细胞数*3%,且在1%的细胞中表达
genesKept <- geneFiltering(exprMat, scenicOptions,
minCountsPerGene = 3 * 0.01 * ncol(exprMat),
minSamples = ncol(exprMat) * 0.01)
exprMat_filtered <- exprMat[genesKept, ]
##计算相关性矩阵
runCorrelation(exprMat_filtered, scenicOptions)
##TF-Targets相关性回归分析
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)
#这一步消耗的计算资源非常大,个人电脑需要几个小时的运行时间
runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)需要注意nParts参数,它的作用是把表达矩阵分成n份分开计算,目的是防止数据量大时内存不够。以上代码运行后,int目录下有不少中间结果产生,简要解释一下:
-
1.2_corrMat.Rds:基因之间的相关性矩阵
-
1.3_GENIE3_weightMatrix_part_1.Rds等:GENIE3的中间结果
-
1.4_GENIE3_linkList.Rds:GENIE3最终结果,是把“1.3_”开头的文件合并在一起。
查看1.4_GENIE3_linkList.Rds文件,head(1.4_GENIE3_linkList
)
> TF Target weight
2003093 XBP1 MZB1 0.1616139
2003094 XBP1 DERL3 0.1489923
2003095 XBP1 PRDX4 0.1393732
1003616 RPS4X RPL34 0.1285798
1254129 RPS4X RPS6 0.1265229
1 RPS4X RPL13 0.1256756</pre>
TF是转录因子名称,Target是潜在靶基因的名字,weight是TF与Target之间的相关性权重。
推断共表达模块
上一步计算了转录因子与每一个基因的相关性,接下来需要过滤低相关性的组合形成共表达基因集(模块)。作者尝试了多种策略(标准)过滤低相关性TF-Target,研究发现没有一种最佳策略,因此他们的建议是6种过滤标准都用。这6种方法分别是:
-
w001:以每个TF为核心保留weight>0.001的基因形成共表达模块;
-
w005:以每个TF为核心保留weight>0.005的基因形成共表达模块;
-
top50:以每个TF为核心保留weight值top50的基因形成共表达模块;
-
top5perTarget:每个基因保留weight值top5的TF得到精简的TF-Target关联表,然后把基因分配给TF构建共表达模块;
-
top10perTarget:每个基因保留weight值top10的TF得到精简的TF-Target关联表,然后把基因分配给TF构建共表达模块;
-
top50perTarget:每个基因保留weight值top50的TF得到精简的TF-Target关联表,然后把基因分配给TF构建共表达模块;
##推断共表达模块
runSCENIC_1_coexNetwork2modules(scenicOptions)
主要运行结果是int目录下的1.6_tfModules_asDF.Rds
> Target TF method corr
1 BCDIN3D ABCF2 w001 1
2 STRADB ABCF2 w001 1
3 HIST1H2BG ABCF2 w001 1
4 GID4 ABCF2 w001 1
5 LPCAT1 ABCF2 w001 -1
6 ZHX2 ABCF2 w001 1</pre>
method是上面提到的6种方法,corr是runCorrelation(exprMat_filtered, scenicOptions)命令得到的,1代表激活,-1代表抑制,0代表中性,SCENIC只会采用corr值为1的数据用于后续分析。
Motif验证共表达模块
经过上述分析每个转录因子都找到了强相关的靶基因,很多基因调控网络分析到此就结束了。SCENIC的创新之处是对此结果提出了质疑,并通过以下步骤修剪共表达模块形成有生物学意义的调控单元(regulons):
-
对每个共表达模块进行motif富集分析,保留显著富集的motif;此项分析依赖gene-motif评分(排行)数据库,其行为基因、列为motif、值为排名,就是我们下载的cisTarget数据库。
-
使用数据库对motif进行TF注释,注释结果分高、低可信度 。数据库直接注释和同源基因推断的TF是高可信结果,使用motif序列相似性注释的TF是低可信结果。
-
用保留的motif对共表达模块内的基因进行打分(同样依据cisTarget数据库),识别显著高分的基因(理解为motif离这些基因的TSS很近);
-
删除共表达模块内与motif评分不高的基因,剩下的基因集作者称为调控单元(regulon)。
##推断转录调控网络(regulon)
runSCENIC_2_createRegulons(scenicOptions)
#以上代码可增加参数coexMethod=c("w001", "w005", "top50", "top5perTarget", "top10perTarget", "top50perTarget"))
#默认6种方法的共表达网络都计算,可以少选几种方法以减少计算量
内存不够的电脑在这一步运行时会报错,重要的结果保存在output文件夹:
Step2_MotifEnrichment.tsv
图片这是各个共表达模块显著富集motif的注释信息,表头官方解释如下:
-
geneSet: Name of the gene set
-
motif: ID of the motif
-
NES: Normalized enrichment score of the motif in the gene-set
-
AUC: Area Under the Curve (used to calculate the NES)
-
TFinDB: Indicates whether the highlightedTFs are included within the high confidence annotation (two asterisks) or low confidence annotation (one asterisk).
-
TF_highConf: Transcription factors annotated to the motif according to ‘motifAnnot_highConfCat’.
-
TF_lowConf: Transcription factors annotated to the motif according to ‘motifAnnot_lowConfCat’.
-
enrichedGenes: Genes that are highly ranked for the given motif.
-
nErnGenes: Number of genes highly ranked
-
rankAtMax: Ranking at the maximum enrichment, used to determine the number of enriched genes.
Step2_MotifEnrichment_preview.html
图片显著富集的motif注释信息,表头含义同上个文件。
Step2_regulonTargetsInfo.tsv
图片最终的regulon文件,将第一个文件的数据整合在了一起,表头含义如下:
-
TF:转录因子名称
-
gene:TF靶基因名称
-
nMotif:靶基因在数据库的motif数量
-
bestMotif:最显著富集的motif名称
-
NES:标准富集分数,分值越高越显著
-
highConfAnnot:是不是高可信注释
-
Genie3Weight:TF与靶基因的相关性权重
请特别注意这个文件,后续分析找到了有价值的regulon,需要回到这个文件找对应的转录因子和靶基因 。SCENIC中regulon的名称有两种,一种是TF名称+extended+靶基因数目,另一种是TF名称+靶基因数目。第一种是转录因子与所有靶基因组成的基因调控网络,第二种是转录因子与高可信靶基因(即highConfAnnot=TRUE的基因)组成的基因调控网络。
Regulon活性评分与可视化
每个Regulon就是一个转录因子及其直接调控靶基因的基因集,SCENIC接下来的工作就是对每个regulon在各个细胞中的活性评分。评分的基础是基因的表达值,分数越高代表基因集的激活程度越高。我们推断regulons虽然只用了1000个随机抽取的细胞,但是regulon评分的时候可以把所有细胞导进来计算。
##==regulon活性评分与可视化==##
runSCENIC_3_scoreCells()是一个多种功能打包的函数,它包含的子功能有:
-
计算regulon在每个细胞中AUC值,最后得到一个以regulon为行细胞为列的矩阵。
结果目录:int/3.4_regulonAUC.Rds
-
计算每个regulon的AUC阈值,细胞中regulonAUC值>阈值,代表此regulon在细胞中处于激活状态,否则代表沉默状态。
结果目录:int/3.5_AUCellThresholds.Rds
-
使用regulonAUC矩阵对细胞进行降维聚类
-
用heatmap图展示regulonAUC矩阵,用t-SNE图分别展示每个regulon的活性分布情况。
结果目录:output/Step3_开头的系列文件
个人觉得分析结果中regulon活性tSNE图更有价值,举例说明如下:
图片左图展示的是某个regulon在所有细胞中AUC值分布情况,红色虚线是自动分配的阈值分界线;左中图蓝色点是regulonAUC值大于阈值的细胞,灰色点是regulonAUC值小于阈值的细胞;右中图是regulonAUC原始值的tSNE图;右图是regulon中转录因子的表达量tSNE图。
Regulon活性二进制转换与可视化
对于细胞类型清晰的数据集而言,构建regulons并对其活性打分足够后续分析。但是,在很多情况下将评分转换为二进制的“开/关”,则既方便解释,又能最大化体现细胞类型的差异。将特定的regulon转换为“0/1”有利于探索和解释关键转录因子。将所有的regulons转换为“0/1”后创建二进制的活性矩阵,则可以用于细胞聚类,对消除技术偏倚特别有用。AUCell会自动计算可能的阈值进行“0/1”转换,作者建议在转换之前手工检查并调整这些阈值。
查看调整阈值的代码(可选)
#使用shiny互动调整阈值
aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_all)
savedSelections <- shiny::runApp(aucellApp)
#保存调整后的阈值
newThresholds <- savedSelections$thresholds
scenicOptions@fileNames$int["aucell_thresholds",1] <- "int/newThresholds.Rds"
saveRDS(newThresholds, file=getIntName(scenicOptions, "aucell_thresholds"))
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
二进制转换及衍生分析
runSCENIC_4_aucell_binarize(scenicOptions, exprMat=exprMat_all)
runSCENIC_4_aucell_binarize()将regulonAUC矩阵转换为二进制矩阵后,会重新降维聚类,输出的结果与runSCENIC_3_scoreCells()类似。
SCENIC结果可视化
runSCENIC_3_scoreCells()和runSCENIC_4_aucell_binarize()运行之后会生成一些可视化的heatmap图与tSNE图,但是他们既不容易与seurat分析的结果联系起来,又不容易调整图形参数和分析内容。我们可以调用SCENIC的分析结果,使用seurat和pheatmap进行可视化。
图片热图图例显示不全,尺寸不可调;中图和右图是runSCENIC_3与runSCENIC_4得到的tSNE图,与seurat的tSNE图很难联系起来。
Seurat可视化SCENIC结果
把SCENIC结果中最重要的regulonAUC矩阵导入Seurat,这样得到的可视化结果更容易与我们之前的分析联系起来。
##导入原始regulonAUC矩阵
AUCmatrix <- readRDS("int/3.4_regulonAUC.Rds")
AUCmatrix <- AUCmatrix@assays@data@listData$AUC
AUCmatrix <- data.frame(t(AUCmatrix), check.names=F)
RegulonName_AUC <- colnames(AUCmatrix)
RegulonName_AUC <- gsub(' \\(','_',RegulonName_AUC)
RegulonName_AUC <- gsub('\\)','',RegulonName_AUC)
colnames(AUCmatrix) <- RegulonName_AUC
scRNAauc <- AddMetaData(scRNA, AUCmatrix)
scRNAauc@assays$integrated <- NULL
saveRDS(scRNAauc,'scRNAauc.rds')
##导入二进制regulonAUC矩阵
BINmatrix <- readRDS("int/4.1_binaryRegulonActivity.Rds")
BINmatrix <- data.frame(t(BINmatrix), check.names=F)
RegulonName_BIN <- colnames(BINmatrix)
RegulonName_BIN <- gsub(' \\(','_',RegulonName_BIN)
RegulonName_BIN <- gsub('\\)','',RegulonName_BIN)
colnames(BINmatrix) <- RegulonName_BIN
scRNAbin <- AddMetaData(scRNA, BINmatrix)
scRNAbin@assays$integrated <- NULL
saveRDS(scRNAbin, 'scRNAbin.rds')
##利用Seurat可视化AUC
dir.create('scenic_seurat')
#FeaturePlot
p1 = FeaturePlot(scRNAauc, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p2 = FeaturePlot(scRNAbin, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p3 = DimPlot(scRNA, reduction = 'tsne', group.by = "celltype_Monaco", label=T)
plotc = p1|p2|p3
ggsave('scenic_seurat/CEBPB_extended_2290g.png', plotc, width=14 ,height=4)
图片
决定单核细胞命运的regulon:转录因子CEBPB主导的调控网络
#RidgePlot&VlnPlot
图片
用山脊图和小提琴图展示CEBPB调控网络的AUC值
pheatmap可视化SCENIC结果****
#RidgePlot&VlnPlot
p1 = RidgePlot(scRNAauc, features = "CEBPB_extended_2290g", group.by="celltype_Monaco") +
theme(legend.position='none')
p2 = VlnPlot(scRNAauc, features = "CEBPB_extended_2290g", pt.size = 0, group.by="celltype_Monaco") +
theme(legend.position='none')
plotc = p1 + p2
ggsave('scenic_seurat/Ridge-Vln_CEBPB_extended_2290g.png', plotc, width=10, height=8)
图片
有兴趣的朋友还可以做个附加题:把CEBPB主导的基因调控网络做GO和KEGG富集分析。