单细胞个性化分析之转录因子篇
转录因子(Transcription Factors，TFs)，是指能够以特定序列与基因专一性结合，从而保证目的基因以特定的强度在特定的时间与空间表达的蛋白质分子。转录因子通过识别特定的DNA序列来控制染色质和转录，以形成指导基因组表达的复杂系统。许多转录因子充当着主调节因子和选择基因的角色，控制着细胞类型的决定、发育模式和特定途径控制(如免疫反应)的过程。
单细胞技术把视角拉到了单个细胞的研究精度，在研究单个细胞基因表达、轨迹发生、相互作用等的同时，细胞转录状态的变化也进入了研究者的目光。在复杂多变的细胞微环境中，受到外界刺激的细胞是如何通过转录因子调节基因表达，从而调整细胞的转录状态以适应新的环境，尤其在肿瘤微环境中转录状态的转变，成为了单细胞数据分析不可或缺的一环。对于单细胞数据转录因子分析，SCENIC[1]无疑是研究者的首选。当然，限于运行速度的关系，python版本的pyscenic也越来越多的投入到实际的运用中，本篇将从SCENIC软件着手，讲解单细胞研究中转录因子带来的生物学意义。
转录因子介绍
基因转录有正调控和负调控之分。如细菌基因的负调控机制是当一种阻遏蛋白(repressor protein)结合在受调控的基因上时，基因不表达；而从靶基因上去除阻遏蛋白后，RNA聚合酶识别受调控基因的启动子，使基因得以表达，这是正调控。这种阻遏蛋白是反式作用因子。而顺式作用因子则指的是基因上与反式作用因子结合的对基因表达起调控作用的基因序列。
图1 TF调控基因表达原理示意图
定义
转录因子(transcription factor)是一群能与基因5`端上有特定序列专一性结合，从而保证目的基因以特定的强度在特定的时间与空间表达的蛋白质分子。
结合位点
转录因子的结合位点（transcription factor binding site，TFBS）是转录因子调节基因表达时，与mRNA结合的区域。按照常识，转录因子（transcription factor，TF）的结合位点一般应该分布在基因的前端，但是，新的研究发现，人21和22号染色体上，只有22％的转录因子结合位点分布在蛋白编码基因的5'端。
图2 转录因子结合位点
分类
真核生物在转录时往往需要多种蛋白质因子的协助。一种蛋白质是不是转录结构的一部分往往是通过体外系统看它是否是转录起始所必须的。一般可将这些转录所需的蛋白质分为三大类：
(1)RNA聚合酶的亚基
RNA聚合酶的亚基，它们是转录必须的，但并不对某一启动子有特异性。
(2)与RNA聚合酶结合形成起始复合物的
某些转录因子能与RNA聚合酶结合形成起始复合物，但不组成游离聚合酶的成分。这些因子可能是所有启动子起始转录所必须的，但亦可能仅是譬如说转录终止所必须的。但是，在这一类因子中，要严格区分开哪些是RNA聚合酶的亚基，哪些仅是辅助因子，是很困难的。
(3)仅与其靶启动子中的特异顺序结合
某些转录因子仅与其靶启动子中的特异序列结合。如果这些序列存在于启动子中，则这些顺式因子是一般转录结构的一部分。如果这些序列仅存在于某些种类的启动子中，则识别这序列的因子是在这些特异启动子上起始转录所必须的。
图3 转录因子结合位点motif
黑腹果蝇的RNA聚合酶需要至少两个转录因子方能起始转录。其中一个是B因子，它与含TATA盒的部位结合。人的因子TFⅡD亦和类似的部位结合。同样，CTF(CAAT结合因子)则与腺病毒的主要晚期启动子中与CAAT盒同源的部位相结合。结合在上游区的另一个转录因子是USF(亦称MLTF)，则可以识别腺病毒晚期启动子中靠近-55的序列。转录因子Sp1则能和GC盒相结合。在SC40启动子中有多个GC盒，位于-70到-110之间，它们均能和Sp1相结合。然而含有GC盒的不同的DNA序列与Sp1的亲和力却各不相同。可见GC盒两侧的序列对Sp1-GC盒的结合情况究竟如何能影响转录。有时候需要几个转录因子才能起始转录。例如胞苷激酶的启动子需要Sp1与GC盒结合和CTF与CAAT盒结合;腺病毒晚期启动子需要TFⅡD与TATA盒结合和USF与其邻近部位相结合。以上所述的因子是一般转录都需要的，似乎并没有什么调节功能。另一些转录因子则可以调控一组特殊基因的转录。热休克基因就是一个很好的例子。真核生物的热休克基因在转录起始点的上游15bp处有一个特殊序列。HSTF因子仅在热休克细胞中有活性。它与包括热休克基因在内的一段DNA相结合，所以这个因子的激活可以引起约20个基因起始转录。在这里，转录因子和RNA聚合酶Ⅱ之间关系很类似细菌的σ因子与核心酶之间的关系。
转录因子是一种具有特殊结构、行使调控基因表达功能的蛋白质分子，也称为反式作用因子。植物中的转录因子分为二种，一种是非特异性转录因子，它们非选择性地调控基因的转录表达，如大麦 (Hordeum vulgare) 中的HvCBF2 (C-repeat/DRE binding factor 2) 。还有一种称为特异型转录因子，它们能够选择性调控某种或某些基因的转录表达。典型的转录因子含有DNA结合区 (DNA-binding domain)、转录调控区 (activation domain)、寡聚化位点(oligomerization site) 以及核定位信号 (nuclear localization signal) 等功能区域。这些功能区域决定转录因子的功能和特性。DNA结合区带共性的结构主要有：1）HTH 和 HLH 结构：由两段α-螺旋夹一段β-折叠构成，α-螺旋与β-折叠之间通过β-转角或成环连接，即螺旋-转角-螺旋结构和螺旋-环-螺旋结构。2）锌指结构： 多见于 TFIII A 和类固醇激素受体中，由一段富含半胱氨酸的多肽链构成。每四个半光氨酸残基或组氨酸残基螯合一分子 Zn2+ ，其余约 12-13 个残基则呈指样突出，刚好能嵌入 DNA 双螺旋的大沟中而与之相结合。3）亮氨酸拉链结构：多见于真核生物 DNA 结合蛋白的 C 端，与癌基因表达调控有关。由两段α - 螺旋平行排列构成，其α - 螺旋中存在每隔 7 个残基规律性排列的亮氨酸残基，亮氨酸侧链交替排列而呈拉链状，两条肽链呈钳状与 DNA 相结合。
转录调控区
同一家族的转录因子之间的区别主要在转录调控区。
转录调控区包括转录激活区 (transcription activation domain) 和转录抑制区 (transcription repression domain) 二种。近年来，转录的激活区被深入研究。它们一般包含DNA结合区之外的30-100个氨基酸残基，有时一个转录因子包含不止一个转录激活区。如控制植物储藏蛋白基因表达的VP1和PvALF转录因子，它们的N-末端酸性氨基酸保守序列都具有转录激活能力，与酵母转录因子GCN4和病毒转录因子的VP16的酸性氨基酸转录激活区有较高同源性。典型的植物转录因子激活区一般富含酸性氨基酸、脯氨酸或谷氨酰胺等，如GBF (G-box binding factor) 含有的GCB盒 (GBF conserved box) 激活结构域。
转录抑制区
也是转录因子调控表达的重要位点，但是对其作用机理研究尚不深入。可能的作用方式有三种：1）与启动子的调控位点结合，阻止其它转录因子的结合；2）作用于其它转录因子，抑制其它因子的作用；3）通过改变DNA的高级结构阻止转录的发生。
转录因子必须在核内作用，才能起到调控表达的目的。因此，转录因子上的核定位序列是其重要的组成部分。一般一个或多个核定位序列在转录因子中不规则分布，同时也存在不含核定位序列的转录因子，它们通过结合到其它转录因子上进入细胞核。核定位序列一般是转录因子中富含精氨酸和赖氨酸残基的区段。目前，水稻中的GT-2、西红柿中的HSFA1-2、玉米的O2和碗豆的PS-IAA4和6等转录因子中的核定位序列都已被鉴定。
绝大多数转录因子结合 DNA前需通过蛋白质-蛋白质相互作用形成二聚体或多聚体。所谓二聚体化就是指两分子单体通过一定的结构域结合成二聚体，它是转录因子结合DNA时最常见的形式。由同种分子形成的二聚体称同二聚体，异种分子间形成的二聚体称异二聚体。这种多聚体的形成是转录因子上的寡聚化位点 (oligomerization site) 相互作用的结果，寡聚化位点的氨基酸序列很保守，大多与DNA结合区相连并形成一定的空间构象。除二聚化或多聚化反应，还有一些调节蛋白不能直接结合DNA，而是通过蛋白质－蛋白质相互作用间接结合DNA，调节基因转录，这样就形成了一个表达调控的复合物。
转录因子的作用
是通过和顺式因子的互作来实现的。这段序列可以和转录因子的DNA结合域实现共价结合，从而对基因的表达起抑制或增强的作用。
SCENIC简介
SCENIC是2017年11月发表在Nature Methods 期刊的一种单细胞转录因子分析方法，也是目前进行单细胞转录因子分析的主流软件，该软件在进行数据分析的同时也能得到可视化结果图。另外，SCENIC是一款开源软件，可以免费下载使用，目前软件有R和python两个版本，每个版本都配备了详细的使用说明（软件官网https://scenic.aertslab.org/）。但有一点需要特别注意，该软件是有物种限制的，目前只能分析人、小鼠和果蝇的数据。
SCENIC分析原理
在输入单细胞基因表达量矩阵后，SCENIC经过以下三个步骤完成转录因子分析：第一步是构建共表达网络、第二步是构建TF-targets网络、第三步是计算Regulons活性，每一个步骤都由一个专门的软件包完成。我们来看一下各步骤详解。这个地方要结合之前的转录因子背景来看。
第一步：GENIE3——共表达网络构建
第一步由GENIE3或GRNBoost软件完成，这里以GENIE3为例介绍。GENIE3 (GEne NetworkInference with Ensemble of trees) ，基于树的基因网络推理，是一种从基因表达数据推断基因调控网络的方法。软件以单细胞基因表达量矩阵为输入文件，以每个目标基因 (gene) 为输出，以转录因子 (TF) 为输入，构建P个随机森林树（P=矩阵中基因数量），并计算每个TF与gene之间的重要性评分 (IM) ，最终可以获得TF-genes共表达模块。最后删除IM低于阈值的基因关系，过滤基因数低于50的模块。如果觉得不太好理解，下图也为大家做了简化~
图4 GENIE3网络构建过程
这里需要注意一点，构建TF-TG共表达网络默认是TF激活基因表达。
第二步：RcisTarget——motif富集及靶基因预测
从第一步获得了TF-genes共表达网络，但这个网络只是基于TF和gene表达量相关性推测的，TF和gene之间是否现实存在调控关系还需要进一步确证。确证的方法主要从TF功能结构入手，从图3可以看出，TF是通过直接与DNA特定位置序列而发挥作用的，因此可以通过反向查看gene上是否存在TF结合的motif序列来验证TF与gene的靶向关系。
这一步可以借助RcisTarget软件完成，该软件运行必备两个数据库：1）gene-motif排名数据库：为每个motif提供所有gene的排名（~分数）；2）motif-TF注释数据库：对每一个motif注释其所对应的TF。由于不同物种基因组不一样，导致每个motif对应靶基因不同，因此针对不同物种需要构建不同的数据库，软件目前配置了人、小鼠、果蝇数据库.
那么具体验证过程，首先基于gene-motif数据库，每个motif对模块中所有基因进行累积，模块中的基因排名越靠前，累积曲线越高，曲线下面积 (AUC) 越大，表明motif在该模块中的富集程度越高，然后对每个模块选取显著富集的motif，并预测其靶基因，最终综合TF-genes模块和靶基因预测结果，构成一个包含了TF和靶基因的基因调控网络模块 (regulons)。
图5 RcisTarget Regulons构建过程
第三步：AUCell——Regulons活性定量
第三步就是Regulons活性定量。这一步由AUCell软件完成，AUCell是一种新的方法，允许在scRNA-seq数据中识别具有活性基因调控网络的细胞。
实际分析过程中，输入到AUCell的是一个基因集，输出的是每个细胞中的基因集“活性” (AUC, Area Under Curve)。在SCENIC中，这些基因集即Regulons中所有基因，针对每个细胞，将细胞中所有基因按照表达量从高到低进行排序，根据Regulons中的基因在序列中的位置，计算累计曲线面积 (AUC) ，即为Regulons在细胞中的活性。
但由于不同regulons包含的基因不同，它们之间的AUC值不具有可比较性，因此基于AUC值在所有细胞中的双峰分布特征，增加了Regulons“on/off”的概念，认为双峰之间的低谷为判断Regulons活性开放的阈值，如果AUC值小于阈值，则判定为该Regulons在该细胞中未开放，即未发挥调控作用。最终获得每个Regulons在每个细胞中的开放性热图。
图6 AUCell 细胞Regulons活性计算过程
结果示例
结果可视化也是SCENIC软件功能之一，默认输出的结果图主要是以下两种，一个是基于AUC值绘制的Regulons活性热图，图中每一行表示一个Regulons，每一列为一个细胞，颜色表示AUC值大小，这类热图主要用于比较同一Regulons在不同细胞中的活性，寻找Regulons特异性活化的细胞亚群；另一种是基于AUC二值矩阵绘制的Regulons开放性热图，这类热图可用于寻找细胞中开放的Regulons，有利于细胞亚群功能判断。
图7 Regulons活性热图（左），Regulons开放性热图（右）
基因网络推断方法之间的优劣势
scRNA-seq 数据集的可用性推动了许多从这些数据进行网络推理的方法的发展，这些方法使用不同类型的模型，包括高斯图模型、信息论方法、随机森林、常微分方程和布尔网络。方法在包含伪时间或估算的去噪信号方面也有所不同。其中一些方法专门对 scRNAseq 数据的统计特性进行建模，而另一些方法则是对现有方法的bulk数据的改进。
文章Identifying strengths and weaknesses of methods for computational network inference from single cell RNA-seq data[2]中比较了11 种基因网络推断方法，根据方法的整体性能和计算要求，将它们分为三大类：top（PIDC、MERLIN、SCENIC、PEARSON）、mid（Inferelator、 SCODE、LEAP、Scribe）和bottom（knnDREMI、SILGGM）。
基因网络推断在文献中的运用
第一篇是发表杂志Cell上的文章Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma[3]，运用到了单细胞空间技术，其实在发现新的细胞亚群TSK中运用SCENIC分析进行佐证。
图8 Clustered heatmap of all recovered transcription factor (TF) regulons from SCENIC analysis (n = 370 regulons) across tumor keratinocyte subpopulations
图9 Heatmap of single-cell regulon scores inferred by SCENIC (g, genes; extended, SCENIC-annotated additional genes).
这里注意，文章中仅仅是对上皮细胞进行subcluster，发现了肿瘤样本中有新的亚群TSK，通过scenic的分析，发现了新的亚群中存在独有的调控网络，AP1 和 ETS 转录因子家族成员为可能控制 TSK 的 TF，至于每个TF控制的基因功能，需要认真研究。
在NC的文章Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma[4]中也运用到了SCENIC
图10 细胞亚群转录活性
研究者定位到了两个细胞亚群 iCAFs and mCAFs，然后针对性的对这两个细胞亚群进行SCENIC分析，取那些在两个细胞亚群有统计学差异的TF的全部细胞的AUC值进行热图可视化，分析两种状态下调控网络的差异。
2020年12月NC的食管鳞状细胞癌微环境文章Immune suppressive landscape in the human esophageal squamous cell carcinoma microenvironment[5]中，也是对DC细胞的亚群进行SCENIC分析，查看每个亚群调控网络的差异。
图11 Heat-map of the t values of AUC scores of expression regulation by transcription factors of the indicated clusters, as estimated using SCENIC
当然例子还有很多，但是核心只有一点，文章对SCENIC主要是用来识别不同条件下相同细胞类型调控网络的变化，或者更细化一点，用来分析相同细胞不同亚群调控网络的差异.
写在后面
细胞的基因调控对于控制细胞的行为非常重要，而细胞本身的调控又受到各种微环境的影响，尤其疾病状态下对细胞状态的转变是如何发生的，以及是否可逆，这些研究对认识生物学功能及疾病治疗都有深刻的借鉴意义。
参考文献
[1]Aibar, S., González-Blas, C., Moerman, T. et al.SCENIC: single-cell regulatory network inference and clustering. Nat Methods 14, 1083–1086 (2017). https://doi.org/10.1038/nmeth.4463
[2]Stone M , Sg M C , Fotuhi S A , et al. Identifying strengths and weaknesses of methods for computational network inference from single cell RNA-seq data. Cold Spring Harbor Laboratory, 2021.
[3] Ji A L , Rubin A J , Thrane K , et al. Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma[J]. Cell, 2020.
[4]Chen Z , Zhou L , Liu L , et al. Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma[J]. Nature Communications.
