CN105825077B

CN105825077B - 一种基于多元格兰杰检验的信号传导网络构建方法

Info

Publication number: CN105825077B
Application number: CN201610137995.5A
Authority: CN
Inventors: 李敏; 郑瑞清; 王建新
Original assignee: Central South University
Current assignee: Central South University
Priority date: 2016-03-14
Filing date: 2016-03-14
Publication date: 2018-10-09
Anticipated expiration: 2036-03-14
Also published as: CN105825077A

Abstract

本发明公开了一种基于多元格兰杰检验的信号传导网络构建方法，对输入的时序基因表达数据进行二范数规范化；对二范数规范化后的t个时间片的基因表达数据，分别构建1～t‑1和2～t时序基因表达矩阵，并建立基于阶马尔卡夫的线性方程，通过奇异值分解计算获得系数矩阵和表达值矩阵，并根据所得的系数矩阵计算残差值Rss*；对每一个基因都将其时序基因表达打乱，并计算打乱后的残差值，获得残差分布Rss,对Rss*在Rss中进行排序，获得系数的显著性α。本发明的方法适用于真实的时序基因表达数据时间片过少的情况，并有效提高了有向推断的准确性。

Description

一种基于多元格兰杰检验的信号传导网络构建方法

技术领域

本发明涉及生物信息学领域，涉及复杂生物网络的构建方法。

背景技术

在生物系统中，常见的以网络形式存在的生物结构包括基因转录调控网络、生物代谢与信号传导网络、蛋白质相互作用网络。其中，基因转录调控网络、生物代谢与信号传导网络都有传递的有向性。基因调控网络是生物体内控制基因表达的机制。研究表明，基因以及基因产物并非是单独起作用的，而是参与在复杂的、相互联系的通路和网络中。这些基因相互作用、相互影响，从而产生了细胞、组织、器官和有机体的机理，形成“基因网络”这样一个复杂系统来推动演化的。代谢网络是一种由代谢通路组成的抽象网络。代谢网络与基因调控网络在生物特性以及网络特性上存在一定的不同。细胞的代谢网络由于代谢的特性，是一种无尺度网络。生物中的信号传导(signal transduction)则是指细胞将一种类型的信号或刺激转换为其他生物信号最终激活细胞反应的过程。同代谢通路一样，信号传导的过程中多个生物分子在酶的作用下发生一连串有序的反应，由此得到了信号传导通路。信号传导网络即是指参与信号传导通路的分子和酶以及其间所发生的生化反应所构成的网络。这些网络中的关系一般都表现为抑制和促进基因表达。磷酸化和去磷酸化，甲基化以及激活、制止作用。随着高通量检测技术的发展，基因表达数据变得更为全面和复杂，从而有力的推动了生物网络的构建。基因表达数据反映的是直接或间接测量得到的基因转录产物mRNA在细胞中的丰度，这些数据可以用于分析基因之间有何相关性，在不同条件下基因的活动是如何受影响的。基于基因表达数据的特点以及生物网络互相作用的特性，可以从局部或者全局的角度去推断基因之间作用的有向性。现有的使用基因表达数据的生物有向网络的构建可以分为两大类：

(1)概率模型：使用基因表达数据计算表达作用关系的概率。其中贝叶斯模型是典型的概率模型。它的基本思想是使用简单的局部概率乘积来近似复杂的高维概率分布。贝叶斯网络引入有向无圈图模型和隐马尔可夫链来描述变量间的联系与相互作用，构建有向网络模型。贝叶斯模型往往需要在没有先验先验概率的情况需要先通过基因表达值的分布，构建先验知识，然后通过计算所得的先验概率计算的到后验概率。在贝叶斯网络的基础上，研究者还提出了动态贝叶斯网络模型(DBN)，这种模型和普通贝叶斯网络模型不同之处在于它考虑了同一基因前后时间点的表达向量，这种模型的优势在于可以将调控的负反馈和延时因素考虑进去，克服了普通贝叶斯网络是一个无环图带来的不足。贝叶斯网络的方向由于需要计算先验概率，往往计算时间较长。

(2)线性模型：线性模型主要是通过构建线性回归的函数，通过对回归方程以及系数的显著性评估来判断基因之间作用的方向。常见的线性模型就是格兰杰检验，主要包括二元和多元格兰杰检验。二元格兰杰检验是对任意两个基因分别作为因变量和自变量构建两个线性回归方程，并对其中因变量的系数做零假设，然后比较这两个零假设的显著性来确定两个的因果关系。多元格兰杰检验是在二元的基础上，将自变量变为多个基因，然后对其中每个自变量做零假设检验。多元模型可以有效的减少二元模型中的间接的因果关系。

目前，虽然已存在一些模型和方法构建有向生物网络，但是这些模型和方法的生物网络构建的精确度还有待于进一步提高。此外，由于生物基因表达数据的时间片特点，一些传统的格兰杰检验方法无法适用。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种基于多元格兰杰检验的信号传导网络构建方法。

为解决上述技术问题，本发明所采用的技术方案是：一种基于多元格兰杰检验的信号传导网络构建方法，包括以下步骤：

1)对输入的时序基因表达数据进行二范数规范化；

2)对二范数规范化后的t个时间片的基因表达数据，分别构建1～t-1和2～t时序基因表达矩阵，并建立基于阶马尔卡夫的线性方程，通过奇异值分解计算获得系数矩阵和表达值矩阵，并根据所得的系数矩阵计算残差值RSS；

3)对每一个基因都将其时序基因表达打乱，并计算打乱后的残差值，获得残差分布RSS^*,对RSS在RSS^*中进行排序，获得系数的显著性α。

二范数规范化的表达式为：g_i，t表示基因i在t时刻的表达值；M表示基因表达数据时间片；x_i，t为二范数规范化后的基因表达值。

系数矩阵和表达值矩阵获取的具体过程为：

1)基于时间平稳性和一阶马尔科夫模型的假设，构建线性回归方程：x_i，t＝∑_{0＜j＜N，i≠j}r_j，t-1x_j，t-1；其中x_j，t-1表示基因j在t-1时刻下的表达值，r_j，t-1代表该表达值对基因i在时刻t-1的影响系数，N表示总的基因个数；

2)根据线性方程获得矩阵方程：

其中r_i，j表示基因i对基因j的影响系数；

3)将步骤2)的矩阵分别用R，T和T^*表示：R*T＝T^*；通过奇异值分解，得到系数矩阵R为：R＝T^*VS^-1U^T；其中，T为表达值矩阵，T＝U*S*V^T，U和V为酉矩阵，S为奇异值分解产生的对角矩阵，s(i,i)为S中的元素。

显著性α的计算过程为：

1)利用表达值矩阵T，根据系数矩阵R按行计算每个基因的残差RSS；

2)将基因i的表达值的顺序打乱，并用打乱后的表达值替换替换前的基因i的表达值；

3)使用奇异值分解获得新的系数矩阵R^*，并重新计算该系数矩阵下的基因残差值RSS^*,并保存；

4)重复步骤1)—2)共k次，根据计算所得的k个残差值估计RSS^*的分布；

5)对RSS在获得的RSS^*分布中以升序进行排名，记录排名值为rank，来获得显著性α:(α＝((k+1-rank))/(k+1))。

与现有技术相比，本发明所具有的有益效果为：本发明的方法适用于真实的时序基因表达数据时间片过少的情况，并有效提高了有向推断的准确性。

附图说明

图1本发明MGT-SM的流程图；

图2为方法MGT-SM和DBN,CGC2SPR,Bivariate Granger方法基于ROC曲线以及其下的面积(AUC)的比较图；

图3为方法MGT-SM和DBN,CGC2SPR,Bivariate Granger方法基于Recall方法评价的比较图。

具体实施方式

一、基因表达数据的预处理

读入时序基因表数据文件，以g_i，t表示基因i在t时刻的表达值，对其在该基因的时序表达组下进行二范数规范化：

M表示基因表达数据时间片个数规范化后的x_i，t在后续步骤中作为基因表达值被使用。

二、构建系数矩阵

基于时间平稳性和一阶马尔科夫模型的假设，构建线性回归方程：

x_i，t＝∑_{0＜j＜N，i≠j}r_j，t-1x_j，t-1 (2)

其中x_j，t-1表示基因j在t-1时刻下的表达值，r_j，t-1代表该表达值对基因i在时刻t-1的影响系数，N表示总的基因个数。根据线性方程获得矩阵方程：

其中r_i，j表示基因i对j的影响系数。对等式中矩阵分别用R，T和T^*表示，

R*T＝T^* (4)

通过奇异值分解可解得系数矩阵R为：

R＝T^*VS^-1U^T (5)

其中T＝U*S*V^TU和V为酉矩阵，S为奇异值分解产生的对角矩阵,

奇异值分解对表达系数矩阵的预测有很好的效果，同时使用奇异值分解，解决了传统格兰杰检验在时间片少的情况下无法使用最小二乘法解得系数的问题。

三、显著性评估

在该步骤中，我们需要对每个系数的显著性进行评估，由于传统的格兰杰检验在T＜np-p时无法使用F检验，因此我们采用蒙特卡罗模拟的方法获得残差的分布，并对残差进行排序。

由于大部分的生物有向网络都是单向边，因此，根据有向边(i,j)和(j,i)显著性α的大小，我们将显著性较大方向的边保存为一个文件，将较小方向的边保存为另一个文件，并将两个文件内的显著性单独做排序。

四、实验验证

为了验证本方法的有效性，我们在乳腺癌的真实数据的4个样本进行了测试，并和传统的二元格兰杰检验和动态贝叶斯(DBN)以及使用岭回归的格兰杰检验CGC2SPR进行比较分析。乳腺癌信号传导网络的数据集来源于Neve等人发表在Cancer Cell的有明显特征的基础乳腺癌亚型的公开数据集，细胞系中包含了EGFR扩增和PTEN、Rb1、Smad4和p53突变。该数据也是人类乳腺癌网络推断挑战赛HPN-DREAM的标准数据集(https:// www.synapse.org/#！Synapse:syn1720047)。该数据集包含4个样本，每个样本有8个时间片。整个乳腺癌信号传导网络包含20个基因和48条有向边。

为了评价推断结果的连续性和准确性，我们采用AUC值和Recall值两个指标进行比较。AUC值是ROC曲线的线下面积，Recall值计算公式为方向预测正确的边数和真实有向边数的比值。AUC值和Recall值的实验结果分别如图2和图3所示。

从图2和图3可以看出，我们的方法在不同样本下，无论基于AUC值还是Recall值均优于其他的方法。由此可见我们提出的方法MGT-SM具有很好的稳定性。有实验表明，DBN在网络规模较大的情况下，时间复杂度呈指数级增长，基本无法计算，而我们的方法MGT-SM受益于奇异值分解和蒙特卡罗模拟，仍能在较短的时间内获得结果。

Claims

1.一种基于多元格兰杰检验的信号传导网络构建方法，其特征在于，包括以下步骤：

a)对输入的时序基因表达数据进行二范数规范化；二范数规范化的表达式为：g_i,t表示基因i在t时刻的表达值；M表示基因表达数据时间片；x_i,t为二范数规范化后的基因表达值；

b)对二范数规范化后的t个时间片的基因表达数据，分别构建1～t-1和2～t时序基因表达矩阵，并建立基于时间平稳性和一阶马尔科夫模型的线性方程，通过奇异值分解计算获得系数矩阵和表达值矩阵，并根据所得的系数矩阵计算残差值RSS；

系数矩阵和表达值矩阵获取的具体过程为：

a)基于时间平稳性和一阶马尔科夫模型的假设，构建线性回归方程：x_i,t＝∑_0<j<N,i≠ _jr_j,t-1x_j,t-1(2)；其中x_j,t-1表示基因j在t-1时刻下的表达值，r_j,t-1代表该表达值对基因i在时刻t-1的影响系数，N表示总的基因个数；

b)根据线性方程获得矩阵方程：

其中r_i,j表示基因i对基因j的影响系数；

c)将步骤b)的矩阵分别用R，T和T^*表示：R*T＝T^*(4)；通过奇异值分解，得到系数矩阵R为：R＝T^*VS^-1U^T(5)；其中，T为表达值矩阵，T＝U*S*V^T，U和V为酉矩阵，S为奇异值分解产生的对角矩阵， if s(i,i)≠0；s(i,i)为S中的元素；

3)对每一个基因都将其时序基因表达打乱，并计算打乱后的残差值，获得残差分布RSS^*,对RSS在RSS^*中进行排序，获得系数的显著性α；显著性α的计算过程为：

i.利用表达值矩阵T，根据系数矩阵R按行计算每个基因的残差RSS；

ii.将基因i的表达值的顺序打乱，并用打乱后的表达值替换替换前的基因i的表达值；

iii.使用奇异值分解获得新的系数矩阵R^*，并重新计算该系数矩阵下的基因残差值RSS^*，并保存；

iv.重复步骤1)—2)共k次，根据计算所得的k个残差值估计RSS^*的分布；

v.对RSS在获得的RSS^*分布中以升序进行排名，记录排名值为rank，来获得显著性α：α＝(k+1-rank)/(k+1)。