CN109360607A - 一种动态基因调控网的网络演化分析方法及装置 - Google Patents

Abstract

本发明提供了一种动态基因调控网的网络演化分析方法及装置，其中的方法对动态基因调控网网络演化的研究分成两个步骤来进行，首先对基因调控网中模体的转换进行了研究，将现有对基因调控网的研究扩展到了动态的领域，能够更好的捕捉到基因调控网的演化模式，其次通过基于隐空间特征的符号判别算法对将来时刻基因调控网快照的连边符号进行了判别，将研究领域扩展到带符号网络。实现了对动态基因调控网的演化分析，改善了预测效果。

Description

一种动态基因调控网的网络演化分析方法及装置

技术领域

本发明涉及计算机科学和生物网络技术领域，具体涉及一种动态基因调控网的网络演化分析方法及装置。

背景技术

在生物体内，基因通过表达和调控相互作用，从而实现它们的生物学功能，并完成复杂的生命活动。基因间的相互作用是一个连续而复杂的动态过程，随着时间、环境的变化而变化。对基因调控网网络演化的研究有许多重要意义，例如可以通过当前的基因调控关系预测未来的基因调控关系，还可以通过基因调控关系的变化预测基因功能的变化，从而探明某些疾病尤其是癌症的发病机制、为疾病的预测和治疗提供依据等。

目前，在基因调控网网络演化的研究中，生物医学领域的研究者和计算机科学的研究者分别从不同的角度对基因调控网的网络演化进行了分析。网络演化与链路预测具有内在的一致性，通过链路预测算法可以对网络演化模型的正确性进行验证。传统的链路预测方法主要分为三类：基于相似性的链路预测、基于最大似然估计的链路预测以及概率模型的链路预测方法。其中，基于相似性的链路预测算法包括共同邻居算法(CN)、AA算法等。基于最大似然估计的链路预测方法通过似然估计值和马尔可夫-蒙特卡洛算法得到两节点之间产生连边的概率，最大似然估计方法在面对有明显层次结构的复杂网络时有较好的效果，但面对非层次网络时效果不如基于相似性的方法。概率模型方法的基本思想是建立一个具有多参数的概率模型，通过调节参数，使模型能够再现该网络的真实连边关系。基于这类思想的经典算法有马尔可夫网络模型(RMN)、朴素贝叶斯等。

本发明申请人在实施本发明的过程中，发现现有方法中，主要存在下述技术问题：基于概率模型的预测方法因计算复杂度过高无法应用于大规模网络，而基于相似度的链路预测方法和基于似然估计的链路预测方法，都是根据网络的某些局部或全局的某些拓扑结构特征来进行预测，而无法对基因调控网进行很好地预测。

由上可知，现有技术的方法对基因调控网的预测效果不佳的技术问题。

发明内容

有鉴于此，本发明提供了一种动态基因调控网的网络演化分析方法及装置，用以解决或者至少部分解决现有技术的方法对基因调控网的预测效果不佳的技术问题。

本发明第一方面提供了一种动态基因调控网的网络演化分析方法，包括：

步骤S1：将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况；

步骤S2：以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻；

步骤S3：对源网络和目标网络的连边分别提取显式特征和隐式特征，其中，源网络为符号已知的网络，目标网络为基因调控网，并基于隐式特征，将无符号的基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

在一种实现方式中，步骤S1具体包括：

步骤S1.1：将网络中所有节点每三个组成一组，进行分组和编号，每组为一个模体，根据三个节点间连边的不同，共有64种基本模体类型；

步骤S1.2：将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM，其中，矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示从t时刻到t+1时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体；

步骤S1.3：将不同时刻的模体转换概率矩阵合成模体转换概率张量(TCT)，其中，TCT＝(TCM₁,TCM₂,…,TCM_T-1)，张量中的元素TCT(i,j,t)＝TCM_t(i,j)。

在一种实现方式中，步骤S2具体包括：

步骤S2.1：对模体转换概率张量进行非负张量分解，分解得到三个因式矩阵A,B,C，其中AB^T表示不同类型模体间的转换关系，C包含上述转换关系在时间维度上的信息，称为时间因式矩阵；

步骤S2.2：采用指数分布模型预测时间因式矩阵的第T行，计算方式如下：

其中，T为常量，表示快照的总数量，t为变量，表示从1到T中的某一个时刻，a表示权重参数，r表示矩阵A、B、C的一列，若共有R列，则r取值为1到R之间的任意整数；

步骤S2.3：根据预测得到的T时刻的时间因式矩阵获得T时刻的模体转换概率矩阵，称为模体转换似然矩阵(TCLM)，计算方式如下：

其中，λ_r为非负张量分解的结果，为已知数，A(i,r)为矩阵A的第i行第r列元素，B(j,r)为矩阵B的第j行第r列元素,C(T,r)为矩阵C的第T行第r列元素。

在一种实现方式中，步骤S3具体包括：

步骤S3.1：对连边进行特征提取并映射到隐空间中，提取的特征包括显式特征和隐式特征，对一个有向边(u,v)，显式特征包括节点的度数、中介中心性，模体个数以及共同邻居，对一个有向边(u,v)，通过非负矩阵三因式分解的方法提取其隐式特征，具体的提取方式如下：

其中，||·||_F为弗罗贝尼乌斯范数，M为基因调控网的规模，U_s,V_s,U_t,V_t是提取出的四个隐空间特征矩阵，U_s,V_s,U_t,V_t为提取出的四个隐空间特征矩阵，U_s的第d行代表源网络第d个节点作为边的出节点的特征向量，V_s的第d行代表源网络

d个节点作为边的入节点的特征向量，U_t第d行代表目标网络第d个节点作为边的出节点的特征向量,V_t的第d行代表目标网络d个节点作为边的入节点的特征向量，α为正则化系数，Σ_k为正则项参数，目标函数的作用为寻找合适的矩阵分解，使A_s≈U_s∑_kV_s ^T且其中，矩阵是源网络和目标网络共有的特征空间；

步骤S3.2：对源网络的显式特征和隐式特征进行训练，得到分类模型；

步骤S3.3：将分类模型应用于目标网络的连边中，获得目标符号网络的边的符号，获得未来时刻的有符号的网络快照。

在一种实现方式中，步骤S1.1具体包括：

将每组模体中三个节点的连边状况与预先存储的基本模体的连边状况进行比对，确定所属的模体类型。

在一种实现方式中，步骤S2.1具体包括：

步骤2.1.1：通过求解目标函数的方式，进行非负张量分解，其中，目标目标函数的形式如下：

s.t.λ_n≥0,a_n≥0,b_n≥0,c_n≥0,n＝1,2,…,r

步骤2.1.2：采用MATLAB工具中预设函数执行上述目标函数，得到分解结果。

基于同样的发明构思，本发明第二方面提供了一种动态基因调控网的网络演化分析装置，包括：

预处理模块，用于将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况；

模体转换模块，用于以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻；

符号判别模块，用于基于隐空间特征，将无符号的基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，并以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

基于同样的发明构思，本发明第三方面提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被执行时实现第一方面的所述的方法。

本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：

在本发明提供的方法，一方面，将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，将基因调控网的演化分析转化为对模体的演化分析，将现有对基因调控网的研究扩展到了动态的领域，能够更好的捕捉到基因调控网的演化模式，并且，通过前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，可以得到较好的预测效果，另一方面，通过基于隐空间特征的符号判别方法对将来时刻基因调控网快照的连边符号进行了判别，将研究领域扩展到带符号网络，可以进一步改善预测效果，解决现有技术的方法对基因调控网的预测效果不佳的技术问题。本发明的方法中基因调控网网络演化的研究，可以为生物学实验节约成本，并为基因组学、药物研发的相关科研人员提供可靠的参考，促进生物医学事业的进步。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例中动态基因调控网的网络演化分析方法的流程图；

图2为图1中所示方法的整体框架图；

图3为一种实施方式中基于模体的基因调控网建模过程示意图；

图4为基于隐空间特征的符号判别过程示意图；

图5为本发明实施例中动态基因调控网的网络演化分析装置的结构框图；

图6为本发明实施例中计算机可读存储介质的结构图。

具体实施方式

本发明实施例提供了一种动态基因调控网的网络演化分析方法及装置，通过基于模体转换概率和隐空间特征的符号判别方法，将以往对基因调控网的研究扩展到了动态的领域，能够更好的捕捉到基因调控网的演化模式。其次通过基于隐空间特征的符号判别算法对将来时刻基因调控网快照的连边符号进行了判别，将研究领域扩展到带符号网络。解决了现有技术的方法对基因调控网的预测效果不佳的技术问题，改善了基因调控网的预测效果。

为了达到上述技术效果，本发明的总体思路如下：

首先将基因调控网映射为模体，对相邻快照间的模体变化进行统计分析，得到未来时刻的无符号网络快照；然后基于隐空间特征的符号判别方法，对源网络和目标网络的连边分别提取显式特征和隐式特征，提取到的特征值构成特征矩阵。在此基础上，以连边符号作为标签，用机器学习的方法训练得到分类模型，并应用于目标网络上，将连边符号进行分类，分为激励边(negative edge)和抑制边(positive edge)。得到未来时刻的带符号的网络快照。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供了一种动态基因调控网的网络演化分析方法，请参见图1，该方法包括：

首先执行步骤S1：将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况。

具体来说，基因调控网的拓扑结构随着时间和外部条件的变化不断变化。以一定的时间间隔对其进行采样，每一次采样得到一个基因调控网在某时刻的静态结构，这个静态结构称之为快照，G_T表示T时刻网络G的快照。一串以一定时间间隔分割的快照序列G＝(G₁,G₂,…,G_t)就是基因调控网，它反映了基因调控网动态演化的过程。时间间隔可以根据实际情况进行设置，例如设置为100秒、200秒等，当时间间隔为100秒时，则记录第0秒、第100秒、第200秒的快照，相邻是指时间上相邻，例如快照1、快照2则是第0秒、第100秒的快照。

然后执行步骤S2：以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻。

具体来说，t时刻的模体转换概率矩阵即为预测结果，在具体的实现过程中，对T+1时刻快照中的节点对(u,v)，对其赋予一个分值，其分值越大，说明这对节点之间存在连边的可能性越高，此时预测结果就是得到T+1时刻快照中节点对的分数。计算方法如下：

其中|tri(u,v)|表示T-1时刻包含节点对(u,v)的模体的总数，T(m)表示包含了边(u,v)的第m个模体从T时刻到T+1时刻的转移概率，W_m是表示第m个模体的影响力的分值。

步骤S3：对源网络和目标网络的连边分别提取显式特征和隐式特征，其中，源网络为符号已知的网络，目标网络为基因调控网，并基于隐式特征，将无符号的动态基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

具体来说，连边符号未知的网络(本实施方式中为基因调控网)称作目标网络，将符号已知的网络称作源网络。显式特征用以表达实例中显而易见的属性。隐式特征不能直接由网络拓扑结构看出，但可以表达源网络和目标网络之间所共有的模式。由于基因调控网的网络演化规律十分复杂，仅使用显式特征无法很好的对基因调控网的符号进行准确的判别，为了更好的利用源网络中已知符号的边所蕴含的信息，本实施还构造了隐空间特征以捕捉蕴藏在拓扑结构之下的源网络和目标网络之间共有的模式，从而进一步改善预测效果。

请参见图2，为动态基因调控演化分析的整体框架图，首先将基因调控网快照进行预处理，即进行模体转换概率统计以及隐空间映射，然后基于模体转换概率的连边预测算法进行网络预测，输出无符号的未来时刻基因调控网快照，接着基于隐空间映射和无符号的未来时刻基因调控网快照进行符号判别，最终获得未来时刻动态基因调控网快照。

基因调控网是由基因表达数据经过推断生成的用来描述基因间调控关系的带符号的有向图。基因表达数据是一个l×m的矩阵，它是表示l个基因在m个不同时刻上的表达强度高低的数据，这一数据是通过生物学实验获得的。根据基因表达数据可以计算基因与基因之间表达强度变化的相关性。基因调控关系分为两类，如果一个基因的表达强度的提高导致另一个基因表达强度提高，称这种调控关系为激励关系；如果一个基因表达强度的提高导致另一个基因的表达强度降低，就称这种调控关系为抑制关系。将这种基因调控关系映射为复杂网络中节点和边的关系，就是基因调控网的基本形式，动态基因调控网快照反映了基因调控网动态演化的过程。

在一种实施方式中，步骤S1具体包括：

步骤S1.1：将网络中所有节点每三个组成一组，进行分组和编号，每组为一个模体，根据三个节点间连边的不同，共有64种基本模体类型；

具体来说，模体是网络中的一种特殊的子图结构，在基因调控网中具有重要的统计学和生物学意义。模体这类子图在网络所有子图中占较大比例，在基因调控网中，某些模体已被确定具有生物学意义。模体种类和结构复杂多样，本发明中研究的模体限于由三个节点构成的子图。在有向网络中，三个节点按照两两之间的连边状况，总共存在64种可能的连边情况。

步骤S1.2：将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM，其中，矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示从t时刻到t+1时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体；

步骤S1.3：将不同时刻的模体转换概率矩阵合成模体转换概率张量(TCT)，其中，TCT＝(TCM₁,TCM₂,…,TCM_T-1)，张量中的元素TCT(i,j,t)＝TCM_t(i,j)。

在一种实施方式中，步骤S2具体包括：

步骤S2.1：对模体转换概率张量进行非负张量分解，分解得到三个因式矩阵A,B,C，其中AB^T表示不同类型模体间的转换关系，C包含上述转换关系在时间维度上的信息，称为时间因式矩阵；

步骤S2.2：采用指数分布模型预测时间因式矩阵的第T行，计算方式如下：

其中，T为常量，表示快照的总数量，t为变量，表示从1到T中的某一个时刻，a表示权重参数，r表示矩阵A、B、C的一列，若共有R列，则r取值为1到R之间的任意整数；

步骤S2.3：根据预测得到的T时刻的时间因式矩阵获得T时刻的模体转换概率矩阵，称为模体转换似然矩阵(TCLM)，计算方式如下：

其中，i,j为变量指针，指模体转换似然矩阵中的每一个元素，取值为1到64的整数。T意义同上，R意义同上，为A、B、C的总列数。λ_r是非负张量分解的结果之一，为已知数。A(i,r)为矩阵A的第i行第r列元素，B(j,r)为矩阵B的第j行第r列元素,C(T,r)为矩阵C的第T行第r列元素。

具体来说，步骤S2.2可以预测第T行，即T时刻的时间因式矩阵，而前T-1行是由非负张量分解已知的，预测得到的第T行，加上前T-1行，共同构成时间因式矩阵。

在一种实施方式中，步骤S3具体包括：

步骤S3.1：对连边进行特征提取并映射到隐空间中，提取的特征包括显式特征和隐式特征，对一个有向边(u,v)，显式特征包括节点的度数、中介中心性，模体个数以及共同邻居，对一个有向边(u,v)，通过非负矩阵三因式分解的方法提取其隐式特征，具体的提取方式如下：

其中，||·||_F为弗罗贝尼乌斯范数，M为基因调控网的规模，U_s,V_s,U_t,V_t是提取出的四个隐空间特征矩阵，U_s,V_s,U_t,V_t为提取出的四个隐空间特征矩阵，U_s的第d行代表源网络第d个节点作为边的出节点的特征向量，V_s的第d行代表源网络d个节点作为边的入节点的特征向量，U_t第d行代表目标网络第d个节点作为边的出节点的特征向量,V_t的第d行代表目标网络d个节点作为边的入节点的特征向量，α为正则化系数，Σ_k为正则项参数，目标函数的作用为寻找合适的矩阵分解，使A_s≈U_s∑_kV_s ^T且其中，矩阵是源网络和目标网络共有的特征空间；

步骤S3.2：对源网络的显式特征和隐式特征进行训练，得到分类模型；

步骤S3.3：将分类模型应用于目标网络的连边中，获得目标符号网络的边的符号，获得未来时刻的有符号的网络快照。

在一种实施方式中，步骤S1.1具体包括：

将每组模体中三个节点的连边状况与预先存储的基本模体的连边状况进行比对，确定所属的模体类型。

具体来说，将所有节点进行三个一组的排列组合，为了避免重复，可以设定每组序号小的节点永远在序号大的节点前面。

在一种实施方式中，步骤S2.1具体包括：

步骤2.1.1：通过求解目标函数的方式，进行非负张量分解，其中，目标目标函数的形式如下：

s.t.λ_n≥0,a_n≥0,b_n≥0,c_n≥0,n＝1,2,…,r

步骤2.1.2：采用MATLAB工具中预设函数执行上述目标函数，得到分解结果。

具体来说，是将三阶张量TCT∈R^I×J×K进行分解，非负张量分解的目的就是寻找向量a_i∈R^I,b_i∈R^J,c_i∈R^K,n＝1,2,…r，使得：

s.t.λ_n≥0,a_n≥0,b_n≥0,c_n≥0,n＝1,2,…,r。

预设函数可以是tensor toolbox工具箱中的cp_emu()函数，通过MATLAB则可以获得分解结果。

为了更清楚地说明本发明实施例提供的方法，下面通过具体的示例对其中的基因调控网以及符号判别过程，予以说明，请参见图3和图4，分为为基于模体的基因调控网建模过程示意图和基于隐空间特征的符号判别过程示意图。

在图3中，表示了快照T、快照T+1、快照T+N等，以及模体间的转换概率，三个节点按照两两之间的连边状况，总共存在64种可能的连边情况。在网络中，任意三个节点都可以映射成一个模体。在一种实例中，若单个基因调控网快照的规模为|V|＝M，则此快照中一共有M×(M-1)×(M-2)个模体。如上节所描述的，网络中的每一个模体都有且仅有一个代表其类型的编号。网络演化的过程可以看作是模体转换的过程，从微观来看，一个模体在下一时刻要么保持原有类型不变，要么转换成另一种类型的模体。通过上述过程，本发明统计一个动态基因调控网络中所有快照中所有模体的转换，进而统计64类模体间相互转换的概率，从而描述网络的演化过程。模体i转换为模体j的概率P(i→j)就是这个网络的一个演化特征。相邻快照间对应的模体是可以由一种类型转换到另一种类型的。由于受到不同生命周期的发展特征、生长环境以及各种内外部因素的影响，模体之间相互转换的概率随时间变化是非线性关系，不能简单地采用线性模型来估计未来时刻的模体转换概率。因此，本发明使用一种基于张量分解的方法进行时间序列预测。张量作为多维数组，是向量和矩阵在高维空间上的推广。用它来描述数据，相比传统的矩阵方法，能够最大程度的保证高维数据原本的结构上不被分解破坏，从而保证预测结果的准确性。

具体来说，本发明实施方式将求解非负张量分解的问题形式表达如下：

s.t.λ_n≥0,a_n≥0,b_n≥0,c_n≥0,n＝1,2,…,r

将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM。矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示t时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体。

然后构建一个模体转换概率张量(TCT)来表示模体转换的时间序列，若一个动态基因调控网中有T个快照，则TCT＝(TCM₁,TCM₂,…,TCM_T-1),张量中的元素TCT(i,j,T)＝TCM_T(i,j)。通过MATLAB的tensor toolbox工具箱中的cp_nmu()函数对TCT张量进行非负CP分解，得到参数λ和三个因式矩阵A,B,C，其中AB^T表达了不同类型模体间的转换关系，而C则包含了该关系在时间维度上的信息，称为时间因式矩阵。使用指数分布模型对时间因式矩阵C的前T-1行进行时间序列预测，得到矩阵C的第T行元素，从而预测得到T-1时刻到T时刻的一种可能的模体转换概率矩阵，称为模体转换似然矩阵(TCLM)，如下所示：

接着，将TCLM按行进行归一化，使矩阵每行之和为1，便得到了T到T+1时刻的模体转换概率矩阵TCM_T。获得了TCM_T，便可以对T+1时刻的链路状况进行预测。连边预测的目的就是给未来时刻的基因调控网快照的每一个节点对(u,v)赋予一个分数score(u,v)。该分数越高意味着该节点对之间存在边的可能性越大。作为一种实现方式，由于本发明所研究的基因调控网是有向网络，所以对一个节点对的两种可能的边的方向——即score(u,v)和score(v,u)分别赋分。从上述描述可知，所有包含了u和v两节点的模体的转换都可以为下一时刻该节点对之间存在边的可能性产生影响，然而不同模体的影响力不是一样的。越有影响力的模体在连边预测中所占的比重越大。为了更为准确地进行预测，本发明将模体的“影响力”定义为两个方面：一是历史快照内该模体中的连边形成频率，二是历史快照内该模体形成闭合的频率。总体来说，该模体越稠密，则说明其内部节点关系越紧密，在链路预测中发挥的作用相比其他稀疏模体更重要。此外，对节点对之间形成连边的概率，即某条历史连边产生的时刻距离待预测时刻越近，则待预测时刻这条边仍然存在的可能性越高。

综上所述，本发明对模体在链路预测中的影响力定义如下所示：

Wi是表示第i个模体的影响力的分值，g(i,t)表示t时刻第i个模体中各节点连边个数，f(i,t)表示t时刻第i个模体是否闭合，即三个节点两两之间是否至少存在一条有向边，闭合时为1，不闭合为0。θ₁,θ₂是控制不同时期的历史快照对Wi贡献的系数，距当前时刻越近，贡献越大。

由此，根据T到T+1时刻的模体转换概率矩阵TCM(T)和模体的影响力W，可以为每一对节点对赋予一个得分，表示T时刻每个边存在的概率高低，如下所示：

其中|tri(u,v)|表示T-1时刻包含节点对(u,v)的模体的总数，T(m)表示包含了边(u,v)的第m个模体从T时刻到T+1时刻的转移概率。由于本发明研究的对象是有向的动态基因调控网，因而不仅要计算有向边(u,v)对应的分数，还需要计算相反方向的边(u,v)对应的分数。

接下来介绍隐空间和符号判别，图4中实心表示抑制边，空心表示激励边。

与其他机器学习方法不同的是，在基因调控网的符号判别中，由于没有“先验”的特征向量可以来对训练集中一条边的符号来进行描述。因此，需要根据源网络和目标网络的拓扑结构进行特征空间的构造。

本发明构造的特征分为两类：(1)显式特征，用以表达实例中显而易见的属性。(2)隐空间特征，不能直接由网络拓扑结构看出，但也表达了源网络和目标网络之间所共有的一些模式。

对一个有向边(u,v),本发明为其定义的显式特征包括节点的度数、中介中心性，模体个数以及共同邻居等。需要注意的是，在为每一个样本定义这些特征的时候，不考虑这条边的符号，因为在目标网络中，对于绝大部分边的符号都是未知的。各个特征的描述具体如下：

(1)节点的度数。对一个有向边(u,v)，通过deg_out(u)和deg_in(v)来分别指代节点u的出度和节点v的入度。节点的度数代表着它与图中其他节点连接的紧密性。

(2)中介中心性。对于一个节点v∈V，其中介中心性为 其中σ_r,s是指节点r到节点s的最短路径的数量，σ_r,s(v)是指节点r到节点s最短路径中经过节点v的数量。对一个有向边(u,v)，采用f_bc(u)和f_bc(v)作为它的两个特征。中介中心性代表着一个点在图中作为中心节点的地位。

(3)模体个数。对一个有向边(u,v)，考虑将包含了(u,v)的模体个数作为其特征。设该模体的第三个节点为w。若该模体中存在有向边(u,w)，则称此边为前向边(F)，若存在有向边(w,v)则称此边为后向边(B)，或者w和u之间没有边存在，则记为N。w和另一个节点v的关系同理。这样对于一条边(u,v)有八个特征，分别为f_FF,f_FB,f_BF,f_BB,f_FN,f_NF,f_BN,f_NB，f_FF表示存在有向边(u,w)和(v,w)的模体的个数。

(4)共同邻居。对一个有向边(u,v)，如果存在另一个节点w，使得w与u和v之间均有边相连。则w为u,v的共同邻居。f_cn(u,v)指有向边(u,v)的两端点的共同邻居的个数。

以上的显式特征是直观的，但基因调控网的网络演化规律十分复杂，仅使用以上特征无法很好的对基因调控网的符号进行准确的判别。为了更好地利用源网络中已知符号的边所蕴含的信息，本发明还构造了隐空间特征以捕捉蕴藏在拓扑结构之下的源网络和目标网络之间共有的模式。

本发明将寻找隐特征空间的问题表示如下：

由于上式所有变量都非负，在求最小值的过程中，∑_k中过大的值将会使U_s,V_s,U_t以及V_t中的某些值趋于0，这会使网络中每个节点的隐空间特征向量难以区分。因此，需要加上一个正则项参数Σ_k。

本发明使用一种迭代更新的算法来求解上式。首先将上式改写成如下便于用代码描述的形式：

其中，tr(·)是指矩阵的迹。

以U_s为例介绍求解上式的方法，由于约束条件中包含U_s≥0，可以使用拉格朗日乘子法解决此问题。本发明引入拉格朗日乘子并使拉格朗日函数最小。设与KKT条件联立可得

基于上式和已知的非负矩阵分解的方法，本发明按以下所示规则迭代更新U_s：

同理可得,V_s,U_t,V_t和Σ_k的迭代规则如下所示：

通过以上算法，可以在若干次迭代后获得使取得最小值的U_s,V_s,U_t,V_t。上述四个矩阵就是本发明要求得的基因调控网络的隐空间特征。对于训练集和测试集的每条边，将两端点的显式特征和隐空间特征向量作为特征、将边的符号作为标签，如图4的过程所示，最后使用LibSVM工具进行训练和预测，得到符号判别的结果，即带符号的目标网络快照。

本发明提供的方法具有如下优点或者有益技术效果：

对动态基因调控网网络演化的研究分成两个步骤来进行，首先对网络中模体的转换进行了研究，将以往对基因调控网的研究扩展到了动态的领域，能够更好的捕捉到基因调控网的演化模式。其次通过基于隐空间特征的符号判别算法对将来时刻基因调控网快照的连边符号进行了判别，将研究领域扩展到带符号网络。对动态基因调控网网络演化的研究，可以为生物学实验节约成本，并为基因组学、药物研发的相关科研人员提供可靠的参考，促进生物医学事业的进步。

基于同一发明构思，本申请还提供了与实施例一中动态基因调控网的网络演化分析方法对应的装置，详见实施例二。

实施例二

本实施例提供动态基因调控网的网络演化分析装置，请参见图5，该装置包括：

预处理模块501，用于将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况；

模体转换模块502，用于以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻；

符号判别模块503，用于基于隐空间特征，将无符号的基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，并以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

在一种实施方式中，预处理模块具体用于执行下述步骤：

步骤S1.1：将网络中所有节点每三个组成一组，进行分组和编号，每组为一个模体，根据三个节点间连边的不同，共有64种可能的模体类型；

步骤S1.2：将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM，其中，矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示从t时刻到t+1时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体；

步骤S1.3：将不同时刻的模体转换概率矩阵合成模体转换概率张量(TCT)，其中，TCT＝(TCM₁,TCM₂,…,TCM_T-1)，张量中的元素TCT(i,j,t)＝TCM_t(i,j)。

由于本发明实施例二所介绍的装置，为实施本发明实施例一中动态基因调控网的网络演化分析方法所采用的装置，故而基于本发明实施例一所介绍的方法，本领域所属人员能够了解该装置的具体结构及变形，故而在此不再赘述。凡是本发明实施例一的方法所采用的装置都属于本发明所欲保护的范围。

施例三

基于同一发明构思，本申请还提供了一种计算机可读存储介质600，请参见图6，其上存储有计算机程序611，该程序被执行时实现实施例一中的方法。

由于本发明实施例三所介绍的计算机可读存储介质，为实施本发明实施例一中动态基因调控网的网络演化分析方法所采用的计算机可读存储介质，故而基于本发明实施例一所介绍的方法，本领域所属人员能够了解该计算机可读存储介质的具体结构及变形，故而在此不再赘述。凡是本发明实施例一的方法所采用的计算机可读存储介质都属于本发明所欲保护的范围。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

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

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

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

Claims (10)

1.一种动态基因调控网的网络演化分析方法，其特征在于，包括：

步骤S1：将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况；

步骤S2：以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻；

步骤S3：对源网络和目标网络的连边分别提取显式特征和隐式特征，其中，源网络为符号已知的网络，目标网络为基因调控网，并基于隐式特征，将无符号的基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

2.如权利要求1所述的方法，其特征在于，步骤S1具体包括：

步骤S1.1：将网络中所有节点每三个组成一组，进行分组和编号，每组为一个模体，根据三个节点间连边的不同，共有64种基本模体类型；

步骤S1.2：将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM，其中，矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示从t时刻到t+1时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体；

步骤S1.3：将不同时刻的模体转换概率矩阵合成模体转换概率张量(TCT)，其中，TCT＝(TCM₁,TCM₂,…,TCM_T-1)，张量中的元素TCT(i,j,t)＝TCM_t(i,j)。

3.如权利要求1所述的方法，其特征在于，步骤S2具体包括：

步骤S2.1：对模体转换概率张量进行非负张量分解，分解得到三个因式矩阵A,B,C，其中AB^T表示不同类型模体间的转换关系，C包含上述转换关系在时间维度上的信息，称为时间因式矩阵；

步骤S2.2：采用指数分布模型预测时间因式矩阵的第T行，计算方式如下：

其中，T为常量，表示快照的总数量，t为变量，表示从1到T中的某一个时刻，a表示权重参数，r表示矩阵A、B、C的一列，若共有R列，则r取值为1到R之间的任意整数；

步骤S2.3：根据预测得到的T时刻的时间因式矩阵获得T时刻的模体转换概率矩阵，称为模体转换似然矩阵(TCLM)，计算方式如下：

其中，λ_r为非负张量分解的结果，为已知数，A(i,r)为矩阵A的第i行第r列元素，B(j,r)为矩阵B的第j行第r列元素,C(T,r)为矩阵C的第T行第r列元素。

4.如权利要求1所述的方法，其特征在于，步骤S3具体包括：

步骤S3.1：对连边进行特征提取并映射到隐空间中，提取的特征包括显式特征和隐式特征，对一个有向边(u,v)，显式特征包括节点的度数、中介中心性，模体个数以及共同邻居，对一个有向边(u,v)，通过非负矩阵三因式分解的方法提取其隐式特征，具体的提取方式如下：

其中，||·||_F为弗罗贝尼乌斯范数，M为基因调控网的规模，U_s,V_s,U_t,V_t是提取出的四个隐空间特征矩阵，U_s,V_s,U_t,V_t为提取出的四个隐空间特征矩阵，U_s的第d行代表源网络第d个节点作为边的出节点的特征向量，V_s的第d行代表源网络d个节点作为边的入节点的特征向量，U_t第d行代表目标网络第d个节点作为边的出节点的特征向量,V_t的第d行代表目标网络d个节点作为边的入节点的特征向量，α为正则化系数，Σ_k为正则项参数，目标函数的作用为寻找合适的矩阵分解，使A_s≈U_s∑_kV_s ^T且其中，矩阵是源网络和目标网络共有的特征空间；

步骤S3.2：对源网络的显式特征和隐式特征进行训练，得到分类模型；

步骤S3.3：将分类模型应用于目标网络的连边中，获得目标符号网络的边的符号，获得未来时刻的有符号的网络快照。

5.如权利要求2所述的方法，其特征在于，步骤S1.1具体包括：

将每组模体中三个节点的连边状况与预先存储的基本模体的连边状况进行比对，确定所属的模体类型。

6.如权利要求4所述的方法，其特征在于，步骤S2.1具体包括：

步骤2.1.1：通过求解目标函数的方式，进行非负张量分解，其中，目标目标函数的形式如下：

s.t.λ_n≥0,a_n≥0,b_n≥0,c_n≥0,n＝1,2,…,r

步骤2.1.2：采用MATLAB工具中预设函数执行上述目标函数，得到分解结果。

7.一种动态基因调控网的网络演化分析装置，其特征在于，包括：

预处理模块，用于将基因调控网以模体的形式进行表示，统计快照间的模体转换概率，并将相邻两个快照间的模体转换概率用矩阵表示，得到模体转换概率矩阵，模体为由三个节点构成的子图，快照为以预设时间间隔对基因调控网络进行采样，得到的预设时刻的静态结构，模体转换概率矩阵中的元素用以表征模体从一个时刻到下一时刻的变化情况；

模体转换模块，用于以包含前T-1个时刻的模体转换概率矩阵的模体转换概率张量为基础，预测t时刻的模体转换概率矩阵，获得无符号的网络快照，其中，T表示快照的总数量，t表示时刻；

符号判别模块，用于基于隐空间特征，将无符号的基因调控网的边通过预设非负矩阵三因式分解的方法映射到隐空间中，并以连边在隐空间位置的坐标为特征、连边的符号为标签，通过机器学习方法进行样本的训练和预测，获得未来时刻的有符号的网络快照。

8.如权利要求7所述的装置，其特征在于，预处理模块具体用于执行下述步骤：

步骤S1.1：将网络中所有节点每三个组成一组，进行分组和编号，每组为一个模体，根据三个节点间连边的不同，共有64种可能的模体类型；

步骤S1.2：将相邻两个快照间模体转换概率用一个64×64的矩阵进行表示，称为模体转换概率矩阵，记作TCM，其中，矩阵中元素的值TCM_t(i,j)＝P(tri_t[i]→tri_t[j])表示从t时刻到t+1时刻编号为i的模体转换为编号为j的模体的概率，其中tri_t[i]表示t时刻编号为i的模体；

步骤S1.3：将不同时刻的模体转换概率矩阵合成模体转换概率张量(TCT)，其中，TCT＝(TCM₁,TCM₂,…,TCM_T-1)，张量中的元素TCT(i,j,t)＝TCM_t(i,j)。

9.如权利要求1所述的装置，其特征在于，模体转换模块具体用于执行下述步骤：

步骤S2.1：对模体转换概率张量进行非负张量分解，分解得到三个因式矩阵A,B,C，其中AB^T表示不同类型模体间的转换关系，C包含上述转换关系在时间维度上的信息，称为时间因式矩阵；

步骤S2.2：采用指数分布模型预测时间因式矩阵的第T行，计算方式如下：

其中，T为常量，表示快照的总数量，t为变量，表示从1到T中的某一个时刻，a表示权重参数，r表示矩阵A、B、C的一列，若共有R列，则r取值为1到R之间的任意整数；

步骤S2.3：根据预测得到的T时刻的时间因式矩阵获得T时刻的模体转换概率矩阵，称为模体转换似然矩阵(TCLM)，计算方式如下：

其中，λ_r为非负张量分解的结果，

A(i,r)为矩阵A的第i行第r列元素，B(j,r)为矩阵B的第j行第r列元素,C(T,r)为矩阵C的第T行第r列元素。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被执行时实现如权利要求1至6中任一项权利要求所述的方法。