CN103268431A - 一种基于学生t分布的癌症亚型生物标志物检测系统 - Google Patents

Abstract

本发明涉及一种基于学生t分布的癌症亚型生物标志物检测系统，包括基因表达谱提取单元，从基因芯片中提取基因样本的表达谱；网络聚类模型建立单元，提取到的基因表达谱输入其中，建立混合模型并对其进行求解，其中混合模型使用多元学生t分布来描述每个成分；癌症亚型发现单元，利用混合模型和聚类方法发现新的癌症亚型；基因间相互关系探索单元，将每个基因看作是网络中的点，通过挖掘网络结构探索不同癌症亚型中基因间的相互关系；生物标志物检测单元，根据不同癌症亚型中基因间的相互关系完成与癌症亚型相关的生物标记物的检测。本发明使用多元学生t分布来描述混合模型中的每个成分，使得聚类、变量选择过程对异常点、噪声具有鲁棒性。

Description

一种基于学生t分布的癌症亚型生物标志物检测系统

技术领域

本发明涉及生物信息学、统计学习和模式识别领域，更具体地，涉及一种基于学生t分布的癌症亚型生物标志物检测系统。

背景技术

随着基因组表达谱数据的日益增加，芯片方法已经成为识别癌症相关基因以及发现新癌症亚型（子类）的关键技术。与病人的年龄、肿瘤大小、类固醇受体状态等临床和病理危险因素相比，对基因的研究能更深入地了解癌症的生理特性，并能更有效地侦测新的癌症亚型。这些癌症亚型可能在基因或蛋白质表达谱、基因调控或蛋白质信号网络等方面存在差异。从基因表达谱中检测这些癌症子类可被看作一个无监督的聚类问题，而寻找用于检测的基因可被看作一个基于高维无类标数据的变量选择问题。

当前，癌症亚型发现所面临的一个挑战是，这些癌症亚型间各种网络或途径水平的差异可能会使传统的基于基因表达谱差异的聚类方法变的不合适。侦测这些网络和途径对于理解基因的集体生物学功能以及它们对病人表型变化的影响有着非常重要的作用。除此之外，人们经常根据每个基因的判别能力，将它们独立的挑选出来。然而，基因通常需要和其它的基因相互作用，共同参与一些生物过程或分子功能。它们中有一些在基因表达谱上没有存在差异，但属于一个子网络，它具有总体的判别能力或是一个特定癌症亚型的有用途径。

现有很多聚类方法被用于分析基因表达数据，对生物样本进行划分。基于模型聚类方法（model-based clustering）具有固定的概率框架，其由于具有好的效果、解释性、并易于实施，而被广泛地应用于生物标志物和癌症亚型的检测。大部分方法通过对似然函数进行惩罚限制得到稀疏解，以达到基因选择的目的。现有较为常用的实现方案主要有以下几种：

（1）使用对角协方差的惩罚基于模型的高斯聚类系统，该系统所采取的算法假设每个聚类具有对角的协方差矩阵，因此基因间是独立的。虽然独立假设在真实生物数据中一般是不成立的，它仍能得到较好的聚类效果。模型中使用混合高斯分布来模拟每个聚类中的样本，并与一些惩罚函数相结合，例如L1惩罚、自适应性L1惩罚和成组L1惩罚。

（2）使用对角协方差的惩罚基于模型的学生t聚类模型系统，该系统借助自助法给出了根据基因对聚类贡献的基因排序方法。

（3）使用无限制协方差的惩罚基于模型的高斯聚类系统，其通过使用无限制协方差矩阵来考虑基因间的各种关系。利用图套索（graphical lasso），使用期望最大化（expectation maximization，EM）算法进行参数估计。其所运用的算法具有更好的癌症亚型和生物标志物检测性能。

（4）基于高斯混合模型的惩罚因子分析系统，其采用的算法介于对角和一般的协方差矩阵之间，类似于因子分析，该系统使用隐变量对协方差矩阵进行建模。与基于无限制协方差矩阵的系统相比，该系统的模型更为复杂，并具有更多的限制。然而，当基因表达谱数据中，协方差矩阵的隐变量诱导假设成立时，该模型具有更好的效果。

但上述系统中存在如下的某些缺点：

由于高斯分布的概率密度函数是基于样本与均值的偏差平方，以上基于高斯分布的系统对异常点较为敏感。而鉴于一些化学或电气原因，基因表达数据中经常包含异常点。另外，使用对角协方差假设的系统忽略了癌症亚型中基因间的相关性。并且，传统的惩罚基于模型聚类算法都以均值为选择基因的标准，而忽略了它们对网络和路径的影响，这对理解基因的集体生物学功能是非常不利的。

发明内容

本发明为克服上述现有技术所述的至少一种缺陷（不足），提供一种对噪声、异常点有良好鲁棒性的基于学生t分布的癌症亚型生物标志物检测系统。

为解决上述技术问题，本发明的技术方案如下：

一种基于学生t分布的癌症亚型生物标志物检测系统，包括：

基因表达谱提取单元，从基因芯片中提取n个独立的包含d个基因的样本的表达谱χ={x₁,…;x_n}，其中x_i=[x_i1 x_i2 … x_id]表示样本i中d个基因的表达谱；

网络聚类模型建立单元，与基因表达谱提取单元连接，提取到的基因表达谱χ={x₁,…,x_n}输入其中，建立混合模型并对该混合模型进行求解，其中，该混合模型使用多元学生t分布来描述混合模型中的每个成分；

癌症亚型发现单元，与网络聚类模型建立单元连接，其利用混合模型和聚类方法发现新的癌症亚型；

基因间相互关系探索单元，与癌症亚型发现单元连接，将每个基因看作是网络中的点，通过挖掘网络结构探索不同癌症亚型中基因间的相互关系；

生物标志物检测单元，根据不同癌症亚型中基因间的相互关系完成与癌症亚型相关的生物标记物的检测。

上述方案中，网络聚类模型建立单元中建立混合模型过程如下：

假设每个样本x_i来自于具有K个成分的混合分布，其概率密度函数如下：

$f(x_i;\mathrm{\Ψ})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\π}}_k{f}_k(x_i;{\mathrm{\θ}}_k),$

其中Ψ={π₁，…，π_k，θ₁，…，θ_k}包括混合模型中所有的参数，f_k表示第k个成分，π_k≥0是第k个成分的混合比例，

θ_k为对应于f_k的未知参数集合;

假设每个成分f_k的参数集合为θ_k={μ_k，Σ_k，v_k}的多元学生t分布T(μ_k，Σ_k，v_k)，其中μ_k=[μ_k1 μ_k2 … μ_kd]为位置参数矩阵，Σ_k=(σ_k，jl)_d×d为尺度参数，而v_k为自由度，多元学生t分布的概率密度函数为:

$f_k(x_i;{\mathrm{\θ}}_k)=\frac{\mathrm{\Γ}\left(\frac{v_k+d}{2}\right){\left|{\mathrm{\Σ}}_k\right|}^{-\frac{1}{2}}}{\mathrm{\Γ}\left(\frac{v_k}{2}\right){\left(\mathrm{\π}{v}_k\right)}^{d/2}{\{1+\mathrm{\δ}(x_i;{\mathrm{\μ}}_k,{\mathrm{\Σ}}_k)/v_k\}}^{(v_k+d)/2}},$

其中Γ(·)为Gamma函数， $\mathrm{\δ}(x_i;{\mathrm{\μ}}_k,{\mathrm{\Σ}}_k)={(x_i-{\mathrm{\μ}}_k)}^{\′}{\mathrm{\Σ}}_k^{-1}(x_i-{\mathrm{\μ}}_k)$ 表示x_i和μ_k间的马氏距离，每个学生t分布的均值和协方差矩阵分别为μ_k和

d表示每个样本包含的基因的个数，j和l是矩阵Σ_k中元素所在的位置，π为圆周率；

记尺度参数矩阵Σ_k的逆为W_k={W_k,jl}_d×d，对每个μ_kj以及W_k={W_k,jl}_d×d的非对角元素进行自适应的L₁惩罚;基于由对数似然函数l(Ψ)和惩罚部分Pen_λ(φ)组成的惩罚对数似然函数，得到模型的目标函数:

$l_{\mathrm{pen}}\left(\mathrm{\Ψ}\right)=l\left(\mathrm{\Ψ}\right)-{\mathrm{pen}}_{\mathrm{\λ}}\left(\text{\Φ}\right)$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\log \left[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\π}}_k{f}_k(x_i;{\mathrm{\θ}}_k)\right]-{\mathrm{\λ}}_1\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{d}{\mathrm{\Σ}}}{w}_{\mathrm{kj}}\left|{\mathrm{\μ}}_{\mathrm{kj}}\right|-{\mathrm{\λ}}_2\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{l=1}{\overset{d}{\mathrm{\Σ}}}\underset{j=1,j\≠l}{\overset{d}{\mathrm{\Σ}}}{v}_{k,\mathrm{jl}}\left|W_{k,\mathrm{jl}}\right|,$

其中，φ={μ₁，…，μ_K，W₁，…，WK}，λ={λ₁λ₂}包括了对应于μ_k和W_k的非负正则化参数λ₁和λ₂，其预先设定并通过外部输入，w_kj、v_k,jl分别为对应于μ_kj和W_k,jl的惩罚的权重。

上述方案中，网络聚类模型建立单元中对混合模型的求解过程如下：

给定K和λ，使用期望最大化算法对目标函数l_pen(Ψ)进行优化：假设每个样本x_i对应于一个无法观察到的指标向量z_i=[z_i1 z_i2 … z_iK]，指定x_i所属的成分；如果x_i来自于成分k，则z_ik=1，否则z_ik=0；引入另一个隐变量u_i=[u_i1 u_i2 … u_iK]，其中，每个u_ik服从Gamma分布；则算法中样本x_i完整的数据变为其对应的惩罚完整数据对数似然函数可被分解为四个部分：

l_c,pen(Ψ)=l₁(π)+l₂(v)+l₃(φ)-pen_λ(φ),

其中

$l_1\left(\mathrm{\π}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{z}_{\mathrm{ik}}\log {\mathrm{\π}}_k,$

$l_2\left(v\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{z}_{\mathrm{ik}}[-\log \mathrm{\Γ}\left(\frac{v_k}{2}\right)+\frac{v_k}{2}\log \left(\frac{v_k}{2}\right)+\frac{v_k}{2}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})-\log u_{\mathrm{ik}}],$

$l_3\left(\mathrm{\Φ\φ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[-\frac{d}{2}\log \left(2\mathrm{\π}\right)+\frac{d}{2}\log u_{\mathrm{ik}}+\frac{1}{2}\log |W_k|-\frac{1}{2}{u}_{\mathrm{ik}}{(x_i-{\mathrm{\μ}}_k)}^{\′}{W}_k(x_i-{\mathrm{\μ}}_k)].$

执行期望最大化算法期望步，在已知第t次迭代的参数估计Ψ^(t)的情况下，计算l_c,pen(Ψ)的期望Q(Ψ;Ψ^(t))，其中先求得E_Ψ(t)(z_ik|x_i)、E_Ψ(t)(u_ik|x_i,z_ik=1)和E_Ψ(t)(logu_ik|x_i,z_ik=1)；

$E_{\mathrm{\Ψ}\left(t\right)}\left(z_{\mathrm{ik}}\right|x_i)=\frac{{\mathrm{\π}}_k^{\left(t\right)}{f}_k(x_i;{\mathrm{\θ}}_k^{\left(t\right)})}{f(x_i;{\mathrm{\Ψ}}^{\left(t\right)})}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}.$

被看作为x_i属于第k个聚类的后验概率；

鉴于Gamma分布是高斯分布的共轭先验，有

$E_{\mathrm{\Ψ}\left(t\right)}\left(u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\frac{v_k^{\left(t\right)}+d}{v_k^{\left(t\right)}+\mathrm{\δ}(x_i;{\mathrm{\μ}}_k^{\left(t\right)},{\mathrm{\Σ}}_k^{\left(t\right)})}\stackrel{\mathrm{\Δ}}{=}{u}_{\mathrm{ik}}^{\left(t\right)},$

$E_{\mathrm{\Ψ}\left(t\right)}\left(\log u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\log u_{\mathrm{ik}}^{\left(t\right)}+\mathrm{\ψ}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right),$

其中 $\mathrm{\ψ}\left(s\right)=\{\∂\mathrm{\Γ}\left(s\right)/\∂s\}/\mathrm{\Γ}\left(s\right)$ 是Digamma函数；

执行期望最大化算法最大化步，得到各个参数的更新：

${\mathrm{\π}}_k^{(t+1)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}/n,k=1,\·\·\·,K.$

v_k在第t+1次迭代的估计满足等式：

$-\mathrm{\ψ}\left(\frac{v_k}{2}\right)+\log \left(\frac{v_k}{2}\right)+1+\frac{1}{n_k^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(k\right)}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})+\mathrm{\ψ}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right)=0,$

其中 $n_k^{\left(t\right)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)};$

在t+1步，定义自适应的权重为：

$w_{\mathrm{kj}}^{\left(t\right)}=\frac{1}{{|\mathrm{\μ}}_{\mathrm{kj}}^{\left(t\right)}|+\∈},v_{k,\mathrm{jl}}^{\left(t\right)}=\frac{1}{\left|W_{k,\mathrm{jl}}^{\left(t\right)}\right|+\∈}\text{.}$

其中，设定∈=0.1，并使用L₁惩罚的结果作为初始估计值和

μ_kj估计的更新方法为：如果

$\left|\frac{1}{w_{\mathrm{kj}}^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(\underset{l=1,l\≠j}{\overset{d}{\mathrm{\Σ}}}(x_{\mathrm{il}}-{\mathrm{\μ}}_{\mathrm{kl}}^{\left(t\right)})W_{k,\mathrm{jl}}^{\left(t\right)}+x_{\mathrm{ij}}{W}_{k,\mathrm{jj}}^{\left(t\right)})\right|\≤{\mathrm{\λ}}_1,$

则 ${\mathrm{\μ}}_{\mathrm{kj}}^{(t+1)}=0;$ 否则

$\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}\right)u_{\mathrm{kj}}^{(t+1)}{W}_{k,\mathrm{jj}}^{\left(t\right)}+{\mathrm{\λ}}_1{w}_{\mathrm{kj}}^{\left(t\right)}\mathrm{sign}\left({\mathrm{\μ}}_{\mathrm{kj}}^{(t+1)}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(x_i{W}_{k,.j}^{\left(t\right)}-{\mathrm{\μ}}_k^{\left(t\right)}{W}_{k,.j}^{\left(t\right)}-{\mathrm{\μ}}_{\mathrm{kj}}^{\left(t\right)}{W}_{k,\mathrm{jj}}^{\left(t\right)}).$

进一步，得到有关W_k的优化问题

$W_k^{(t+1)}=\arg \underset{W_k}{\max }\log \left[\det \left(W_k\right)\right]-\mathrm{tr}\left({\tilde{S}}_k{W}_k\right)-\underset{j\≠l}{\mathrm{\Σ}}{\mathrm{\λ}}_{k,\mathrm{jl}}^{\′}\left|W_{k,\mathrm{jl}}\right|,$

其中

${\tilde{S}}_k=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}{(x_i-{\mathrm{\μ}}_k^{\left(t\right)})}^{\′}(x_i-{\mathrm{\μ}}_k^{\left(t\right)})}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}},$ ${\mathrm{\λ}}_{k,\mathrm{jl}}^{\′}=\frac{{2\mathrm{\λ}}_2{v}_{k,\mathrm{jl}}^{\left(t\right)}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}}.$

该优化问题使用图套索求解。

上述方案中，网络聚类模型建立单元基于证据近似权重准则进行混合模型的选择，具体为：估计三个参数，聚类的个数K、惩罚参数λ₁和λ₂，使用证据近似权重准则进行模型选择，

AWE=-2l(Ψ)+2d_e(3/2+log(n))，

其中d_e=K-1+K+Kd-q₁+Ks+Kd(d-l)/2-q₂/2是混合模型中有效参数的个数，q₁=#{(k,j）：μ_kj=0)，q₂=#{(k,j,l)：W_k，lj=0}，通过网格搜索法寻找使AWE最小的最优组合(K，λ_λ，λ₂)，n为样本个数，d为每个样本的基因的个数，i(Ψ)为对数似然函数。

上述方案中，癌症亚型发现单元发现癌症亚型的过程如下:

通过网络聚类模型建立单元对混合模型中参数进行估计之后，聚类被定义为服从相似分布的样本，样本的分布由后验概率τ_ik决定;给定样本x_i，将其聚类为使后验概率最大的类，即y_i=arg max_k τ_ik;获取属于相同类的样本组成的集合判断为新的癌症亚型。

上述方案中，基因间相互关系探索单元探索不同癌症亚型中基因间的相互关系的过程为:

将不同癌症亚型的每个基因看作是网络中的点，由尺度矩阵的逆W_k来决定点与点间的边，若W_k，jl≠0，则对应聚类k的网络中，j和l相连，否则j和l不相连;该网络中单连通的点集以及它们之间的边构成一系列的子网络;在每个聚类中不同的子网络表示不同癌症亚型中基因间的关系。

上述方案中，生物标志物检测单元检测癌症亚型相关的生物标记物的过程如下:

从基因间相互关系探索单元中获得一系列的子网络;

由于噪声基因和有信息基因间是不相关的，如果一个子网络中存在一个基于均值的有判别力基因，判断它是一个子网络生物标志物，从而识别出在基因表达谱上不存在差异的基因，所述不存在差异的基因与有判别能力的基因相关，共同完成一项生物功能；

最后，在剩下的子网络中，若其在各个聚类对应的W₁，…，W_K上存在差异，则判断该子网络为生物标志物；

将检测到的生物标志物数据输出。

与现有技术相比，本发明技术方案的有益效果是：

（1）利用本发明提供的系统可以检测到与癌症亚型相关的生物标记物，简单方便。

（2）本发明所提供的系统使用多元学生t分布来描述混合模型中的每个成分，使得系统不会受极端或异常基因影响，从而使得聚类、变量选择过程对异常点具有鲁棒性，能更好的处理基因表达谱中包含的噪音，利用本发明的系统能够发现真实的癌症亚型以及检测在生物化学或生物医学研究中得到验证的与癌症亚型相关生物标志物。

（3）本发明所提供的系统使用无限制协方差矩阵来考虑基因间的关系,这更有利于仿真基因间存在的相互联系，通过对均值和协方差的逆进行自适应的L1惩罚来实现生物标志物检测以及基因间的网络重构，并进一步借助图套索，使用期望最大化算法来拟合模型。

（4）本发明所提供的系统提供了基于网络的基因选择准则，其能够识别出具有较低判别能力，但和有判别能力的基因相关或具有聚类特定的网络结构，从而能够更加准确地识别出与癌症亚型相关的生物标记物。

附图说明

图1为本发明的一种基于学生t分布的癌症亚型生物标志物检测系统实施例的架构图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；

对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

如图1所示，为本发明的一种基于学生t分布的癌症亚型生物标志物检测系统实施例的架构图。参见图1，本实施例的一种基于学生t分布的癌症亚型生物标志物检测系统，包括：

基因表达谱提取单元100，从基因芯片中提取n个独立的包含d个基因的样本的表达谱X={x₁,…，x_n}，其中x_i=[x_i1 x_i2 … x_id]表示样本i中d个基因的表达谱；每个基因均预先经标准化处理，样本均值为0，样本方差为1。

网络聚类模型建立单元101，与基因表达谱提取单元100连接，提取到的基因表达谱X={x₁，…，x_n}输入其中，建立混合模型并对该混合模型进行求解，其中，该混合模型使用多元学生t分布来描述混合模型中的每个成分；

癌症亚型发现单元102，与网络聚类模型建立单元101连接，其利用混合模型和聚类方法发现新的癌症亚型；

基因间相互关系探索单元103，与癌症亚型发现单元102连接，将每个基因看作是网络中的点，通过挖掘网络结构探索不同癌症亚型中基因间的相互关系；

生物标志物检测单元104，与基因间相互关系探索单元103连接，根据不同癌症亚型中基因间的相互关系完成与癌症亚型相关的生物标记物的检测。

其中，网络聚类模型建立单元101中建立混合模型过程如下：

假设每个样本x_i来自于具有K个成分的混合分布，其概率密度函数如下：

$f(x_i;\mathrm{\Ψ})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\π}}_k{f}_k(x_i;{\mathrm{\θ}}_k),,---\left(1\right)$

其中Ψ={π₁,…,π_k,θ₁,…,θ_k}包括混合模型中所有的参数，f_k表示第k个成分，π_k≥0是第k个成分的混合比例，

θ_k为对应于f_k的未知参数集合；

假设每个成分f_k的参数集合为θ_k={μ_k,Σ_k,v_k}的多元学生t分布T(μ_k,Σ_k,v_k)，其中μ_k=[μ_k1 μ_k2 … μ_kd]为位置参数矩阵，Σ_k=(σ_k,jl)_d×d为尺度参数，而v_k为自由度，多元学生t分布的概率密度函数为：

$f_k(x_i;{\mathrm{\θ}}_k)=\frac{\mathrm{\Γ}\left(\frac{v_k+d}{2}\right){\left|{\mathrm{\Σ}}_k\right|}^{-\frac{1}{2}}}{\mathrm{\Γ}\left(\frac{v_k}{2}\right){\left(\mathrm{\π}{v}_k\right)}^{d/2}{\{1+\mathrm{\δ}(x_i;{\mathrm{\μ}}_k,{\mathrm{\Σ}}_k)/v_k\}}^{(v_k+d)/2}},---\left(2\right)$

其中Γ(·)为Gamma函数， $\mathrm{\δ}(x_i;{\mathrm{\μ}}_k,{\mathrm{\Σ}}_k)={(x_i-{\mathrm{\μ}}_k)}^{\′}{\mathrm{\Σ}}_k^{-1}(x_i-{\mathrm{\μ}}_k)$ 表示x_i和μ_k间的马氏距离，每个学生t分布的均值和协方差矩阵分别为μ_k和

d表示每个样本包含的基因的个数，j和l是矩阵Σ_k中元素所在的位置，π为圆周率；此处系统的网络聚类模型建立单元101使用学生t分布模拟数据；

记尺度参数矩阵Σ_k的逆为W_k={W_k,jl}_d×d，对每个μ_kj以及W_k={W_k,jl}_d×d的非对角元素进行自适应的L₁惩罚；基于由对数似然函数l(Ψ)和惩罚部分pen_λ(φ)组成的惩罚对数似然函数，得到模型的目标函数：

$l_{\mathrm{pen}}\left(\mathrm{\Ψ}\right)=l\left(\mathrm{\Ψ}\right)-{\mathrm{pen}}_{\mathrm{\λ}}\left(\mathrm{\Φ}\right)$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\log \left[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\π}}_k{f}_k(x_i;{\mathrm{\θ}}_k)\right]-{\mathrm{\λ}}_1\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{j=1}{\overset{d}{\mathrm{\Σ}}}{w}_{\mathrm{kj}}\left|{\mathrm{\μ}}_{\mathrm{kj}}\right|-{\mathrm{\λ}}_2\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{l=1}{\overset{d}{\mathrm{\Σ}}}\underset{j=1,j\≠l}{\overset{d}{\mathrm{\Σ}}}{v}_{k,\mathrm{jl}}\left|W_{k,\mathrm{jl}}\right|,,---\left(3\right)$

其中，φ={μ₁,…,μ_K,W₁,…,W_K},λ={λ₁,λ₂}包括了对应于λ_k和W_k的非负正则化参数λ₁和λ₂，其预先设定并通过外部输入，w_kj、v_k,jl分别为对应于μ_kj和W_k,jl的惩罚的权重。

其中，网络聚类模型建立单元101中对混合模型的求解过程如下：

给定K和λ，使用期望最大化算法对目标函数l_pen(Ψ)进行优化：假设每个样本x_i对应于一个无法观察到的指标向量z_i=[z_i1 z_i2 … z_iK]，指定x_i所属的成分；如果x_i来自于成分k，则z_ik=1，否则z_ik=0；引入另一个隐变量u_i=[u_i1 u_i2 … u_iK]，其中，每个u_ik服从Gamma分布；则算法中样本x_i完整的数据变为其对应的惩罚完整数据对数似然函数可被分解为四个部分：

l_c,pen(Ψ)=l₁(π)+l₂(v)+l₃(φ)-pen_λ(φ),   （4）

其中

$l_1\left(\mathrm{\π}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{z}_{\mathrm{ik}}\log {\mathrm{\π}}_k,---\left(5\right)$

$l_2\left(v\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{z}_{\mathrm{ik}}[-\log \mathrm{\Γ}\left(\frac{v_k}{2}\right)+\frac{v_k}{2}\log \left(\frac{v_k}{2}\right)+\frac{v_k}{2}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})-\log u_{\mathrm{ik}}],---\left(6\right)$

$l_3\left(\mathrm{\Φ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[-\frac{d}{2}\log \left(2\mathrm{\π}\right)+\frac{d}{2}\log u_{\mathrm{ik}}+\frac{1}{2}\log |W_k|-\frac{1}{2}{u}_{\mathrm{ik}}{(x_i-{\mathrm{\μ}}_k)}^{\′}{W}_k(x_i-{\mathrm{\μ}}_k)].---\left(7\right)$

执行期望最大化算法期望步，在已知第t次迭代的参数估计Ψ^(t)的情况下，计算l_c,pen(Ψ)的期望Q(Ψ;Ψ^(t))，其中先求得E_Ψ(t)(z_ik|x_i)、E_Ψ(t)(u_ik|x_i,z_ik=1)和E_Ψ(t)(log u_ik|x_i,z_ik=1)；

$E_{\mathrm{\Ψ}\left(t\right)}\left(z_{\mathrm{ik}}\right|x_i)=\frac{{\mathrm{\π}}_k^{\left(t\right)}{f}_k(x_i;{\mathrm{\θ}}_k^{\left(t\right)})}{f(x_i;{\mathrm{\Ψ}}^{\left(t\right)})}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}.---\left(8\right)$

被看作为x_i属于第k个聚类的后验概率；

鉴于Gamma分布是高斯分布的共轭先验，有

$E_{\mathrm{\Ψ}\left(t\right)}\left(u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\frac{v_k^{\left(t\right)}+d}{v_k^{\left(t\right)}+\mathrm{\δ}(x_i;{\mathrm{\μ}}_k^{\left(t\right)},{\mathrm{\Σ}}_k^{\left(t\right)})}\stackrel{\mathrm{\Δ}}{=}{u}_{\mathrm{ik}}^{\left(t\right)},---\left(9\right)$

$E_{\mathrm{\Ψ}\left(t\right)}\left(\log u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\log u_{\mathrm{ik}}^{\left(t\right)}+\mathrm{\ψ}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right),---\left(10\right)$

其中 $\mathrm{\ψ}\left(s\right)=\{\∂\mathrm{\Γ}\left(s\right)/\∂s\}/\mathrm{\Γ}\left(s\right)$ 是Digamma函数；

执行期望最大化算法最大化步，得到各个参数的更新：

${\mathrm{\π}}_k^{(t+1)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}/n,k=1,\·\·\·,K.---\left(11\right)$

v_k在第t+1次迭代的估计满足等式：

$-\mathrm{\ψ}\left(\frac{v_k}{2}\right)+\log \left(\frac{v_k}{2}\right)+1+\frac{1}{n_k^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(k\right)}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})+\mathrm{\ψ}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right)=0,$

                                                            (12)

其中 $n_k^{\left(t\right)}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)};$

在t+1步，定义自适应的权重为：

$w_{\mathrm{kj}}^{\left(t\right)}=\frac{1}{{|\mathrm{\μ}}_{\mathrm{kj}}^{\left(t\right)}|+\∈},v_{k,\mathrm{jl}}^{\left(t\right)}=\frac{1}{\left|W_{k,\mathrm{jl}}^{\left(t\right)}\right|+\∈}.---\left(13\right)$

其中，设定∈=0.1，并使用L₁惩罚的结果作为初始估计值

和

μ_kj估计的更新方法为：如果

$\left|\frac{1}{w_{\mathrm{kj}}^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(\underset{l=1,l\≠j}{\overset{d}{\mathrm{\Σ}}}(x_{\mathrm{il}}-{\mathrm{\μ}}_{\mathrm{kl}}^{\left(t\right)})W_{k,\mathrm{jl}}^{\left(t\right)}+x_{\mathrm{ij}}{W}_{k,\mathrm{jj}}^{\left(t\right)})\right|\≤{\mathrm{\λ}}_1,---\left(14\right)$

则 ${\mathrm{\μ}}_{\mathrm{kj}}^{(t+1)}=0;$ 否则

$\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}\right){\mathrm{\μ}}_{\mathrm{kj}}^{(t+1)}{W}_{k,\mathrm{jj}}^{\left(t\right)}+{\mathrm{\λ}}_1{w}_{\mathrm{kj}}^{\left(t\right)}\mathrm{sign}\left({\mathrm{\μ}}_{\mathrm{kj}}^{(t+1)}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(x_i{W}_{k,.j}^{\left(t\right)}-{\mathrm{\μ}}_k^{\left(t\right)}{W}_{k,.j}^{\left(t\right)}-{\mathrm{\μ}}_{\mathrm{kj}}^{\left(t\right)}{W}_{k,\mathrm{jj}}^{\left(t\right)}).$

                                                        (15)

进一步，得到有关W_k的优化问题

$W_k^{(t+1)}=\arg \underset{W_k}{\max }\log \left[\det \left(W_k\right)\right]-\mathrm{tr}\left({\tilde{S}}_k{W}_k\right)-\underset{j\≠l}{\mathrm{\Σ}}{\mathrm{\λ}}_{k,\mathrm{jl}}^{\′}\left|W_{k,\mathrm{jl}}\right|,---\left(16\right)$

其中

${\tilde{S}}_k=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}{(x_i-{\mathrm{\μ}}_k^{\left(t\right)})}^{\′}(x_i-{\mathrm{\μ}}_k^{\left(t\right)})}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}},$ ${\mathrm{\λ}}_{k,\mathrm{jl}}^{\′}=\frac{{2\mathrm{\λ}}_2{v}_{k,\mathrm{jl}}^{\left(t\right)}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\τ}}_{\mathrm{ik}}^{\left(t\right)}}.$

该优化问题使用图套索求解。

其中，网络聚类模型建立单元101基于证据近似权重（approximate weight ofevidence，AWE）准则进行混合模型的选择，具体为：估计三个参数，聚类的个数K、惩罚参数λ₁和λ₂，使用证据近似权重准则进行模型选择，

AWE=-2l(Ψ)+2d_e(3/2+log(n)),   （17）

其中d_e=k-1+K+Kf-q₁+kd+Kd(d-1)/2-q₂/2是混合模型中有效参数的个数，q₁=#{(k,j)：μ_kj=0}，q₂=#{(k，j，l):W_k，lj=0}，通过网格搜索法寻找使AWE最小的最优组合(K，λ₁，λ₂)，n为样本个数，d为每个样本的基因的个数，l(Ψ)为对数似然函数。

其中，癌症亚型发现单元102发现癌症亚型的过程如下：

通过网络聚类模型建立单元101对混合模型中参数进行估计之后，聚类被定义为服从相似分布的样本，样本的分布由后验概率τ_ik决定；给定样本x_i，将其聚类为使其后验概率最大的类，即y_i=arg max_k τ_ik；获取属于相同类的样本组成的集合判断为新的癌症亚型。

其中，基因间相互关系探索单元103探索不同癌症亚型中基因间的相互关系的过程为：

将不同癌症亚型的每个基因看作是网络中的点，由尺度矩阵的逆W_k来决定点与点间的边，若W_k,jl≠0，则对应聚类k的网络中，j和l相连，否则j和l不相连；该网络中单连通的点集以及它们之间的边构成一系列的子网络；在每个聚类中不同的子网络表示不同癌症亚型中基因间的关系。

其中，生物标志物检测单元104检测癌症亚型相关的生物标记物的过程如下：

从基因间相互关系探索单元中获得一系列的子网络；

由于噪声基因和有信息基因间是不相关的，如果一个子网络中存在一个基于均值的有判别力基因，判断它是一个子网络生物标志物，从而识别出在基因表达谱上不存在差异的基因，所述不存在差异的基因与有判别能力的基因相关，共同完成一项生物功能；

最后，在剩下的子网络中，若其在各个聚类对应的W₁,…,W_K上存在差异，则判断该子网络为生物标志物；

将检测到的生物标志物数据输出。

本发明所提供的系统鉴于基因间网络或路径、和高维基因表达数据集中包含的异常点所带来的挑战，以及以上方法应用于基因表达数据集时存在的弊端，本发明的系统使用无限制协方差的惩罚基于模型的学生t聚类，用于与癌症相关生物标志物的检测，并探索基因间的关系。使用多元学生t分布来模拟样本的生成过程，使得系统不会受极端或异常基因影响。使用无限制的协方差代替对角的协方差，用于考虑基因间的关系从而发现那些在网络结构上存在差异的癌症子类。通过对均值和协方差矩阵的逆进行自适应的L1惩罚，达到生物标志物检测以及网络重构的目的。此外，本发明的系统还提供了一种新的基因选择准则，用于检测如下的生物标志物：在不同的癌症亚型中具有不同均值的基因、在基因表达谱上不存在差异但与其它有判别信息的基因相互作用共同完成一项生物功能、以及具有不同网络结构的基因。

在上述方案中，本发明的系统主要是基于如下算法构建整个架构，以完成于癌症亚型相关的生物标记物的检测，具体算法为：

1）使用K均值聚类初始化每个样本属于的癌症子类；

2）分别使用常数、样本概率、样本均值、和对应于样本协方差矩阵的对角阵的逆来初始化v_k、π_k、μ_k和W_k。

3）根据上述式（8）和（9）计算τ_ik和u_ik；

4）根据上述式（13）计算权重w_kj和v_k,jl；

5）根据上述式（11）、（12）、（14）‐（16）更新π_k、v_k、μ_k和W_k；

6）重复步骤(3)‐(5)直到|l_pen(Ψ^(t+1))-l_pen(Ψ^(t))|<10^-6；

7）若对于k′≠k，有

则样本x_i属于y_i=k的癌症子类，其中t_e为最终迭代的次数；

8）由尺度矩阵的逆

得到的与癌症亚型相关的子网络;

9）基于

和

使用基于网络的准则得到子网络生物标志物。

相同或相似的标号对应相同或相似的部件；

附图中描述位置关系的用于仅用于示例性说明，不能理解为对本专利的限制；

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (7)

1.一种基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，包括：

基因表达谱提取单元，从基因芯片中提取n个独立的包含d个基因的样本的表达谱χ={x₁,…,x_n}，其中x_i=[x_i1 x_i2 … x_id]表示样本i中d个基因的表达谱；

网络聚类模型建立单元，与基因表达谱提取单元连接，提取到的基因表达谱χ={x₁,…,x_n}输入其中，建立混合模型并对该混合模型进行求解，其中，该混合模型使用多元学生t分布来描述混合模型中的每个成分；

癌症亚型发现单元，与网络聚类模型建立单元连接，其利用混合模型和聚类方法发现新的癌症亚型；

基因间相互关系探索单元，与癌症亚型发现单元连接，将每个基因看作是网络中的点，通过挖掘网络结构探索不同癌症亚型中基因间的相互关系；

生物标志物检测单元，根据不同癌症亚型中基因间的相互关系完成与癌症亚型相关的生物标记物的检测。

2.根据权利要求1所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，网络聚类模型建立单元中建立混合模型过程如下：

假设每个样本x_i来自于具有K个成分的混合分布，其概率密度函数如下：

$f(x_i;\mathrm{\&Psi;})=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_k{f}_k(x_i;{\mathrm{\&theta;}}_k),$

其中Ψ={π₁，…，π_k，θ₁，…,θ_k}包括混合模型中所有的参数，f_k表示第k个成分，π_k≥0是第k个成分的混合比例，

θ_k为对应于f_k的未知参数集合；

假设每个成分f_k的参数集合为θ_k={μ_k,Σ_k,v_k}的多元学生t分布T(μ_k,Σ_k,v_k)，其中μ_k=[μ_k1 μ_k2 … μ_kd]为位置参数矩阵，Σ_k=(σ_k,jl)_d×d为尺度参数，而v_k为自由度，多元学生t分布的概率密度函数为：

$f_k(x_i;{\mathrm{\&theta;}}_k)=\frac{\mathrm{\&Gamma;}\left(\frac{v_k+d}{2}\right){\left|{\mathrm{\&Sigma;}}_k\right|}^{-\frac{1}{2}}}{\mathrm{\&Gamma;}\left(\frac{v_k}{2}\right){\left(\mathrm{\&pi;}{v}_k\right)}^{d/2}{\{1+\mathrm{\&delta;}(x_i;{\mathrm{\&mu;}}_k,{\mathrm{\&Sigma;}}_k)/v_k\}}^{(v_k+d)/2}},$

其中Γ(·)为Gamma函数， $\mathrm{\&delta;}(x_i;{\mathrm{\&mu;}}_k,{\mathrm{\&Sigma;}}_k)={(x_i-{\mathrm{\&mu;}}_k)}^{\&prime;}{\mathrm{\&Sigma;}}_k^{-1}(x_i-{\mathrm{\&mu;}}_k)$ 表示x_i和μ_k间的马氏距离，每个学生t分布的均值和协方差矩阵分别为μ_k和d表示每个样本包含的基因的个数，j和l是矩阵Σ_k中元素所在的位置，π为圆周率;

记尺度参数矩阵Σ_k的逆为W_k={W_k,jl}_d×d，对每个μ_kj以及W_k={W_k，jl}_d×d的非对角元素进行自适应的L₁惩罚;基于由对数似然函数l(Ψ)和惩罚部分pen_λ(φ)组成的惩罚对数似然函数，得到模型的目标函数:

$l_{\mathrm{pen}}\left(\mathrm{\&Psi;}\right)=l\left(\mathrm{\&Psi;}\right)-{\mathrm{pen}}_{\mathrm{\&lambda;}}\left(\mathrm{\&Phi;}\right)$

$=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\log \left[\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_k{f}_k(x_i;{\mathrm{\&theta;}}_k)\right]-{\mathrm{\&lambda;}}_1\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{d}{\mathrm{\&Sigma;}}}{w}_{\mathrm{kj}}\left|{\mathrm{\&mu;}}_{\mathrm{kj}}\right|-{\mathrm{\&lambda;}}_2\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{l=1}{\overset{d}{\mathrm{\&Sigma;}}}\underset{j=1,j\&NotEqual;l}{\overset{d}{\mathrm{\&Sigma;}}}{v}_{k,\mathrm{jl}}\left|W_{k,\mathrm{jl}}\right|,$

其中，φ={μ₁,…,μ_K,W₁,…,W_K}，λ={λ₁,λ₂}包括了对应于μ_k和W_k的非负正则化参数λ₁和λ₂，其预先设定并通过外部输入，w_kj、v_k,jl分别为对应于μ_kj和W_k,jl的惩罚的权重。

3.根据权利要求2所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，网络聚类模型建立单元中对混合模型的求解过程如下：

给定K和λ，使用期望最大化算法对目标函数l_pen(Ψ)进行优化：假设每个样本x_i对应于一个无法观察到的指标向量z_i=[z_i1 z_i2 … z_iK]，指定x_i所属的成分；如果x_i来自于成分k，则z_ik=1，否则z_ik=0；引入另一个隐变量u_i=[u_i1 u_i2 … u_iK]，其中，每个u_ik服从Gamma分布；则算法中样本x_i完整的数据变为

其对应的惩罚完整数据对数似然函数可被分解为四个部分：

l_c,pen(Ψ)=l₁(π)+l₂(v)+l₃(φ)-pen_λ(φ),

其中

$l_1\left(\mathrm{\&pi;}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}}\log {\mathrm{\&pi;}}_k,$

$l_2\left(v\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{z}_{\mathrm{ik}}[-\log \mathrm{\&Gamma;}\left(\frac{v_k}{2}\right)+\frac{v_k}{2}\log \left(\frac{v_k}{2}\right)+\frac{v_k}{2}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})-\log u_{\mathrm{ik}}],$

$l_3\left(\mathrm{\&phi;}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}[-\frac{d}{2}\log \left(2\mathrm{\&pi;}\right)+\frac{d}{2}\log u_{\mathrm{ik}}+\frac{1}{2}\log |W_k|-\frac{1}{2}{u}_{\mathrm{ik}}{(x_i-{\mathrm{\&mu;}}_k)}^{\&prime;}{W}_k(x_i-{\mathrm{\&mu;}}_k)].$

执行期望最大化算法期望步，在已知第t次迭代的参数估计Ψ^(t)的情况下，计算l_c,pen(Ψ)的期望Q(Ψ;Ψ^(t))，其中先求得E_Ψ(t)(z_ik|x_i)、E_Ψ(t)(u_ik|x_i,z_ik=1)和E_Ψ(t)(logu_ik|x_i,z_ik=1)；

$E_{\mathrm{\&Psi;}\left(t\right)}\left(z_{\mathrm{ik}}\right|x_i)=\frac{{\mathrm{\&pi;}}_k^{\left(t\right)}{f}_k(x_i;{\mathrm{\&theta;}}_k^{\left(t\right)})}{f(x_i;{\mathrm{\&Psi;}}^{\left(t\right)})}\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}.$

被看作为x_i属于第k个聚类的后验概率；

鉴于Gamma分布是高斯分布的共轭先验，有

$E_{\mathrm{\&Psi;}\left(t\right)}\left(u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\frac{v_k^{\left(t\right)}+d}{v_k^{\left(t\right)}+\mathrm{\&delta;}(x_i;{\mathrm{\&mu;}}_k^{\left(t\right)},{\mathrm{\&Sigma;}}_k^{\left(t\right)})}\stackrel{\mathrm{\&Delta;}}{=}{u}_{\mathrm{ik}}^{\left(t\right)},$

$E_{\mathrm{\&Psi;}\left(t\right)}\left(\log u_{\mathrm{ik}}\right|x_i,z_{\mathrm{ik}}=1)=\log u_{\mathrm{ik}}^{\left(t\right)}+\mathrm{\&psi;}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right),$

其中 $\mathrm{\&psi;}\left(s\right)=\{\&PartialD;\mathrm{\&Gamma;}\left(s\right)/\&PartialD;s\}/\mathrm{\&Gamma;}\left(s\right)$ 是Digamma函数；

执行期望最大化算法最大化步，得到各个参数的更新：

${\mathrm{\&pi;}}_k^{(t+1)}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}/n,k=1,\&CenterDot;\&CenterDot;\&CenterDot;,K.$

v_k在第t+1次迭代的估计满足等式：

$-\mathrm{\&psi;}\left(\frac{v_k}{2}\right)+\log \left(\frac{v_k}{2}\right)+1+\frac{1}{n_k^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(k\right)}(\log u_{\mathrm{ik}}-u_{\mathrm{ik}})+\mathrm{\&psi;}\left(\frac{v_k^{\left(t\right)}+d}{2}\right)-\log \left(\frac{v_k^{\left(t\right)}+d}{2}\right)=0,$ 其中 $n_k^{\left(t\right)}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)};$

在t+1步，定义自适应的权重为：

$w_{\mathrm{kj}}^{\left(r\right)}=\frac{1}{\left|{\mathrm{\&mu;}}_{\mathrm{kj}}^{\left(t\right)}\right|+\&Element;},$ $v_{k,\mathrm{jl}}^{\left(t\right)}=\frac{1}{\left|W_{k,\mathrm{jl}}^{\left(t\right)}\right|+\&Element;}.$

其中，设定∈=0.1，并使用L₁惩罚的结果作为初始估计值

和

μ_kj估计的更新方法为：如果

$\left|\frac{1}{w_{\mathrm{kj}}^{\left(t\right)}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(\underset{l=1,l\&NotEqual;j}{\overset{d}{\mathrm{\&Sigma;}}}(x_{\mathrm{il}}-{\mathrm{\&mu;}}_{\mathrm{kl}}^{\left(t\right)})W_{k,\mathrm{jl}}^{\left(t\right)}+x_{\mathrm{ij}}{W}_{k,\mathrm{jj}}^{\left(t\right)})\right|\&le;{\mathrm{\&lambda;}}_1,$

则 ${\mathrm{\&mu;}}_{\mathrm{kj}}^{(t+1)}=0;$ 否则

$\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}\right){\mathrm{\&mu;}}_{\mathrm{kj}}^{(t+1)}{W}_{k,\mathrm{jj}}^{\left(t\right)}+{\mathrm{\&lambda;}}_1{w}_{\mathrm{kj}}^{\left(t\right)}\mathrm{sign}\left({\mathrm{\&mu;}}_{\mathrm{kj}}^{(t+1)}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}(x_i{W}_{k,.j}^{\left(t\right)}-{\mathrm{\&mu;}}_k^{\left(t\right)}{W}_{k,.j}^{\left(t\right)}-{\mathrm{\&mu;}}_{\mathrm{kj}}^{\left(t\right)}{W}_{k,\mathrm{jj}}^{\left(t\right)}).$

进一步，得到有关W_k的优化问题

$W_k^{(t+1)}=\arg \underset{W_k}{\max }\log \left[\det \left(W_k\right)\right]-\mathrm{tr}\left({\tilde{S}}_k{W}_k\right)-\underset{j\&NotEqual;l}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{k,\mathrm{jl}}^{\&prime;}\left|W_{k,\mathrm{jl}}\right|,$

其中

${\tilde{S}}_k=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}{u}_{\mathrm{ik}}^{\left(t\right)}{(x_i-{\mathrm{\&mu;}}_k^{\left(t\right)})}^{\&prime;}(x_i-{\mathrm{\&mu;}}_k^{\left(t\right)})}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}},$ ${\mathrm{\&lambda;}}_{k,\mathrm{jl}}^{\&prime;}=\frac{{2\mathrm{\&lambda;}}_2{v}_{k,\mathrm{jl}}^{\left(t\right)}}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&tau;}}_{\mathrm{ik}}^{\left(t\right)}}.$

该优化问题使用图套索求解。

4.根据权利要求3所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，网络聚类模型建立单元基于证据近似权重准则进行混合模型的选择，具体为:估计三个参数，聚类的个数K、惩罚参数λ₁和λ₂，使用证据近似权重准则进行模型选择，

AWE=-2l(Ψ)+2d_e(3/2+log(n))，

其中d_e=K-1+K+Kd-q₁+Kd+Kd(d-l)/2-q₂/2是混合模型中有效参数的个数，q₁=#{(k,j):μ_kj=0}，q₂=#{(k,j，l):W_k,lj=0}，通过网格搜索法寻找使AWE最小的最优组合(K，λ₁，λ₂)，n为样本个数，d为每个样本的基因的个数，l(Ψ)为对数似然函数。

5.根据权利要求4所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，癌症亚型发现单元发现癌症亚型的过程如下:

通过网络聚类模型建立单元对混合模型中参数进行估计之后，聚类被定义为服从相似分布的样本，样本的分布由后验概率τ_ik决定;给定样本x_i，将其聚类为使其后验概率最大的类，即y_i=arg max_k τ_ik;获取属于相同类的样本组成的集合判断为新的癌症亚型。

6.恨掂权利要求5所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，基因间相互关系探索单元探索不同癌症亚型中基因间的相互关系的过程为:

将不同癌症亚型的每个基因看作是网络中的点，由尺度矩阵的逆W_k来决定点与点间的边，若W_k，jl≠0，则对应聚类k的网络中，j和l相连，否则j和l不相连;该网络中单连通的点集以及它们之间的边构成一系列的子网络;在每个聚类中不同的子网络表示不同癌症亚型中基因间的关系。

7.根据权利要求6所述的基于学生t分布的癌症亚型生物标志物检测系统，其特征在于，生物标志物检测单元检测癌症亚型相关的生物标记物的过程如下:

从基因间相互关系探索单元中获得一系列的子网络;

由于噪声基因和有信息基因间是不相关的，如果一个子网络中存在一个基于均值的有判别力基因，判断它是一个子网络生物标志物，从而识别出在基因表达谱上不存在差异的基因，所述不存在差异的基因与有判别能力的基因相关，共同完成一项生物功能;

最后，在剩下的子网络中，若其在各个聚类对应的W₁，…，W_K上存在差异，则判断该子网络为生物标志物;

将检测到的生物标志物数据输出。