基因转录调控网络推导方法——伪逆矩阵模型
在基因调控网络推导中，使用到的
基因芯片数据通常具有
样本个数少（通常小于10）而
基因数目大（通常大于1000）的局限性，也就是说，实验样本个数远远小于基因的个数。另外，调控矩阵具有较强的
稀疏性，即每个基因只被少量的转录因子调控，而每个转录因子只调控少量的基因。
伪逆矩阵法就是为了解决基因调控网络线性模型中基因个数远远大于样本个数的问题而发展的一种有效的基因调控网络推导方法。Andrecut和Kauffman给出了推导基因调控网络的伪逆矩阵模型，并用贪婪算法进行了计算，得到了较好的效果。
在建模过程中，考虑含有N个基因在M次干扰实验中的基因调控网络，每一次干扰实验可以用一个N维动力系统建模
d
x
n
m
d
t
=
∑
i
=
1
N
w
n
i
x
i
m
−
μ
n
x
n
m
+
u
n
m
,
n
=
1
,
2
,
.
.
.
,
N
m
=
1
,
2
,
.
.
.
,
M
(
1
)
\frac{\mathrm{d} x_{nm}}{\mathrm{d} t} = \sum_{i=1}^{N}w_{ni}x_{im}-\mu _nx_{nm+u_{nm}}, n=1,2,...,N m=1,2,...,M (1)
dtdxnm=i=1∑Nwnixim−μnxnm+unm,n=1,2,...,Nm=1,2,...,M(1)
式中，
x
n
m
x_{nm}
xnm是基因n在实验m中的转录活性，即mRNA的浓度；
w
n
i
w_{ni}
wni是基因i对基因n的调控强度，
w
n
i
>
0
w_{ni}>0
wni>0表示激活，
w
n
i
<
0
w_{ni}<0
wni<0表示抑制；
μ
n
\mu _n
μn是由基因n转录生成的mRNA的降解率，
u
n
m
u_{nm}
unm是基因n在实验m中的噪声。
如果系统具有稳定性或处于平衡态，则稳态条件为：
d
x
n
m
d
t
=
0
⇔
∑
i
=
1
N
w
n
i
x
i
m
−
μ
n
x
n
m
+
u
n
m
=
0
(
2
)
\frac{\mathrm{d} x_{nm}}{\mathrm{d} t} = 0 \Leftrightarrow \sum_{i=1}^{N}w_{ni}x_{im}-\mu _nx_{nm}+u_{nm}=0 (2)
dtdxnm=0⇔i=1∑Nwnixim−μnxnm+unm=0(2)
设
X
=
[
x
n
m
]
N
×
M
X=[x_{nm}]_{N\times M}
X=[xnm]N×M为N×M阶矩阵，
G
=
[
w
n
n
]
N
×
N
−
d
i
a
g
(
μ
1
,
μ
2
,
.
.
.
,
μ
N
)
G=[w_{nn}]_{N\times N}-diag(\mu _1, \mu _2, ..., \mu _N)
G=[wnn]N×N−diag(μ1,μ2,...,μN)表示基因相互作用的N×N矩阵，
U
=
[
u
n
m
]
N
×
M
U=[u_{nm}]_{N\times M}
U=[unm]N×M为与转录干扰有关的N×M阶矩阵，则稳态条件等价于一个线性方程系统：
G
X
+
U
=
0
（
3
）
GX+U=0 （3）
GX+U=0（3）
对方程取转置，可得
X
T
G
T
+
U
T
=
0
(
4
)
X^T G^T+U^T = 0 (4)
XTGT+UT=0(4)
这等价于N个线性方程：
X
T
g
i
+
u
j
=
0
,
j
=
1
,
2
,
.
.
.
,
N
(
5
)
X^T g_i + u_j = 0, j= 1,2, ..., N(5)
XTgi+uj=0,j=1,2,...,N(5)
式中，
g
i
g_i
gi和
u
j
u_j
uj分别是矩阵
X
T
X^T
XT、
U
T
U^T
UT的第j列。既然M＜N，则该方程是欠定的。
一般方法是考虑等价的
l
2
l_2
l2优化问题：
g
j
=
a
r
g
min
g
j
∣
∣
g
j
∣
∣
2
(
5
)
g_j = arg\min_{g_j} ||g_j||_2 (5)
gj=arggjmin∣∣gj∣∣2(5)
s
.
t
.
X
T
g
j
+
u
j
=
0
s.t.X^Tg_j+u_j=0
s.t.XTgj+uj=0
式中，
∣
∣
g
j
∣
∣
2
||g_j||_2
∣∣gj∣∣2为
g
j
g_j
gj的欧几里得范数。这种情况可以通过最小化
l
2
l_2
l2范数
∣
∣
g
j
∣
∣
2
||g_j||_2
∣∣gj∣∣2得到唯一解
g
j
=
−
(
X
T
)
†
u
j
(
6
)
g_j=-(X^T)^\dagger u_j (6)
gj=−(XT)†uj(6)
式中，
(
X
T
)
†
=
X
(
X
T
X
)
−
1
(X^T)^\dagger =X(X^TX)^{-1}
(XT)†=X(XTX)−1为
X
T
X^T
XT的伪逆。但是所有解
g
j
g_j
gj的参数都是非零的，所以结果将违背基因调控网络的稀疏性。
为了量化稀疏性，考虑
l
0
l_0
l0范数
∣
∣
g
j
∣
∣
0
||g_j||_0
∣∣gj∣∣0，它能够方便地计算向量
g
j
g_j
gj中非零项的个数，所以，为了达到最稀疏的要求需要解
l
0
l_0
l0优化问题：
g
j
=
a
r
g
min
g
j
∣
∣
g
j
∣
∣
0
(
7
)
g_j = arg\min_{g_j} ||g_j||_0 (7)
gj=arggjmin∣∣gj∣∣0(7)
s
.
t
.
X
T
g
j
+
u
j
=
0
s.t.X^Tg_j+u_j=0
s.t.XTgj+uj=0
Andrecutd等用贪婪迭代试探法结合稀疏模型对微阵列芯片数据进行分析，数值实验结果表明该算法是简单有效的，能够对稀疏模型问题给出比较好的近似解，满足了基因调控网络的稀疏性要求。
补充：欠定方程与超定方程
-
欠定方程组：方程个数小于未知量个数的方程组
对于方程组 R a = y Ra = y Ra=y，R为n×m的矩阵，且n＜m。则方程组有无穷多组解，此时称方程组为欠定方程组。
内点法和梯度投影法是目前解欠定方程组的常用方法 -
超定方程组：方程个数大于未知量个数的方程组。
对于方程组 R a = y Ra = y Ra=y，R为n×m矩阵，如果R列满秩，且n＞m（n为方程组行数，m未未知量个数）
超定方程一般是不存在解的矛盾方程。例如，如果给定的三点不在一条直线上，我们将无法得到这样一条直线，使得这条直线同时经过给定这三个点。也就是说给定的条件（限制）过于严格，导致解不存在。
在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。形象地说，就是在无法完全满足给定的这些条件的情况下，求一个最接近的解。