基于单细胞多组学数据无监督构建基因调控网络

在单细胞分辨率下识别基因调控网络（GRNs，gene regulatory networks）一直是一个巨大的挑战，而单细胞多组学数据的出现为构建GRNs提供了机会。

来自：Unsupervised construction of gene regulatory network based on single-cell multi-omics data of colorectal cancer

前置内容

DNA没有功能。 它包含一组指令（instructions），这些指令必须先转化为RNA，然后再转化为蛋白质。在大多数情况下，我们可以将 RNA 视为 DNA 和蛋白质之间的信使（messenger）。DNA由包含如何制造蛋白质指令的基因组成。蛋白质负责在细胞中执行生物功能，例如代谢葡萄糖为细胞产生能量。一般来说，人体内的每种蛋白质都由一个基因编码。

• 注意：基因不是蛋白质，而是DNA序列的一部分，可以编码蛋白质。基因是生物体内遗传信息的基本单位，包含了一系列的DNA序列，其中编码了蛋白质所需的氨基酸序列信息。基因通过转录作用，将DNA序列转录成RNA分子，然后再通过翻译作用，将RNA分子翻译成蛋白质。因此，基因可以被视为控制生物体内蛋白质合成的遗传信息的存储和传递的基本单位。

调控因子是指能够调控基因表达的蛋白质或其他分子。它们作为转录因子或其他类型的调控分子，可以结合到基因组的某些区域，例如启动子或增强子，来激活或抑制基因的转录。调控因子在细胞分化、发育和应对环境变化等生物学过程中发挥重要作用。在复杂的生物系统中，调控因子通常以复杂的网络形式相互作用，共同调控基因表达。

背景介绍

结直肠癌是威胁人类生命健康的主要癌症之一，在所有癌症中排名前五，发病率和死亡率仍在增长。然而，阐明生物分子之间复杂的调控关系对结直肠癌的治疗具有重要的研究意义。

基因调控网络GRN由调控因子（regulatory factors）和靶基因（target genes）组成。最常见的调控因子是转录因子（transcription factors），它可以通过控制复杂的细胞内相互作用来决定细胞的表型和命运。一旦调节过程发生变化，就会导致疾病出现。例如，如果AREG在GRN中高水平表达，则该组织很有可能是结直肠癌样本。因此，准确构建GRN可以有效地提高我们对驱动细胞类型和特定基因表达的调控机制的理解。

---

AREG是基因，虽然AREG编码的蛋白质可以作为一种信号分子参与调控，但在基因级别上，AREG本身是一个基因而不是调控因子。

---

在单细胞多组学数据中，同一组样本可以获得不同分子过程的全基因组数据，如转录组、表观基因组等组学数据。这些数据为单细胞分辨率下的疾病研究提供了多种模态信号，可以更准确、系统地分析临床疾病的发病机制，识别疾病的重要治疗靶点。

为了更好理解基因调控机制，近年来开发了多种算法得到GRNs。目前算法大多基于单一组学数据，然而现有研究表明，整合多种组学技术的数据可以显著提高医疗机构预测患者临床结果的准确性。基于多组学识别GRNs的研究比较少，比如LinkedSOMs，这些方法往往复杂度过高。因此，开发一种基于多组学数据的低复杂度策略对于探索GRNs至关重要。

方法

数据收集和组织

作者使用的所有数据集都是从具有相同基因集的公开数据库中下载的。scRNA-seq数据（GEND000035，2018年11月29日公开）来自Gene Expression Nebulas database，其中包含1150例人类结直肠癌癌症样本。scATAC-seq数据（GSE201336，2022年4月28日公开）来自Gene Expression Omnibus database，包含6例人类结直肠癌癌症样本。

对于scATAC-seq数据，Trimmomatic被用于移除低质量碱基，FastQC用于质量控制。应用Bowtie2将reads与参考基因组（hg19）进行对齐，所有样本的对齐率>92%。BedTools根据染色体的位置计算基因符号，scATAC-seq数据维度为8097×6，scRNA-seq数据维为8097×1150。其中，8097代表对齐的基因。

细胞类型标注

scATAC-seq提供了染色质可及性的信息，并揭示了单个细胞中基因的转录活性。直观地说，基因的转录活性和表达值之间的分布是一致的。在这项工作中，scATAC-seq和scRNA-seq从左到右连接，通过连接获得的多组学数据的维数为8097×1156。首先，去除在少于三个样本上表达的基因，和在少于200基因上表达的样本，记为$X\in R^{M\times N}$，其中$M$和$N$代表基因和样本。$x_{ji}$代表第$j$个基因处第$i$个样本的表达值。使用下式进行标准化：$$log(10000\times \frac{{\sum }_{j=1}^M{L}_{(x_{ji}\ne 0)}}{M})$$SingleR用于对标准化后的数据进行细胞类型注释。注释结果显示，大部分细胞是上皮细胞（epithelial cells），少数是星形胶质细胞（astrocytes）。上皮细胞用于随后的分析，星形胶质细胞被过滤掉，t-SNE可视化如图1A所示。

• 图1A：样本的t-sne可视化。
• 图1B：所有样本的CNV，颜色的变化表示CNV的degree，方框表示CNV更显著的区域。
• 图1C：截取了五类样本的CNV得分分布，以及得分>500的样本进行下游分析。

CNV（copy number variation）是由基因组重排引起的，导致DNA片段拷贝的增加或丢失，这种现象在癌症中普遍存在。因此，使用inferCNV分析上皮细胞的CNV，分析结果如图1B所示。图中的黑框显示大多数细胞的拷贝数增加或减少，但少数细胞没有显著的CNV。因此，为了挑选具有显著拷贝的细胞，作者使用k-means将上皮细胞分为五个组，分析每组细胞的CNV分数。并过滤得分<500的样本。具体计算如下：

• 数据被缩放为$[-1,1]$的$X^{\prime }$，其中$x_{i_{max}}$和$x_{i_{min}}$分别是第$i$个样本中所有基因表达的最小值和最大值，$$X_i^{\prime }=2\times \frac{x_{ji}-x_{i_{min}}}{x_{i_{max}}-x_{i_{min}}}-1,i=\mathrm{1...}N,j=\mathrm{1...}M$$
• 计算CNV：$$CN{V}_i=sum(X_i^{\prime }{)}^2,i=\mathrm{1...}N$$其中$CN{V}_i$代表第$i$个样本的CNV得分。这五组数据的CNV得分分布如图1C所示。其中，CNV得分<500的8个样本被认为是正常样本，其余497个样本被视为具有更显著拷贝数的样本。这497个样本的表达矩阵被提取用于下游GRN分析。

构建GRNs

构建GRN有两个主要步骤。首先，使用无监督人工神经网络自组织映射（SOM）将所有基因映射到不同的神经元（参考机器学习笔记本第三十三课一些经典的优化策略与神经网络变种），其中每个神经元代表一个基因集。SOM可以保存输入空间的拓扑结构，并将相似的基因映射在一起，从而确保这组基因的功能高度相似，这对于后续构建调控网络至关重要。其次，pySCENIC用于计算转录因子并调节每组基因。pySCENIC的创新之处在于引入转录因子基序序列，以验证通过统计方法推断的基因共表达网络，从而识别出由转录因子主导的高度可靠的GRN。

---

在TF-靶基因（转录因子-靶基因）图中，比如：

从数据上，它们都来自于scRNA-seq，三角形是一个节点，其实属于scRNA-seq中的某个基因，但在这个图中，它代表该基因编码的蛋白，这个蛋白是对靶基因有调控关系的转录因子。

---