风景:单细胞调节网络推理和聚类

（一个)第一步，利用GENIE3 (Random Forest)或GRNBoost (Gradient Boosting)推断转录因子与候选靶基因之间的共表达模块。每个模块由一个转录因子及其预测靶标组成，完全基于共表达。（b)第二步，用RcisTarget对每个共表达模块进行分析，找出富集的基序;只保留TF基序丰富的模块和目标。每个TF及其潜在的直接目标是一个规则。（c)第三步，使用AUCell评估每个细胞中每个调控子的活性，AUCell计算恢复曲线下的面积(Area Under the recovery Curve)，整合一个调控子中所有基因的表达等级。AUCell分数用于生成规则活动矩阵。该矩阵可以通过为每个规则子设置AUC阈值来进行二值化，这将确定规则子在哪些细胞中“开启”。（d)调节活性矩阵可用于细胞聚类(例如t-SNE)，从而根据调节子网的共享活动识别细胞类型和状态。

（一个)用AUCell对小鼠脑数据集进行评分的多个基因集的AUC分布。AUC代表每个细胞中调控子或基因标记的活性。具有“活动”规则的单元格的选择是基于数据集中所有单元格的AUC分布。理想的情况是，一个调控子或基因信号仅在一部分细胞中活跃，将返回双峰分布(如神经元或少突胶质细胞)，或长尾分布(如小胶质细胞)。相反，正态分布更可能发生在非差异表达的基因集上。这种情况在这里通过随机基因集(例如，从数据集中随机取的基因名称)和类似管家的基因集(在大多数细胞中检测到的基因)来说明。AUCell自动探索AUC分数的分布，并为每个基因集计算几个可能的阈值:(1)密度曲线的拐点，对于双峰分布(蓝色)的理想情况，这通常是一个很好的选择;(2)全局分布(灰色/绿色)子分布的异常值(调整两个或三个正态分布的混合物，红色或粉红色)。与这些分布相关的阈值在直方图上以虚线表示;每个基因集的选定阈值用较粗的连续线突出显示。请注意，当前版本中的阈值选择并不详尽，我们强烈建议检查AUC直方图，并在需要时手动调整阈值。 We also recommend being cautious about gene-sets with few genes (10-15) and thresholds that are set extremely low. (b基于表达的t-SNEs (Zeisel等人的小鼠大脑数据集)根据给定基因集的每个细胞的AUC着色。当单元格AUC大于分配阈值时，使用粉红色/红色阴影，否则使用蓝色阴影。（c)使用Zeisel等人提供的细胞类型作为正确标签计算灵敏度和特异性，以及自动AUCell分配阈值。用Scran归一化后的AUC(左)和基因集表达的平均值²⁶(右)。

基因集的来源: Cahoy等。⁴⁵(超过1000个基因的基因签名，上一行)，Lein等。⁴⁶(星形胶质细胞和神经元少于100个基因的基因签名)，Lavin等。⁴⁷(小胶质细胞:小胶质细胞与其他组织内巨噬细胞)。

（一个)比较由SCENIC调控子活性和单独由TF表达引起的细胞聚类。左:聚类二进制规则活动矩阵。右:基于92个tf归一化表达的聚类(簇内大小因子归一化与scran，热图颜色:中位数以基因为中心)。（b)几个关键规则的AUC直方图。AUC允许分裂具有高活性和低活性调控子的细胞群。（c)表达式矩阵上的t-SNE(与SCENIC相同的输入:UMI计数，没有进一步的规范化)和(d)的t-SNE。这两种t-SNEs都是基于pca的，并根据每个细胞中检测到的基因数量(表达量大于0)着色。基于SCENIC的聚类有效地纠正了簇内偏差，而神经元(表达的基因较多)和胶质细胞(表达的基因较少)之间的真正生物学差异不受影响。^12日,48（e) SCENIC与其他鉴别细胞类型相关tf的方法的比较(详见方法)。右边的条形图显示了每种方法识别的tf数量(白色)和验证集中的tf数量(彩色)。与差异表达分析相比，SCENIC保留了更多的转录因子。（f) SCENIC可检测到的tf比例。Zhou等在蛋白水平上比较小鼠脑内检测到的tf的维恩图。⁴⁹，在Zeisel等人的scRNA-seq数据集中，⁴⁹以及RcisTarget数据库中可用的tf(即已知motif)。

（一个)根据小鼠脑数据推断出的与小胶质细胞相关的调控，可以根据相关TF (iRegulon中构建的网络)的结合基序进行总结。先前发表的小胶质细胞特征中包含的基因(Lavin等人)。⁴⁷)用较大的字体表示;节点的颜色表示调节器的数量(浅:少，深:多)。预测的小胶质细胞网络包含许多众所周知的小胶质细胞命运和/或小胶质细胞激活的调节因子，包括PU.1、Nfkb、Irf和AP-1/Maf。（b)当我们将预测的小胶质细胞网络与先前发表的小鼠阿尔茨海默病模型中小胶质细胞“激活”的基因特征进行比较时，我们发现小胶质细胞网络在阿尔茨海默病的进展过程中被强烈激活，而神经元网络被下调，这表明小胶质细胞网络捕获了相关的调控程序。所示的图是GSEA分析的结果，该分析针对阿尔茨海默病(AD)小鼠模型中与每种野生型细胞类型(调控基因的结合)相关的网络，以及基于基因表达的排序。数据集由Gjoneska等人提供。⁵⁰: p25诱导后2周和6周，CK-p25 AD小鼠模型与CK同窝对照相比海马转录变化。

（一个(Zeisel等人)在100个随机的小鼠大脑数据集上运行的SCENIC结果与在整个数据集上运行的结果相似(左列:与出版物中提供的细胞类型匹配的细胞状态，每个状态的相似相关规则，以及目标的显著重叠)。然后用AUCell在所有细胞上对100个细胞推断的GRN进行评估，以确认该网络可推广到未包含在GRN推断中的细胞。在右列，同样的方法，但模拟稀疏数据集(UMI计数矩阵除以3并截断，结果中位数为1121个检测到的基因)。在稀疏数据集中没有检测到许多相关的tf(因此会错过相关的规则)，但SCENIC仍然能够识别主要的细胞类型。（b)在小鼠大脑数据集上运行10次SCENIC，评估SCENIC结果的稳定性(Zeisel等人)。上图:二进制规则活性矩阵和t-SNEs(由作者的细胞类型标记着色)。10个二进制矩阵的聚合说明了运行结果的稳定性，并且大多数顶级调节器在10/10运行中被发现。（c) t-SNE对在100个细胞的10个随机子集上运行SCENIC所产生的AUC调控活性的影响。

（一个)二元正则活动矩阵。（b)三个基于grn的人类中间神经元亚群，具有定义这些组的主要DNA基序和转录因子(VIP中间神经元的NFIX和PV中间神经元的ESRRG)，以及在每个组中确定的已知标记(注意tf本身也在各自的集群中上调)。下框:中间神经元亚型标记的表达。^51、52在第一个图中，中间神经元根据z得分最高的标记着色。其余的图是根据基因表达来着色的(灰色:无表达，暗红色:高表达，黄色:中等表达)。

基于合并表达矩阵(Z-score归一化)的人类和小鼠脑细胞的分层聚类。聚类先按种类，再按细胞类型对细胞进行分组。缩略图显示图2 c，其中SCENIC生成单元类型的主要聚类。

（一个)少突胶质细胞瘤数据集上批量效应去除方法的比较。原始表达矩阵(第一行)上的t-SNEs和扩散图，经过原始肿瘤用Combat或Limma(第2-3行)校正后，或来自SCENIC的二元活性矩阵(第4行)上的t-SNEs和扩散图。细胞根据原始肿瘤或GRN活性着色(红色:星形细胞样调节，绿色:少突胶质细胞样调节，蓝色:与细胞周期或干细胞相关的调节)。（b)少突胶质细胞瘤数据集的简化二进制规则活动矩阵(SCENIC的输出)。高亮显示的调控(彩色TF名称)分别是少突胶质细胞或星形胶质细胞的特征。

二元活性矩阵突出转录因子组及其基序。在得到的t-SNE中，细胞根据三组选定的调控子(红、绿、蓝)的平均二进制活性着色，这对应于分化轨迹中的三种主要状态:OPC网络，由Bcl6和共调节因子驱动;由Bach2等因素驱动的中间网络;成熟的少突胶质细胞网络，包括Etv1和Nfel2l2等转录因子。Sox10被发现是所有少突胶质细胞亚型的调节因子。在一组成熟的少突胶质细胞(蓝色细胞)中，检测到两个稍微落在分化轨迹之外的外群:具有神经元特性的少突胶质细胞(B)和具有AP-1激活特征的少突胶质细胞(C)。在t-SNE的每个簇旁边，显示了富集的氧化石墨烯项。在方框中，t-SNEs以另一种方案着色，显示其他少突胶质细胞网络和状态。请注意，尽管是处于分化中的细胞，但主导网络是相当离散的。这表明在这些主要状态之间的转换必须迅速发生，因为只有少数细胞被发现处于转换状态。

（一个) t-SNE对二元活度矩阵的影响。（b)黑色素瘤数据集的二进制规则活动矩阵。热图上方的颜色条表示肿瘤的起源;与细胞周期相关的规则(绿色)，MITF^低侵袭性(粉色)和MITF^高，通常被称为增殖(蓝色)状态放大。（c)三个最主要的网络的详细信息。“Confirmed by ChIP-seq”:打勾表示该调控子在ChIP-seq数据集中呈现相同转录因子的富集靶标。（d) TF表达与调控蛋白活性的比较。对于四种转录因子:AUC值的直方图，以及选择的截止点(橙色虚线)。在第二列中，AUC值超过截止值的单元格以蓝色显示。在这些细胞中，调控子被认为是活跃的(即二进制活动矩阵中的“1”)。在第三列中，实际AUC值用于为单元格上色。第四列显示了转录因子本身的表达。TFs的判别能力远低于法规的判别能力。（e)已知黑色素瘤标志物的表达。注意，MITF^低簇显示WNT5A、LOXL2和ZEB1表达上调(均为已知的侵袭状态标志物)^53、54)，并与先前发表的侵袭性基因特征显著相关(图S19)。然而，与“经典的”侵袭性细胞状态不同，这种MITF^低state保留SOX10表达式。(f)黑色素瘤体积特征与单细胞状态的比较。左图:来自侵袭性和增殖性黑色素瘤状态(例如Hoek和Verfaillie)的大量特征在基于单细胞状态的基因排名中各自的上调或下调侧显着富集。在GSEA中，排名(x轴)是基于mitf高和mitf低状态之间的对比。右图:类似的GSEA分析，但现在只显示NES分数。这个分析是倒数的前一个，其中排名是基于散装样品的对比，测试的签名是来自单细胞mitf -高和mitf -低状态之间的差异表达。

注意，在MITF中的细胞^高state也有高活性的STAT和IRF下游靶标。这在大量样本中很难检测到，因为恶性细胞与肿瘤浸润淋巴细胞(TIL)的复杂混合物中STAT和IRF也起着重要作用⁵⁵。在这里，我们发现MITF^高细胞本身具有比MITF更高的STAT活性^低细胞(我们从分析中排除了所有良性细胞，包括免疫细胞)。这对于解释和预测免疫治疗耐药性具有重要意义，因为这些具有高STAT和IRF活性的癌细胞可能对免疫治疗最敏感。事实上，最近的一项研究发现JAK-STAT-IRF轴是免疫治疗中两个主要靶点表达的驱动因素:PD-L1和PD-L2;一方面抑制抗肿瘤免疫反应，但另一方面增加抗pd (L)1免疫疗法的反应⁵⁵。同样值得注意的是，14例肿瘤活检中有2例的mitf低“侵袭性”状态，均从辅助淋巴结切除。与AP-1和TEAD因子驱动的体外侵袭状态不同，这种状态具有不同的转录因子，包括NFATC2和NFIB，我们证实这些转录因子在早期转移性黑色素瘤细胞中表达(即通过免疫组织化学在前哨淋巴结的初始小肿瘤中表达)。利用NFATC2敲除后的基因表达分析补充Fig.13)，我们发现NFATC2是AP-1和TEAD靶基因的转录抑制因子。因此，这些观察结果表明，NFATC2可能作为一个转录断裂，细胞需要克服它才能切换到成熟的侵袭性细胞状态。NFATC2本身就是一个JUN目标⁵⁶，并可能构成负反馈机制。此前在乳腺癌中也观察到类似的NFATC2抑制功能⁵⁷。我们认为这是在本数据集中AP-1和TEAD未被检测为调控的(生物学)原因。值得注意的是，对于TEAD，还有一个原因是它不能作为调控子被检测到，因为我们的SCENIC实验选择了与其靶标共表达的tf，而TEAD在蛋白质水平上受到调控。

（一个)根据与细胞周期相关的基因集的z分数来鉴定细胞周期细胞。（b)显示细胞周期基因集在少突胶质细胞瘤数据集上的z分数的热图。热图顶部的蓝/红条突出了通过三种方法选择的循环细胞(红色):(1)SCENIC的E2F1调控的AUC分数;(2) Tirosh et al. G1/S评分和Tirosh et al. G2/M评分;Tirosh等人采用允许的临界值(如果细胞的G1/S和G2/M评分高于细胞群中平均值的两倍，则将其归类为循环细胞);Tirosh等人采用了更严格的临界值(如果细胞的G1/S和G2/M评分高于细胞群中平均值的四倍，则将其归类为循环细胞);(3)基于GO基因集的Z-score。（c, d)不同方法识别循环细胞的能力比较:(c)回收的细胞数，按真阳性率(TPR，也称为灵敏度或召回率)分类。彩色条表示通过该方法鉴定的细胞周期细胞的数量，白色条表示在所选簇(细胞周期细胞数量最多的簇)中包含的非细胞周期细胞的数量。（d召回率vs精度(风景区_一个:高可信度单元，SCENIC_B:较低的信心/规律活动)。

补充图3我。在这里，我们显示了不同阶段黑色素瘤进展的额外活检(RGP:径向生长期，VGP:垂直生长期，SLN:前哨淋巴结(小转移)，转移(完全)转移)。

NFATC2和NFIB最强烈的阳性信号可见于前哨淋巴结。

（一个) NFIB和NFATC2在cosmos黑色素瘤细胞系中的Z-score归一化表达。根据类似MITF-low状态的关键标记的表达，选择A375进行敲除。（b)在A375黑色素瘤细胞系中敲低NFATC2。NFATC2敲除后差异表达基因的GSEA图:预测的NFATC2靶点在NFATC2敲除后显著上调。(氟)ChIP-seq信号在选定调控中的富集。c-d: MITF和STAT1 ChIP-seq信号在MITF、STAT和NFATC2(即各自调控中基因的调控区域)的预测靶细胞上的聚集图。e-f:调控子上的ChIP-seq信号(在TSS周围10kb的窗口内预测的CRMs)与启动子近端区域(随机选择的基因和调控子外共表达的基因)中TF基序的发生率的比较。与对照相比，调节子中基因的富集证实，与单个/替代方法(例如仅共表达或基序分析)相比，SCENIC增加了寻找直接靶标的特异性。

使用SCENIC软件对小鼠视网膜细胞的Drop-seq数据进行分析，方法是在约11K细胞的子集上运行GENIE3，并对所有细胞的GRN结果进行评估。（一个) tSNE根据Macosko等人注释的预期细胞类型着色。(b，c)。主要网络根据活动规律，以红绿蓝配色方案表示。在这些tSNEs中，显示为灰色的细胞没有考虑到颜色。在iRegulon中发现的重要基序对应于在m神经胶质细胞中识别的规则和相应的GO术语。已确定的Muller胶质细胞的主调控因子(如Sox8/9、Hes1、Rax、Nr1h4、Srf和Nr2e1)已在文献中得到证实^{58 - 61}，说明即使在稀疏数据上也可以推断出正确的网络。子采样方法在稀疏数据集上特别有趣，比如由Drop-seq或10x产生的数据集，因为最好的质量单元可以用来推断GRN，然后这个高质量的GRN可以在所有单元上得分。

（一个) GRNBoost和GENIE3的性能比较。左作为GRNBoost的生物学验证，我们设计了一个基于基因集的精度和召回基准。我们从Zeisel等人的小鼠脑表达数据中询问GRNBoost和GENIE3是否预测了一系列转录因子(Olig1, Sox10, Rel, Lef1, Neurod1, Dlx1)的相似靶基因集。利用GENIE3推断的网络，我们为6个tf中的每一个构建了目标基因的排序列表。对于每个排名列表，我们使用GOrilla⁶²查询每个TF目标列表的top富集GO项。对于6个最热门的GO术语，我们都咨询了QuickGO⁶³获得关联蛋白注释基因符号集(按小鼠分类单元过滤)，仅保留表达矩阵中存在的基因，获得6个基因符号“主列表”。这些主列表最后用于计算GENIE3和GRNBoost对上述每种tf预测的目标基因集的精度和召回分数。正确的为了对GRNBoost的时间性能进行基准测试，并作为概念验证，我们使用了Chromium Megacell演示数据集，其中包含来自胚胎小鼠大脑的100多万个细胞。我们随机取了3000、10000和100000个细胞的子集来推断grn。（b)跨多次运行的稳定性评估。每个散点图比较每个TF在两个独立运行(具有相同或不同数量的细胞)中目标的排名。