去年我们在《生信技能树》公众号带领大家一起学习过：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 分析选择潜在的直接结合靶点，我也是写了两个教程：

• 基因集的转录因子富集分析
• 为什么是AUC值而不是GSEA来挑选转录因子呢

再比如那个*AUCell* 包，我分享了：使用AUCell包的AUCell_calcAUC函数计算每个细胞的每个基因集的活性程度

安装它们的话基本上复制粘贴下面的代码即可：

Python我不怎么使用，所以 Python (pySCENIC).  先略过。虽然这次略过了，但其实是躲不过去的，因为R里面的计算速度真心很慢，后期我们会补上这个 Python (pySCENIC). 教程哈。

另外，运行单细胞转录因子分析之SCENIC流程还需要下载配套数据库，不同物种不一样， 在 https://resources.aertslab.org/cistarget/ 查看自己的物种，按需下载：

每个文件都是1G，下载好了存放在共享文件夹，可以多台电脑传输使用，没有必要每次都下载。而且可以看到，参考基因组版本并不是最新哦，人类的话这里使用hg19，小鼠使用mm9，都是很久以前的参考基因组啦。

SCENIC分析流程

下面才进入正餐！

输入数据的准备：表达矩阵加上表型信息

跟 CellPhoneDB运行需要的输入数据一样，也是表达矩阵加上表型信息即可，官网给出的示例数据是基于loom文件，实际上没有这个必要哈，最重要的是 表达矩阵加上表型信息 ：

可以看到是，这个官方例子是从loom文件里面提取的小鼠的表达矩阵，200个细胞分成5个亚群，基因数量就862个。这个数据量可以说是非常小啦，所以运行这个流程会很快哈。

但是呢，我们自己的数据一般是在Seurat的对象里面，很容易制作这样的的两个信息：表达矩阵加上表型信息。后面我也有一个示例来讲解。

数据库文件的准备

不同物种不一样， 在 https://resources.aertslab.org/cistarget/ 查看自己的物种，按需下载，因为示例数据是小鼠表达矩阵，所以这里我们下载小鼠的，并且构建了一个 cisTarget_databases 文件夹 来存放下载好的数据库文件。

需要注意的是，nCores=10 在部分电脑上面不适用哦，主要是没办法开启并行计算。会报错如下：

解决方案是修改代码为  nCores=1 即可。

这个时候，该流程会自动新建两个文件夹，一个是 int文件夹，一个是output文件夹。

运行SCENIC流程

首先是最耗费时间的步骤，Co-expression network

因为这个示例数据是小鼠表达矩阵200个细胞分成5个亚群，基因数量就862个。这个数据量可以说是非常小啦，所以运行这个流程会很快哈。但是实际情况下，这个步骤至少耗费几个小时，甚至好几天，建议挂在后台慢慢运行。

可以看到，我们的小鼠表达矩阵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个函数，有了前面的数据和数据库挖掘，完成起来非常快，如下：

这4个函数分别完成下面的4个分析环节：

1. GENIE3/GRNBoost：基于共表达情况鉴定每个TF的潜在靶点；
2. RcisTarget：基于DNA-motif 分析选择潜在的直接结合靶点；
1. 基因集的转录因子富集分析
2. 为什么是AUC值而不是GSEA来挑选转录因子呢
3. AUCell：分析每个细胞的regulons活性；
4. 细胞聚类：基于regulons的活性鉴定稳定的细胞状态并对结果进行探索

前面的3个步骤依次使用一个R包，所以原则上需要独立去学习每个R包的用法。大家千万不要妄想就看几次教程就学会了，一定要亲自实践。代码如下：

运行的log日志是：

到此为止，全部的SCENIC流程运行完毕，输出数据都是在output文件夹，剩下来的就是解读它了。不过，因为是测试数据，小鼠的表达矩阵，200个细胞分成5个亚群，基因数量就862个。所以这个结果解读的意义也不大，我们还是直接来一个实战吧！

实战（以Seurat的pbmc3K数据集为例）

下面的代码复制粘贴即可运行，超级简单，如果是你自己的数据，你只需同样的模式做出来 exprMat 表达矩阵，和cellInfo的临床表型，就可以走这个SCENIC流程的4个步骤啦。

因为我们这个是实战案例，表达矩阵很大，接近3000个细胞，是全部的人类基因，所以耗费了一个晚上才完成这个流程，运行的log日志如下：

可以看到，对于我们这个真实数据，就是PBMC3K的，也只有19个regulon被挑选出来了，涉及到11个TF基因

作者推荐的运算结果保存是：

实际上我们也用不上哈！

输出结果的解读

首先看看转录因子富集结果：

每个基因的motif数量：

可视化IRF7基因的motif序列特征：

这个时候的IRF7基因有 56 个motif，如下所示：

如果加上活性单元（regulon）的限定后：

就只有12个啦，不过我们需要的并不是这些结果啦。

如果要理解（regulon），需要看我分享的：使用AUCell包的AUCell_calcAUC函数计算每个细胞的每个基因集的活性程度

可以使用的结果：

其中output文件夹本来就已经自动绘制了大量的图表供使用，而图表对应的数据就存储在 loomFile 里面，可以使用下面的代码重新获取：

可以可视化其中一些TF的AUC值，也可以根据这些TF的AUC值对细胞亚群进行重新降维聚类分群。可以很容易看到血液里面的不同细胞亚群的特异性的转录调控因子：

如果使用这些转录调控因子进行 降维聚类分群 ，可以得到：

更多文章实例图表可以看：SCENIC转录因子分析结果的解读 ，这里面我埋下了两个伏笔，都是关于R里面的这个单细胞转录因子分析之SCENIC流程运行超级慢的问题，仅仅是接近3000个细胞就耗费了一个晚上才完成这个流程，现在的单细胞研究，动辄是几万个细胞，这个流程要是搞几个星期就不得了了。我后面会继续讲解关于这个问题的两个解决方案，第一个是对细胞亚群里面的单细胞进行抽样，第二个是 Python (pySCENIC). 教程，开启多线程！