代码R版本，在这里R版本scenic

代码解析转自：单细胞SCENIC分析原理和流程介绍

1. 输入

输入：SCENIC需要输入的是单细胞RNA-seq表达矩阵—— 每列对应于样品（细胞），每行对应一个基因。基因ID应该是gene-symbol并存储为rownames （尤其是基因名字部分是为了与RcisTarget数据库兼容）；表达数据是Gene的reads count。根据作者的测试，提供原始的或Normalized UMI count，无论是否log转换，或使用TPM值，结果相差不大。(Overall, SCENIC is quite robust to this choice, we have applied SCENIC to datasets using raw (logged) UMI counts, normalized UMI counts, and TPM and they all provided reliable results (see Aibar et al. (2017)).)

另外，运行单细胞转录因子分析之SCENIC流程还需要下载配套数据库，不同物种不一样， 在 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)
  # 
}

2. SCENIC分析流程

SCENIC主要包含三个步骤：

• GENIE3(随机森林)/GRNBoost(Gradient Boosting)：基于共表达情况鉴定每个TF的潜在靶点，推断转录因子与候选靶基因之间的共表达模块。每个模块包含一个转录因子及其靶基因。（R语言版SCENIC这一步用的随机森林的算法，Python版用的是Gradient Boosting算法，两种算法不太一样，但算出来的结果是类似的）
• RcisTarget：由于GENIE3只是推断共表达，因此会有假阳性和间接targets。使用RcisTarget基于DNA-motif分析识别具有正确上游调控子且显著富集的motif（转录因子直接结合的motif），修剪掉缺乏motif支持的间接靶标。修剪后的每个TF和其潜在的直接作用的target genes被称为为一个regulon。（这一步是SCENIC和其他大多数共表达算法的重要区别）
• AUCell：分析每个细胞的regulons活性并进行打分，打分的基础是基因的表达值。对于regulon来说，比较细胞间的AUCell得分可以鉴定出哪种细胞具有显著更高的regulon活性。打分值可以进一步转化为二进制regulon活性矩阵（binarized activity matrix），这将最大化细胞类型的差异，确定regulon在哪些细胞中处于“开放”状态。

image

2.1 上游分析

image

在分析前首先要对initializeScenic()进行设置，后续的分析基于这个设定。

1. runCorrelation()
   共表达分析的结果中既有正向调控也有负向调控，GENIE3无法区分。因此需要相关性矩阵帮助筛选共表达模块中和TF正相关的基因。
2. runGenie3()
   参考数据库找出输入基因中哪些是TF，计算每个TF和各个基因之间的相关性权重。权重其实也就是TF对gene表达量的贡献。
3. runSCENIC_1_coexNetwork2modules()
   得到上面的TF和gene的权重矩阵以后，就可以生成以TF为核心的geneset。（比如TF1它可能对gene1、gene5、gene12...的预测都比较好，就可以得到以TF1为核心的geneset）。随后使用以下6种方法过滤掉低相关性的TF-genes共表达基因集，得到以TF为核心的共表达基因集。这一步运行得到的结果中包含了一列corr，是runCorrelation()得到的结果。1代表激活，-1代表抑制，0代表中性，SCENIC只会采用corr值为1的数据用于后续分析，以得到正调控共表达模块。

4. runSCENIC_2_createRegulons()
   通过RcisTarget数据库对得到的共表达模块进行修剪。通俗的讲就是在计算出的TF-gene对中，结合数据库查看该gene上游序列是否存在该TF结合的motif。从而剔除TF非直接调控基因的共表达模块，保留Motif分析共表达模块内与TF有直接调控关系的基因，得到Regulon。（RcisTarget数据库中的数据目前只支持人，鼠，果蝇三个物种。）
5. runSCENIC_3_scoreCells()
   每个Regulon就是一个转录因子及其直接调控靶基因的基因集，SCENIC接下来的工作就是对每个regulon在各个细胞中的活性评分，得到每个基因集在每个细胞的AUC score矩阵（AUC代表来与细胞内其他基因相比，特征基因中代表基因的比例及其相对表达值）。评分是基于recovery analysis，根据基因的表达值进行，分数越高代表基因集的激活程度越高。
6. runSCENIC_4_aucell_binarize()
   对于细胞类型清晰的数据集而言，构建regulons并对其活性打分足够后续分析。但是，在很多情况下将评分转换为二进制的“开/关”(on|off)，则既方便解释，又能最大化体现细胞类型的差异。将特定的regulon转换为“0/1”有利于探索和解释关键转录因子。将所有的regulons转换为“0/1”后创建二进制的活性矩阵，则可以用于细胞聚类，由于regulon是基于整体评分的，对消除技术偏倚如个别基因的dropout特别有用。AUCell会自动计算可能的阈值进行“0/1”转换，作者建议在转换之前手工检查并调整这些阈值。

⚠️在内存不足的情况下，可以使用抽样的方法，选取部分基因来得到Regulons，在runSCENIC_3_scoreCells()和runSCENIC_4_aucell_binarize()这两步时，再把表达矩阵替换成全部基因的表达矩阵来计算评分。

2.2 下游分析

三个方向：
1、降维聚类发现新亚群（cell type/state由转录调控网络的差异决定）
2、case-control之间的regulons差异分析
3、寻找cell type/state特异性的regulon/TF

代码参考

单细胞转录组高级分析二：转录调控网络分析-腾讯云开发者社区
SCENIC单细胞转录因子分析
单细胞转录调控网络分析工具（TF，motifs）-SCENIC
R-SCENIC 单细胞转录因子调控网络分析
SCENIC运行
SCENIC单细胞转录因子分析介绍及操作 - 简书 (jianshu.com)
R包SCENIC：分析流程 - 简书 (jianshu.com)
2019-08-28第七天：关于SCENIC的使用教程 - 简书 (jianshu.com)
使用SCENIC分析singlecell-RNAseq实战演练 - 简书 (jianshu.com)

3、正式代码

一、SCENIC分析准备

二、共表达网络计算

runGenie3产生文件如下

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种方法分别是：

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

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文件夹：

runSCENIC_2_createRegulons产生文件如下：

Step2_MotifEnrichment.tsv

image.png

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

image.png

Step2_regulonTargetsInfo.tsv

image.png

• 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活性tSNE图更有价值，举例说明如下：

image.png

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进行可视化。

image.png

四、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)

image.png

决定单核细胞命运的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值

image.png

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)

image.png

有兴趣的朋友还可以做个附加题：把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)

复制

获取帮助