CN108090326A - 单细胞网络调控关系的构建方法 - Google Patents

Abstract

本发明公开的单细胞网络调控关系的构建方法，包括：对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值；根据剔除了异常值的样本空间，得到功能单元相关性矩阵；对所述功能单元相关性矩阵进行稀疏化处理，得到稀疏化处理过的功能单元相关性矩阵；根据所述稀疏化处理过的功能单元相关性矩阵，进行网络调控关系计算；根据所述网络调控关系，划分调控功能模块；以hdf5格式存储用于识别调控功能模块的拓扑相似性矩阵。

Description

单细胞网络调控关系的构建方法

技术领域

本发明涉及生物技术领域，特别涉及一种单细胞网络调控关系的构建方法。

背景技术

细胞内转录本，因组织特异性，时空变化，微环境改变而引起的调控关系的改变。为了识别同质型细胞内基因间的调控关系、预测功能调控模块，通常通过构建基于表达调控数据的网络调控关系来实现。

目前在基于表达调控数据的网络调控关系构建的方法中，主要以布尔网络构建表达相关关系矩阵和加权共表达网络为主。这些方法的第一步往往都是需要计算相关关系矩阵，但其计算的时间复杂度和空间复杂度均为O(n²)，即其与数据的维度密切相关，呈指数级别增长。随着基因功能研究的深入，表达调控数据的维度不断增大，导致计算量过大。

同时，现有的技术往往还具有以下缺点：

1.相关矩阵的计算耗时，存储驻留内存，耗费大量的资源；

2.模块识别基于分层聚类算法，距离矩阵计算耗时，划分功能模块效率低；

3.无法支持较大维度数据的计算，在维度超过五万的情况下，目前一般的计算节点都无法满足；

4.对样本群体本身没有有效的划分，模型的鲁棒性不足。

发明内容

为解决以上问题，本发明提供一种单细胞网络调控关系的构建方法，包括步骤：

对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值；

根据剔除了异常值的样本空间，得到功能单元相关性矩阵；

对所述功能单元相关性矩阵进行稀疏化处理，得到稀疏化处理过的功能单元相关性矩阵；

根据所述稀疏化处理过的功能单元相关性矩阵，进行网络调控关系计算；

根据所述网络调控关系，划分调控功能模块；

以hdf5格式存储用于识别调控功能模块的拓扑相似性矩阵。

优选的，所述对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值，具体实施为：

对表达调控数据的样本空间进行主成分分析，对样本空间进行降维；

根据降维后样本空间的马式距离，剔除异常值，其中，剔除异常值的方法为：

记D为降维后样本空间的马式距离；

假定D²服从χ²分布，将满足D²>χ² _0.05的样本值作为异常值进行剔除。

优选的，所述根据剔除了异常值的样本空间，得到功能单元相关性矩阵，具体实施为：

将剔除了异常值的样本空间划分为不同的区块，其中，划分方法为，将不同的样本的第一个基因或功能单元划分为第一个区块，将不同的样本的第i个基因或功能单元划分为第i个区块；

计算功能单元相关性矩阵，将矩阵的自相关矩阵的求解变为相关系数矩阵的求解；

以hdf5格式存储求解得到的相关系数矩阵的上三角部分。

优选的，所述对所述功能单元相关性矩阵进行稀疏化处理，具体实施为：

通过对功能单元进行聚类，得到对功能单元的第一划分，具体实施为：

对于给定的p个基因或功能单元的样本空间X＝{x₁，x₂,…,x_n}，x_i为单个样本的全部p个基因的表达水平，求解其对应的相关系数矩阵S_p×p；

执行以下步骤：

遍历全部样本，for i＝1to n:Ci＝x_i，其中，n为样本的数量，x_i为第i个样本，Ci为中间变量；

Cn1,Cn2＝arg max_(Cu,Cv)S(Cu,Cv)，其中S(Cu,Cv)为Cu和Cv的相关系数矩阵，Cn1,Cn2为待求变量；

合并相似样本，Cj＝Cn1∪Cn2，其中，Cj为待求变量；

C:＝C\{Cn1,Cn2}∪Cj，得到C个基因类别，完成对功能单元相关性矩阵的稀疏化处理。

优选的，所述进行网络调控关系计算，具体实施为：

对所述C个基因类别，

计算邻接矩阵：

a_ij＝|cor(x_i,x_j)|^β

其中，对每个类别，xi和xj∈ci；

计算拓扑相似性(Topological overlap similarity，TOM)矩阵：

其中，ki为调控网络节点中心性的度量，k_i＝∑_ja_ij。

优选的，所述根据所述网络调控关系，划分调控功能模块，具体实施为：

根据每个c_i下的拓扑相似性矩阵，通过对功能单元进行聚类，得到对功能单元的第二划分，具体实施为：

执行以下步骤：

for i＝1to m:Ei＝Yi，其中，m为c_i下的样本的数量，Yi为c_i下的第i个样本，Ei为中间变量；

En1,En2＝arg max_(Eu,Ev)TOM(Eu,Ev)，其中TOM(Eu,Ev)为Eu和Ev的拓扑相似性矩阵，En1,En2为待求变量；

Ej＝En1∪En2，其中，Ej为待求变量；

E:＝E\{En1,En2}∪Ej，得到E个功能单元的聚类，视为c_i下的调控功能模块。

本发明的一些有益效果可以包括：

通过本发明的方法，能够更为快速和准确的识别同质型细胞内，基因间的调控关系，有效地预测功能调控模块。解决了数据规模较大的情况下，即高维度数据(维度大于2万)情况下，计算的效率与存储问题。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1为本发明实施例中一种单细胞网络调控关系的构建方法的流程图；

图2为本发明实施例中将剔除了异常值的样本空间划分为不同的区块的示意图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

图1为本发明实施例中一种单细胞网络调控关系的构建方法的流程图，如图1所示，包括步骤：

步骤S101、对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值；

步骤S102、根据剔除了异常值的样本空间，得到功能单元相关性矩阵；

步骤S103、对所述功能单元相关性矩阵进行稀疏化处理，得到稀疏化处理过的功能单元相关性矩阵；

步骤S104、根据所述稀疏化处理过的功能单元相关性矩阵，进行网络调控关系计算；

步骤S105、根据所述网络调控关系，划分调控功能模块；

步骤S106、以hdf5格式存储用于识别调控功能模块的拓扑相似性矩阵。

通过本发明的方法，通过剔除异常值、稀疏化处理以及使用hdf5格式进行压缩存储，能够更为快速和准确的识别同质型细胞内，基因间的调控关系，有效地预测功能调控模块。解决了数据规模较大的情况下，即高维度数据情况下，计算的效率与存储问题。

表达调控数据往往具有极强的时空性，即不同品种，不同的发育阶段和组织均有其特异性的表达，因此，在进行调控网络构建时，需要尽可能的剔除异常值，避免样本群体混杂而导致网络关系识别错误。在本发明的一个实施例中，所述对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值，具体实施为：

对表达调控数据的样本空间进行主成分分析，对样本空间进行降维；

根据降维后样本空间的马式距离，剔除异常值，其中，剔除异常值的方法为：

记D为降维后样本空间的马式距离；

假定D²服从χ²分布，将满足D²>χ² _0.05的样本值作为异常值进行剔除。

实施为，首先利用主成分分析方法，对数据降维，采用如下(matlab伪代码如下),对数据X(n*p矩阵，n为样本数目，p为维度)：

(1)X＝mean(X)，对数据中心化处理

(2)X＝std(X)，数据标准化

(3)t0＝init_vector()，初始化处理处理

迭代求解：

(4)pvector＝norm(X’*t0/(t0’*t0))

(5)tnew＝X*pvector

(6)若tnew-t0<1e^-8则结束，否则t0＝tnew返回(4)

然后根据降维后样本间的计算马式距离，剔除奇异样本，方法如下：

(1)D＝sqrt((T-u)’*cov(T)^-1*(T-u))

其中，u为均值向量，此时，D²近似服从卡方分布,利用D²>χ² _0.05剔除样本空间中的异常值。采用此方法，能够有效的降维，减少存储和计算消耗，尤其在求解马氏距离时，若维度较大，则计算T的协方差的逆矩阵会大大增加时间和存储的消耗。

为了降低对存储空间的需求，在本发明的一个实施例中，所述根据剔除了异常值的样本空间，得到功能单元相关性矩阵，具体实施为：

将剔除了异常值的样本空间划分为不同的区块，其中，划分方法为，将不同的样本的第一个基因或功能单元划分为第一个区块，将不同的样本的第i个基因或功能单元划分为第i个区块，划分方式如图2所示，为本发明实施例中将剔除了异常值的样本空间划分为不同的区块的示意图；

计算功能单元相关性矩阵，将矩阵的自相关矩阵的求解变为相关系数矩阵的求解；

由于求解得到的相关系数矩阵为对称矩阵，所以可以以hdf5格式存储求解得到的相关系数矩阵的上三角部分。

常用的稀疏化处理方法，往往通过相关关系的统计检验，给定显著性水平，从而实现稀疏化，但这这往往会导致给定阈值的不同，而导致无法对2个群体进行有效的比较。为了解决这一问题，通过对功能单元合并聚类的方法，先对功能模块做一个初步的划分。在本发明的一个实施例中，所述对所述功能单元相关性矩阵进行稀疏化处理，具体实施为：

通过对功能单元进行聚类，得到对功能单元的第一划分，具体实施为：

对于给定的p个基因或功能单元的样本空间X＝{x₁，x₂,…,x_n}，x_i为单个样本的全部p个基因的表达水平，求解其对应的相关系数矩阵S_p×p；

执行以下步骤：

遍历全部样本，for i＝1to n:Ci＝x_i，其中，n为样本的数量，x_i为第i个样本，Ci为中间变量；

Cn1,Cn2＝arg max_(Cu,Cv)S(Cu,Cv)，其中S(Cu,Cv)为Cu和Cv的相关系数矩阵，Cn1,Cn2为待求变量；

合并相似样本，Cj＝Cn1∪Cn2，其中，Cj为待求变量；

C:＝C\{Cn1,Cn2}∪Cj，得到C个基因类别，完成对功能单元相关性矩阵的稀疏化处理。

为了减少计算量，针对每个子分类(小矩阵)，对数据在局部水平上进行鉴定功能模块。在本发明的一个实施例中，所述进行网络调控关系计算，具体实施为：

对所述C个基因类别，

计算邻接矩阵：

a_ij＝|cor(x_i,x_j)|^β

其中，对每个类别，xi和xj∈ci；

计算拓扑相似性(Topological overlap similarity，TOM)矩阵：

其中，ki为调控网络节点中心性的度量，k_i＝∑_ja_ij。

在本发明的一个实施例中，所述根据所述网络调控关系，划分调控功能模块，具体实施为：

根据每个c_i下的拓扑相似性矩阵，通过对功能单元进行聚类，得到对功能单元的第二划分，具体实施为：

执行以下步骤：

for i＝1to m:Ei＝Yi，其中，m为c_i下的样本的数量，Yi为c_i下的第i个样本，Ei为中间变量；

En1,En2＝arg max_(Eu,Ev)TOM(Eu,Ev)，其中TOM(Eu,Ev)为Eu和Ev的拓扑相似性矩阵，En1,En2为待求变量；

Ej＝En1∪En2，其中，Ej为待求变量；

E:＝E\{En1,En2}∪Ej，得到E个功能单元的聚类，视为c_i下的调控功能模块。

通过本发明的方法，能够更为快速和准确的识别同质型细胞内，基因间的调控关系，有效地预测功能调控模块。解决了数据规模较大的情况下，即高维度数据(维度大于2万)情况下，计算的效率与存储问题，本研究方法能够支撑上十万甚至上亿级别变量数目的计算，可以同时针对编码基因的多种转录本，表观修饰区域及其他功能元件进行处理，且处理速度更快，同时占用的内存小，数据通过hdf5存储具有高效的读写能力，同时能满足超大维度数据的计算。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种单细胞网络调控关系的构建方法，其特征在于，包括步骤：

对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值；

根据剔除了异常值的样本空间，得到功能单元相关性矩阵；

对所述功能单元相关性矩阵进行稀疏化处理，得到稀疏化处理过的功能单元相关性矩阵；

根据所述稀疏化处理过的功能单元相关性矩阵，进行网络调控关系计算；

根据所述网络调控关系，划分调控功能模块；

以hdf5格式存储用于识别调控功能模块的拓扑相似性矩阵。

2.根据权利要求1所述的方法，其特征在于，所述对表达调控数据的样本空间进行聚类，并剔除样本空间中的异常值，具体实施为：

对表达调控数据的样本空间进行主成分分析，对样本空间进行降维；

根据降维后样本空间的马式距离，剔除异常值，其中，剔除异常值的方法为：

记D为降维后样本空间的马式距离；

假定D²服从χ²分布，将满足D²>χ² _0.05的样本值作为异常值进行剔除。

3.根据权利要求1所述的方法，其特征在于，所述根据剔除了异常值的样本空间，得到功能单元相关性矩阵，具体实施为：

将剔除了异常值的样本空间划分为不同的区块，其中，划分方法为，将不同的样本的第一个基因或功能单元划分为第一个区块，将不同的样本的第i个基因或功能单元划分为第i个区块；

计算功能单元相关性矩阵，将矩阵的自相关矩阵的求解变为相关系数矩阵的求解；

以hdf5格式存储求解得到的相关系数矩阵的上三角部分。

4.根据权利要求3所述的方法，其特征在于，所述对所述功能单元相关性矩阵进行稀疏化处理，具体实施为：

通过对功能单元进行聚类，得到对功能单元的第一划分，具体实施为：

对于给定的p个基因或功能单元的样本空间X＝{x₁，x₂,…,x_n}，x_i为单个样本的全部p个基因的表达水平，求解其对应的相关系数矩阵S_p×p；

执行以下步骤：

遍历全部样本，for i＝1 to n:Ci＝x_i，其中，n为样本的数量，x_i为第i个样本，Ci为中间变量；

Cn1,Cn2＝arg max_(Cu,Cv)S(Cu,Cv)，其中S(Cu,Cv)为Cu和Cv的相关系数矩阵，Cn1,Cn2为待求变量；

合并相似样本，Cj＝Cn1∪Cn2，其中，Cj为待求变量；

C:＝C\{Cn1,Cn2}∪Cj，得到C个基因类别，完成对功能单元相关性矩阵的稀疏化处理。

5.根据权利要求4所述的方法，其特征在于，所述进行网络调控关系计算，具体实施为：

对所述C个基因类别，

计算邻接矩阵：

a_ij＝|cor(x_i,x_j)|^β

其中，对每个类别，xi和xj∈ci；

计算拓扑相似性矩阵：

<mrow> <msub> <mi>TOM</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mo>&amp;Sigma;</mo> <mi>u</mi> </msub> <msub> <mi>a</mi> <mrow> <mi>i</mi> <mi>u</mi> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mi>a</mi> <mrow> <mi>u</mi> <mi>j</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>a</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <mi>min</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>a</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> </mrow>

其中，ki为调控网络节点中心性的度量，k_i＝∑_ja_ij。

6.根据权利要求5所述的方法，其特征在于，所述根据所述网络调控关系，划分调控功能模块，具体实施为：

根据每个c_i下的拓扑相似性矩阵，通过对功能单元进行聚类，得到对功能单元的第二划分，具体实施为：

执行以下步骤：

for i＝1 to m:Ei＝Yi，其中，m为c_i下的样本的数量，Yi为c_i下的第i个样本，Ei为中间变量；

En1,En2＝arg max_(Eu,Ev)TOM(Eu,Ev)，其中TOM(Eu,Ev)为Eu和Ev的拓扑相似性矩阵，En1,En2为待求变量；

Ej＝En1∪En2，其中，Ej为待求变量；

E:＝E\{En1,En2}∪Ej，得到E个功能单元的聚类，视为c_i下的调控功能模块。