单细胞转录因子分析之SCENIC流程
去年我们在《生信技能树》公众号带领大家一起学习过： SCENIC转录因子分析结果的解读 ，提到了在做单细胞转录因子分析，首选的工具就是SCENIC流程，其工作流程 两次发表在nature系列杂志 足以说明它的优秀 :
-
SCENIC : single-cell regulatory network inference and clustering（2017年的nature methods）
- https://www.nature.com/articles/nmeth.4463
-
A scalable SCENIC workflow for single-cell gene regulatory network analysis（2020年的nature protocls）
- https://www.nature.com/articles/s41596-020-0336-2
SCENIC (Single-Cell rEgulatory Network Inference and Clustering) is a computational method to infer Gene Regulatory Networks and cell types from single-cell RNA-seq data. 官网教程非常清晰：
- Introduction and setup 安装
- running SCENIC 使用
提供了 (R / Python)两个版本的运行方式 ，SCENIC is implemented in R (this package and tutorial) and Python (pySCENIC).
安装SCENIC流程
其中R语言版本的 SCENIC流程 依赖于3个R包：
- GENIE3 to infer the co-expression network (faster alternative: GRNBoost2)
- RcisTarget for the analysis of transcription factor binding motifs
- AUCell to identify cells with active gene sets (gene-network) in scRNA-seq data
所以，实际上学习这个SCENIC流程，必须要 先学习这3个R包的 。
比如*RcisTarget* 包，基于DNA-motif 分析选择潜在的直接结合靶点，我也是写了两个教程：
再比如那个*AUCell* 包，我分享了： 使用AUCell包的AUCell_calcAUC函数计算每个细胞的每个基因集的活性程度
安装它们的话基本上复制粘贴下面的代码即可：
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::version()
# If your bioconductor version is previous to 4.0, see the section bellow
## Required
BiocManager::install(c("AUCell", "RcisTarget"),ask = F,update = F)
BiocManager::install(c("GENIE3"),ask = F,update = F) # 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"),ask = F,update = F)
# For various visualizations and perform t-SNEs:
BiocManager::install(c("DT", "NMF", "ComplexHeatmap", "R2HTML", "Rtsne"),ask = F,update = F)
# To support paralell execution (not available in Windows):
BiocManager::install(c("doMC", "doRNG"),ask = F,update = F)
# To export/visualize in http://scope.aertslab.org
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCENIC")
packageVersion("SCENIC")Python我不怎么使用，所以 Python (pySCENIC). 先略过。虽然这次略过了，但其实是躲不过去的，因为R里面的计算速度真心很慢，后期我们会补上这个 Python (pySCENIC). 教程哈。
另外，运行单细胞转录因子分析之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)
}每个文件都是1G，下载好了存放在共享文件夹， 可以多台电脑传输使用 ，没有必要每次都下载。而且可以看到，参考基因组版本并不是最新哦，人类的话这里使用hg19，小鼠使用mm9，都是很久以前的参考基因组啦。
SCENIC分析流程
下面才进入正餐！
输入数据的准备：表达矩阵加上表型信息
跟 CellPhoneDB运行需要的输入数据一样，也是 表达矩阵加上表型信息 即可，官网给出的示例数据是基于loom文件，实际上没有这个必要哈，最重要的是 表达矩阵加上表型信息 ：
rm(list = ls())
library(SCENIC)
## Load data
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
loom <- open_loom(loomPath)
exprMat <- get_dgem(loom)
cellInfo <- get_cell_annotation(loom)
close_loom(loom)
dim(exprMat)
exprMat[1:4,1:4]
head(cellInfo)
table(cellInfo$CellType)可以看到是，这个官方例子是从loom文件里面提取的 小鼠的表达矩阵 ，200个细胞分成5个亚群，基因数量就862个。这个数据量可以说是非常小啦，所以运行这个流程会很快哈。
> dim(exprMat)
[1] 862 200
> exprMat[1:4,1:4]
1772066100_D04 1772063062_G01 1772060224_F07 1772071035_G09
Arhgap18 0 0 0 0
Apln 0 0 0 0
Cnn3 0 0 0 0
Eya1 0 0 0 4
> head(cellInfo)
CellType nGene nUMI
1772066100_D04 interneurons 170 509
1772063062_G01 oligodendrocytes 152 443
1772060224_F07 microglia 218 737
1772071035_G09 pyramidal CA1 265 1068
1772067066_E12 oligodendrocytes 81 273
1772066100_B01 pyramidal CA1 108 191
> table(cellInfo$CellType)
astrocytes_ependymal interneurons microglia oligodendrocytes pyramidal CA1
29 43 14 55 59
> 但是呢， 我们自己的数据一般是在Seurat的对象里面 ，很容易制作这样的的两个信息： 表达矩阵加上表型信息 。后面我也有一个示例来讲解。
数据库文件的准备
不同物种不一样， 在 https://resources.aertslab.org/cistarget/ 查看自己的物种，按需下载，因为示例数据是小鼠表达矩阵，所以这里我们下载小鼠的，并且构建了一个 cisTarget_databases 文件夹 来存放下载好的数据库文件。
### Initialize settings
library(SCENIC)
# 保证 cisTarget_databases 文件夹下面有下载好2个1G的文件
scenicOptions <- initializeScenic(org="mgi", dbDir="cisTarget_databases", nCores=10)
# scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 需要注意的是，nCores=10 在部分电脑上面不适用哦，主要是 没办法开启 并行计算 。会报错如下：
Using 10 cores.
Warning in .AUCell_calcAUC(geneSets = geneSets, rankings = rankings, nCores = nCores, :
Using only the first 58.99 genes (aucMaxRank) to calculate the AUC.
Error in .AUCell_calcAUC(geneSets = geneSets, rankings = rankings, nCores = nCores, :
Valid 'mctype': 'snow' or 'doMC'解决方案是修改代码为 nCores=1 即可。
这个时候，该流程会 自动新建两个文件夹 ，一个是 int文件夹，一个是output文件夹。
运行SCENIC流程
首先是最耗费时间的步骤，Co-expression network
因为这个示例数据是小鼠表达矩阵200个细胞分成5个亚群，基因数量就862个。这个数据量可以说是非常小啦，所以运行这个流程会很快哈。但是 实际情况下 ，这个步骤至少耗费几个小时，甚至好几天，建议挂在后台慢慢运行。
### Co-expression network
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
exprMat_filtered[1:4,1:4]
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions)可以看到，我们的小鼠表达矩阵200个细胞分成5个亚群，基因数量就862个，这里面呢， 属于转录因子的居然仅仅是8个基因 ，所以函数会 提醒一下 ：Only 8 (0%) of the 1721 TFs in the database were found in the dataset. Do they use the same gene IDs?
其中这个 runGenie3 流程是限速步骤，可以转移到 服务器 上面运行，或者走对应的 Python (pySCENIC). 版本。
然后是SCENIC流程的4个步骤
SCENIC流程其实是 包装好的4个函数 ，有了前面的数据和数据库挖掘，完成起来非常快，如下：
runSCENIC_1_coexNetwork2modules(scenicOptions)
runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
runSCENIC_4_aucell_binarize(scenicOptions)这4个函数分别完成下面的 4个分析环节 ：
- GENIE3/GRNBoost：基于共表达情况鉴定每个TF的潜在靶点；
- RcisTarget：基于DNA-motif 分析选择潜在的直接结合靶点；
- AUCell：分析每个细胞的regulons活性；
- 细胞聚类：基于regulons的活性鉴定稳定的细胞状态并对结果进行探索
前面的3个步骤依次使用一个R包，所以 原则上需要独立去学习每个R包的用法 。大家千万不要妄想就看几次教程就学会了，一定要亲自实践。代码如下：
### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
tsneAUC(scenicOptions, aucType="AUC") # choose settings运行的log日志是：
> scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
17:23 Creating TF modules
Number of links between TFs and targets (weight>=0.001): 6139
nTFs 8
nTargets 770
nGeneSets 47
nLinks 17299
> scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
17:23 Step 2. Identifying regulons
tfModulesSummary:
top5perTarget 8
17:23 RcisTarget: Calculating AUC
Scoring database: [Source file: mm9-tss-centered-10kb-7species.mc9nr.feather]
17:23 RcisTarget: Adding motif annotation
Using BiocParallel...
Number of motifs in the initial enrichment: 867
Number of motifs annotated to the matching TF: 126
17:23 RcisTarget: Prunning targets
17:24 Number of motifs that support the regulons: 126
Preview of motif enrichment saved as: output/Step2_MotifEnrichment_preview.html
Warning messages:
1: In RcisTarget::importRankings(dbFilePath, columns = randomCol) :
The following columns are missing from the database:
2: In importRankings(rnkName, columns = allGenes) :
The following columns are missing from the database:
3: In importRankings(dbNames[motifDbName], columns = allGenes) :
The following columns are missing from the database:
4: In .addSignificantGenes(resultsTable = resultsTable, geneSets = geneSets, :
The rakings provided only include a subset of genes/regions included in the whole database.
> scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
17:24 Step 3. Analyzing the network activity in each individual cell
Number of regulons to evaluate on cells: 14
Biggest (non-extended) regulons:
Tef (405g)
Sox9 (150g)
Irf1 (104g)
Sox8 (97g)
Sox10 (88g)
Dlx5 (35g)
Stat6 (27g)
Quantiles for the number of genes detected by cell:
(Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).
min 1% 5% 10% 50% 100%
46.00 58.99 77.00 84.80 154.00 342.00
Warning in .AUCell_calcAUC(geneSets = geneSets, rankings = rankings, nCores = nCores, :
Using only the first 58.99 genes (aucMaxRank) to calculate the AUC.
17:24 Finished running AUCell.
17:24 Plotting heatmap...
17:24 Plotting t-SNEs...
> scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
Binary regulon activity: 7 TF regulons x 200 cells.
(13 regulons including 'extended' versions)
7 regulons are active in more than 1% (2) cells.
> scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
Binary regulon activity: 7 TF regulons x 200 cells.
(13 regulons including 'extended' versions)
7 regulons are active in more than 1% (2) cells.
> tsneAUC(scenicOptions, aucType="AUC") # choose settings
[1] "int/tSNE_AUC_50pcs_30perpl.Rds"到此为止，全部的SCENIC流程运行完毕，输出数据都是在output文件夹，剩下来的就是解读它了。不过，因为是测试数据， 小鼠的表达矩阵 ，200个细胞分成5个亚群，基因数量就862个。所以这个结果解读的意义也不大，我们还是直接来一个实战吧！
实战（以Seurat的pbmc3K数据集为例）
下面的 代码复制粘贴即可运行 ，超级简单，如果是你自己的数据，你只需 同样的模式 做出来 exprMat 表达矩阵，和cellInfo的临床表型，就可以走这个SCENIC流程的4个步骤啦。
rm(list = ls())
library(Seurat)
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
AvailableData()
# InstallData("pbmc3k") # (89.4 MB)
data("pbmc3k")
exprMat <- as.matrix(pbmc3k@assays$RNA@data)
dim(exprMat)
exprMat[1:4,1:4]
cellInfo <- pbmc3k@meta.data[,c(4,2,3)]
colnames(cellInfo)=c('CellType', 'nGene' ,'nUMI')
head(cellInfo)
table(cellInfo$CellType)
### Initialize settings
library(SCENIC)
# 保证cisTarget_databases 文件夹下面有下载好2个1G的文件
scenicOptions <- initializeScenic(org="hgnc",
dbDir="cisTarget_databases", nCores=1)
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
### Co-expression network
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
exprMat_filtered[1:4,1:4]
dim(exprMat_filtered)
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions)
### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions,
coexMethod=c("top5perTarget")) # Toy run settings
library(doParallel)
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log )
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
tsneAUC(scenicOptions, aucType="AUC") # choose settings因为我们这个是实战案例，表达矩阵很大，接近3000个细胞，是全部的人类基因，所以耗费了 一个晚上 才完成这个流程，运行的log日志如下：
> runGenie3(exprMat_filtered_log, scenicOptions)
Using 480 TFs as potential regulators...
Running GENIE3 part 1
Running GENIE3 part 10
Running GENIE3 part 2
Running GENIE3 part 3
Running GENIE3 part 4
Running GENIE3 part 5
Running GENIE3 part 6
Running GENIE3 part 7
Running GENIE3 part 8
Running GENIE3 part 9
Finished running GENIE3.
> ### Build and score the GRN
> exprMat_log <- log2(exprMat+1)
> scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
> scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
07:33 Creating TF modules
Number of links between TFs and targets (weight>=0.001): 1773984
nTFs 480
nTargets 5318
nGeneSets 3835
nLinks 2516422
> scenicOptions <- runSCENIC_2_createRegulons(scenicOptions,
+ coexMethod=c("top5perTarget")) # Toy run settings
07:34 Step 2. Identifying regulons
tfModulesSummary:
top5perTarget 59
07:34 RcisTarget: Calculating AUC
Scoring database: [Source file: hg19-tss-centered-10kb-7species.mc9nr.feather]
07:35 RcisTarget: Adding motif annotation
Using BiocParallel...
Number of motifs in the initial enrichment: 18696
Number of motifs annotated to the matching TF: 334
07:35 RcisTarget: Prunning targets
07:37 Number of motifs that support the regulons: 334
Preview of motif enrichment saved as: output/Step2_MotifEnrichment_preview.html
> library(doParallel)
Loading required package: iterators
> scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log )
07:37 Step 3. Analyzing the network activity in each individual cell
Number of regulons to evaluate on cells: 20
Biggest (non-extended) regulons:
SPI1 (538g)
CEBPB (47g)
CEBPD (41g)
SPIB (39g)
TBX21 (27g)
IRF8 (17g)
STAT1 (14g)
IRF7 (12g)
Quantiles for the number of genes detected by cell:
(Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).
min 1% 5% 10% 50% 100%
212.00 325.00 434.95 539.90 816.00 3400.00
07:38 Finished running AUCell.
07:38 Plotting heatmap...
07:38 Plotting t-SNEs...
> scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
Binary regulon activity: 11 TF regulons x 1678 cells.
(19 regulons including 'extended' versions)
11 regulons are active in more than 1% (16.78) cells.
> tsneAUC(scenicOptions, aucType="AUC") # choose settings
[1] "int/tSNE_AUC_50pcs_30perpl.Rds"可以看到，对于我们这个真实数据，就是PBMC3K的，也只有19个regulon被挑选出来了，涉及到11个TF基因
作者推荐的运算结果保存是：
export2loom(scenicOptions, exprMat)
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 实际上我们也用不上哈！
输出结果的解读
首先看看转录因子富集结果：
rm(list = ls())
library(Seurat)
library(SCENIC)
library(doParallel)
scenicOptions=readRDS(file="int/scenicOptions.Rds")
### Exploring output
# Check files in folder 'output'
# Browse the output .loom file @ http://scope.aertslab.org
# output/Step2_MotifEnrichment_preview.html in detail/subset:
motifEnrichment_selfMotifs_wGenes <- loadInt(scenicOptions, "motifEnrichment_selfMotifs_wGenes")
as.data.frame(sort(table(motifEnrichment_selfMotifs_wGenes$highlightedTFs),decreasing = T))每个基因的motif数量：
> as.data.frame(sort(table(motifEnrichment_selfMotifs_wGenes$highlightedTFs),decreasing = T))
Var1 Freq
1 SPI1 61
2 IRF7 59
3 TBX21 47
4 STAT1 33
5 SPIB 27
6 MAFB 26
7 IRF8 23
8 CEBPD 21
9 FOS 10
10 POU2AF1 5
11 CEBPB 3
12 TCF7 3
13 MAX 2
14 FOSB 1
15 LEF1 1
16 LYL1 1可视化IRF7基因的motif序列特征：
tableSubset <- motifEnrichment_selfMotifs_wGenes[highlightedTFs=="IRF7"]
viewMotifs(tableSubset) 这个时候的IRF7基因有 56 个motif，如下所示：
如果加上活性单元（regulon）的限定后：
regulonTargetsInfo <- loadInt(scenicOptions, "regulonTargetsInfo")
tableSubset <- regulonTargetsInfo[TF=="IRF7" & highConfAnnot==TRUE]
viewMotifs(tableSubset) 就只有12个啦，不过我们需要的并不是这些结果啦。
如果要理解（regulon），需要看我分享的： 使用AUCell包的AUCell_calcAUC函数计算每个细胞的每个基因集的活性程度
可以使用的结果：
其中output文件夹本来就已经自动绘制了大量的图表供使用，而图表对应的数据就存储在 loomFile 里面，可以使用下面的代码重新获取：
rm(list = ls())
library(Seurat)
library(SCENIC)
library(doParallel)
library(SCopeLoomR)
scenicOptions=readRDS(file="int/scenicOptions.Rds")
scenicLoomPath <- getOutName(scenicOptions, "loomFile")
loom <- open_loom(scenicLoomPath)