SCENIC全称Single-Cell rEgulatory Network Inference and Clustering，即转录因子分析，用于构建GRN，gene regulatory network基因调控网络。2017年于Nat Methods首次发表R包版本(https://doi:10.1038/nmeth.4463.)，后又于2020年NATURE PROTOCOLS发表python版本(
https://doi.org/10.1038/s41596-020-0336-2
)，分析速度得到大幅度提升，而分析方法基本相同。
这篇笔记首先简单了解下分析步骤与原理，然后参考文章提供代码进行实操
1、TF转录因子
#
转录因子(Transcription factor)是指能够识别、结合在某基因上游特异核苷酸序列上的蛋白质;
这类蛋白通过介导RNA聚合酶与DNA模板的结合，从而调控下游靶基因的转录；也可以和其它转录因子形成转录因子复合体来影响基因的转录。
转录因子特异识别、结合的DNA序列称为转录因子结合位点(Transcription factor binding site,TFBS)，长度在5~20bp范围。由于同一个TF可以调控多个，而其具体结合到每个靶基因的TFBS不完全相同，但具有一定的保守性。
基于此，使用motif来概括TFBS的特征序列，具体是使用PWM(Postion Weight Matrix)表示。如下例图，即有的位点是某一种特定碱基，而有的位点是多种可能。
输入
：单细胞测序数据；
输出
：单细胞表达矩阵
这一步主要就是读入单细胞数据，过滤低质量的细胞、低表达的基因；然后将单细胞count表达矩阵(行名是细胞，列名是基因)存储为csv格式（或者loom）;
文章建议使用Scanpy预处理，因为更适合大的单细胞数据集。但个人还是更习惯于Seurat的流程。
2.2 网络推断Network inference
#
输入
：单细胞表达矩阵，转录因子列表；
输出
：TF-gene.csv文件。第1列是转录因子，第2列是共表达的基因，第3列是共表达强度(weight)
针对转录因子列表里的每一个转录因子，探索在单细胞表达矩阵中哪些基因与之存在共表达关系。
通过基于树(tree-based)的回归模型，根据转录因子的表达去预测基因的表达情况。具体算法可选择GRNBoost2或者GENIE3，前者的计算效率会更高。
2.3 得到TF-regulon
#
输入
：上一步得到的TF-gene.csv，motif-gene注释数据，motif-TF注释数据，单细胞表达数据
输出
：motif富集结果或者regulon数据
分为如下几个子步骤
(1) Module generation
#
首先根据第二步得到的TF-gene.csv文件，计算中TF -module(one Translate factor and its target genes)。scenic提供下面三种思路计算得到。
Based on a percentile score; the default is to select sets based on the 75th or 90th percentile of a factor’s targets.
Based on the top N targets for each factor, with a default of 50 targets.
Based on the top N regulators for a target; the default settings select sets of 5, 10, and 50 regulators.
按照如上默认值设置的话，有 2 + 1 + 3 = 6种思路的module结果。所以在最终得到的结果里，可能同一个TF有多个module。
根据调控基因与转录因子的相关性方向可分为 activating module(ρ ≥ +0.03)与 repressing module(ρ ≤ −0.03)。由于负相关模块(TF介导对靶基因的转录抑制)的数量一般较少，默认后续分析只保留正相关的activating module进行后续分析。
最后可限制module所包含的最少基因数，默认为20。
(2) Motif enrichment
#
在提供的motif-gene注释文件中，提供了每个motif的基因排名，即哪些基因的上游区域含有最类似该motif的区域。
基于此计算出每个module中所富集的motif有哪些。然后根据motif-TF注释数据，判断module 所富集到的motif是否为该module TF的motif。再然后，通过module的motif预测模块中真正的靶基因，去除其中的假阳性或者间接靶点。
最后通过合并同一个TF的所有保留基因，作为一个Regulon(one TF and its target gene)
2.4 Regulon活性评分
#
基于上一步得到的所有Regulon(可认为是Gene Set)，结合单细胞表达矩阵去评价每一个细胞对于每一个Regulon的表达程度评分。
不同细胞对于同一Regulon的表达活性是可以比较的，但同一细胞的不同Regulon的表达与否是不好比较的。
文章提供一种分析思路，将Regulon的活性得分二值化，即是否激活(On/Off)。这主要通过对那些具有 bimodal distribution双峰分布曲线的Regulon进行混合高斯分布拟合，确定合适的分割阈值，从而将细胞划分为两种状态。
原先参考代码是pySCENIC文章里提供的示范代码；如今版本以及数据库均已更新，现主要参考官方手册：
https://pyscenic.readthedocs.io/en/latest/index.html
3、pySCENIC代码实操
#
原文的预处理步骤是使用Scanpy进行单细胞预处理的，我还是比较熟悉Seurat流程，因此对文章给的代码有所删改。
3.1 准备conda环境
#
conda create -y -n pyscenic python=3.10
conda activate pyscenic
pip install pyscenic
3.2 下载参考数据#
（1）motif–gene注释数据(genome ranking database)
https://resources.aertslab.org/cistarget/databases/ 提供了不同物种以及不同基因组版本的注释数据，可根据需要选择；
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
##下载单细胞示例数据
# wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.tar.gz
# tar xvf pbmc_10k_v3_filtered_feature_bc_matrix.tar.gz
library(Seurat)
sce = Read10X(data.dir = "filtered_feature_bc_matrix")
sce = CreateSeuratObject(sce)
## 省略质控相关步骤
sce_count = GetAssayData(sce[["RNA"]], slot = "counts")
sce_count[1:4,1:4] #行名基因，列名细胞
write.csv(sce_count, file="sce_count.csv")
Note：
上面输出、保存的表达矩阵行名是基因，列名是细胞；
而SCENIC分析需要的表达矩阵格式为：行名是细胞，列名是基因；
因此可在之后pyscenic分析时添加参数 –transpose，进行转置。
3.4 scenic分析(命令行)#
分别对应上述的2.2~2.4简介
## expression matrix
f_ex_matrix_csv=/path/to/sce_count.csv
## referrence databases
f_db_names=$(ls /path/to/*feather | tr "\n" " ")
f_motif_path=/path/to/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl
f_tfs=/path/to/hs_hgnc_tfs.txt
(1) pyscenic grn#
pyscenic grn \
--output adj.tsv \
--method grnboost2 \
--transpose \
${f_ex_matrix_csv} \
${f_tfs} \
--num_workers 2
head adj.tsv
(2) pyscenic ctx#
pyscenic ctx adj.tsv \
${f_db_names} \
--annotations_fname ${f_motif_path} \
--expression_mtx_fname ${f_ex_matrix_csv} \
--transpose \
--output sce_regulon.gmt \
--mask_dropouts \
--mode "dask_multiprocessing" \
--num_workers 20
less -SN sce_regulon.gmt
#每一行是一个regulon：第一列是转录因子，第二列是说明性文字，第三列往后是target gene
(3) pyscenic aucell#
pyscenic aucell \
${f_ex_matrix_csv} \
sce_regulon.gmt \
--transpose \
--output sce_regulon_AUC.csv \
--num_workers 20
less -SN regulon_AUC.csv
将AUCell结果二值化：on/off。需要进入python环境
from pyscenic import binarization
import pandas as pd
auc_mtx = pd.read_csv("sce_regulon_AUC.csv", index_col="Regulon")
# time consuming
auc_binary = binarization.binarize(auc_mtx, num_workers=20)
# 结果返回一个tuple，包含两个元素：第一个是二值化后的表达矩阵；第二个是每个转录因子二值化的阈值
auc_binary[0]
auc_binary[1]
Note :
（1）输入单细胞表达矩阵也支持loom格式；输出结果中也支持不同类型的数据格式，具体可参看命令帮助文档 。
（2）官方手册中也提供了多种不同使用方法，例如在python中的分析流程；docker环境中的使用方法。
