1 简介：
转录因子 (transcription factors, TFs) 是与特定DNA序列结合 (TFBS/motif) ，调控DNA转录过程的一类蛋白质。转录因子在生命活动过程中，参与调节基因组DNA开放性、募集RNA聚合酶进行转录过程、募集辅助因子调节特定的转录阶段，调控生命体的免疫、发育等进程。本文我们一起学习单细胞转录组中转录因子的表达及调控活性。
SCENIC（single-cell regulatory network inference and clustering）是一个基于共表达和motif分析，计算单细胞转录组数据基因调控网络重建以及细胞状态鉴定的方法。2017年11月发表在Nature Methods 期刊上，是目前进行单细胞转录因子分析的主流软件，提供了R和python两个版本，（软件官网https://scenic.aertslab.org/），本文我们先学习R版本。SCENIC通过使用生物学驱动的features自动清除肿瘤样本特异性等批次效应。当前支持分析人、小鼠和果蝇的数据。
2 原理
基因调控网络（gene regulatory network）包括 转录因子（transcription factor，TF）、共调因子（cofactor）、靶基因（target gene）三部分组成。输入单细胞基因表达量矩阵后，SCENIC经过三个步骤完成转录因子分析：
- 构建共表达网络、
- 构建TF-targets网络、
- 计算Regulons活性，
三个步骤各由以下一个独立的软件完成：
2.1 GENIE3 - 构建共表达网络
GENIE3 (GEne NetworkInference with Ensemble of trees) ，基于树的基因网络推理，从基因表达数据推断基因调控网络。软件以单细胞基因表达量矩阵和TF表达作为输入文件，构建P个随机森林树（P=矩阵中基因数量），并计算每个TF与gene之间的重要性评分 (IM) ，IM代表了TF（input gene）在预测靶标时的权重，并测量它们各自的相关性以预测每个靶基因的表达，最终可以获得TF-genes共表达模块。最后删除IM低于阈值（0.001或0.005）的基因关系，过滤基因数低于50的模块。
2.2 RcisTarget-motif富集及靶基因预测
TF是通过直接与DNA结合而发挥作用的，结合的特异性序列，我们称之为motif，因此我们可以通过反向查看gene上是否存在TF结合的motif序列来验证TF与gene的靶向关系。RcisTarget可以从一个基因列表识别富集的TF-binding motifs和候选转录因子。
RcisTarget软件运行必备两个数据库：（1）gene-motif 排名数据库：为每个motif提供所有gene的排名（分数）；（2）motif-TF注释数据库：对每一个motif注释其所对应的TF。软件目前提供了人、小鼠、果蝇数据库，其他物种需要自己构建数据库。
RcisTarget首先基于gene-motif数据库，遍历每个motif对模块中所有基因进行累积，模块中的基因排名越靠前，累积曲线越高，曲线下面积 (AUC) 越大，表明motif在该模块中的富集程度越高，然后对每个模块选取显著富集的motif，并预测其靶基因，最终综合TF-genes模块和靶基因预测结果，构成一个包含了TF和靶基因的基因调控网络模块 (regulons)。进而对第一步预测出的TF-gene共表达网络进行验证。
2.3 AUCell-Regulons活性定量
第三步是对第二步得到的Regulons活性定量。通过比较所有细胞间的AUCell（area under the recovery curve）打分值，我们可以识别哪些细胞具有更显著高的regulon活性。
AUCell软件提供了一种新的方法，允许在scRNA-seq数据中识别具有活性基因调控网络的细胞。每一个regulons作为一个基因集输入到AUCell，这些基因集即Regulons中所有基因，针对每个细胞，将细胞中所有基因按照表达量从高到低进行排序，根据Regulons中的基因在序列中的位置，计算累计曲线面积 (AUC) ，即为Regulons在细胞中的活性。
但由于不同regulons包含的基因不同，它们之间的AUC值不具有可比较性，因此基于AUC值在所有细胞中的双峰分布特征，增加了Regulons“on/off”的概念，认为双峰之间的低谷为判断Regulons活性开放的阈值，如果AUC值小于阈值，则判定为该Regulons在该细胞中未开放，即未发挥调控作用。最终获得每个Regulons在每个细胞中的开放性热图。
3 实战代码
#### 引入需要的R包
library(foreach)
library(SCENIC)
library(Seurat)
library(GENIE3)
library(AUCell)
library(RcisTarget)
library(loomR)
library(stringr)
#### 1. 导入数据
setwd("/datafile2/scRNA/SCENIC/") # 设置工作路径
adata = readRDS("/datafile2/scRNA/SCENIC/pbmc_final.rds") # suerat分析的中间文件
exprMat = adata@assays$RNA@counts
exprMat = as.matrix(exprMat)
cellinfo = data.frame(adata@meta.data)
cellinfo = cellinfo[,c("sample","seurat_clusters","celltype","group")]
saveRDS(cellinfo,file="./int/cellinfo.rds")
# initializeScenic初始化，导入评分数据库
scenicOptions <- initializeScenic(org="mgi",dbDir = "./db/mmu/",nCores=20,dbs=c("mm9-500bp-upstream-7species.mc9nr.feather","mm9-tss-centered-10kb-7species.mc9nr.feather")) # 如果是人，则org="hgnc"
scenicOptions@inputDatasetInfo$cellinfo <- "./int/cellinfo.rds"
saveRDS(scenicOptions,file="./int/scenicOptions.rds")
3.1 共表达网络计算
# 1 转录调控网络推断
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
# 2 计算相关性矩阵
runCorrelation(exprMat_filtered, scenicOptions)
# 3 TF-Targets相关性回归分析。
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions)
# 输出文件说明：
# ...corrMat.rds 基因间相关性矩阵。
# ...GENIE3_weightMatrix_part_1.rds GENIE3 中间结果
# ...GENIE3_linkList.rds 最终结果
3.2 推断共表达模块
GENIE3计算了转录因子与每个基因的相关性,接下来需要过滤掉低相关性的组合，形成共表达基因集，有六种阈值可选
(1)w001、w005 分别以每个TF为核心保留weight>0.001和0.005的基因形成共表达模块。
(2)top50以每个TF为核心保留weight值top50的基因形成共表达模块。
(3)top5perTarget，top10perTarget，top50perTarget 为每个基因保留weight值top5，10，50的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块。
runSCENIC_1_coexNetwork2modules(scenicOptions)
3.3 motif 验证共表达模块
# 推断转录调控单元regulon
runSCENIC_2_createRegulons(scenicOptions)
# coexMethod过滤方法，默认6种方法都做，如果为提高运行速度，可以适当选择不同的阈值进行过滤
3.4 regulon 评分与可视化
#regulons计算AUC值并进行下游分析
exprMat_all = exprMat
exprMat_all <- log2(exprMat_all+1)
runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)
3.5 生成二进制的regulon活性矩阵。
runSCENIC_4_aucell_binarize(scenicOptions)
3.6 数据可视化
运行完以上的步骤后，我们发现SCENIC已经为我们绘制了非常多的图片，比如step3有所有regulon在细胞的AUCscore热图Step3_RegulonActivity_heatmap.pdf；step4生成多个热图Step4_BinaryRegulonActivity_Heatmap*等，此外我们如果对其他数据感兴趣，也可以进行提取：
# 1 转录因子富集结果
scenicOptions=readRDS(file="int/scenicOptions.Rds")
motifEnrichment_selfMotifs_wGenes <- loadInt(scenicOptions, "motifEnrichment_selfMotifs_wGenes")
gene_motif = as.data.frame(sort(table(motifEnrichment_selfMotifs_wGenes$highlightedTFs),decreasing = T)) # 基因的motif数量
# 可视化APT基因的motif
tableSubset <- motifEnrichment_selfMotifs_wGenes[highlightedTFs=="APT"]
viewMotifs(tableSubset)
# 2 loomFile提取AUC值,拿到数据后可以进行一些常规画图。
scenicOptions=readRDS(file="int/scenicOptions.Rds")
scenicLoomPath <- getOutName(scenicOptions, "loomFile")
loom <- open_loom(scenicLoomPath)
regulons_incidMat <- get_regulons(loom)
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonsAUC <- get_regulons_AUC(loom)
regulonsAucThresholds <- get_regulon_thresholds(loom)
embeddings <- get_embeddings(loom)
# 3 细胞类型对应的特异性活性单元
regulonAUC <-loadInt(scenicOptions, "aucell_regulonAUC")
cellInfo1 <- readRDS("./int/cellinfo.rds")
colnames(cellInfo1) = c("sample","cluster","celltype","Group")
rss <- calcRSS(AUC=getAUC(regulonAUC),
cellAnnotation=cellInfo1[colnames(regulonAUC), "celltype"])
rssPlot <-plotRSS(rss)#大小 rss评分，颜色 Z-score
ggsave("plotRSS_plot.pdf",rssPlot$plot)
