基于Dirichlet過程無限混合模型的基因表達數據聚類算法

張 林，劉 輝

（中國礦業大學 信息與電氣工程學院，江蘇 徐州 221116）

0 引言

近年來基因芯片、基因表達數據的出現導致基因表達數據的爆炸性增長，有關DNA基因表達數據的生物信息學研究方法發展迅速。基于高通量的DNA基因表達數據實現相關預測已成為生物醫學實驗檢測技術的重要補充。如何從具備高維度、小樣本、高冗余特征的基因表達數據中利用計算機分析工具得到有用的信息，已成為基因表達數據分析的主要內容。聚類分析作為一種重要的數據分析方法，在基因表達數據的分析中已有廣泛的應用，如通過對樣本的聚類分析自動對不同的疾病亞型或實驗條件實現區分，通過對基因聚類發現未知基因的功能，等等。

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

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

1 基因表達數據的Dirichlet無限混合模型

1.1 DNA基因表達數據

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

1.2 基因表達數據的Dirichlet過程混合模型

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

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

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

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

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

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

2 參數估計

根據Bayesian理論:

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

2.1 Dirichlet過程先驗分布

2.1.1 參數S和G的先驗分布

（1）參數S的先驗分布

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

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

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

（2）G的先驗分布

G的先驗分布，即（μk,）的先驗分布。假設μk是p維均值向量，是p×p協方差矩陣。因為μk,（j=1,…,K）相互獨立，在給定的條件下μk服從正態分布。即

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

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

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

2.1.2 參數π的先驗分布

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

2.2 參數的滿條件分布

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

為描述方便，記

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

2.2.2的后驗分布

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

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

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

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

3 仿真結果與分析

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

3.1 仿真數據

首先構造一組包含K=5個成分模型、N=400個樣本的模擬數據集加以檢驗，該模擬數據集基于一個四維正態混合模型產生。其中，各成分模型的權重π={0.18,0.10,0.18,0.25,0.29}，對應的各成分正態成分模型的均值由正態先驗分布隨機產生：μ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}，對應的各成分正態混合模型的方差由InverseWishart分布隨機產生：σ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過程無限混合模型對模擬數據展開聚類分析，經過200次MCMC Gibbs采樣，估計出模型中的各參數。采樣過程中聚類個數K、Dirichlet過程的集中度參數α及聚類個數K的直方圖如圖1中所示。可以看出，對上述模擬數據，通過Dirichlet過程無限混合模型聚類分析估計出的數據中潛在的成分模型的個數，完全符合產生該模擬數據集的條件。

3.2 基因表達數據

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

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

對上述基因表達數據，同樣建立Dirichlet過程無限混合模型進行聚類分析，并估計出模型中的各參數。采樣過程中聚類個數K、Dirichlet過程的集中度參數α及聚類個數K的直方圖如圖2中所示。可以看出，對上述基因表達數據，通過建立Dirichlet過程無限混合模型聚類分析估計出的數據中潛在的成分模型的個數，也是符合該測試數據集的生物學描述的。

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

4 結論

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

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

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

[3]姚宗靜.基于Dirichlet過程的非參數貝葉斯分析[D].西南交通大學,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.