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富集分析。
网友评论