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

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


奇异值分解是线性代数中的一种矩阵分解方法,在信号处理和统计学等领域应用广泛。目前,奇异值分解已被广泛应用在大规模基因表达数据分析。例如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)​=−λi​xi​(t)+∑j=1n​Wij​xi​(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​={10​i=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)

执行结果:
【转录调控网络】典型的基因转录调控网络推导方法——奇异值分解
这里的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
上一篇:矩阵的旋转


下一篇:多元最短路-Floyd算法