利用单细胞转录和基因敲除数据推断基因调控网络的方法

摘要

本发明公开了一种利用单细胞转录和基因敲除数据推断基因调控网络的方法，通过分析基因敲除前后的稳态表达数据对基因进行分类，作为先验知识以降低时空复杂度并提升推断准确度；利用单细胞转录数据计算基因间分布距离，结合基因分类的结果为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子；对基因敲除数据进行相应的算法分析以去除部分假阳性判断，弥补分析动态数据的缺陷；本发明有效解决了分析时间序列单细胞数据中高计算复杂度的问题，提升了推断基因调控网络的准确度。

说明书

技术领域

本发明涉及基因调控网络研究与分析领域，特别涉及一种利用单细胞转录和基因敲除数据推断基因调控网络的方法。

背景技术

基因调控网络作为一种解释分析基因数据的工具，能够揭示基因、蛋白质和小分子之间的调控关系、了解生物细胞内的生理活动和功能、通路中的相互作用以及如何使机体发生变化。从单细胞转录的角度研究基因调控网络能够揭示细胞详细的表达动态和功能关系，阐明不同关键生理过程中细胞间变异的功能作用。但单细胞表达数据是具有高时间分辨率的时间过程数据，给推断基因调控网络算法带来了很高的时空复杂度的同时，也丧失了对稳态基因表达水平的分析，这将会降低最终实验结果的精确度。

发明内容

为解决上述现有技术的不足，本发明提供一种利用单细胞转录和基因敲除数据推断基因调控网络的方法。

为达到上述目的，本发明的技术方案为：

一种利用单细胞转录和基因敲除数据推断基因调控网络的方法，包括：利用基因敲除技术分别收集多个细胞中的稳态基因转录表达数据，刺激各个细胞后，在多个时间点收集单个细胞的转录表达数据；分析基因敲除前后的稳态表达数据对基因进行分类；利用单细胞转录数据计算基因间分布距离，结合基因分类的先验知识为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子；对基因敲除数据进行相应的算法分析去除部分假阳性判断；根据所推断的基因间关系确定基因调控网络。

其中，利用基因敲除技术分别收集多个细胞中的稳态基因转录表达数据，刺激各个细胞后，在多个时间点收集单个细胞的转录表达数据的步骤包括：

利用GNW模拟基因敲除实验，得到野生型和敲除型生物体的稳态基因表达数据，其中野生型表达数据来自于没有经历任何突变的生物体的原始菌株。敲除数据来自于原始菌株敲除或削弱一个或多个基因后得到的生物体的菌株；

在多个时间点从多个单细胞中采集的基因表达值构成时间序列单细胞表达数据。

其中，分析基因敲除前后的稳态表达数据对基因进行分类，包括步骤：

利用基因敲除技术依此扰动各基因，并收集相应时刻的所有基因的稳态表达值，将基因分为四类，分别为URGRRGNRGISG，具体的方法如下：

(a)不受调控的调控基因(URG)：扰动其他所有基因对当前基因均无影响，但该基因的扰动将影响其他基因的表达水平；

(b)受调控的调控基因(RRG)：扰动其他任一基因时当前基因的表达水平有所波动，且该基因的扰动将影响其他基因的表达水平；

(c)非调控基因(NRG)：扰动其他任一基因时当前基因的表达水平有所波动，但该基因的扰动不影响其他任一基因的表达水平；

(d)独立基因(ISG)：扰动其他所有基因对当前基因均无影响，且该基因的扰动不影响其他任一基因的表达水平。

其中，利用单细胞转录数据计算基因间分布距离，结合基因分类的先验知识为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子，包括步骤：

利用时间序列单细胞表达数据计算基因间分布距离；

用给定时间窗口的其他基因的表达分布来“预测”下一个时间窗口目标基因的表达分布，即建立每个基因多时间点的多维回归模型；

模型中包含基因间作用关系因子，即为所求的解向量，利用带惩罚项的最小二乘法求得解向量，较大的因子表示对应的调控推断具有更高的置信度。

其中，利用单细胞转录数据计算基因间分布距离，包括步骤：

时间序列单细胞表达数据包含多个基因、多个时间点、多个细胞。令g为基因数，n为所测量时间点的数量，C_t为第t个时间点样本中的细胞数(t>矩阵元素是基因j的转录表达值，即在第k个时间点，基因j的mRNA分子在第i个细胞中的数量。

利用单细胞基因表达数据集中所包含的信息，特别是基因表达分布的变化进行GRN推断，首先要计算基因在两个时间点的表达分布距离来量化每个个体基因表达的时间变化情况。基因j在t时刻的距离量化如下：

DD_j，t＝max|F_t+1(A_j)-F_t(A_j)|

其中，F_t(A_j)表示A_j在时间t上的累积分布函数，即为基因j从0时刻到>j，t表示基因j在t时刻与t+1时刻的表达变化水平。

其中，结合基因分类的先验知识为每个基因建立多时间点的多维回归模型，包括步骤：

用给定时间窗口的其他基因的表达分布来“预测”下一个时间窗口目标基因的表达分布，以解释某基因受其他基因的调控情况。基因j在t+1时刻的线性关系如下所示：

DD_j，t+1＝α_1，jDD_1，t+α_2，jDD_2，t+…+α_g，jDD_g，t

其中，α_g，j表示基因g对基因j的调控关系作用因子；将基因j在所有时刻的线性关系表示出来得到以下矩阵：

其中，g表示网络中的基因数，n表示所测量时间点的数量，DD_j，n-1表示第n-1>

矩阵中的α向量即为所需求解的各基因间作用关系因子，利用基因分类的结果先确定部分α向量的值，再通过数学计算求得所有的解向量。

其中，对基因敲除数据进行相应的算法分析去除部分假阳性判断，包括步骤：

分析基因敲除数据去除初步GRN推断中误判的两基因间的直接调控关系；

分析基因敲除数据去除初步GRN推断中将间接调控误判为直接调控的推断；

其中，分析基因敲除数据去除初步GRN推断中误判的两基因间的直接调控关系，包括步骤：

分别在野生型和敲除后的菌株中采集稳态的基因表达值。若后者的除敲除基因以外的其他基因的表达水平较前者有较大变化，则说明所敲除的基因对该基因有调控作用，反之，没有调控作用。因此，利用这样的方法可以去除部分初步推断的假阳性结果，具体步骤如下：

(1)利用基因芯片技术采集野生菌株所有基因的表达水平；

(2)利用基因敲除技术依此敲除每个基因，并同时采集该菌株所有基因的表达水平；

(3)对比敲除前后各基因表达水平的变化情况，采用一定的度量方式来判断其是否有变化，从而判断当前基因是否受敲除基因的调控作用。

本文采取两者差的绝对值来来作为度量，假设GK为敲除菌株的表达值，GW 为野生菌株的表达值，则：GK_i,j-GW_j>α时，认为基因j的表达水平有变化，则基因j受到基因i的调控。其中GK_i,j表示敲除基因i时基因j的表达值，>

其中，分析基因敲除数据去除初步GRN推断中将间接调控误判为直接调控的推断，包括步骤：

首先需要确定所需判断的调控关系，即调控网络中存在间接调控关系的两基因间可能存在或不存在的直接调控关系，称之为“不确定调控”。实现这一步骤首先要确定基因调控网络的上限GU和下限GL集合，上限集合可由前期推断的GRN得到，下限集合则是将所有“不确定调控”去掉。通过后续步骤不断更新上下限集合，直至两集合全等。

通过断开“不确定调控”的间接调控路径后扰动调控基因，观察当前GRN 中基因表达水平的变化，来确定两基因间是否真正存在直接调控关系。为了实现这一过程，需要找到能够断开间接路径的最佳基因组。

其中，找到能够断开间接路径的最佳基因组，包括步骤：

将以下规则确定的集合成为“边分离”：

1.S₁(i，j)＝GU中i的子代∩GU中j的祖先

2.S₂(i，j)＝GU中i的后代∩GU中j的亲本

3.S₃(i，j)＝GU中i的后代∩GU中j的祖先

当网络中存在多条间接调控时，依此找到他们的边分离，再统计共同的集合数量最大的集合，即为所要敲除的最佳基因组。

其中，根据所推断的基因间关系确定基因调控网络，包括步骤：

根据真实网络中的调控数量确定一个参数因子；

根据参数因子，按照所推断的关系因子表中的概率值确定最终的基因调控网络。

相对于现有技术，本发明的有益效果为：本发明所述的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法，通过分析基因敲除前后的稳态表达数据对基因进行分类，作为先验知识以降低时空复杂度并提升推断准确度；利用单细胞转录数据计算基因间分布距离，结合基因分类的结果为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子；对基因敲除数据进行相应的算法分析以去除部分假阳性判断，弥补分析动态数据的缺陷；本发明有效解决了分析时间序列单细胞数据中高计算复杂度的问题，提升了推断基因调控网络的准确度。

附图说明

图1是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法的流程示意图；

图2是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法中,所采用的数据--时间序列单细胞表达谱的示意图。

图3是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法中，某一个基因在多个时间点的单细胞表达值。

图4是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法中，若干基因在若干时间点的基因表达变化水平折线图。

图5是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法中，某一个基因调控网络的最终推断结果示意图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步更详细的描述。显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

参阅图1-5，图1是本发明提供的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法的流程示意图。该方法的步骤包括：

S110：利用基因敲除技术分别收集多个细胞中的稳态基因转录表达数据，刺激各个细胞后，在多个时间点收集单个细胞的转录表达数据。

所属步骤S110中包括：

1.利用GNW模拟基因敲除实验，得到野生型和敲除型生物体的稳态基因表达数据，其中野生型表达数据来自于没有经历任何突变的生物体的原始菌株。敲除数据来自于原始菌株敲除或削弱一个或多个基因后得到的生物体的菌株；

2.在多个时间点从多个单细胞中采集的基因表达值构成时间序列单细胞表达数据。

S120：分析基因敲除前后的稳态表达数据对基因进行分类。

利用基因敲除技术依此扰动各基因，并收集相应时刻的所有基因的稳态表达值，将基因分为四类，分别为URGRRGNRGISG，具体的方法如下：

(a)不受调控的调控基因(URG)：扰动其他所有基因对当前基因均无影响，但该基因的扰动将影响其他基因的表达水平；

(b)受调控的调控基因(RRG)：扰动其他任一基因时当前基因的表达水平有所波动，且该基因的扰动将影响其他基因的表达水平；

(c)非调控基因(NRG)：扰动其他任一基因时当前基因的表达水平有所波动，但该基因的扰动不影响其他任一基因的表达水平；

(d)独立基因(ISG)：扰动其他所有基因对当前基因均无影响，且该基因的扰动不影响其他任一基因的表达水平。

S130：利用单细胞转录数据计算基因间分布距离，结合基因分类的先验知识为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子。

利用单细胞转录数据计算基因间分布距离，包括步骤：

时间序列单细胞表达数据包含多个基因、多个时间点、多个细胞。令g为基因数，n为所测量时间点的数量，C_t为第t个时间点样本中的细胞数(t>矩阵元素是基因j的转录表达值，即在第k个时间点，基因j的mRNA分子在第i个细胞中的数量。

利用单细胞基因表达数据集中所包含的信息，特别是基因表达分布的变化进行GRN推断，首先要计算基因在两个时间点的表达分布距离来量化每个个体基因表达的时间变化情况。用公式(1)对基因j在t时刻的分布距离进行量化。

DD_j,t＝max|F_t+1(A_j)-F_t(A_j)|>

其中，F_t(A_j)表示A_j在时间t上的累积分布函数，即为基因j从0时刻到>j，t表示基因j在t时刻与t+1时刻的表达变化水平。

利用公式(1)得到的分布距离，用给定时间窗口的其他基因的表达分布来“预测”下一个时间窗口目标基因的表达分布，以解释某基因受其他基因的调控情况。基因j在t+1时刻的分布距离表示为公式(2)中的线性关系。

DD_j,t+1＝α_1,jDD_1,t+α_2,jDD_2,t+…+α_g，jDD_g，t>

其中，α_g，j表示基因g对基因j的调控关系作用因子。将基因j在所有时刻的线性关系表示出来得到公式(3)中的矩阵。矩阵中的α向量即为所需求解的各基因间作用关系因子，利用基因分类的结果先确定部分α向量的值，再通过数学计算求得所有的解向量。

其中，g表示网络中的基因数，n表示所测量时间点的数量，DD_j，n-1表示第n-1>

S140：对基因敲除数据进行相应的算法分析去除部分假阳性判断。

所属步骤S140中包括：

1.分析基因敲除数据去除初步GRN推断中误判的两基因间的直接调控关系。

分别在野生型和敲除后的菌株中采集稳态的基因表达值。若后者的除敲除基因以外的其他基因的表达水平较前者有较大变化，则说明所敲除的基因对该基因有调控作用，反之，没有调控作用。因此，利用这样的方法可以去除部分初步推断的假阳性结果，具体步骤如下：

(1)利用基因芯片技术采集野生菌株所有基因的表达水平；

(2)利用基因敲除技术依此敲除每个基因，并同时采集该菌株所有基因的表达水平；

(3)对比敲除前后各基因表达水平的变化情况，采用一定的度量方式来判断其是否有变化，从而判断当前基因是否受敲除基因的调控作用。

本文采取两者差的绝对值来来作为度量，假设GK为敲除菌株的表达值，GW 为野生菌株的表达值，则：GK_i,j-GW_j>α时，认为基因j的表达水平有变化，则基因j受到基因i的调控。其中GK_i,j表示敲除基因i时基因j的表达值，>

2.分析基因敲除数据去除初步GRN推断中将间接调控误判为直接调控的推断。

首先需要确定所需判断的调控关系，即调控网络中存在间接调控关系的两基因间可能存在或不存在的直接调控关系，称之为“不确定调控”。实现这一步骤首先要确定基因调控网络的上限GU和下限GL集合，上限集合可由前期推断的GRN得到，下限集合则是将所有“不确定调控”去掉。通过后续步骤不断更新上下限集合，直至两集合全等。

通过断开“不确定调控”的间接调控路径后扰动调控基因，观察当前GRN 中基因表达水平的变化，来确定两基因间是否真正存在直接调控关系。为了实现这一过程，需要找到能够断开间接路径的最佳基因组。

将以下规则确定的集合成为“边分离”：

1.S₁(i，j)＝GU中i的子代∩GU中j的祖先

2.S₂(i，j)＝GU中i的后代∩GU中j的亲本

3.S₃(i，j)＝GU中i的后代∩GU中j的祖先

当网络中存在多条间接调控时，依此找到他们的边分离，再统计共同的集合数量最大的集合，即为所要敲除的最佳基因组。

S150：根据所推断的基因间关系确定基因调控网络，包括步骤：

根据真实网络中的调控数量确定一个参数因子；

根据参数因子，按照所推断的关系因子表中的概率值确定最终的基因调控网络。

区别于现有技术，本发明所述的一种利用单细胞转录和基因敲除数据推断基因调控网络的方法，通过分析基因敲除前后的稳态表达数据对基因进行分类，作为先验知识以降低时空复杂度并提升推断准确度；利用单细胞转录数据计算基因间分布距离，结合基因分类的结果为每个基因建立多时间点的多维回归模型，采用数学方法计算模型中的关系因子；对基因敲除数据进行相应的算法分析以去除部分假阳性判断，弥补分析动态数据的缺陷；本发明有效解决了分析时间序列单细胞数据中高计算复杂度的问题，提升了推断基因调控网络的准确度。

以上仅为本发明的实施方式，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。