CN105160208A - 用于疾病亚型问题的基于网络的聚类方法 - Google Patents

Abstract

本发明涉及逆向研究疾病亚型领域，旨在提供用于疾病亚型问题的基于网络的聚类方法。该用于疾病亚型问题的基于网络的聚类方法包括过程：获得O-G矩阵以及基因调控网络；选取适用于具体问题的基于网络的距离定义，构建距离矩阵；运用k-medoids算法对O-G矩阵进行聚类分析，聚类时距离的选择用基于网络的距离。本发明对于特定的基因网络，将有更好的组间相似性，能更有效地还原三种亚型；当有大量的基因需要测定其表达数据时，本发明通过优先精确测量信息基因的表达数据，不会大大地削弱对于疾病亚型的鉴定效果。

Description

用于疾病亚型问题的基于网络的聚类方法

技术领域

本发明是关于逆向研究疾病亚型领域，特别涉及用于疾病亚型问题的基于网络的聚类方法。

背景技术

对于由基因变异导致的疾病的研究一直是一个非常热门的议题。在这些疾病中，很多疾病都有不同的亚型。所谓亚型(subtype)，是同一个疾病下的不同的种型。它们可能由不同的病因引起，并且有不同的临床表征。例如HIV有1型和2型，肿瘤亚型有卵巢癌、肺癌、子宫癌等。

对于疾病亚型的很多研究，现阶段还是集中在对于其病理的正向研究。而运用逆向工程技术(reverseengineering)，逆向研究疾病亚型也逐渐成为一个热门的话题。“逆向工程技术”是一个研究主体系统的过程。它通过研究主体系统来鉴定系统的各个成分以及它们之间的相互关联，并用另一种方式在更加抽象，更加上层的水平上对该系统进行代表。逆向工程技术在疾病亚型鉴定与分类方面研究上的一个非常重要的应用，就是通过已有的临床信息，例如基因表达数据(geneexpressiondata)等，运用包括聚类分析(clusteranalysis)在内的各种手段，反向研究并鉴定疾病的亚型。然而，由于基因的表达之间并非是独立的，而是会通过各种关系相互作用。因此，最终获取的基因表达数据，也应该是相互关联的。而在以往的研究中，人们只是使用传统的聚类方法，并没有考虑这种基因表达数据间的相互关联。因此，将有关基因的作用关系的信息加入聚类分析中是一个自然、新颖的想法并且值得一试。

发明内容

本发明的主要目的在于克服现有技术中的不足，提供能更好的将疾病亚型进行分类，更好的还原真实的疾病亚型的基于网络的聚类方法。为解决上述技术问题，本发明的解决方案是：

提供用于疾病亚型问题的基于网络的聚类方法，具体包括下述过程：

(1)获得O-G矩阵以及基因调控网络；

(2)选取适用于具体问题的基于网络的距离定义，构建距离矩阵；

(3)运用k-medoids算法对O-G矩阵进行聚类分析；聚类时距离的选择用基于网络的距离；

(4)得出最终关于疾病亚型的分类；

所述过程(1)具体包括下述步骤：

步骤A：根据基因调控网络(即基因-蛋白质调控网络，是一个细胞中DNA片段集合通过相互间的各种非间接作用，比如RNA作用以及蛋白质表达作用，来影响其mRNA以及蛋白质表达水平的相互关系)的特性(例如网络的平均出度、入度等参数)，构建随机的有向图来代表基因调控网络G(V，E)；其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i(转录因子i)；每条有向边e_ji∈E代表着“转录因子j调控基因i的转录”这种调控关系；

步骤B：根据产生的基因调控网络G(V，E)，对每个基因i建立激活函数f_i(·)，具体建立方式为：

对于任意的基因i∈V，i＝1，2，K，n，我们从G(V，E)中找出所有与i相邻且以i为有向边终点的点，构成影响因子集合{q₁，q₂，K，q_sn}；其中，q₁表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，q₂表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，q_1n表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，n表示基因调控网络中基因的数量；

确定解离常数k_ij，且k_ij从定义在[0.01，1]区间上的均匀分布中选取；

确定希尔系数n_ij，且n_ij服从[1，10]区间中的高斯分布函数

确定相对活性α_i，且α_i从定义在[0，1]区间上的均匀分布上采样；

步骤C：确定无噪声动态基因调控模型，即确定公式(2.1)的各个参数；

$\left\{\begin{array}{c}{F}_i^{mRNA}(x,y)=\frac{{\mathrm{dx}}_i}{dt}=m_i\·f_i\left(y\right)-{\mathrm{\λ}}_i^{RNA}\·x_i\\ F_i^{\mathrm{Pr}ot}(x,y)=\frac{{\mathrm{dy}}_i}{dt}=r_i\·x_i-{\mathrm{\λ}}_i^{\mathrm{Pr}ot}\·y_i\end{array}\right.---\left(2.1\right)$

式(2.1)中，x_i表示基因i的浓度；y_i表示蛋白质i的浓度；表示mRNA_i的浓度变化率；表示蛋白质i的浓度变化率；m_i表示基因i的最大转录速率；r_i表示mRNA_i的翻译速率；表示mRNA_i的降解速率；表示蛋白质i的降解速率；f_i(·)表示基因i的激活函数；

确定公式(2.1)中各个参数的具体方式为：mRNA的半衰期以及蛋白质的半衰期(以分钟为单位)从定义在[5，50]区间上的高斯分布上采样；

根据公式(2.9)，获得mRNA以及蛋白质的降解速率，最大转录速率m_i以及翻译速率r_i服从[0.01，0.011]区间上的均匀分布；

$\left\{\begin{array}{c}{\mathrm{\λ}}_i^{RNA}=\frac{ln2}{T_i^{RNA}}\\ {\mathrm{\λ}}_i^{\mathrm{Pr}ot}=\frac{1n2}{T_i^{\mathrm{Pr}ot}}\end{array}\right.---\left(2.9\right)$

式(2.9)中，表示mRNA_i的降解速率；表示蛋白质i的降解速率；mRNA的半衰期以及蛋白质的半衰期(以分钟为单位)；

步骤D：在获得了基因调控网络以及无噪声动态基因调控模型之后，选定mRNA浓度x(x₁，x₂，K，x_n)以及蛋白质浓度y(y₁，y₂，K，y_n)的初始值(可以令各个x_i和y_i服从[0，1]区间上的均匀分布，并随机选取作为初始值)，然后求解公式(2.1)，得到最终的基因表达数据；

所述过程(2)具体是指：根据过程(1)所获得的基因网络的拓扑关系G(V，E)，定义三种基于网络的距离，用于比较x₁(x₁₁，x₁₂，K，x_1n)与x₂(x₂₁，x₂₂，K，x_2n)的差别；其中x₁(x₁₁，x₁₂，K，x_1n)、x₂(x₂₁，x₂₂，K，x_2n)分别表示两个被试者P₁和P₂的mRNA浓度；

令G(V，E)代表该基因调控网络，其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i(转录因子i)；它关联的x_i表示该基因转录的mRNA_i浓度；令每条有向边e_ji∈E代表着“转录因子j调控基因i的转录”这种调控关系；记T_i表示与节点i相连的边数(即节点i的度)，I_i表示节点i的入度，O_i表示节点i的出度；

其中，基于网络的Jaccard距离定义为：

$EJ=\frac{2\underset{i=1}{\overset{n}{\Σ}}{T}_i\left|x_{1i}-x_{2i}\right|}{\underset{i=1}{\overset{n}{\Σ}}{T}_i\left(x_{1i}+x_{2i}\right)+\underset{i=1}{\overset{n}{\Σ}}{T}_i\left|x_{1i}-x_{2i}\right|}$ 公式(3.10)；

其中，令G(V，E)代表该基因调控网络，其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i(转录因子i)；它关联的x_i表示该基因转录的mRNA_i浓度；T_i表示与节点i相连的边数(即节点i的度)，I_i表示节点i的入度，O_i表示节点i的出度；x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA_i浓度；n表示基因调控网络中基因的数量；

基于网络的Euclidean距离：

$EE=\sqrt{\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\mathrm{\δ}}_{ij}(x_{1i}-x_{2i})(x_{1j}-x_{2j})}$ 公式(3.11)；

式中，

其中，x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA_i浓度；x_1j指被试者P1的mRNA_j浓度；x_2j指被试者P2的mRNA_j浓度；n表示基因调控网络中基因的数量；

基于网络的Pearson距离：

$EP=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{1i}-\stackrel{\‾}{x_1})}^2}\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{2i}-\stackrel{\‾}{x_2})}^2}}$ 公式(3.12)；

其中，x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA₁浓度；n表示基因调控网络中基因的数量；

I_i表示节点i的入度；这里的x_i1指被试者Pi的mRNA₁浓度；这里的x_i2指被试者Pi的mRNA₂浓度；

所述过程(3)具体是指：将过程(2)中定义的距离引入聚类分析中，使用k-medoids聚类分析方法，对过程(1)所获得的基因表达数据进行聚类；

假设有n个被试者，我们将n个被试者划分为k类，K-medoids聚类算法是，基于网络的Pearson距离具体的算法具体方法如下：

(a)从n个数据对象中任意选取k个数据对象作为medoids-聚类的中心，

(b)选定基于网络的Person距离，即：

$EP=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{1i}-\stackrel{\‾}{x_1})}^2}\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{2i}-\stackrel{\‾}{x_2})}^2}},$

然后分别计算余下的数据对象到各个聚类中心的距离，并将余下的数据对象分配到离自己最近的聚类中，最终得到k组划分，G₁，G₂，...，G_k；

(c)数据对象分配完成后，顺序选取一个数据对象来代替原来的聚类中心，并计算代替后的优化目标函数 $f=\underset{i=1}{\overset{k}{\Σ}}\underset{js.t.x_j\∈G_i}{\Σ}d\left(x_j^{\left(i\right)},x_i^{*}\right);$

其中，d(x₁，x₂)定义如下：

$d(x_1,x_2)=1-\frac{{\mathrm{\Σ}}_{i=1}^n{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{\left(x_{1i}-\stackrel{\‾}{x_1}\right)}^2\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{\left(x_{2i}-\stackrel{\‾}{x_2}\right)}^2}}};$

同理定义d(x_i，x_j)和其中，为从x₁，x₂，...，x_n中选取的k个聚类中心；表示x_j∈G_i；

再选择f最小的数据对象来代替聚类中心，这样K个mediods就改变了；

(d)与前一次的聚类中心相比较，如果发生变化转到方法(b)，如果不发生变化转到方法(e)；

(e)将聚类的结果输出；

所述过程(4)具体是指：根据过程(3)的聚类结果，得出最终关于疾病亚型的分类。

与现有技术相比，本发明的有益效果是：

对于特定的基因网络，基于网络的聚类方法将有更好的组间相似性，更有效地还原三种亚型。此外，当有大量的基因需要测定其表达数据时，现有的方法可能无法同时对所有的基因进行精确的测量。此时，我们提出的“基于网络的聚类”法使得我们通过优先精确测量信息基因的表达数据，并不会大大地削弱对于疾病亚型的鉴定效果。

附图说明

图1为本发明的操作流程图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

现在，我们假定一共有32位被试者P₁、P₂、...、P₃₂，其中被试者P₁、P₂、...、P₈为正常状态病人，被试者P₉、P₁₀、...、P₁₆为患有基因疾病亚型D₁的病人，被试者P₁₇、P₁₈、...、P₂₄患有基因疾病亚型D₂，被试者P₂₅、P₂₆、...、P₃₂患有基因疾病亚型D₃。D₁、D₂、D₃中的每一种亚型都代表着某些基因表达的失常。为了模拟这个表达失常过程，对于某一种亚型，我们从整个基因调控网络中随机选取一定的节点(也就是基因)，对其最大转录速率m_i进行扰动。对于不同的亚型，我们选取不同的基因进行扰动。我们希望做的是通过对32位被试者最后的mRNA浓度向量进行聚类分析，试图分出对照组与三种疾病亚型。

步骤A：我们根据基因调控网络的某些特性(例如网络的平均出度、入度等参数)来构建随机的有向图来代表基因调控网络构建基因调控网络。假设我们要产生由n个基因组成的基因调控网络，根据基因调控网络的特性，我们将产生一张平均入度为2，且分布满足幂定理分布(powerlawdistribution)的随机有向网络G(V，E)，其中|V|＝n。此外，图中不允许有自环的出现。

步骤B：根据我们产生的基因调控网络G(V，E)，对每个基因i建立激活函数f_i(·)。建立过程如下：(1)对于任意的基因i∈V，i＝1，2，K，n，我们从G(V，E)中找出所有与i相邻且以i为有向边终点的点，构成影响因子集合{q₁，q₂，K，q_sn}。(2)确定解离常数k_ij：k_ij将从定义在[0.01，1]区间上的均匀分布中选取。(3)确定希尔系数n_ij：n_ij服从[1，10]区间中的高斯分布函数(4)确定相对活性α_i：α_i将从定义在[0，1]区间上的均匀分布上采样。

步骤C：确定无噪声动态基因调控模型，也就是方程

$\left\{\begin{array}{c}{F}_i^{mRNA}(x,y)=\frac{{\mathrm{dx}}_i}{dt}=m_i\·f_i\left(y\right)-{\mathrm{\λ}}_i^{RNA}\·x_i\\ F_i^{\mathrm{Pr}ot}(x,y)=\frac{{\mathrm{dy}}_i}{dt}=r_i\·x_i-{\mathrm{\λ}}_i^{\mathrm{Pr}ot}\·y_i\end{array}\right.,$ (2.1)的各个参数。

式中：x_i——基因i的浓度；

y_i——蛋白质i的浓度；

——mRNA_i的浓度变化率；

——蛋白质i的浓度变化率；

m_i——基因i的最大转录速率；

r_i——mRNA_i的翻译速率；

——mRNA_i的降解速率；

——蛋白质i的降解速率；

f_i(·)——基因i的激活函数。

首先，mRNA的半衰期以及蛋白质的半衰期(以分钟为单位)将从定义在[5，50]区间上的高斯分布上采样。根据公式：

$\left\{\begin{array}{c}{\mathrm{\λ}}_i^{RNA}=\frac{ln2}{T_i^{RNA}}\\ {\mathrm{\λ}}_i^{\mathrm{Pr}ot}=\frac{1n2}{T_i^{\mathrm{Pr}ot}}\end{array}\right.---\left(2.9\right)$

我们便可以获得mRNA以及蛋白质的降解速率。此外，最大转录速率m_i以及翻译速率r_i服从[0.01，0.011]区间上的均匀分布。

步骤D：在获得了基因调控网络以及动态基因模型之后，求解方程(2.1)，并做出各个mRNA浓度随时间变化的时间序列图。为此，我们首先要选定mRNA浓度x(x₁，x₂，K，x_n)以及蛋白质浓度y(y₁，y₂，K，y_n)的初始值。我们的做法是令各个x_i和y_i服从[0，1]区间上的均匀分布，并随机选取作为初始值。在获得了各个mRNA以及蛋白质浓度的初始值后，我们便可以通过方程(2.1)求解最终的基因表达数据。我们一共取n＝10，20，50个基因，计算t＝0到t＝500之间的时间序列图。

步骤F：我们用一个仿真给出一个观察。我们首先选定一个基因调控网络。接着，我们随机产生四组不同的初值x_i(x_i1，x_i2，K，x_in)，n＝20，i＝1，2，3，4以及y_i(y_i1，y_i2，K，y_in)，n＝20，i＝1，2，3，4作为mRNA及蛋白质初始浓度。最后，我们用方程(2.1)计算并绘制四种不同初值情形下的mRNA浓度时间序列图，并给出四种状态下mRNA浓度终态的数据。

步骤G：我们使用基因调控模型(2.1)进行仿真，并最终得到一个O-G矩阵(32×20)。然后我们运用基于网络的Pearson距离结合k-medoids算法对被试者进行分类具体做法为将32个被试者划分为4类，K-medoids聚类算法是，基于网络的Pearson距离具体的算法过程如下：

A：从32数据对象中任意选取4个数据对象作为medoids-聚类的中心；

B：选定基于网络的Person距离， $EP=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{1i}-\stackrel{\‾}{x_1})}^2}\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{2i}-\stackrel{\‾}{x_2})}^2}}$ 分别计算余下的数据对象到各个聚类中心的距离，并将余下的数据对象分配到离自己最近的聚类中，最终得到4组划分；

C：数据对象分配完成后，顺序选取一个数据对象来代替原来的聚类中心，并计算代替后的优化目标函数f，选择f最小的数据对象来代替聚类中心，这样4个mediods就改变了；

D：与前一次的聚类中心相比较，如果发生变化转到(B)，如果不发生变化转到(E)；

E：将聚类的结果输出，并将分类的结果与真实情况进行对比(真实情况即为正常类型与每种异常各有8位被试者)。

步骤H：综上所述，我们提出“基于网络的聚类(Network-basedClustering)”方法，得出最终关于疾病亚型的分类。

最后，需要注意的是，以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有很多变形。本领域的普通技术人员能从本发明公开的内容中直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (1)

1.用于疾病亚型问题的基于网络的聚类方法，其特征在于，具体包括下述过程：

(1)获得O-G矩阵以及基因调控网络；

(2)选取适用于具体问题的基于网络的距离定义，构建距离矩阵；

(3)运用k-medoids算法对O-G矩阵进行聚类分析；聚类时距离的选择用基于网络的距离；

(4)得出最终关于疾病亚型的分类；

所述过程(1)具体包括下述步骤：

步骤A：根据基因调控网络的特性，构建随机的有向图来代表基因调控网络G(V，E)；其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i；每条有向边e_ji∈E代表着“转录因子j调控基因i的转录”这种调控关系；

步骤B：根据产生的基因调控网络G(V，E)，对每个基因i建立激活函数f_i(·)，具体建立方式为：

对于任意的基因i∈V，i＝1，2，K，n，我们从G(V，E)中找出所有与i相邻且以i为有向边终点的点，构成影响因子集合{q₁，q₂，K，q_sn}；其中，q₁表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，q₂表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，q_sn表示与i相邻且以i为有向边终点的某基因中对基因i起影响作用的因子，n表示基因调控网络中基因的数量；

确定解离常数k_ij，且k_ij从定义在[0.01，1]区间上的均匀分布中选取；

确定希尔系数n_ij，且n_ij服从[1，10]区间中的高斯分布函数

确定相对活性α_i，且α_i从定义在[0，1]区间上的均匀分布上采样；

步骤C：确定无噪声动态基因调控模型，即确定公式(2.1)的各个参数；

$\left\{\begin{array}{c}{F}_i^{mRNA}(x,y)=\frac{{\mathrm{dx}}_i}{dt}=m_i\·f_i\left(y\right)-{\mathrm{\λ}}_i^{RNA}\·x_i\\ F_i^{\mathrm{Pr}ot}(x,y)=\frac{{\mathrm{dy}}_i}{dt}=r_i\·x_i-{\mathrm{\λ}}_i^{\mathrm{Pr}ot}\·y_i\end{array}\right.---\left(2.1\right)$

式(2.1)中，x_i表示基因i的浓度；y_i表示蛋白质i的浓度；F_i ^mRNA表示mRNA_i的浓度变化率；F_i ^Prot表示蛋白质i的浓度变化率；m_i表示基因i的最大转录速率；r_i表示mRNA_i的翻译速率；表示mRNA_i的降解速率；表示蛋白质i的降解速率；f_i(·)表示基因i的激活函数；

确定公式(2.1)中各个参数的具体方式为：mRNA的半衰期T_i ^RNA以及蛋白质的半衰期T_i ^Prot从定义在[5，50]区间上的高斯分布上采样；

根据公式(2.9)，获得mRNA以及蛋白质的降解速率，最大转录速率m_i以及翻译速率r_i服从[0.01，0.011]区间上的均匀分布；

$\left\{\begin{array}{c}{{\mathrm{\λ}}_i}^{\mathrm{RNA}}=\frac{\ln 2}{{T_i}^{\mathrm{RNA}}}\\ {{\mathrm{\λ}}_i}^{\mathrm{Rrot}}=\frac{\ln 2}{{T_i}^{\mathrm{Rrot}}}\end{array}\right.---\left(2.9\right)$

式(2.9)中，表示mRNA_i的降解速率；表示蛋白质i的降解速率；mRNA的半衰期T_i ^RNA以及蛋白质的半衰期T_i ^Prot；

步骤D：在获得了基因调控网络以及无噪声动态基因调控模型之后，选定mRNA浓度x(x₁，x₂，K，x_n)以及蛋白质浓度y(y₁，y₂，K，y_n)的初始值，然后求解公式(2.1)，得到最终的基因表达数据；

所述过程(2)具体是指：根据过程(1)所获得的基因网络的拓扑关系G(V，E)，定义三种基于网络的距离，用于比较x₁(x₁₁，x₁₂，K，x_1n)与x₂(x₂₁，x₂₂，K，x_2n)的差别；其中x₁(x₁₁，x₁₂，K，x_1n)、x₂(x₂₁，x₂₂，K，x_2n)分别表示两个被试者P₁和P₂的mRNA浓度；

令G(V，E)代表该基因调控网络，其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i；它关联的x_i表示该基因转录的mRNA_i浓度；令每条有向边e_ji∈E代表着“转录因子j调控基因i的转录”这种调控关系；记T_i表示与节点i相连的边数，I_i表示节点i的入度，O_i表示节点i的出度；

其中，基于网络的Jaccard距离定义为：

$EJ=\frac{2\underset{i=1}{\overset{n}{\Σ}}{T}_i\left|x_{1i}-x_{2i}\right|}{\underset{i=1}{\overset{n}{\Σ}}{T}_i\left(x_{1i}+x_{2i}\right)+\underset{i=1}{\overset{n}{\Σ}}{T}_i\left|x_{1i}-x_{2i}\right|}$ 公式(3.10)；

其中，令G(V，E)代表该基因调控网络，其中每个顶点i∈V代表基因i及其产生的mRNA_i和蛋白质i；它关联的x_i表示该基因转录的mRNA_i浓度；T_i表示与节点i相连的边数，I_i表示节点i的入度，O_i表示节点i的出度；x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA_i浓度；n表示基因调控网络中基因的数量；

基于网络的Euclidean距离：

$EE=\sqrt{\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\mathrm{\δ}}_{ij}(x_{1i}-x_{2i})(x_{1j}-x_{2j})}$ 公式(3.11)；

式中，

其中，x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA_i浓度；x_1j指被试者P1的mRNA_j浓度；x_2j指被试者P2的mRNA_j浓度；n表示基因调控网络中基因的数量；

基于网络的Pearson距离：

$EP=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{1i}-\stackrel{\‾}{x_1})}^2}\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{2i}-\stackrel{\‾}{x_2})}^2}}$ 公式(3.12)；

其中，x_1i指被试者P1的mRNA_i浓度；x_2i指被试者P2的mRNA₁浓度；n表示基因调控网络中基因的数量；

I_i表示节点i的入度；这里的x_i1指被试者Pi的mRNA₁浓度；这里的x_i2指被试者Pi的mRNA₂浓度；

所述过程(3)具体是指：将过程(2)中定义的距离引入聚类分析中，使用k-medoids聚类分析方法，对过程(1)所获得的基因表达数据进行聚类；

假设有n个被试者，我们将n个被试者划分为k类，K-medoids聚类算法是，基于网络的Pearson距离具体的算法具体方法如下：

(a)从n个数据对象中任意选取k个数据对象作为medoids-聚类的中心，

(b)选定基于网络的Person距离，即：

$EP=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{1i}-\stackrel{\‾}{x_1})}^2}\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{(x_{2i}-\stackrel{\‾}{x_2})}^2}},$

然后分别计算余下的数据对象到各个聚类中心的距离，并将余下的数据对象分配到离自己最近的聚类中，最终得到k组划分，G₁，G₂，…，G_k；

(c)数据对象分配完成后，顺序选取一个数据对象来代替原来的聚类中心，并计算代替后的优化目标函数 $f=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\underset{\mathrm{j\; s}.t.x_j\∈G_i}{\mathrm{\Σ}}d(x_j^{\left(i\right)},x_i^{*});$

其中，d(x₁，x₂)定义如下：

$d(x_1,x_2)=1-\frac{{\mathrm{\Σ}}_{i=1}^n{I}_i(x_{1i}-\stackrel{\‾}{x_1})(x_{2i}-\stackrel{\‾}{x_2})}{\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{\left(x_{1i}-\stackrel{\‾}{x_1}\right)}^2\sqrt{{\mathrm{\Σ}}_{i=1}^n{I}_i{\left(x_{2i}-\stackrel{\‾}{x_2}\right)}^2}}};$

同理定义d(x_i，x_j)和其中，为从x₁，x₂，…，x_n中选取的k个聚类中心；表示x_j∈G_i；

再选择f最小的数据对象来代替聚类中心，这样K个mediods就改变了；

(d)与前一次的聚类中心相比较，如果发生变化转到方法(b)，如果不发生变化转到方法(e)；

(e)将聚类的结果输出；

所述过程(4)具体是指：根据过程(3)的聚类结果，得出最终关于疾病亚型的分类。