基因调控网络的边预测


扩展功能
	加入收藏夹

	复制引文信息

	加入引用管理器

	Email Alert

	RSS
本文作者相关文章
	黎妍

	张晓飞

	易鸣

	刘妍岩

基因调控网络的边预测

黎妍¹, 张晓飞², 易鸣³, 刘妍岩¹

1 武汉大学数学与统计学学院武汉 430072;
2 华中师范大学数学与统计学学院武汉 430079;
3 中国科学院武汉物理与数学研究所武汉 430071

收稿日期: 2014-12-10; 修订日期: 2015-06-20

基金项目: 国家自然科学基金(11275259,91330113)资助

作者简介: 易鸣,yiming@wipm.ac.cn

摘要: 为构建基因调控网络,提出了一个基于基因表达水平和网络反传递的算法.该算法用网络反传递思想来分析由传统相关性计算方法产生的间接效果,并考虑了调控网络的稀疏性,在模型算法中加入了控制网络稀疏性的l₁范数惩罚项.在大肠杆菌实验数据上测试该算法,这种方法提高了相关性分析对调控网络中边的预测能力,皮尔逊相关系数提高了6.42%,斯皮尔曼相关系数提高了5.92%,互信息提高了9.35%.总的来说,这个模型为修饰大量系统的相关性数据提供一种新思路,可以应用到网络边的预测和推断生物网络的控制动力学中.

关键词: 基因调控网络直接相关间接相关皮尔逊相关系数斯皮尔曼等级相关系数互信息

Link Prediction for the Gene Regulatory Network

Li Yan¹, Zhang Xiaofei², Yi Ming³, Liu Yanyan¹

1 School of Mathematics and Statistics, Wuhan University, Wuhan 430072;
2 College of Mathematics and Statistics, Central China Normal University, Wuhan 430079;
3 Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071

Received: 2014-12-10; Revised: 2015-06-20

Abstract: In order to build gene regulatory network, we proposed an algorithm based on gene expression levels and network anti-delivery. The algorithm used the network anti-delivery idea to analyze indirect effect produced from the traditional correlated calculation methods, and added an norm penalty term for controlling network sparsity after considering the sparsity of regulatory networks. We use the algorithm on the E. coli experimental data. This method improves the edge predictive ability of the correlation analysis on the regulatory network, Pearson correlation coefficient increased by 6.42%, Spearman correlation coefficient increased by 5.92%, mutual information improves 9.35%. Overall, this model provides a new idea on modifying a large number of system-related data, and can be applied to network edge prediction and the control dynamics of biological network inference.

Key words: Gene regulatory networks Direct correlation Indirect correlation Pearson correlation Spearman rank correlation Mutual information

1 引言

最近几年,网络科学在各种背景下在被广泛使用,包括分子和细胞生物学、社会科学、信息科学、数据挖掘等. 网络可以高效的表示变量之间的相互依存关系^[1]. 基因调控网络是由一组基因、蛋白质、小分子以及它们之间的相互调控作用所构成的一种生化网络. 网络中,节点代表基因, 两个节点间的连边表示这两个基因在物理或功能上具有某种相关性^[2]. 这种相关性通常用基因表达水平的相关系数来预测,但是节点间相关性的传递会产生一些间接的相互作用,并且随着网络的结构不断增大,节点之间的二阶传递、三阶传递或更高阶转递作用导致了更多间接边的产生. 传统测量方法比如皮尔逊相关系数、斯皮尔曼相关系数、互信息等并不能区分这种非直接的相互作用^{[3, 4, 5]}, 所以我们需要研究出一种方法能够从生物实验数据中推断出基因调控的真实网络结构. 目前已经有一些方法被用来推断网络中的直接作用关系. 例如,利用矩阵特征分解和无限级数求和的方法预测网络中真实存在的相关关系^[3], 或者是利用具有全局概念的概率转移矩阵方法消除网络中的间接相关关系^[4]. 然而这些方法都具有一定的缺陷性. 基因之间的相互作用往往具有局部的特性, 所以考虑全局特性会导致误差增大和结果不可信,并且这些方法都涉及矩阵求逆, 因而会有矩阵不可逆的问题和时间复杂度相对较高的问题存在. 这些缺陷在一定程度上限制了方法的应用.

本文在分析了基因相互作用性质之后,提出了从基因表达相关性数据中推断直接调控关系的模型. 模型假设成对基因的调控作用由他们之间的一阶传递、二阶传递和三阶传递作用的加权平均和决定,而认为更高阶传递作用的影响较为微弱且有产生间接相互作用的较高可能性,因而可以将其忽略不予考虑. 本文模型中还考虑了基因调控网络的稀疏性^[6],即一个基因通常只与一个或少数几个基因发生相互作用, 所以在模型中加入了控制基因调控网络稀疏性的调控因子,以更符合现有的生物现象. 将模型应用到文献[4]中 DREAM5 的大肠杆菌数据集中得到了较好的结果, 从而证明了模型能提高基因调控网络中边的预测能力.

2 方法

2.1 符号说明

基因之间调控网络矩阵分为观察矩阵 $G=[G_{ij}]$ 和直接矩阵 $S=[S_{ij}]$ . $G_{ij}$ 表示矩阵中节点 $i$ 和节点 $j$ 之间的相关性,包含直接相关性和间接相关性, $S_{ij}$ 表示节点 $i$ 和节点 $j$ 之间的直接相关性. ${\| }\cdot{\| }_{F}$ 表示矩阵的 $F$ 范数^[7], ${\| }\cdot{\| }_{l_{1}}$ 表示矩阵的 $l_{1}$ 范数^[8], $\lambda$ 为稀疏调节参数, $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ 为模型权重参数.

2.2 模型构建

给出了模型的理论框架图. 图中节点 1 和节点 2、节点 2 和节点 3 之间均有很强的相互作用,那么根据相关系数计算得到的观察矩阵 $G$ 中节点 1 和节点 3 有相关性的可能性就很大,所以 $G$ 中节点 1 和节点 3 之间有一条连边, 但是事实上他们之间没有直接相互作用. 在节点 2 和节点 3 之间的相互作用关系中, 2→3 为一阶相互作用,2→4→3 为二阶相互作用,2→4→5→3 为三阶相互作用. 我们认为这三种传递关系的影响并不一定是相等的,所以我们在模型中给不同的传递关系加上了权重因子 $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ . 消除观察矩阵 $G$ 中所有非直接的边就可以得到直接矩阵 $S$ . 本文模型假设成对基因的调控作用由他们之间的一阶传递、二阶传递和三阶传递作用的加权平均和决定,而认为更高阶传递作用的影响较为微弱且有产生间接相互作用的较高可能性,因而可以将其忽略不予考虑. 另外,模型利用直接矩阵 $S$ 的 $l_{1}$ 范数来控制网络的稀疏性, 并且加上了控制参数 $\lambda$ . 如果已知直接相关矩阵 $S$ , 那么基于马尔科夫链^[9]概率转移的性质,我们可以用直接矩阵 $S$ 的前三阶传递作用的加权平均和构造 $G$ 的近似模型. 相反的, 如果已知观察矩阵 $G$ ,我们可以设计一个迭代算法求出直接矩阵 $S$ 的最优值, $S$ 即为我们的得到的基因调控的真实网络.

图 1 网络净化模拟图

传递作用

$G=S+S^{2}+S^{3}+\cdots \cong \alpha_{1}S+\alpha_{2}S^{2} +\alpha_{3}S^{3}.$

模型

$\min\limits_{\alpha_{1},\alpha_{2},\alpha_{3},S} \bigg\| G-\sum\limits_{k=1}^{3}\alpha_{k}S^{k} \bigg\| _{F}^{2}+\lambda{\| S\| }_{l_{1}}.$

2.3 参数估计

因为矩阵 $S$ 中的元素和参数 $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ 均要求为非负的,所以我们在求解模型时增加了控制非负性的矩阵 $\Phi=[\Phi_{ij}]$ 和向量 $\beta=[\beta_{1} \quad \beta_{2} \quad \beta_{3}]$ ,其中 $\Phi_{ij}$ 控制 $S_{ij}\geqslant0$ , $\beta$ 用来限制权重向量 ${\alpha}$ 的元素都不小于 0. 加上控制项之后,模型的目标函数为

$\begin{equation} L(S,\Phi,\beta)=\min_{\alpha_{1},\alpha_{2},\alpha_{3},S}\bigg\| G -\sum_{k=1}^{3}\alpha_{k}S^{k}\bigg\|_{F}^{2}+\lambda{\| S\| }_{l_{1}}+{\Phi}S +\sum_{m=1}^{3}\beta_{m}\alpha_{m} \end{equation}$

(2.1)

模型计算时,矩阵 $S$ 的更新公式为

$\begin{equation} S\leftarrow\frac{1}{2}S+\frac{1}{2}S{\cdot}\frac{2(\alpha_{1}G+\alpha_{2} (S^{\rm T}G+G{S^{\rm T}})+\alpha_{3}(S^{\rm T}G{S^{\rm T}}+G(S^{2})^{{\rm T}} +(S^{2})^{{\rm T}}G))} {2(\alpha_{1}\tilde{G}+\alpha_{2}(S^{\rm T}\tilde{G} +\tilde{G}{S^{\rm T}})+\alpha_{3}(S^{\rm T}\tilde{G}{S^{\rm T}}+ \tilde{G}(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}\tilde{G}))+\lambda}, \end{equation}$

(2.2)

其中 $\tilde{G}=\sum\limits_{k=1}^{3}\alpha_{k}S^{k}$ . 在公式的右边,矩阵 $S$ 与后一项的乘积是点乘,表示矩阵对应元素相乘,最后一项中除法为点除, 表示矩阵对应元素相除. $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ 的迭代更新公式为

$\begin{equation} \alpha_{i}\leftarrow\frac{1}{2}\alpha_{i}+\frac{1}{2}\alpha_{i} \frac{{\rm tr}((S^{i})^{\rm T}G)}{{\rm tr}((S^{i})^{\rm T}\tilde{G})} ,\qquad i=1,2,3. \end{equation}$

(2.3)

2.4 算法分析

2.4.1 算法流程

(1) 确定常数 $\lambda$ 取值,初始化变量 $S=G$ , $\alpha={\alpha}_{0}$ ;

(2) 第 $i$ 步迭代,计算 $L_{i}=\bigg\| G-\sum\limits_{k=1}^{3}{\alpha}_{k}S_{i}^{k}\bigg\|_{F}^{2} +\lambda{\| S_{i}\| }_{l_{1}}$ ;

(3) 第 $i+1$ 步,根据迭代更新公式更新 $S_{i}$ 为 $S_{i+1}$ , $\alpha_{i}$ 为 $\alpha_{i+1}$ 直到达到停止准则;

(4) 计算 $L_{i+1}= \bigg\| G-\sum\limits_{k=1}^{3}\alpha_{k}S_{i+1}^{k}\bigg\|_{F}^{2} +\lambda{\| S_{i+1}\| }_{l_{1}}$ ;

(5) 计算 $odd=\frac{|L_{i}-L_{i+1}|}{L_{i}}$ ,直到达到停止准则.

为了提高运算效率又不丧失结果精度,我们在用迭代更新公式更新矩阵 $S$ 和向量 $\alpha$ 时,限制最大迭代次数为 50,容忍阈值为 0.01. 当累计迭代次数达 50 次或连续两次目标函数的变化率小于 0.01 时,迭代就停止.

2.4.2 算法推导

步骤 1 计算目标函数关于矩阵 $S$ 的梯度

$\begin{eqnarray} \frac{\partial(L(S,\Phi,\beta))}{\partial(S)} &=&2[\alpha_{1}\tilde{G}+\alpha_{2}(S^{\rm T}\tilde{G}+\tilde{G}{S^{\rm T}}) +\alpha_{3}(S^{\rm T}\tilde{G}{S^{\rm T}}+\tilde{G}(S^{2})^{{\rm T}} +(S^{2})^{{\rm T}}\tilde{G})\nonumber\\ &&-(\alpha_{1}G+\alpha_{2}(S^{\rm T}G+G{S^{\rm T}})+ \alpha_{3}(S^{\rm T}G{S^{\rm T}}+G(S^{2})^{{\rm T}})+(S^{2})^{{\rm T}}G)]+\lambda+\Phi. \nonumber\\ \end{eqnarray}$

(2.4)

根据 KKT 条件^[10], ${\Phi}S=0$ ,可以得到以下关于 $S$ 的方程

$\begin{eqnarray} &&2[\alpha_{1}G+\alpha_{2}(S^{\rm T}G+G{S^{\rm T}})+\alpha_{3}(S^{\rm T}G{S^{\rm T}}+G(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}G)]S \nonumber\\ &=&[2(\alpha_{1}\tilde{G}+\alpha_{2}(S^{\rm T}\tilde{G}+\tilde{G}{S^{\rm T}})+\alpha_{3}(S^{\rm T}\tilde{G}{S^{\rm T}} +\tilde{G}(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}\tilde{G}))+\lambda]S. \end{eqnarray}$

(2.5)

然后很容易得到如下的更新准则^[11]

$\begin{equation} S{\leftarrow}S\cdot\frac{2(\alpha_{1}G+\alpha_{2}(S^{\rm T}G +G{S^{\rm T}})+\alpha_{3}(S^{\rm T}G{S^{\rm T}}+G(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}G))}{2(\alpha_{1}\tilde{G}+\alpha_{2}(S^{\rm T}\tilde{G}+\tilde{G}{S^{\rm T}})+\alpha_{3}(S^{\rm T}\tilde{G}{S^{\rm T}}+\tilde{G}(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}\tilde{G}))+\lambda}. \end{equation}$

(2.6)

在实际操作中,变换为以下的更新准则计算速度更快

$\begin{equation} S\leftarrow{\frac{1}{2}}S+\frac{1}{2}S{\cdot}\frac{2(\alpha_{1}G+\alpha_{2}(S^{\rm T}G+G{S^{\rm T}})+\alpha_{3}(S^{\rm T}G{S^{\rm T}}+G(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}G))} {2(\alpha_{1}\tilde{G}+\alpha_{2}(S^{\rm T}\tilde{G}+\tilde{G}{S^{\rm T}})+\alpha_{3}(S^{\rm T}\tilde{G}{S^{\rm T}}+\tilde{G}(S^{2})^{{\rm T}}+(S^{2})^{{\rm T}}\tilde{G}))+\lambda}. \end{equation}$

(2.7)

步骤 2 计算目标函数关于 $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ 的梯度

$\begin{equation} \frac{\partial(L(S,\Phi,\beta))}{\partial(\alpha_{i})}= -2{\rm tr}[(S^{i})^{{\rm T}}(G-\tilde{G})]+\beta_{i},\qquad i=1,2,3. \end{equation}$

(2.8)

同样可以得到 $\alpha_{1}\mbox{、}\ \alpha_{2}\mbox{、}\ \alpha_{3}$ 的迭代更新公式

$\begin{equation} \alpha_{i}\leftarrow\frac{1}{2}\alpha_{i}+\frac{1}{2}\alpha_{i} \frac{{\rm tr}((S^{i})^{\rm T}G)}{{\rm tr}((S^{i})^{\rm T}\tilde{G})},\qquad i=1,2,3. \end{equation}$

(2.9)

3 结果分析

3.1 数据

我们将模型应用到文献[]中 DREAM5 的大肠杆菌数据集中.大肠杆菌数据集包含有 4511 个基因在 805 组不同的试验条件下的表达水平,并且给出了 141 个转录因子. 在建模分析时为了和 DREAM5 的协议保持一致,我们并不必要构建全部 $4511\times 4511$ 的相关性矩阵,只需要构建 141 个转录因子和 4511 个目标基因之间的 $141\times 4511$ 的相关性矩阵. 最后我们将其余所有节点的对角线元素设定为 1 得到方阵 $G$ . 我们构建观察矩阵 $G$ ,分别基于在网络边的预测中常用的三种方法: (1) Pearson correlations^[12]; (2) Spearman rank correlations^[13]; (3) Mutual information^[14]. 模型中不区分正相关性和负相关性, 即矩阵 $G$ 中的元素取为其对应基因表达量的相关系数的绝对值.

(1) Pearson correlations : 我们定义 $M=[M_{ni}]$ , $M_{ni}$ 代表第 $i=1,\cdots ,4511$ 个基因在第 $n=1,\cdots ,805$ 组试验条件下的表达水平, $M_{i}$ 为表达矩阵 $M$ 的第 $i$ 列. 对每一个转录因子 $j$ , 我们可以计算其与目标基因 $i$ 之间的 Pearson correlations

$\begin{equation} G_{ij}=\frac{E(M_{i}-EM_{i})(M_{j}-EM_{j})}{\sqrt{D(M_{i})}\sqrt{D(M_{j})}}. \end{equation}$

(3.1)

(2) Spearman rank correlation: 我们将表达水平向量 $M_{i}$ 和 $M_{j}$ 按次序编号, 转化等级变量 $m_{i}$ 和 $m_{j}$ ,那么转录因子与目标基因之间的 Spearman rank correlations 就为等级变量 $m_{i}$ 和 $m_{j}$ 的 Pearson correlations.

(3) Mutual information: $M_{i}$ 和 $M_{j}$ 之间的互信息定义为

$\begin{equation} G_{ij}=\sum_{x_{i},x_{j}}P(x_{i},x_{j})\log\frac{P(x_{i},x_{j})}{P(x_{i})P(x_{j})}. \end{equation}$

(3.2)

3.2 评估方法

我们应用模型到观察矩阵 $G$ 得到对应的直接矩阵 $S$ ,然后比较观察矩阵 $G$ 和相应直接矩阵 $S$ 的预测效果. 为了证实我们的预测,我们的比对标准是 DREAM5 中给出的基因相互作用标准对照表,它包括 2066 个已经被验证确定存在的基因调控相互作用对. 网络中 141 个转录因子和 4511 个目标基因之间总共有 $141\times 511=636,$ 051条可能的边,将这些边按概率大小降序排列并取其中前 $n$ ( $n$ 通常取值 100,000)条边(参见文献[]),得到 $L_{G}$ 和 $L_{S}$ , 分别代表矩阵 $G$ 和 $S$ 对网络中边的预测. 为了定量评估网络构建的准确性,我们使用真阳率( $TPR$ ) 和假阳率( $FPR$ ) 作为评价指标. 为此首先需要计算以下 2 个值.

(1) 真阳性 $TP$ (true positive): 构建的网络中正确的边数.

(2) 假阳性 $FP$ (false positive): 构建的网络中错误的边数.

AUROC曲线^[15]: 我们定义真阳率为

$\begin{equation} TPR(n)=\frac{TP(n)}{P}, \end{equation}$

(3.3)

其中 $TP(n)$ 代表 $L_{G}$ 和 $L_{S}$ 中预测正确的边数, $P$ 是标准比对表中正确的边数. 类似的,定义假阳率为

$\begin{equation} FPR(n)=\frac{FP(n)}{Q}, \end{equation}$

(3.4)

其中 $FP(n)$ 代表 $L_{G}$ 和 $L_{S}$ 中预测错误的边数, $Q$ 是标准比对表中错误的边数. AUROC 值即为 AUROC 曲线与 $X$ 轴围成区域的面积. AUROC 值越高说明预测效果越好.

3.3 正则化参数

$\lambda$ 和初值

$\alpha_{0}$ 的确定

为了得到模型的全局最优解,避免结果陷入局部最小值,我们用两种方法来求解模型, 然后比较这两种方法的结果,选取最高的 AUROC 值对应的参数 $\alpha$ 和 $\lambda$ 为模型的最优参数值. 第一种方法是固定 $\lambda$ 的值,分别取 $\lambda$ 为 $0\mbox{、}\ 0.0001\mbox{、}\ 0.001\mbox{、}\ 0.005\mbox{、}\ 0.01\mbox{、}\ 0.05\mbox{、}\ 0.1\mbox{、}\ 0.5\mbox{、}$ 10 这 9 个值, 取初值 $\alpha_{0} =[0.3333~ 0.3333~ 0.3333].$ 然后求出 9 组相对应的 $\alpha$ 和 $\lambda$ 的值,最后取最高的 AUROC 值对应的参数 $\alpha$ 和 $\lambda$ 为模型的最优参数值. 第二种方法是分别固定 $\alpha$ 为 $[0.95~ 0.05~ 0]\mbox{、}\ [0.85~ 0.1~ 0.05]\mbox{、}\ [0.75~ 0.25~ 0]\mbox{、}\ [0.75~ 0.15~ 0.1]\mbox{、}\ [0.65~ 0.35~ 0]\mbox{、}\ [0.65~ 0.25~ 0.1]\mbox{、}$ $[0.55~ 0.45~ 0]\mbox{、}\ [0.5~ 0.5~ 0]$ 这 8 个向量,然后求出最高的 AUROC 值对应的 $\alpha$ 和 $\lambda$ 为模型的最优参数值.

3.4 实验结果分析

通过计算我们发现,皮尔逊相关系数提高了 $\frac{0.5995-0.5935}{0.5935-0.5}=6.42\%$ , 斯皮尔曼相关系数提高了 $\frac{0.5895-0.5845}{0.584-0.5}=5.92\%$ ,互信息提高了 $\frac{0.6064-0.5973}{0.5973-0.5}=9.35\%$ . 如下表所示.

表1 观察矩阵

$G$ 和直接矩阵

$S$ 在三种相关系数下求得的 AUROC 值

表1 观察矩阵

$G$ 和直接矩阵

$S$ 在三种相关系数下求得的 AUROC 值

对 Pearson 相关系数,第一种方法的计算结果如下表.

表2 Pearson 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

表2 Pearson 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

从中我们发现第一种方法得到的结果都不理想,于是再按第二种方法计算. 实验结果发现当 $\lambda =0,$ $\alpha =[0.95~ 0.05~ 0]$ 时 AUROC 的值达到最高 0.5995.

对 Spearman 相关系数,第一种方法的计算结果如下表.

表3 Spearman 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

表3 Spearman 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

并且按第二种方法将固定 $\alpha$ 取值之后,得到的结果都没有上述表格中当 $\lambda$ 的值取 0.005、 $\alpha =[0.96~ 0.03~ 0.01]$ 时的结果好.

对 Mutual 相关系数,第一种方法的计算结果如下表.

表4 Mutual 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

表4 Mutual 相关系数直接矩阵在各

$\lambda$ 初值条件下得到的 AUROC 值

从中我们发现第一种方法的最优值为 $\lambda =0.001$ , ${\rm AUROC}=0.6059.$ 我们再按第二种方法进行计算. 实验结果发现当 $\lambda =0.0001$ , $\alpha =[0.55~ 0.45~ 0]$ 时 AUROC 的值达到最高为 0.6064.

4 小结

我们证明了本文的算法对基因调控网络是有效的,其时间复杂度^[16]远小于 $n^{3}$ , 并且不需要对网络的拓扑结构作任何假设,也没有用到任何专业领域的知识, 所以这个方法具有一定的普适性. 算法通过消除网络中间接作用的影响, 把原始的相关性数据转化为更能反映基因之间相互作用的直接矩阵 $S$ . 因此,本文中的算法可以提高各类观测网络的预测能力,消除传递作用产生的间接影响力. 此外,本文提出的算法有几个优于其他网络推断方法的特点. 第一点,我们不要求模型中矩阵 $G$ 是可逆的,而且在对模型的目标函数求梯度的过程中没有使用任何近似处理,计算结果非常可信. 第二点,我们在模型中添加了优化网络 $S$ 的稀疏性参数,使得迭代得到的直接网络 $S$ 更符合生物网络现象, 即一个基因通常只和少数几个基因之间有相互作用关系. 本文的模型有一个主要缺陷: 算法的目标函数不是凸函数^[17],这样得到的解经常是局部极小解. 为了避免迭代结果是一个局部最小解,建议多次重复整个计算,然后选择目标函数值最小的结果.

参考文献

[1]	方锦清,汪小帆,郑志刚,等. 一门崭新的交叉科学:网络科学(上). 物理学进展,2007, 27(3):239-343
[2]	王沛, 吕金虎. 基因调控网络的控制:机遇与挑战. 自动化学报,2013, 39(12):1969-1979
[3]	Feizi S, Marbach D, Médard M, Kellis M. Network deconvolution as a general method to distinguish direct dependencies in networks. Nature Biotechnology, 2013, 31:726-733
[4]	Barzel B, Barab\'asi A L. Network link prediction by global silencing of indirect correlations. Nat Biotechnol, 2013, 31:720-725
[5]	Marbach D, Costello J, Kuffner R, et al. Wisdom of crowds for robust gene network inference. Nature Methods, 2012, 9:796-804
[6]	朱陈平,张永梅,刘小廷,等. 复杂网络稀疏性的统计物理研究综述. 上海理工大学学报,2011, 33(5):425-432
[7]	李华云. F范数及矩阵分解实例研究. 现代情报,2008, 28(10):223-225
[8]	胡正平,路亮,许成谦. 基于L1范数稀疏距离测度学习的单类分类算法. 电子学报, 2012, 40(1):134-140
[9]	夏乐天,朱元甡,沈永梅. 加权马尔可夫链在降水状况预测中的应用. 水利水电科技进展, 2007, 26(6):20-23
[10]	文波,单甘霖,段修生. 基于KKT条件与壳向量的增量学习算法研究. 计算机科学, 2013, 40(3):255-258
[11]	Lee D D, Seung H S. Algorithms for non-negative matrix factorization. Adv Neural Inf Process Syst, 2001, 13:556-562
[12]	张宇镭,党琰,贺平安. 利用Pearson相关系数定量分析生物亲缘关系. 计算机工程与应用, 2006, 41(33):79-82
[13]	杨遵庆. 等级相关系数方法的应用. 北京商学院学报(社会科学版), 1985,(2):31-39
[14]	丁晶,王文圣,赵永龙. 以互信息为基础的广义相关系数. 四川大学学报(工程科学版), 2002, 34(3):1-5
[15]	宋花玲. ROC曲线的评价研究及应用[D]. 上海:第二军医大学, 2006
[16]	刘怀愚,朱昌杰,李王景. 时间复杂度的几种计算方法. 电脑知识与技术, 2011, 7(7):4636-4638
[17]	刘建国,葛仁东,夏尊铨,郭强. 非凸函数极小问题的BFGS算法. 运筹与管理,2004, 13(2):62-65