文献介绍

文献题目： SCENIC: 单细胞调控网络推断和聚类
研究团队： Stein Aerts（比利时鲁汶大学）
发表时间： 2017-10-09
发表期刊： Nature Methods
影响因子： 25.0（2017年）
DOI： 10.1038/nmeth.4463

摘要

作者提出了 SCENIC，这是一种从单细胞 RNA-seq 数据中同时进行基因调控网络重建和细胞状态识别的计算方法（http://scenic.aertslab.org）。通过对来自肿瘤和大脑的单细胞数据集的综合分析，作者证明了顺式调控（cis-regulatory）分析可用于指导转录因子和细胞状态的识别。SCENIC 为驱动细胞异质性的机制提供了关键的生物学见解。

研究结果

细胞的转录状态源于其内在的基因调控网络（GRN），在这个网络中，有限数量的转录因子（TFs）和辅助因子（cofactors）相互调控，并调控其下游的靶基因（target）。近年来，单细胞转录组分析技术的进步为高分辨率识别转录状态及状态间的转换（例如在分化过程中）提供了令人兴奋的机遇。针对单细胞 RNA-seq 优化的统计技术和生物信息学方法已经带来了新的生物学见解，但目前尚不清楚是否能够确定稳定细胞状态背后的特异且稳健的 GRNs。鉴于在单细胞水平上，由于转录爆发和其他来源导致的基因表达随机变异，基因表达可能与转录因子输入的动态部分脱节，这一任务确实可能具有挑战性。已有一些方法可以从单细胞 RNA-seq 数据中推断共表达网络，但这些方法并未利用调控序列分析来预测转录因子与靶基因之间的相互作用。

作者认为，将顺式调控序列与单细胞基因表达联系起来可以克服数据丢失和技术变异，从而优化细胞状态的发现和表征。为此，作者开发了单细胞调控网络推断与聚类（single-cell regulatory network inference and clustering, SCENIC）方法，以绘制 GRNs，并通过评估每个细胞中 GRNs 的活性来识别稳定的细胞状态。SCENIC 的工作流程包括三个步骤（Fig. 1a, Supplementary Fig. 1 and see Online Methods）。第一步，使用 GENIE3 识别与转录因子共表达的基因集（Supplementary Fig. 1a）。由于 GENIE3 模块仅基于共表达，它们可能包含许多假阳性和间接靶标。为了识别假定的直接结合靶标，每个共表达模块都使用 RcisTarget 进行顺式调控 motif 分析（Supplementary Fig. 1b and see Online Methods）。仅保留具有正确上游调控因子显著基序富集的模块，并对其进行修剪以去除缺乏 motif 支持的间接靶标。作者将这些处理后的模块称为调控子（regulons）。

Figure 1: SCENIC 的工作流及其在鼠脑中的应用

(a) 在 SCENIC 工作流程中，首先使用 GENIE3 或 GRNBoost 推断转录因子（TFs）与候选靶基因之间的共表达模块。接着，RcisTarget 识别那些在靶基因中显著富集调控因子结合 motif 的模块，并创建仅包含直接靶基因的调控子（regulons）。AUCell 对每个细胞中每个调控子的活性进行评分，从而生成二值化的活性矩阵。细胞状态的预测基于调控子网络的共享活性。
(b) SCENIC 在小鼠大脑数据上的结果。聚类标签与参考文献中使用的标签一致；主调控因子（master regulators）的颜色与其控制的细胞类型相匹配。
(c) 通过文献确认的转录因子（A）或具有小鼠基因组信息库（Mouse Genome Informatics）中脑表型的转录因子（B）；其对应的富集 DNA-binding motifs 以 Logo 图形式展示。
(d) 基于二值化调控子活性矩阵的 t-SNE 图。每个细胞被分配为最活跃 GRN 的颜色。
(e) 不同聚类方法在该数据集上的准确性。

作为 SCENIC 的一部分，作者开发了 AUCell 算法来对每个细胞中每个调控子（regulon）的活性进行评分（Supplementary Figs. 1c and 2, and see Online Methods）。对于给定的调控子，通过比较不同细胞的 AUCell 评分，可以识别哪些细胞具有显著更高的子网络活性。生成的二值化活性矩阵降低了维度，这对下游分析非常有用。例如，基于该矩阵的聚类能够根据调控子网络的共享活性识别细胞类型和状态。由于调控子是作为一个整体进行评分，而不是使用单个基因的表达，这种方法对数据丢失具有鲁棒性（Supplementary Fig. 3）。

为了评估 SCENIC 的性能，作者将其应用于一个包含已知细胞类型的成年小鼠大脑 scRNA-seq 数据集（Fig. 1b–e）。该分析从 1,046 个初始共表达模块中筛选出 151 个调控子，这些调控子具有对应转录因子的显著富集 motifs（占初始转录因子的 7%）。对每个细胞的调控子活性进行评分后，揭示了预期的细胞类型（Fig. 1d,e），并列出了每种细胞类型的潜在主调控因子（e.g., the microglia network in Supplementary Fig. 4）。基于细胞类型的聚类（总体灵敏度为 0.88，特异性为 0.99，调整兰德指数（ARI）> 0.80）比许多专门用于单细胞聚类的方法更为准确。

为了评估 SCENIC 的鲁棒性，作者重新分析了小鼠大脑数据：包括完整数据集、随机选择 100 个细胞的子集以模拟小数据集、以及仅使用三分之一的测序读长以模拟低覆盖度数据。SCENIC 成功识别了仅由少数细胞代表的细胞类型（例如，小胶质细胞、星形胶质细胞或中间神经元中的 2 至 6 个细胞；Supplementary Fig. 5）。此外，预测的转录因子与细胞类型的关联与先前已知的功能一致（Fig. 1c），且其准确性优于标准分析流程（Supplementary Fig. 3e）。

为了验证在小鼠中间神经元中识别到的 Dlx1/2 网络，作者分析了人类大脑的 snRNA-seq 数据集（Supplementary Fig. 6）。在人类数据中，SCENIC 同样识别出一个由 DLX1/2 强烈驱动的中间神经元群，其识别 motif 与小鼠一致，并发现了一组保守的靶基因，包括 DLX1 本身（Fig. 2a,b）。接下来，作者将这一跨物种分析扩展到其他细胞类型。与基于标准化表达的标准聚类方法（其聚类结果主要由物种驱动；Supplementary Fig. 7）相比，SCENIC 分析有效地按细胞类型对细胞进行了分组（Fig. 2c）。这表明网络活性的评分具有鲁棒性，可用于克服批次效应或技术差异（Supplementary Fig. 3d）。

Figure 2: 跨物种比较神经元网络和细胞类型

(a) 从小鼠和人类大脑 scRNA-seq 数据中推断出的 DLX1/2 调控子。红色标记的基因在 GeneMANIA 中与 Dlx1/2 存在关联。
(b) 人类和小鼠 Dlx1/2 调控子在小鼠和人类单细胞数据上的互作活性。在每个 SCENIC t-SNE 图中，细胞根据相应的二值化调控子活性着色。插图展示了调控子的 AUCell 评分分布。
(c) 基于 GRN 活性的人类和小鼠大脑单细胞 RNA-seq 数据的联合聚类。着色的转录因子名称表示在人类和小鼠 SCENIC 分析中均被识别的调控子。

作者还将 SCENIC 应用于少突胶质细胞瘤（来自 6 个肿瘤的 4,043 个细胞）和黑色素瘤（来自 14 个病变的 1,252 个细胞）的 scRNA-seq 数据集，以识别复杂的细胞状态。由于肿瘤特异性突变和复杂的基因组异常，癌细胞状态的识别比正常细胞状态更具挑战性。标准聚类方法根据细胞的来源肿瘤进行分组（Fig. 3a,b），但 SCENIC 揭示了不同的结果。对于少突胶质细胞瘤，SCENIC 在多个肿瘤中识别出三种癌细胞状态（Fig. 3c–e），每种状态均由预期的转录因子驱动——包括少突胶质细胞样状态中的 SOX10/4/8、OLIG1/2 和 ASCL1；星形胶质细胞样状态中的 SOX9、NFIB、AP-1；以及循环细胞中的 E2F 和 FOXM1。

Figure 3: SCENIC 克服了肿瘤效应并揭示了癌症中的相关细胞状态和 GRNS

(a,b) 基于表达矩阵的 t-SNE 图，按来源肿瘤着色。
(c,d and f,g) 应用 SCENIC 后基于二值化活性矩阵（e,h）的 t-SNE 图。在 d 和 g 中，细胞根据 GRN 活性着色。
(i) 使用 NFATC2、NFIB、ZEB1 和 EPHA2 抗体对 25 例人类黑色素瘤进行免疫组织化学（IHC）分析。热图显示了每个样本中标记物阳性细胞的百分比。右侧展示了前哨淋巴结中 NFIB 的 IHC 代表性示例（更多图像见 Supplementary Fig. 13）。NA, not applicable。
(j) MITF 和 STAT1 ChIP-seq 信号在预测靶区域及随机选择的含有 MITF/STAT motif 的基因组区域上的聚合图，后者作为对照。

此外，将扩散映射应用于二值化的 SCENIC 矩阵（Supplementary Fig. 8）重建了一条从干细胞样到少突胶质细胞样和星形胶质细胞样分支的分化轨迹。需要注意的是，与正常的少突胶质细胞分化相比，这条路径代表了一种不同的“轨迹”（关于 5,069 个少突胶质细胞的 SCENIC 分析见 Supplementary Fig. 9）。作者在黑色素瘤数据中观察到了类似的肿瘤效应校正，SCENIC 识别出跨肿瘤的细胞群（Supplementary Fig. 10），包括由与少突胶质细胞瘤相似的转录因子驱动的循环细胞群（例如 E2F1/2/8 和 MYBL2；Fig. 3f–h and Supplementary Fig. 10）。与需要预先指定批次效应来源的专用批次效应去除方法（如 Combat 和 Limma；Supplementary Fig. 11）不同，SCENIC 通过使用生物学驱动的特征自动消除了肿瘤效应。

黑色素瘤细胞主要分为两组，一组对应于 MITF 高表达状态（典型的增殖状态），以 MITF 和 STAT/IRF 为关键调控因子；另一组对应于 MITF 低表达状态，其 WNT5A、LOXL2 和 ZEB1 表达上调——这些是已知的侵袭状态标志物（Supplementary Fig. 10e,f）。SCENIC 在 MITF 低表达状态中识别出两个新的转录因子：NFATC2（预测有 114 个靶基因）和 NFIB（预测有 15 个靶基因）。NFATC2 是 JNK/MAPK 通路中的转录抑制因子，参与黑色素瘤的去分化和免疫逃逸。而 NFIB 则与毛囊和黑色素细胞干细胞的干细胞行为相关，并在小细胞肺癌的转移进展中起重要作用。

为了进一步探索 NFATC2 和 NFIB 在 MITF 低表达状态中的潜在作用，作者对 25 个不同进展阶段的黑色素瘤标本进行了免疫组织化学分析。作者发现，在前哨淋巴结中 NFIB 和 NFATC2 的表达最高。这与 ZEB1 的表达共定位，表明这些标志物的表达与早期转移事件之间存在关联（Fig. 3i and Supplementary Fig. 12）。当作者在 A375（一种高表达 NFATC2 和 NFIB 的黑色素瘤细胞系；Supplementary Fig. 13）中使用 siRNA 敲低 NFATC2 时，作者发现 NFATC2 调控子中的基因显著上调（see Online Methods）。这与 NFATC2 先前被确立为抑制因子的角色一致。此外，参与细胞粘附和细胞外基质调控的基因以及多个先前发表的代表黑色素瘤侵袭状态的基因特征也出现上调（Supplementary Table 1），这表明 NFATC2 可能在疾病进展中确实发挥了重要作用。作为黑色素瘤调控子的第二次验证，作者使用 ChIP-seq 数据确认了 MITF 和 STAT 的预测靶基因（Fig. 3j）。

随着单细胞数据集规模的增加，作者提出了两种互补的方法来扩展网络推断。第一种方法是从子样本数据集中推断 GRN，并在评分步骤中使用 AUCell 对所有细胞进行评分。作者在一个包含超过 40,000 个小鼠视网膜单细胞的数据集上展示了这种方法（Supplementary Fig. 14）。第二种方法旨在使用更高效的机器学习和大数据处理解决方案。作者在 Apache Spark 上使用 Scala 实现了 GRNBoost，这是 GENIE3 的一个新变体，用梯度提升替代了随机森林回归。这一实现大大减少了推断 GRN 所需的时间（Supplementary Fig. 15），并将为在非常大的数据集上进行网络推断铺平道路，例如即将到来的人类细胞图谱。

SCENIC 是一种适用于单细胞 RNA-seq 数据分析的通用方法，它利用转录因子和顺式调控序列来指导细胞状态的发现。作者的结果表明，GRN 是识别细胞状态的稳健指南，而单细胞 RNA-seq 数据非常适合追踪由特定转录因子组合驱动的细胞类型特异性转录组的基因调控程序。

--------------- 结束 ---------------

注：本文为个人学习笔记，仅供大家参考学习，不得用于任何商业目的。如有侵权，请联系作者删除。