单细胞转录组高级分析二：转录调控网络分析

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！（扫描文末二维码）

单细胞转录组基础分析专题

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分析流程

官方介绍的主要分析有四步：

1. GENIE3/GRNBoost：基于共表达情况鉴定每个TF的潜在靶点；

2. RcisTarget：基于DNA-motif 分析选择潜在的直接结合靶点；

3. AUCell：分析每个细胞的regulons活性；

4. 细胞聚类：基于regulons的活性鉴定稳定的细胞状态并对结果进行探索
   

分析流程可以细化为以下步骤：

SCENIC安装

详细的安装教程参阅官网：

https://rawcdn.githack.com/aertslab/SCENIC/701cc7cc4ac762b91479b3bd2eaf5ad5661dd8c2/inst/doc/SCENIC_Setup.html

安装代码如下

在安装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([e-mail])
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

TF是转录因子名称，Target是潜在靶基因的名字，weight是TF与Target之间的相关性权重。

推断共表达模块

上一步计算了转录因子与每一个基因的相关性，接下来需要过滤低相关性的组合形成共表达基因集（模块）。作者尝试了多种策略（标准）过滤低相关性TF-Target，研究发现没有一种最佳策略，因此他们的建议是6种过滤标准都用。这6种方法分别是：

1. w001：以每个TF为核心保留weight>0.001的基因形成共表达模块；

2. w005：以每个TF为核心保留weight>0.005的基因形成共表达模块；

3. top50：以每个TF为核心保留weight值top50的基因形成共表达模块；

4. top5perTarget：每个基因保留weight值top5的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；

5. top10perTarget：每个基因保留weight值top10的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；

6. 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

method是上面提到的6种方法，corr是runCorrelation(exprMat_filtered, scenicOptions)命令得到的，1代表激活，-1代表抑制，0代表中性，SCENIC只会采用corr值为1的数据用于后续分析。

Motif验证共表达模块

经过上述分析每个转录因子都找到了强相关的靶基因，很多基因调控网络分析到此就结束了。SCENIC的创新之处是对此结果提出了质疑，并通过以下步骤修剪共表达模块形成有生物学意义的调控单元（regulons）：

1. 对每个共表达模块进行motif富集分析，保留显著富集的motif；此项分析依赖gene-motif评分（排行）数据库，其行为基因、列为motif、值为排名，就是我们下载的cisTarget数据库。

2. 使用数据库对motif进行TF注释，注释结果分高、低可信度 。数据库直接注释和同源基因推断的TF是高可信结果，使用motif序列相似性注释的TF是低可信结果。

3. 用保留的motif对共表达模块内的基因进行打分（同样依据cisTarget数据库），识别显著高分的基因（理解为motif离这些基因的TSS很近）；

4. 删除共表达模块内与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活性评分与可视化==##
##regulons计算AUC值并进行下游分析
exprMat_all <- as.matrix(scRNA@assays$RNA@counts)
exprMat_all <- log2(exprMat_all+1)
runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)

runSCENIC_3_scoreCells()是一个多种功能打包的函数，它包含的子功能有：

1. 计算regulon在每个细胞中AUC值，最后得到一个以regulon为行细胞为列的矩阵。

   结果目录：int/3.4_regulonAUC.Rds

2. 计算每个regulon的AUC阈值，细胞中regulonAUC值>阈值，代表此regulon在细胞中处于激活状态，否则代表沉默状态。

   结果目录：int/3.5_AUCellThresholds.Rds

3. 使用regulonAUC矩阵对细胞进行降维聚类

4. 用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
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调控网络的AUC值

pheatmap可视化SCENIC结果

library(pheatmap)
cellInfo <- readRDS("int/cellInfo.Rds")
celltype = subset(cellInfo,select = 'celltype')
AUCmatrix <- t(AUCmatrix)
BINmatrix <- t(BINmatrix)
#挑选部分感兴趣的regulons
my.regulons <- c('ETS1_2372g','ETV7_981g','IRF7_239g','XBP1_854g','ATF4_37g',
                 'KLF13_78g','ATF6_129g','CREB3L2_619g','TAGLN2_13g',
                 'STAT1_extended_1808g','CEBPB_extended_2290g','IRF5_extended_422g',
                 'SPI1_1606g','HMGA1_14g','SPIB_1866g','IRF8_348g','BCL11A_136g',
                 'EBF1_40g','MAF_45g','BATF_131g','FOXP3_55g','TBX21_388g',
                 'EOMES_extended_101g','TCF7_extended_31g','LEF1_extended_49g')
myAUCmatrix <- AUCmatrix[rownames(AUCmatrix)%in%my.regulons,]
myBINmatrix <- BINmatrix[rownames(BINmatrix)%in%my.regulons,]
#使用regulon原始AUC值绘制热图
pheatmap(myAUCmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myAUCmatrix_heatmap.png',
         width = 6, height = 5)
#使用regulon二进制AUC值绘制热图
pheatmap(myBINmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myBINmatrix_heatmap.png',
         color = colorRampPalette(colors = c("white","black"))(100),
         width = 6, height = 5)

有兴趣的朋友还可以做个附加题：把CEBPB主导的基因调控网络做GO和KEGG富集分析。

后记

写这篇教程我花了一周的业余时间，这几天我都在忙些什么呢？考虑到服务器的配置和使用门槛较高，大家很可能在个人电脑上运行教程的内容，我所有的代码都会在自己的笔记本上验证一遍。此前文章的代码都没有问题，这次运行SCENIC遇到过不去的坎了。SCENIC有几步的计算量非常大，我的电脑运行一次要10几个小时，连续尝试了几个晚上都报错。后来搞清楚根本原因是内存不够，最后用服务器才顺利跑完了流程。如果你也遇到以下警告和错误提示信息，请放弃分析或减少分析的细胞数。

Warning message:
In mclapply(argsList, FUN, mc.preschedule = preschedule, mc.set.seed = set.seed,  :
 scheduled cores 3, 6, 8 did not deliver results, all values of the jobs will be affected

Error in signif(trhAssignment, 3)

获取帮助

本教程的目的在于把单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。