張 林，劉 輝

（中國礦業(yè)大學(xué) 信息與電氣工程學(xué)院，江蘇 徐州 221116）

0 引言

近年來基因芯片、基因表達(dá)數(shù)據(jù)的出現(xiàn)導(dǎo)致基因表達(dá)數(shù)據(jù)的爆炸性增長，有關(guān)DNA基因表達(dá)數(shù)據(jù)的生物信息學(xué)研究方法發(fā)展迅速。基于高通量的DNA基因表達(dá)數(shù)據(jù)實(shí)現(xiàn)相關(guān)預(yù)測已成為生物醫(yī)學(xué)實(shí)驗(yàn)檢測技術(shù)的重要補(bǔ)充。如何從具備高維度、小樣本、高冗余特征的基因表達(dá)數(shù)據(jù)中利用計(jì)算機(jī)分析工具得到有用的信息，已成為基因表達(dá)數(shù)據(jù)分析的主要內(nèi)容。聚類分析作為一種重要的數(shù)據(jù)分析方法，在基因表達(dá)數(shù)據(jù)的分析中已有廣泛的應(yīng)用，如通過對樣本的聚類分析自動對不同的疾病亞型或?qū)嶒?yàn)條件實(shí)現(xiàn)區(qū)分，通過對基因聚類發(fā)現(xiàn)未知基因的功能，等等。

大多基于模型的聚類算法在假設(shè)給定聚類數(shù)的前提下，根據(jù)待聚類樣本的屬性，建立有限混合模型對基因表達(dá)數(shù)據(jù)展開研究，模型中聚類數(shù)的確定問題通常通過模型選擇問題加以解決，因而聚類的準(zhǔn)確性和泛化性受到模型選擇準(zhǔn)則的影響。作為無限混合模型核心的Dirichlet過程，則被廣泛應(yīng)用于解決傳統(tǒng)的有限混合模型中子模型個數(shù)不確定的難題。

本文擬基于DNA基因表達(dá)數(shù)據(jù)建立Dirichlet過程無限混合模型展開聚類分析，其中的聚類數(shù)將由模型和數(shù)據(jù)自主計(jì)算得出[1～3]，無需獨(dú)立確定，因而更為靈活。

1 基因表達(dá)數(shù)據(jù)的Dirichlet無限混合模型

1.1 DNA基因表達(dá)數(shù)據(jù)

DNA基因表達(dá)數(shù)據(jù)通過荷載成千上萬個基因片段，實(shí)現(xiàn)高通量的生物學(xué)檢測，使得從整個人類基因組研究基因的表達(dá)與調(diào)控成為可能。但DNA基因表達(dá)數(shù)據(jù)實(shí)驗(yàn)中如熒光標(biāo)記效率、掃描參數(shù)設(shè)置、空間位置差異等各種變異都是基因表達(dá)水平原始數(shù)據(jù)中噪聲產(chǎn)生的來源。因此，DNA基因表達(dá)數(shù)據(jù)需要經(jīng)過標(biāo)準(zhǔn)化，以消除由于系統(tǒng)變異引起的誤差。經(jīng)過預(yù)處理的DNA基因表達(dá)數(shù)據(jù)通常服從正態(tài)分布p(μ,σ2)。

1.2 基因表達(dá)數(shù)據(jù)的Dirichlet過程混合模型

用X={x1,x2,…xN}表示DNA基因表達(dá)數(shù)據(jù)，N表示數(shù)據(jù)中包含的樣本個數(shù)，xi={xi1,…xiG}T表示第i個樣本，則xig表示第i個樣本第g個基因的表達(dá)水平，各樣本間相互獨(dú)立。X可由K個正態(tài)模型混合而成（K未知）。為求解K，本文定義一隱變量s={s1,…,sN}(si∈{1,…,K})表示樣本的聚類標(biāo)簽，si=k表示第i個樣本經(jīng)過聚類分析屬于第k類。用p(?)表示模型中各成分模型，各自遵照不同的分布參數(shù)θk，θk={μk,σk2}，μk表示第k個成分模型的均值，σk2表示方差。πk表示第k個成分模型的混合系數(shù)，滿足 πk≥0,k={1,…,K}并且假設(shè) Θ={π1,…,πK;θ1,…,θK;K}，則 Θ 表示模型中所有待求參數(shù)，則基于DNA基因表達(dá)數(shù)據(jù)可建立如式（1）所示的正態(tài)混合模型。

因此，本文的目標(biāo)即估計(jì)上述正態(tài)混合模型中的參數(shù)Θ。目前，求解此類問題的方法大概有兩種，一是期望值最大化算法（EM：Expectation Maximization），二是Bayesian隨機(jī)采樣算法。EM算法主要用于在極大似然準(zhǔn)則下估計(jì)模型參數(shù)，一直是有限混合模型參數(shù)估計(jì)問題的標(biāo)準(zhǔn)算法之一，但該算法容易陷入局部極值點(diǎn)。本文基于后種方法求解，采用Dirichlet過程作為先驗(yàn)分布，建立DNA基因表達(dá)數(shù)據(jù)的Dirichlet過程無限混合模型，利用基于Gibbs采樣的MCMC方法估計(jì)出Θ的后驗(yàn)分布[4]。

Dirichlet過程基于Dirichlet分布生成，作為分布上的分布，是Dirichlet分布在連續(xù)空間上的擴(kuò)展。通常，Dirichlet過程表示為

其中，G0是基分布；α(α＞0)是集中度參數(shù)，表示G逼近G0的程度；G表示基于Dirichlet過程在基分布和集中度參數(shù)基礎(chǔ)上產(chǎn)生的某隨機(jī)分布。α越大，G越接近G0。

因此，本文基于DNA基因表達(dá)數(shù)據(jù)建立Dirichlet過程無限混合模型[5～7]，如式(3)所示。

其中，p(xi|θsi)取正態(tài)分布。

2 參數(shù)估計(jì)

根據(jù)Bayesian理論:

因此，式(3)中的未知參數(shù)的估計(jì)可通過計(jì)算其先驗(yàn)分布及似然函數(shù)來實(shí)現(xiàn)。本文首先確定各參數(shù)的先驗(yàn)分布，在此基礎(chǔ)上確定各參數(shù)的滿條件分布，最后通過MCMC方法估計(jì)出參數(shù)Θ。

2.1 Dirichlet過程先驗(yàn)分布

2.1.1 參數(shù)S和G的先驗(yàn)分布

（1）參數(shù)S的先驗(yàn)分布

可用于描述Dirichlet過程的模型有多種，本文基于經(jīng)典的Blackwell-MacQueen Urn抽球模型，計(jì)算si的條件概率。其過程描述如下：某罐子中裝有K種不同顏色的球。初始，罐中有紅色球α1個，有黃色球α2個，…等，假設(shè)從罐中隨機(jī)取出一個球（取N次），每次取完之后將兩個相同顏色的球放回罐中。當(dāng)趨近于無窮時，罐中各種顏色的球所占的比例α1,…,αK將遵循Dir(α1,…,αK)。若將該罐中球的顏色的種類K擴(kuò)展到無限集，就得到Dirichlet過程[1,8,9]。即

因此，Si的條件概率可采用式（5）計(jì)算。

其中，δsj(?)表示第j類中的樣本個數(shù)。

（2）G的先驗(yàn)分布

G的先驗(yàn)分布，即（μk,）的先驗(yàn)分布。假設(shè)μk是p維均值向量，是p×p協(xié)方差矩陣。因?yàn)棣蘫,（j=1,…,K）相互獨(dú)立，在給定的條件下μk服從正態(tài)分布。即

服從p維Inverse Wishart分布，表示為

其中，δ表示形狀參數(shù)，δ=n-p+1，n表示自由度；其均值為Q/(δ-2)。

因此，（μk,σ2k）的共軛先驗(yàn)分布為normal-inverse Wishart分布，所以G的先驗(yàn)分布即選擇normal-inverse Wishart。

2.1.2 參數(shù)π的先驗(yàn)分布

π和（μ,σ2）相互獨(dú)立，并且π的先驗(yàn)分布為Dirichlet分布：π～Dir(α1,…,αK)。

2.2 參數(shù)的滿條件分布

本文利用Gibbs采樣MCMC方法估計(jì)DNA基因表達(dá)數(shù)據(jù)的Dirichlet過程無限混合模型中的參數(shù)。由上述參數(shù)的先驗(yàn)分布，估計(jì)模型中各參數(shù)的滿條件分布。

為描述方便，記

2.2.1μk的后驗(yàn)分布p(μk|μ-k,σ2,π,s,X)(k=1,…,K)

2.2.2的后驗(yàn)分布

2.2.3 πk的后驗(yàn)分布p(πk|π-k,μ,σ2,s,X)

若 給 定 s，則 π 與（μ,σ2,X）相互獨(dú)立 ，所以p(πk|π-k,μ,σ2,s,X)=p(πk|π-k,s)∝p(π|s)。對于先驗(yàn)分布為Dirichlet分布而言，后驗(yàn)分布也一定是Dirichlet分布。因此

2.2.4si的后驗(yàn)分布p(si|s-i,μ,σ2,s,X)(i=1,…,N)

給定X和μ,σ2,π，si相互獨(dú)立，所以p(si|s-i,μ,σ2,s,π

3 仿真結(jié)果與分析

為了驗(yàn)證本文所述算法的有效性，本文分別對仿真數(shù)據(jù)和IRIS測試數(shù)據(jù)集建立Dirichlet過程無限混合模型，展開聚類分析。

3.1 仿真數(shù)據(jù)

首先構(gòu)造一組包含K=5個成分模型、N=400個樣本的模擬數(shù)據(jù)集加以檢驗(yàn)，該模擬數(shù)據(jù)集基于一個四維正態(tài)混合模型產(chǎn)生。其中，各成分模型的權(quán)重π={0.18,0.10,0.18,0.25,0.29}，對應(yīng)的各成分正態(tài)成分模型的均值由正態(tài)先驗(yàn)分布隨機(jī)產(chǎn)生：μ1={19.7,6.5,5.6,28.6}，μ2={1.2,7.2,21.6,20.6}，μ3={7.4,0.2,1.1,1.0}，μ4={24.9,26.4,14.8,4.6}，μ5={20.4,9.2,9.2,1.1}，對應(yīng)的各成分正態(tài)混合模型的方差由InverseWishart分布隨機(jī)產(chǎn)生：σ21={3.5,7.9,4.7,3.4}，σ22={6.4,4.8,4.7,1.7}，σ23={2.8,3.4,3.7,7.4}，σ24={9.0,6.1,4.8,2.7}，σ25={9.0,4.3,5.9,6.1}。

通過Dirichlet過程無限混合模型對模擬數(shù)據(jù)展開聚類分析，經(jīng)過200次MCMC Gibbs采樣，估計(jì)出模型中的各參數(shù)。采樣過程中聚類個數(shù)K、Dirichlet過程的集中度參數(shù)α及聚類個數(shù)K的直方圖如圖1中所示。可以看出，對上述模擬數(shù)據(jù)，通過Dirichlet過程無限混合模型聚類分析估計(jì)出的數(shù)據(jù)中潛在的成分模型的個數(shù)，完全符合產(chǎn)生該模擬數(shù)據(jù)集的條件。

3.2 基因表達(dá)數(shù)據(jù)

圖1 聚類個數(shù)K、集中度參數(shù)α及K直方圖

本文以著名的白血病樣本基因表達(dá)數(shù)據(jù)集作為測試數(shù)據(jù)集。該數(shù)據(jù)集于1999年由Golub收集，共包含72個急性白血病樣本的7219個基因的表達(dá)水平，其中，急性淋巴瘤白血病（ALL）47例，急性骨髓瘤白血病（AML）25例。

對上述基因表達(dá)數(shù)據(jù)，同樣建立Dirichlet過程無限混合模型進(jìn)行聚類分析，并估計(jì)出模型中的各參數(shù)。采樣過程中聚類個數(shù)K、Dirichlet過程的集中度參數(shù)α及聚類個數(shù)K的直方圖如圖2中所示?？梢钥闯觯瑢ι鲜龌虮磉_(dá)數(shù)據(jù)，通過建立Dirichlet過程無限混合模型聚類分析估計(jì)出的數(shù)據(jù)中潛在的成分模型的個數(shù)，也是符合該測試數(shù)據(jù)集的生物學(xué)描述的。

圖2 聚類個數(shù)K、集中度參數(shù)α及K直方圖

4 結(jié)論

本文提出建立Dirichlet過程無限混合模型進(jìn)行DNA基因表達(dá)數(shù)據(jù)的聚類分析，該模型無需預(yù)先設(shè)定成分模型的個數(shù)，因而具有更好的靈活性和適應(yīng)性，有利于我們挖掘數(shù)據(jù)中的各種有用信息。模擬測試數(shù)據(jù)集和Golub急性白血病DNA基因表達(dá)測試數(shù)據(jù)集的聚類分析結(jié)果表明了該方法在無監(jiān)督聚類方法中的優(yōu)越性。通過建立Dirichlet過程無限混合模型開展的聚類分析算法，能夠正確地估計(jì)出DNA基因表達(dá)數(shù)據(jù)中隱含的成分模型的個數(shù)。

[1]Teh,Y.Dirichlet Process[EB/OL].http://www.gatsby.ucl.ac.uk/～ywtehe/search/npbayes/dp.pdf,2011.

[2]徐謙,周俊生,陳家駿.Dirichlet過程及其在自然語言處理中的應(yīng)用[J].中文信息學(xué)報(bào),2009，23(005).

[3]姚宗靜.基于Dirichlet過程的非參數(shù)貝葉斯分析[D].西南交通大學(xué),2007.

[4]Neal,R.Markov Chain Sampling Methods for Dirichlet Process Mix?ture Models[J].Journal of Computational and Graphical Statistics,2000,9(2).

[5]Gelman,A.,Carlin,J.,Stern,H.,Rubin,D.BayesianDataAnalysis[M].New York:CRC Press,2004.

[6]Kim,S.Mahlet G.Tadesse.Marina Vannucci.Bayesian Variable Selec?tion in Clustering Via dirichlet Process Mixture Models[J].Biometrika,2006,93(4).

[7]Dahl,D.Model-based Clustering for Expression Data Via a Dirichlet Process Mixture Model[J].Bayesian Inference for Gene Expression and Proteomics,2006.

[8]Teh,Y.,et al.Hierarchical Dirichlet Processes[J].Journal of the Ameri?can Statistical Association,2006,101(476).

[9]Rasmussen,C.,Z.Ghahramani.Infinite Mixtures of Gaussian Process Experts[C].Advances in Neural Information Processing Systems 14:Proceedings of the 2002 Conference,2002.