基于回归理论恢复基因调控网络

Restoration of Gene Regulatory Network Based on Regression Theory

摘要: 基因之间的调控关系隐含在基因表达数据里，需要分析该数据从而揭示基因调控网络的拓扑结构。由于静态基因表达数据的样本较少，因此本文提出基于距离相关性扩充样本数据量的方法。接着，本文提出恢复基因调控网络拓扑结构的方法，基于距离样本数据根据回归理论建立基因调控网络线性回归模型，对模型应用最小二乘估计和假设检验判断基因之间是否存在调控关系。此外，提出可以控制假阳性的方法，利用统计检验控制错误发现率提高模型预测的准确性。最后，在DREAM3数据集上验证方法的可行性。

Abstract: The regulatory relationship between genes is implicit in the gene expression data, which needs to be analyzed to reveal the topology of the gene regulatory network. Since the small sample size of static gene expression data, this paper proposes a method to expand the sample data size based on distance correlation. Then, this paper proposes a method to restore the topology of gene regulatory network. Based on the distance sample data, a linear regression model of gene regulatory network is established according to the regression theory. The least square estimation and hypothesis test-ing are applied to the model to determine whether there is a regulatory relationship between genes. In addition, a method to control false positives is proposed. Statistical test is used to control the false discovery rate to improve the accuracy of model prediction. Finally, the feasibility of the method is verified on the DREAM3 dataset.

1. 引言

基因调控网络(Gene Regulatory Network, GRN)是由细胞内基因之间相互作用构成的网络。节点代表基因，节点间的连边代表基因之间的调控关系。它是一个抽象的概念，但可以利用数学表达式来建立基因调控的过程，从而帮助人们理解细胞内基因之间的关系。重构基因调控网络属于逆向工程问题，目前该领域的学者提出了很多用于研究基因调控网络的数学模型 [1] 。例如最简单的、离散的布尔网络(Boolean Networks)模型 [2] ；一种概率描述的贝叶斯网络(Bayesian Networks)模型 [3] 以及微分方程(Differential Equation)模型 [4] 等。

基因表达数据反映的是直接或间接测量得到的基因转录产物 mRNA在细胞中的丰度，它包含的基因活动信息能够反映当前细胞的状态，采用反向工程的方法分析基因表达数据可以还原基因的调控网络。基因表达数据可以分为时间序列基因表达数据和静态基因表达数据，它们都是基因表达水平的数据矩阵。时间序列基因表达数据记录了生命活动周期或不同状态下连续变化的基因表达水平，反映了基因调控网络的动态变化，可以通过相应的算法推断基因之间的调控关系 [5] ；静态基因表达数据记录了稳定状态下细胞内的基因表达水平，包含了基因间的调控信息，直接反映了基因之间的调控关系，可以通过基于信息论等方法建立基因调控网络 [6] [7] 。静态基因表达数据具有高维度、小样本的特点，即与大量的基因相比，样本数量较少。对于这一困难多种方法被提出验证。Fujita等对于样本数据量低于基因数量的问题，提出SVAR方法建模基因表达调控网络 [8] ；Barrera等利用互信息降维的方法，通过计算基因之间的互信息，最终选择基因之间互信息较高的基因恢复基因调控网络 [9] ；Székely等扩充了样本数据量，通过距离相关性概念将m个样本数据扩充为m²个样本数据。通过增加样本数据 [10] ，可以提高模型的预测性能。

真实生物体基因间的调控关系错综复杂，模型识别真正的边显得十分重要。而模型在恢复基因调控网络拓扑结构过程中可能会出现真实基因中没有调控关系但被模型错误的预测为有调控关系的情况，即预测结果中出现比较多假阳性边的情况。因此需要降低假阳性边的数量，提高模型预测的准确率。Fujita等提出通过统计检验控制错误发现率的方法 [8] 。罗霄提出通过DSWLasso逐步删除恢复的基因调控网络中的弱调控关系的方法 [11] 。

本文在此基础上进行研究，将基因调控网络的恢复问题转化为识别方程参数的问题。通过假设检验判断基因间是否存在调控关系，此时将基因调控网络拓扑结构恢复问题转化为分类问题进行解决。并采取适当的方法减少结果中出现的假阳性边的数量。

2. 方法

回归分析是用来确定两种或两种以上变量之间相互依赖的定量关系的一种统计方法。本文所研究的基因调控网络包含多个自变量，并且因变量和自变量之间呈线性关系，因此可以利用多元线性回归模型判断基因调控网络中基因间的调控关系。

2.1. 基因调控网络线性回归模型

基于距离样本数据采用回归理论建立基因调控网络线性回归模型，以此判断基因i与其他基因之间是否存在调控关系，模型表示为：

$X_i=X_{\backslash i}{B}_i+{\epsilon }_i$ (1)

其中 $X_{\backslash i}=\left(X_1,X_2,\cdots ,X_{i-1},X_{i+1},\cdots ,X_N\right)$ 表示除第i个基因外，其余 $N-1$个基因的距离样本向量， ${\epsilon }_i$ 代表误差，满足相互独立且分布服从 $N\left(0,{\sigma }^2\right)$ 。

样本数据代入(1)中，则(1)可以改写为：

$Y_i=\left(\begin{array}{c}{x}_{i1}\\ x_{i2}\\ ⋮\\ x_{iM}\end{array}\right)=\left(\begin{array}{ccccccc}{x}_{11}& x_{21}& \cdots & x_{i-1,1}& x_{i+1,1}& \cdots & x_{N1}\\ x_{12}& x_{22}& \cdots & x_{i-1,2}& x_{i+1,2}& \cdots & x_{N2}\\ ⋮& ⋮& \cdots & ⋮& ⋮& \cdots & ⋮\\ x_{1M}& x_{2M}& \cdots & x_{i-1,M}& x_{i+1,M}& \cdots & x_{NM}\end{array}\right)\left(\begin{array}{c}\begin{array}{l}{b}_{i1}\\ b_{i2}\end{array}\\ ⋮\\ b_{i,i-1}\\ b_{i,i+1}\\ ⋮\\ b_{iN}\end{array}\right)+\left(\begin{array}{c}{\epsilon }_{i1}\\ ⋮\\ {\epsilon }_{i,i-1}\\ {\epsilon }_{i,i+1}\\ ⋮\\ {\epsilon }_{iN}\end{array}\right)=X{B}_i+E_i$ (2)

恢复基因调控网络拓扑重构解释为从可测变量 $X$ 和 $Y_i$ 中估计 $B_i$ ，则转化为求解多元线性回归方程(2)的参数 $B_i$ 的问题。在这里采用对每个基因进行单独回归分析的方法，用最小化误差的平方和估计每个基因的 $B_i$ ，得到预测的邻接矩阵 $\hat{B}$ 。损失函数利用均方误差的形式，同时它也是目标函数：

$L\left(B_i\right)=\frac{1}{2}{\left(X{B}_i-Y_i\right)}^{\text{T}}\left(X{B}_i-Y_i\right)$ (3)

为确定模型的参数 $B_i$ ，同时使预测值 ${\hat{Y}}_i$ 尽可能地接近真实值 $Y_i$ ，等价于求 $L\left(B_i\right)$ 的最小值。由于目标函数(3)中包含要求解的参数 $B_i$ ，因此在求解(3)最小值的过程中便可以将模型(2)中的参数 $B_i$ 求解出来，在这里我们选取最小二乘法求解。用最小二乘法估计 $B_i$ ，只需要令其对 $B_i$ 的一阶偏导数为0，偏导数为0的点是(3)的最小值点，求得：

${\hat{B}}_i={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}{Y}_i,i=1,\cdots ,N$ (4)

(4)就是 $B_i$ 的最小二乘估计，该方法适用于 $X^{\text{T}}X$ 的逆矩阵存在 $\left(M>N\right)$ 时的情况。此时 ${\hat{B}}_i$ 是 $B_i$ 的无偏估计，无偏性保证了最小二乘估计具有较好的结果。

2.2. 基因表达数据处理

基因之间的调控关系隐含在基因表达数据里，需要分析该数据从而揭示基因调控网络的拓扑结构。而现实难以收集到大量的基因表达数据，但对于(4)式，当 $M<N$ 时， $X^{\text{T}}X$ 不可逆，此时最小二乘估计失效，因此需要扩充样本量。本文根据Székely的方法引入距离相关性来扩充样本量。假设有N个基因，实验收集到这N个基因在m个样本上的表达水平，记为 $G=\left(G_1,G_2,\cdots ,G_N\right)$ ，其中 $G_i={\left(g_{i1},g_{i2},\cdots ,g_{im}\right)}^{\text{T}},i=1,\cdots ,N$ 。根据距离相关性和 $G$ 可以生成一个矩阵 $J$ ， $J$ 包含m²个元素，各元素定义为：

$j_{sp}={\Vert g_{is}-g_{ip}\Vert }_2,i=1,2,\cdots ,N$

其中 ${\Vert \cdot \Vert }_2$ 表示 $l_2$ 范数。对矩阵 $J$ 中心化可得对称矩阵 $D$ ， $D$ 中元素定义为：

$d_{sp}=j_{sp}-\frac{1}{m}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{j}_{sk}-}\frac{1}{m}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{j}_{kp}}+\frac{1}{m^2}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{j}_{kl}}}$

将 $D$ 中下三角元素重新组合为一列，可得到一个包含 $M=\frac{m\left(m+1\right)}{2}$ 个元素的列向量 $X_i$ ， $X_i$ 就是对

原始基因表达数据 $G_i$ 的扩充，由原来的m个样本扩充为M个样本，因此达到了扩充样本基因表达数据的目的，称 $X_{}=\left(X_1,X_2\cdots ,X_N\right)$ 为距离样本数据。

2.3. 显著性检验

2.3.1. 基因调控网络线性回归模型的可行性

本文采用多元线性回归方程预测基因i和其它基因间的调控关系，建立的恢复模型为(2)。但基因i和其它基因之间是否为线性关系是未知的，因此需要判断线性关系是否成立。通过对模型(2)进行整体检验来判断回归系数是否不全为0，如此便可说明其他基因是否对基因i有显著影响。如果所有回归系数为0，则代表它们对基因i没有线性调控关系，因此这相当于检验原假设 $H_0^i$ 是否成立：

$H_0^i:b_{i1}=b_{i2}=\cdots =b_{i,i-1}=b_{i,i+1}=\cdots =b_{iN}=0,i=1,\cdots ,N$

$H_0^i$ 对应的备择假设为：

$H_1^i:b_{ij}\ne 0,j=1,\cdots ,N-1$

如果有充分的理由说明 $H_0^i$ 不成立，那么意味着基因i与其它 $N-1$ 个基因之间的线性关系成立。对模型整体检验选取 $F$ 检验。在 $H_0^i$ 为真时构造检验统计量：

$S_R={{\displaystyle \sum \left(X_{\backslash i}{\hat{B}}_i-{\bar{Y}}_i\right)}}^2$

$S_e={{\displaystyle \sum \left(Y_i-X_{\backslash i}{\hat{B}}_i{}_i\right)}}^2$

$F=\frac{S_R/\left(N-1\right)}{S_e/\left(M-N-2\right)}\sim F\left(N-1,M-N-1\right)$

其中 $F\left(N-1,M-N-1\right)$ 表示自由度为 $\left(N-1,M-N-1\right)$ 的 $F$ 分布，取显著性水平为 $\alpha \left(0<\alpha <1\right)$ 。我们选定 $\alpha =0.05$ ，当 $F\ge F_{1-\alpha }\left(N-1,M-N-1\right)$ 时拒绝原假设 $H_0^i$ ，认为基因i与其它 $N-1$ 个基因之间有线性关系，说明利用基因调控网络线性回归模型恢复拓扑结构是可行的。

2.3.2. 基因调控关系的判定

模型整体检验说明基因i与其它 $N-1$ 个基因之间的线性关系成立，但是不能确定某个基因j对基因i是否有调控关系。因此需要对回归系数进行检验，即讨论基因与基因之间是否有调控关系。

本节采用的基因调控矩阵元素为0和1，其中0代表两个基因之间没有调控关系，1表示两个基因之间有调控关系。利用模型(2)判断哪些基因之间存在调控关系，等价于寻找非0的 $b_{ij}$ 。因此为检验基因i与其他基因之间是否存在调控关系，建立如下假设：

$H_0{}^{ij}:b_{ij}=0,H_1{}^{ij}:b_{ij}=1$

如果有充分的理由说明 $H_0^{ij}$ 不成立，那么就可以接受 $H_1^{ij}$ ，认为基因i与基因j之间存在调控关系。选取的假设检验方法为t检验。已知 ${\hat{b}}_{ij}\sim N\left(b_{ij},u_{jj}{\sigma }^2\right)$ ，当假设 $H_0^{ij}$ 为真时构造检验统计量：

$t_{ij}=\frac{{\hat{b}}_{ij}}{\sqrt{u_{jj}{\hat{\sigma }}^2}}\sim t\left(m-N-1\right)$

其中 $t\left(m-N-1\right)$ 表示自由度为 $\left(M-N-1\right)$ 的t分布， $u_{jj}$ 为 ${\left(X^{\text{T}}X\right)}^{-1}$ 的第j个对角线元素， ${\hat{\sigma }}^2=\frac{S_e}{M-N-1}$ 为 ${\sigma }^2$ 的无偏估计。取显著性水平为 $\alpha \left(0<\alpha <1\right)$ 。我们选定 $\alpha =0.05$ ，当 $\left|t_{ij}\right|\ge t_{1-\frac{\alpha }{2}}\left(M-N-1\right)$ 时拒绝原假设 $H_0^{ij}$ ，接受备择假设 $H_1^{ij}$ ，认为 $b_{ij}=1$ ，说明基因j对基因i的作用显著，即确定基因j对基因i有调控关系；反之，当 $\left|t_{ij}\right|<t_{1-\frac{\alpha }{2}}\left(M-N-1\right)$ 时不拒绝原假设 $H_0^{ij}$ ，则认为 $b_{ij}=0$ ，说明基因j对基因i的作用不显著，即确定基因j对基因i没有调控关系。

通过对模型(2)基于距离样本数据利用最小二乘估计得到初始的预测矩阵后，接着对回归系数进行假设检验便可以确定 $\hat{B}$ 中元素为1的位置， $\hat{B}$ 中元素为1意味着基因之间有调控关系，此时模型只保留了对基因i有显著作用的调控基因。因此可根据构建的模型确定具有调控关系的边，即完成了恢复基因调控网络的拓扑结构。

对模型恢复的基因调控网络拓扑结构进行整体评估，定义评价指标：

$E=\frac{FN}{FN+TP}$ 为模型未识别率，其中 $FN$ 表示基因间有调控关系而模型没有识别出的个数， $FN+TP$ 表示基因间真实又调控关系的个数；

$FPR=\frac{FP}{FP+TN}$ 为假阳性率，其中 $FP$ 表示基因间没有调控关系而模型识别为有调控关系的个数， $FP+TN$ 表示基因间没有调控关系的个数；

$FDR=\frac{FP}{FP+TP}$ 为误报率，其中 $FP+TP$ 表示模型识别基因间有调控关系的个数。

2.4. 控制假阳性结果的数量

t检验的拒绝域为 $W=\{\left|t\right|\ge t_{1-\frac{\alpha }{2}}\left(M-N-1\right)\}$ ，而由于随机性等因素的干扰，在进行判断时可能出现两

种错误，假阴性错误和假阳性错误。若出现原假设成立，但计算得到的 $b_{ij}$ 落入拒绝域 $W$ 中的错误，称为假阳性错误(第Ⅰ类错误)。在Neyman-Pearson假设检验中，可以要求假阳性率在可接受的范围内，设置的检验水平 $\alpha =0.05$ 就是限制犯假阳性错误的概率，通过 $\alpha$ 达到对单重假设检验的统计推断的错误控制。在单个假设检验中假阳性率是可控的，但是对于整体多个假设而言，假阳性结果会随着检验次数的增加而增大，出现超出预设范围的情况。因此需要对回归系数整体进行多重检验校正，减少出现假阳性结果的次数。选择控制FDR的BH(Benjamini-Hochberg)法，该方法可以将假阳性率控制在统计显著水平 $\alpha$ 以下。调整 $p$ 值控制FDR，过程如下：

对于假设检验 $\left\{H_0^{ij}\right\}i,j=1,\cdots ,N,i\ne j$ ， $H_0^{ij}$ 表示第 $ij$ 个检验的原假设。将所有检验的 $p$ 值从小到大排序生成顺序数 $\left\{p_k\right\},k=1,\cdots ,v$ 。令 $q_{\max }=p_{\max }$ ，定义

$q_k=\min \left\{p_k\frac{v}{k},p_{k-1}\right\},k=1,\cdots ,v-1$

其中v代表检验次数。当 $q_k<\alpha$ 时，拒绝前k个假设，认为基因之间有调控关系。

3. 模型应用及恢复结果

3.1. DREAM数据集上的实验

通过以上分析可以根据模型(2)推断基因调控网络，为了验证其可行性，本节在DREAM(Dialogue for Reverse Engineering Assessments and Method)提供的静态基因表达数据上进行实验。采用DREAM3数据集，数据集可在https://gnw.sourceforge.net/dreamchallenge.html下载。DREAM3数据集提供了静态基因表达数据和确定的基因调控网络，本节实验选取数据集InSilicoSize10-Ecoli1-null-mutants，它包含10个基因，将基因标号为 $G1,\cdots ,G10$ 。该数据集对应的网络结构如图1所示。其中红线代表基因之间是抑制作用，蓝线代表基因之间是激活作用。本文不区分基因之间的抑制和激活作用，主要关注基因之间是否存在调控关系。

Figure 1. Topology of gene regulatory network

图1. 基因调控网络拓扑结构

由于实验会受到随机性的影响，为了确保模型预测的稳定性，结果为多次实验的平均。每次实验对距离样本采用随机抽样用于模型(2)来构建基因调控网络，每次抽样循环多次并取平均值作为最终结果。

3.2. 实验结果及讨论

3.2.1. 基因调控网络线性回归模型的可行性

DREAM3提供了确定的基因调控网络拓扑图，本文只考虑基因之间是否存在调控关系，因此调控矩阵中的元素只有0，1两种。数据集InSilicoSize10-Ecoli1-null-mutants中仅包含11个基因表达数据，首先对其进行扩充得到包含66个样本数据距离样本，从中随机选取不同长度的样本用于恢复拓扑结构。

对每个基因在不同样本数据时的模型(2)进行检验以判断利用其恢复基因调控网络拓扑结构是否可行，表1包含了在不同数据下对每一个模型进行整体检验的F值。

在 $\alpha =0.05$ 情况下，由表1以及Friedman检验表可知，对于G8，G9，G10基因，由于 $F<F_{1-\alpha }\left(9,M-11\right)$ ，因此不拒绝原假设 $H_0^i,i=8,9,10$ ，认为回归方程的所有系数为0，与真实的基因调控网络一致；对于另外7个基因，由于 $F\ge F_{1-\alpha }\left(9,M-11\right)$ ，因此拒绝原假设 $H_0^i$ ，认为每一个基因 $i,\left(i=1,\cdots ,7\right)$ 与其它基因之间有线性关系，说明利用线性回归模型预测基因调控矩阵是可行的。

Table 1. F-value of the overall model test

表1. 模型整体检验的F值

3.2.2. 基因调控关系的判定及控制假阳性结果的数量

对于原假设 $H_0{}^{ij}:b_{ij}=0$ 和备择假设 $H_1{}^{ij}:b_{ij}=1$ ，采用t检验可以判定矩阵 $\hat{B}$ 中位置为1的元素，即可以判定基因之间的调控关系。

图2显示曲线呈下降趋势，模型的未识别率E在逐渐降低，说明随着数据量M的增加，模型识别到的连接个数逐渐增加。但初始样本数据会限制模型对基因调控网络拓扑结构的恢复能力，当 $M=55$ 时， $E=0.4802$ ，此后距离样本再增加不会明显提高模型的恢复能力。

Figure 2. For the gene regulatory network linear regression model (2), the unrecognized rate E with its standard deviation under different sample data size M

图2. 对于基因调控网络线性回归模型(2)，不同样本数据量M下的未识别率E及其标准差

由图3可知， $FP=5$ ，说明模型(2)在恢复基因间的调控关系时出现较多的假阳性边，即出现较多实际基因没有调控关系但模型错误的判断为有调控关系的结果。图4是在 $M=40$ 时对模型(2)应用控制FDR的BH法后恢复的基因调控网络中假阳性边的情况。由图可知 $FP=1$ ，此时假阳性率低于0.05。与图3相比，模型预测的假阳性由5条减少为1条，达到了控制模型识别假阳性边的个数的目的。因此对于在构造基因调控网络出现假阳性边较多的情况，可以采用控制FDR的BH法来减少出现模型恢复过程中出现假阳性边的次数。

不失一般性，在不同数据量时均可对整体多个假设采用FDR校正。由图5可得，对于不同的样本数据量M，FDR校正后假阳性率较之前会降低，达到了对构建的基因调控网络假阳性边数量控制的目的，假阳性率始终在(0, 0.05)的范围之内。

Figure 3. For the gene regulatory network linear regression model (2), the restored gene regulatory relationship when M = 40. The solid line represents the true gene regulatory relationship, and the dashed line indicates the edges where there is no regulatory relationship between genes but the model identifies as having a regulatory relationship

图3. 对于基因调控网络线性回归模型(2)，M = 40时模型恢复的基因调控关系。实线表示真实的基因调控关系，虚线表示基因间不存在调控关系但模型识别为有调控关系的边

Figure 4. For the gene regulatory network linear regression model (2), the restored gene regulatory relationship when M = 40 after FDR correction. The solid line represents the true gene regulatory relationship, and the dashed line indicates the edges where there is no regulatory relationship between genes but the model identifies as having a regulatory relationship

图4. 对于基因调控网络线性回归模型(2)，经过FDR校正后在M = 40时模型恢复的基因调控关系。实线表示真实的基因调控关系，虚线表示基因间不存在调控关系但模型识别为有调控关系的边

Figure 5. For the gene regulatory network linear regression model (2), the probability and standard deviation of false positive error under different sample data size M. The red line and blue line respectively represent FPR values before and after FDR correction

图5. 对于基因调控网络线性回归模型(2)，不同样本数据量M下犯假阳性错误的概率及其标准差。红线蓝线分别代表FDR校正前和校正后的FPR的值

由于假阳性错误和假阴性错误在相同的样本数据量下不能同时增大或减小，只有增加样本数据量才可以使两种错误同时减少，因此采用FDR校正降低模型识别假阳边的数量后，不可避免的会增加假阴性边的个数。但由图6可知随着数据量的增加，FDR校正前后模型出现假阴性错误的差距越来越小。

Figure 6. For the gene regulatory network linear regression model (2), the unrecognized rate E with its standard deviation under different sample data size M. The red and black lines represent the unrecognized rate E before and after correction, respectively

图6. 对于基因调控网络线性回归模型(2)，不同样本数据量M下的未识别率E。红色线和黑色线分别代表校正前和校正后的未识别率E

4. 结论

基因调控网络存在于生物体内，是由细胞内基因之间相互作用构成的网络，可以利用数学表达式建立基因调控的过程，这可以帮助我们了解细胞内基因之间的关系。静态基因表达数据记录了稳定状态下细胞内的基因表达水平，可以反映基因之间的调控关系，但是其样本数较少，本文对此介绍了扩充样本数据产生距离样本的方法，并基于距离样本根据回归理论建立了基因调控网络线性回归模型，通过最小二乘估计和假设检验来确定哪些基因之间存在调控关系。接着应用控制FDR的BH法有效地减少了模型识别假阳性结果的数量，提高模型预测的准确率。最后选取DREAM3数据集验证模型的可行性，并在M = 40时直观展现了该方法降低假阳性边的有效性。

基金项目

本文工作由国家自然科学基金(No. 11502062)支持。

NOTES

^*通讯作者。