【转录调控网络】典型的基因转录调控网络推导方法——奇异值分解

╰+哭是因爲堅強的太久メ 2022-09-15 14:55 115阅读 0赞

### 基因转录调控网络推导方法——奇异值分解 ###

`奇异值分解`是线性代数中的一种矩阵分解方法，在信号处理和统计学等领域应用广泛。目前，奇异值分解已被广泛应用在大规模基因表达数据分析。例如Yeung等利用奇异值分解法在稀疏条件下重建了基因调控网络，实验证明该方法具有较高的准确性，而且具有较低的计算复杂度。

为了减少复杂度，这种方法**只考虑了稳态情况下的调控关系**，并用一个`线性微分方程组`(方程1)来表示：  
 d x i ( t ) d t = − λ i x i ( t ) + ∑ j = 1 n W i j x i ( t ) + b ( t ) + ξ i ( t ) , i = 1 , 2 , . . . n \\frac\{\\mathrm\{d\} x\_i(t)\}\{\\mathrm\{d\} t\} = -\\lambda\_i x\_i(t)+\\sum\\nolimits\_\{j=1\}^\{n\}\{W\_\{ij\} x\_i(t)\}+b(t)+\\xi \_i(t), i=1,2,...n dtdxi(t)=−λixi(t)\+∑j=1nWijxi(t)\+b(t)\+ξi(t),i=1,2,...n  
上式中， x i ( t ) ( i = 1 , 2 , . . . , n ) x\_i(t)(i=1,2,...,n) xi(t)(i=1,2,...,n)反应了基因 i i i的表达水平； λ i \\lambda\_i λi是基因 i i i的自我降解率； b ( t ) b(t) b(t)是外部因素; ξ i ( t ) \\xi\_i(t) ξi(t)是噪声； W i j W\_\{ij\} Wij是调控矩阵元素，表示基因 j j j对基因 i i i的调控类型和调控强度，正数表示激活，负数表示降解，零表示没有调控关系。

在外部刺激条件 ( b 1 , b 2 , . . . , b n ) T (b\_1,b\_2,...,b\_n)^T (b1,b2,...,bn)T干扰下，n个不同基因的表达量分别记为 ( x 1 , x 2 , . . . , x n ) T (x\_1,x\_2,...,x\_n)^T (x1,x2,...,xn)T：  
 X = ( x 1 ( 1 ) x 1 ( 2 ) . . . x 1 ( m ) x 2 ( 1 ) x 2 ( 2 ) . . . x 2 ( m ) . . . . . . . . . . . . x n ( 1 ) x n ( 2 ) . . . x n ( m ) ) X = \\begin\{pmatrix\} x\_1(1) &x\_1(2) &... &x\_1(m) \\\\ x\_2(1) &x\_2(2) &... &x\_2(m) \\\\ ... &... &... &... \\\\ x\_n(1) &x\_n(2) &... &x\_n(m) \\end\{pmatrix\} X=⎝⎜⎜⎛x1(1)x2(1)...xn(1)x1(2)x2(2)...xn(2)............x1(m)x2(m)...xn(m)⎠⎟⎟⎞  
上式中， x i ( j ) x\_i(j) xi(j)表示在第 j j j个样本（干扰）中第 i i i个基因的表达。

将(方程1)写成矩阵的形式(方程2)：  
 X ˇ = J X + B \\check\{X\} = JX+B Xˇ=JX\+B  
这里忽略了噪声，自我降解率 ξ i ( t ) \\xi\_i(t) ξi(t)被 J i j J\_\{ij\} Jij合并用来简化公式。 X ˇ = d X ( t ) d t \\check\{X\} = \\frac\{\\mathrm\{d\} X(t)\}\{\\mathrm\{d\} t\} Xˇ=dtdX(t)， B = \[ b i ( j ) \] , i = 1 , 2 , . . . , n , j = 1 , 2 , . . . , m , J i j = W i j − δ i j λ i B=\[b\_i(j)\], i=1,2,...,n,j=1,2,...,m,J\_\{ij\}=W\_\{ij\}-\\delta \_\{ij\}\\lambda \_i B=\[bi(j)\],i=1,2,...,n,j=1,2,...,m,Jij=Wij−δijλi为函数，即（方程3）：  
 δ i j = \{ 1 i = j 0 i ≠ j \\delta \_\{ij\}=\\begin\{cases\} 1 & i = j \\\\ 0 & i \\ne j \\end\{cases\} δij=\{ 10i=ji=j

通过奇异值分解，我们可以把 X T X^T XT分解为 X T = U E V T X^T=UEV^T XT=UEVT，其中 U U U是一个左特征向量的 m × n m\\times n m×n阶酉矩阵， E = d i a g ( e 1 , e 2 , . . . , e n ) E=diag(e\_1,e\_2,...,e\_n) E=diag(e1,e2,...,en)是包含n个特征值的n×n阶主对角线矩阵， V T V^T VT是右特征向量的n×n阶酉矩阵，我们可以得到方程的解(方程4)：  
 J ^ = ( X ˇ − B ) U E − 1 V T \\hat\{J\} = (\\check\{X\}-B ) UE^\{-1\}V^T J^=(Xˇ−B)UE−1VT  
式中， E − 1 = d i a g ( 1 / e i ) E^\{-1\}=diag(1/e\_i) E−1=diag(1/ei)，如果 e i = 0 e\_i=0 ei=0，则令 1 / e i = 0 1/e\_i=0 1/ei=0。所以可以得到方程的一般相互关系矩阵解 J = \[ J i j \] n × n J=\[J\_\{ij\}\]\_\{n\\times n\} J=\[Jij\]n×n，即（方程5）：  
 J = ( X ˇ − B ) U E − 1 V T + Y V T = J ^ + Y V T J=(\\check\{X\}-B)UE^\{-1\}V^T+YV^T=\\hat\{J\}+YV^T J=(Xˇ−B)UE−1VT\+YVT=J^\+YVT  
式中， Y = \[ y i j \] n × n Y=\[y\_\{ij\}\]\_\{n\\times n\} Y=\[yij\]n×n是一个n×n阶矩阵，若 e i ≠ 0 e\_i \\ne 0 ei=0，则 y i j = 0 y\_\{ij\}=0 yij=0；若 e i = 0 e\_i = 0 ei=0，则 y i j y\_\{ij\} yij为任意常数。每个J代表一个依赖于Y与芯片数据的网络。

由于奇异值分解会导致非唯一解，这就需要一定的约束从所有解之中分离出真正的解。由于真正的基因调控网络通常是稀疏的，就是说每个基因一般只与所有基因中的一小部分基因相互作用，因此稀疏性是合理的约束条件。Yeung等基于此点利用优化方法使J的 L 1 L\_1 L1值最小，即（方程6）：  
 min ⁡ Y ∣ J ^ + Y V T ∣ \\min\_\{Y\}|\\hat\{J\}+YV^T| Ymin∣J^\+YVT∣

Yeung等在文献中通过多个实验证明，奇异值分解方法推出的基因调控网络具有较高的准确性，但是需要注意的是，这种方法只适合基因调控网络的线性模型。

代码实现：  
numpy当中封装了现成的SVD分解方法。我们直接调用np.linalg.svd接口即能完成矩阵的分解：

import numpy as np
    from sklearn.datasets import load_breast_cancer
    cancer = load_breast_cancer()
    data, label = cancer.data, cancer.target
    
    U, Sigma, V = np.linalg.svd(data)

执行结果：  
![SVD][]  
这里的Sigma返回的是一个向量，代替了对角矩阵，节省了存储开销。可以**通过找出最小的K，使得K个奇异值占据整体奇异值95%以上的和**。这里可以看到，我们选出了5个奇异值就占据所有奇异值和的99%以上：

总结：  
SVD和PCA一样底层都是基于矩阵的线性操作完成的，通过SVD的性质，我们可以对原数据进行压缩和转化。基于这一点，衍生出了许多的算法和应用场景，其中最经典的要属**推荐系统中的协同过滤**了。

--------------------

**矩阵近似值**  
**奇异值分解在统计中的主要应用为主成分分析（PCA）**，它是一种数据分析方法，用来找出大量数据中所隐含的“模式”，它可以用在模式识别，数据压缩等方面。PCA算法的作用是把数据集映射到低维空间中去。  
数据集的特征值（在SVD中用奇异值表征）按照重要性排列，**降维的过程就是舍弃不重要的特征向量的过程，而剩下的特征向量张成空间为降维后的空间**。

**正交矩阵**  
**正交矩阵是实数特殊化的酉矩阵**，因此总是正规矩阵。尽管我们在这里只考虑实数矩阵，  
这个定义可用于其元素来自任何域的矩阵。正交矩阵毕竟是从内积自然引出的，对于复  
数的矩阵这导致了归一要求。注意正交矩阵的定义：n阶‘实矩阵’ A称为正交矩阵，如果：A×A′=E（E为单位矩阵，  
A’表示“矩阵A的转置矩阵”。） 若A为正交阵，则下列诸条件是等价的:

1.  A 是正交矩阵
2.  A×A′=E（E为单位矩阵）
3.  A′是正交矩阵
4.  A的各行是单位向量且两两正交
5.  A的各列是单位向量且两两正交
6.  (Ax,Ay)=(x,y) x,y∈R

[SVD]: /images/20220828/55c1778a7c9a4c2c881f0acf011326e2.png