scRNA-Seq | SCENIC 推断基因调控网络和细胞类型
转录因子 (transcription factors, TFs) 是直接作用于基因组，与特定DNA序列结合 (TFBS/motif) ，调控DNA转录过程的一类蛋白质。转录因子可以调节基因组DNA开放性、募集RNA聚合酶进行转录过程、募集辅助因子调节特定的转录阶段，调控诸多生命进程，诸如免疫反应、发育模式等。因此，分析转录因子表达及其调控活性对于解析复杂生命活动具有重要意义。
了解转录因子对于了解细胞的功能及生命活动有重要作用。传统的转录因子分析大多是对已知转录因子的表达和未知蛋白的转录因子预测分析，涉及活性分析的极少，而转录因子活性是其发挥作用的关键指标。
因此，今天想给大家介绍的就是一款可以分析转录因子活性的软件SCENIC[2]，这款软件是基于单细胞转录组数据开发的，可以解析单个细胞中转录因子活性。对于单细胞转录组而言，转录因子活性差异不仅为细胞异质性研究带来了新的思路，同时还可以从转录因子入手快速解析关键的分析调控机制。
SCENIC是2017年11月发表在 Nature Methods 期刊的一种单细胞转录因子分析方法，也是目前进行单细胞转录因子分析的主流软件，该软件在进行数据分析的同时也能得到可视化结果图。另外，SCENIC是一款开源软件，可以免费下载使用，目前软件有R和python两个版本，每个版本都配备了详细的使用说明（软件官网https://scenic.aertslab.org/）。但有一点需要特别注意，该软件是有物种限制的，目前只能分析人、小鼠和果蝇的数据，具体限制原因，我们在后面的分析原理中揭晓~
一、原理介绍
SCENIC就是一个常见的、基于单细胞转录组数据分析转录因子活性、基因调控网路的工具。SCENIC的分析主要分为三步，
- 第一步通过基因之间的共表达，找到可能的转录因子；
- 第二步是进行转录因子-motif的富集分析并找到对应的靶基因（调控组regulon）；
- 第三步是对调控组的活跃程度进行评分；
第一步由 GENIE3 包或 GRNBoost 实现。
GENIE3 (GEne NetworkInference with Ensemble of trees) ，基于树的基因网络推理，是一种从基因表达数据推断基因调控网络的方法[3]。软件以单细胞基因表达量矩阵为输入文件，以每个目标基因 (gene) 为输出，以转录因子 (TF) 为输入，构建P个随机森林树（P=矩阵中基因数量），并计算每个TF与gene之间的重要性评分 (IM) ，最终可以获得TF-genes共表达模块。最后删除IM低于阈值的基因关系，过滤基因数低于50的模块。
GENIE3网络构建过程
GENIE3用转录因子的表达量，通过训练随机森林（random forest）模型来预测各基因的表达量，从而得到转录因子在预测每个基因转录时的权重。这个权重反映了转录因子对于预测基因转录水平的相关性。相关性越高，则代表基因更有可能是该转录因子的靶基因。
随机森林是由多个决策树形成的分类器，它通过有放回的抽样训练出多个决策树，再以决策树结果中的众数为最终的结果。更详细的解释，可见：https://zhuanlan.zhihu.com/p/57965634
GENIE3的输入为基因表达矩阵，可以是UMI、TPM，或者FPKM/RPKM。而GENIE3的输出为基因、可能参与该基因的转录因子，以及它们的该转录因子的重要性（importance measure, IM），即其在预测基因转录水平时的权重。只有当权重高于0.001时，该转录因子才被认为是可能参与该基因调控的转录因子。
因为随机森林需要进行多次抽样，训练出多个决策树，当数据量很大时，这一步非常花时间，因此针对较大的数据，第一步可以用GRNBoost，它使用了梯度提升算法，在训练新的决策树时，会提高上一个决策树出错的样本比例，以针对模型预测欠缺的地方进行优化。
第二步由 RcisTarget 包实现。
它的主要作用在于通过一个基因列表，找到富集的转录因子及转录因子结合模序（motif），即可能的转录因子结合位点的模板序列。
它通过两步进行。
- 首先，它找到基因列表里基因的转录起始位点（transcript ion start site, TSS）。并找出转录起始位点周围高频出现的DNA motif。它会搜寻一个包含了跨物种基因组范围内各motif信息的数据库，筛选出和目标转录因子相关联的、标准化富集指数（normalised enrichment score, NES）高于3.0的motif。
- 接着，针对每一组motif和基因列表，RcisTarget会预测可能的目标基因。目标基因为基因列表中预测结果排序靠前的基因。所有motif共有的目标基因加上与之对应的转录因子，即为调控组regulon。
RcisTarget Regulons构建过程
第三步由AUCell实现，它能找到每个细胞中一组基因的活跃程度。
在这，SCENIC通过AUCell计算regulon的活跃程度。AUCell计算曲线下的面积（area under recovery curve, AUC），依据每个基因的表达水平，来计算出regular的活跃程度。各基因根据基因表达水平在x轴排序来绘制曲线。因此AUC反映了在每个细胞中给定的一组基因相对于其它基因的表达水平。
通过AUCell，我们能得到一个矩阵，包含每组regulon在每个细胞内相对于其它基因的表达水平，即它们的活跃程度。通过这一个矩阵，我们可以对细胞进行聚类，也可以看不同细胞类型中都有什么regulon是活跃的。
总之，SCENIC在R 中实现基于三个R包：
1.GENIE3：推断基因共表达网络
2.RcisTarget：用于分析转录因子结合motif
3.AUCell：用于鉴定scRNA-seq数据中具有活性基因集（基因网络）的细胞
a 使用GENIE3或GRNBoost推断转录因子与候选靶基因之间的共表达模块。RcisTarget可识别那些调节子的结合基序在目标基因中显着富集的模块。并创建仅具有直接target的调节单元。AUCell对每个细胞中每个调节单元的活性进行评分，从而产生活性矩阵。细胞状态基于调节子网络的共有的活性。
b SCENIC在小鼠大脑数据上的结果; 聚类标签的颜色对应主调节子。
c 显示了文献（A）证实的转录因子或具有MGI（B）的脑表型的转录因子，以及丰富的DNA基序。
d 活性矩阵上得到的t-SNE。每个细胞标记了最活跃的基因调节网络的颜色。
e 此数据集上不同聚类方法的准确性。
原文：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5937676/
以上是SCENIC的算法概述
二、pySCENIC python版本分析
适用背景
单细胞转录组调控网络分析是单细胞转录组分析内容的高级分析之一，本文将介绍SCENIC/pySCENIC的流程，具体原理和内容不展开，主要展示代码复现流程。R的SCENIC基于AUCell，RcisTarget和GENIE3三个包进行分析，所以要先安装这些依赖包，而pySCENIC则已经封装好，直接用pip安装即可。只用SCENIC或pySCENIC也可以单独完成分析，但R语言运行起来很慢，pySCENIC可以有效提升分析速度，还用SCENIC是因为可视化用R语言会简单一些。
1. 软件安装
创建一个环境，安装Python和R
conda create -n pyscenic
conda activate pyscenic
conda install r-base=4.2.1
conda install python=3.9.13
### Python安装pyscenic
pip install pyscenic -i https://mirrors.aliyun.com/pypi/simple/
pip install pandas -i https://mirrors.aliyun.com/pypi/simple/
pip install scanpy -i https://mirrors.aliyun.com/pypi/simple/
pip install opencv-python -i https://mirrors.aliyun.com/pypi/simple/
pip install loompy -i https://mirrors.aliyun.com/pypi/simple/
R安装SCENIC
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("AUCell", "RcisTarget","GENIE3"))
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github("aertslab/SCENIC")
这样就配置好基本环境了，接下来是安装R包，确保以下需要library的包就可以。
library(Seurat) #
library(SCopeLoomR) #
library(AUCell) #
library(SCENIC) #
library(dplyr) #
library(KernSmooth) #
library(RColorBrewer) #
library(plotly) #
library(BiocParallel) #
library(optparse) #
可以先用BiocManager::install() 安装，如果安装不了，谷歌搜一下conda install [package name]（百度搜不出来，不推荐）就能找到conda 安装这个包的官网，虽然conda很万能，但是慢啊，BiocManager::install会快很多。当然还有一些包是装不了，需要用devtools安装，例如SCopeLoomR，那哪些包要用devtools装呢？一般BiocManager::install和conda装不了的包就极大概率要用devtools装了，所以直接搜这个包的GitHub或官网，例如谷歌搜SCopeLoomR GitHub，GitHub里面说明需要这样安装：
devtools::install_github("aertslab/SCopeLoomR")
谈谈SCopeLoomR安装过程遇到的一些坑
本来想偷懒直接用别人的库路径，library的时候结果报错HDF5 library version mismatched error，R还abort突然中断了，就是说我的HDF5库版本跟我这个包安装时的版本不匹配，只能重装。用devtools重装又各种报错，缺失某个包，缺失某个动态库……反正各种缺失，没办法，只能缺什么装什么了，而且装devtools的时候就很容易报错，建议先用conda安装以下包和动态库再安装SCopeLoomR：
conda install -c anaconda git
conda install -c anaconda hdf5
conda install -c conda-forge libgit2
conda install -c conda-forge r-ragg
conda install -c conda-forge r-xml
再在R里运行：
devtools::install_github("aertslab/SCopeLoomR")
2.以人为例分析
2.1 从Seurat 对象生成loom文件
只需要使用build_loom函数，在get_count_from_seurat.R文件后面添加几行代码即可。
get_count_from_seurat.R文件内容如下：
library(optparse)
op_list <- list(
make_option(c("-i", "--inrds"), type = "character", default = NULL, action = "store", help = "The input of Seurat RDS",metavar="rds"),
make_option(c("-d", "--ident"), type = "character", default = NULL, action = "store", help = "The sample Ident of Seurat object",metavar="idents"),
make_option(c("-s", "--size"), type = "integer", default = NULL, action = "store", help = "The sample size of Seurat object",metavar="size"),
make_option(c("-l", "--label"), type = "character", default = "out", action = "store", help = "The label of output file",metavar="label"),
make_option(c("-a", "--assay"), type = "character", default = "Spatial", action = "store", help = "The assay of input file",metavar="assay")
)
parser <- OptionParser(option_list = op_list)
opt = parse_args(parser)
assay <- opt$assay
library(Seurat)
obj <- readRDS(opt$inrds)
if (!is.null(opt$ident)) {
Idents(obj) <- opt$ident
size=opt$size
if (!is.null(size)) {
obj <- subset(x = obj, downsample = opt$size)
}
saveRDS(obj,"subset.rds")
}
if (is.null(opt$label)) {
label1 <- 'out'
}else{
label1 <- opt$label
}
library(SCopeLoomR)
outloom <- paste0(label1,".loom")
build_loom(file.name = outloom,dgem = obj@assays[[assay]]@counts)
write.table(obj@meta.data,'metadata_subset.xls',sep='\t',quote=F)
使用方法：
-i：输入Seurat对象的RDS文件
-d：随机取样分组的列名，例如groups，如果不赋值则表示不随机取样，使用全部细胞
-s：随机取样的大小，例如20，因为这里用的是pbmc_small，所以设置小一点，实际情况可能需要设置大一点
-l：输出文件的标签，默认为out,
-a：使用的Seurat对象assay，默认为Spatial。
Rscript get_count_from_seurat.R \
-i test.rds \
-d groups \
-s 20 \
-l out \
-a Spatial
需要注意的是，新脚本添加了assay这个参数，默认为Spatial，因为最近在用跑空间组数据，如果要应用到单细胞需要改为RNA。
2.2 运行SCENIC的python版本，pyscenic
pyscenic的配置文件在官方提供的网站可以找到，需要三个文件：
-
1、TF列表文件：此处 或 此处
人的列表文件：https://github.com/aertslab/pySCENIC/blob/master/resources/hs_hgnc_tfs.txt -
2、数据库feather文件，需注意版本问题，一个是 pySCENIC < 0.12.0 和ctxcore < 0.2.0软件版本需要进到old文件夹，8月份最近的更新，另一个则是基因组版本，以及数据库来源的版本： 链接
人的feather文件：
https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-10species.mc9nr.feather
-
3、motif2tf数据库tbl文件： 链接
人的文件：
https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
以上文件选取仅供参考，实际情况可能需要选择不同文件。SCENIC的标准流程分为3步，
第一步利用GENIE3构建转录因子与基因的调控网络，
第二步利用RcisTarget验证转录因子与基因的调控网络的真实性；
第三步计算AUC曲线筛选调控网络。
pyscenic_from_loom.sh文件代码如下：
#default value
input_loom=out.loom
n_workers=20
#help function
function usage() {
echo -e "OPTIONS:\n-i|--input_loom:\t input loom file"
echo -e "-n|--n_workers:\t working core number"
echo -e "-h|--help:\t Usage information"
exit 1
}
#get value
while getopts :i:n:h opt
do
case "$opt" in
i) input_loom="$OPTARG" ;;
n) n_workers="$OPTARG" ;;
h) usage ;;
:) echo "This option -$OPTARG requires an argument."
exit 1 ;;
?) echo "-$OPTARG is not an option"
exit 2 ;;
esac
done
#需要更改路径
tfs=hs_hgnc_tfs.txt
feather=hg19-tss-centered-10kb-7species.mc9nr.feather
tbl=motifs-v9-nr.hgnc-m0.001-o0.0.tbl
pyscenic=/app/bin/pyscenic
# grn
$pyscenic grn \
--num_workers $n_workers \
--output grn.tsv \
--method grnboost2 \
$input_loom $tfs
# cistarget
$pyscenic ctx \
grn.tsv $feather \
--annotations_fname $tbl \
--expression_mtx_fname $input_loom \
--mode "dask_multiprocessing" \
--output ctx.csv \
--num_workers $n_workers \
--mask_dropouts
# AUCell
$pyscenic aucell \
$input_loom \
ctx.csv \
--output aucell.loom \
--num_workers $n_workers
使用方法：
-i，输入的loom文件，例如第二步生成的out.loom
-n，运行线程数，设置越大越快，根据实际情况设置
sh pyscenic_from_loom.sh -i out.loom -n 10
2.3 计算 RSS
此步骤的作用是通过计算RSS值找到细胞类型的特异TF。
calcRSS_by_scenic.R文件代码如下：
library(optparse)
op_list <- list(
make_option(c("-l", "--input_loom"), type = "character", default = NULL, action = "store", help = "The input of aucell loom file",metavar="rds"),
make_option(c("-m", "--input_meta"), type = "character", default = NULL, action = "store", help = "The metadata of Seurat object",metavar="idents"),
make_option(c("-c", "--celltype"), type = "character", default = NULL, action = "store", help = "The colname of metadata to calculate RSS",metavar="lab
el")
)
parser <- OptionParser(option_list = op_list)
opt = parse_args(parser)
library(Seurat)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
loom <- open_loom(opt$input_loom)
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')
regulonAucThresholds <- get_regulon_thresholds(loom)
close_loom(loom)
meta <- read.table(opt$input_meta,sep='\t',header=T,stringsAsFactor=F)
cellinfo <- meta[,c(opt$celltype,"nFeature_RNA","nCount_RNA")]
colnames(cellinfo)=c('celltype', 'nGene' ,'nUMI')
cellTypes <- as.data.frame(subset(cellinfo,select = 'celltype'))
selectedResolution <- "celltype"
sub_regulonAUC <- regulonAUC
rss <- calcRSS(AUC=getAUC(sub_regulonAUC),
cellAnnotation=cellTypes[colnames(sub_regulonAUC),
selectedResolution])
rss=na.omit(rss)
try({
rssPlot <- plotRSS(rss)
save(regulonAUC,rssPlot,regulons,file='regulon_RSS.Rdata')
})
saveRDS(rss,paste0(celltype,"_rss.rds"))
使用示例：
-i，第三步得到aucell.loom文件
-m，第一步得到的metadata_subset.xls文件
-c，计算RSS的分组列
Rscript calcRSS_by_scenic.R -l aucell.loom -m metadata_subset.xls -c groups
运行后生成regulon_RSS.Rdata和groups_rss.rds文件
Ref:
https://zhuanlan.zhihu.com/p/434003188
https://www.jianshu.com/p/d6feefe70513/
SCENIC | 从单细胞数据推断基因调控网络和细胞类型_生信宝典的博客-CSDN博客
SCENIC: 单细胞RNA-seq数据推断基因调控网络和细胞功能聚类 - 简书 (jianshu.com)
https://www.jianshu.com/p/c47332b880aa?v=1671031123892
[scRNA-seq]单细胞转录因子分析——SCENIC实操示例scenic分析小L的生信笔记的博客-CSDN博客
单细胞转录因子分析之SCENIC流程 - 知乎 (zhihu.com)
单细胞SCENIC分析原理和流程介绍 - 简书 (jianshu.com)
SCENIC单细胞转录因子调控网络分析_基因cell信息 (sohu.com)
单细胞SCENIC分析——寻找驱动基因 - 腾讯云开发者社区-腾讯云 (tencent.com)
Seurat_Satija - 简书 (jianshu.io)
https://blog.csdn.net/weixin_40594350/article/details/125493430
https://www.jianshu.com/p/bd0df23fbd60
