采用粒子群優(yōu)化的基因轉(zhuǎn)錄差異分析模型

馮偉興王科俊

(哈爾濱工程大學(xué)模式識別與智能系統(tǒng)研究所，哈爾濱 150001)

采用粒子群優(yōu)化的基因轉(zhuǎn)錄差異分析模型

馮偉興*王科俊

(哈爾濱工程大學(xué)模式識別與智能系統(tǒng)研究所，哈爾濱 150001)

高通測序技術(shù)可以更加全面地檢測實(shí)驗樣本中基因轉(zhuǎn)錄水平，這使得對樣本間基因轉(zhuǎn)錄差異的分析精度越來越準(zhǔn)確。根據(jù)通過DNA高通測序技術(shù)獲得的基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白個數(shù)，提出了兩個樣本間基因轉(zhuǎn)錄差異分析模型。該模型在考慮同一基因間轉(zhuǎn)錄差異的同時，引入了全基因的轉(zhuǎn)錄分布特性以提高模型分析精度。模型采用預(yù)測概率最大化統(tǒng)計算法進(jìn)行參數(shù)求取。在概率最大化步驟中，由于難以得到模型參數(shù)解析解和提高算法效率，采用粒子群優(yōu)化算法直接進(jìn)行求取模型參數(shù)數(shù)值解。乳腺癌實(shí)驗數(shù)據(jù)測試表明，該模型可有效分析不同樣本間基因轉(zhuǎn)錄差異。

基因;轉(zhuǎn)錄差異;粒子群;分析模型

Abstract:Gene transcriptional level in an experiment can be more widely measured with high-throughput sequencing technology，which can obviously improve the analysis accuracy of gene transcriptional change between two experimental samples.Using numbers of PolⅡproteins contained inside each gene transcription region checked with DNA sequencing technology，a novel model was proposed to analyze gene transcription change considering gene transcription levels and the whole genome transcriptional distribution factor as well.Expectation-maximization(EM)statistics algorithm was utilized to resolve the model parameters.For obtaining analytical solution of model parameters and improving the efficiency of algorithm，particle swarm optimization(PSO)was used in maximization step to directly search parameters optimization numerical values.The test to analyze gene transcription change between normal and breast cancer samples verified the effectiveness of the proposed model.

Key words:gene;transcription difference;particle swarm;analysis model

引言

國際基因組測序計劃的順利完成，形成了由幾億個堿基代碼組成的DNA序列集，并在其上發(fā)現(xiàn)了25 000個以上的人體基因。在此基礎(chǔ)上，隨著研究的逐步深入，發(fā)現(xiàn)在眾多調(diào)控因子高度非線性的綜合調(diào)控下，以各種功能蛋白為橋梁，基因間形成著錯綜復(fù)雜的調(diào)控關(guān)系。一個基因的表達(dá)發(fā)生變異，就有可能導(dǎo)致整個細(xì)胞功能異常，甚至癌病變。因此，在正常和異常兩個樣本間高精度地分析出表達(dá)發(fā)生變異的基因，可以為相關(guān)疾病診治提供有意義的參考依據(jù)。

隨著計算機(jī)、控制和信息技術(shù)的飛速發(fā)展，人們獲取細(xì)胞樣本內(nèi)基因表達(dá)數(shù)據(jù)的方法越來越成熟［1］。20世紀(jì)末出現(xiàn)并逐步發(fā)展起來的基因芯片技術(shù)，提供了測量樣本中基因表達(dá)程度的有效方法，基于通過基因芯片技術(shù)獲得的基因表達(dá)數(shù)據(jù)，Kerr和Wolfinger設(shè)計了方差分析模型進(jìn)行不同樣本間基因表達(dá)差異測試［2-3］。Dudoit則采用統(tǒng)計學(xué)中的T測試法組合比較多樣本間的基因表達(dá)差異［4］。由于基因芯片僅能做到大體上的全基因表達(dá)測量，因此采用簡單的統(tǒng)計學(xué)方法進(jìn)行樣本間基因表達(dá)差異分析是足夠的。但隨著近幾年高通測序技術(shù)的出現(xiàn)和成熟，真正意義的全基因表達(dá)測量成為可能。與已有的不同樣本間基因表達(dá)差異分析方法相比，基于高通測序技術(shù)，不僅可以參考樣本間基因表達(dá)測得數(shù)據(jù)，還可借鑒樣本中基因表達(dá)的分布特性，從而可以更準(zhǔn)確地對基因表達(dá)差異進(jìn)行分析。

相比較而言，DNA高通測序技術(shù)已經(jīng)比較成熟，可用來獲取DNA上所結(jié)合轉(zhuǎn)錄因子、PolⅡ等蛋白信息，該信息因可被用于分析基因前轉(zhuǎn)錄調(diào)控機(jī)制，并間接反映基因表達(dá)水平及其基因間調(diào)控關(guān)系而得到廣泛應(yīng)用［4-7］。在此基礎(chǔ)上，RNA高通測序技術(shù)也已逐漸成熟。該技術(shù)可直接反映基因表達(dá)水平并可用于基因變異分析。為了更全面地獲取有用信息，相比于一個基因芯片僅能獲得幾十萬條有用信息，目前一次高通測序周期可同時處理30億個堿基并生成接近1億條讀數(shù)。據(jù)估計，2009年年底該數(shù)字將達(dá)到1 000億個堿基和30億條讀數(shù)。面對更加全面的數(shù)據(jù)，通過設(shè)計新方法對其進(jìn)行有效的生物信息處理，無疑將得到與原有技術(shù)相比精度更高的分析成果［8-11］。由于DNA高通測序技術(shù)發(fā)展得比較成熟，測得數(shù)據(jù)精度也較高，本文利用該技術(shù)測得的數(shù)據(jù)進(jìn)行不同樣本間基因表達(dá)差異分析。具體就是利用DNA高通測序技術(shù)測量基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白數(shù)目，進(jìn)而分析不同樣本間基因轉(zhuǎn)錄差異。盡管基因轉(zhuǎn)錄差異是基因表達(dá)差異的有效間接反映，但對于分析不同樣本間致病基因及其致病機(jī)理仍然具有很大的參考價值?；诖?，本文利用DNA高通測序技術(shù)測得的數(shù)據(jù)提出一種分析不同樣本間基因轉(zhuǎn)錄差異的分析模型，該模型不僅參考樣本間基因轉(zhuǎn)錄測得數(shù)據(jù)，還借鑒了樣本中基因轉(zhuǎn)錄的分布特性。鑒于模型的復(fù)雜性，在模型參數(shù)的求取過程中采用智能尋優(yōu)算法。

1 基因轉(zhuǎn)錄差異分析模型

在基因轉(zhuǎn)錄過程中，PolⅡ蛋白結(jié)合在基因轉(zhuǎn)錄區(qū)內(nèi)并對基因進(jìn)行轉(zhuǎn)錄，因此，基于基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白數(shù)量可以進(jìn)行不同樣本間基因轉(zhuǎn)錄差異分析［12］。

1.1 基因轉(zhuǎn)錄數(shù)據(jù)獲取

利用DNA高通測序中的 ChIP-seq技術(shù)可以對基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白數(shù)量進(jìn)行測量和統(tǒng)計，以直接反映基因轉(zhuǎn)錄水平。該技術(shù)首先利用超聲波將DNA鏈降解為DNA片段，然后利用特制的抗體俘獲結(jié)合在DNA片段上的PolⅡ蛋白，再利用沉淀技術(shù)(IP)將含有抗體的DNA片段濾出，隨后通過測序技術(shù)(seq)對所有濾出的DNA片段測序并通過序列比對映射回DNA上，最后根據(jù)基因轉(zhuǎn)錄區(qū)在DNA上的位置定義即可實(shí)現(xiàn)對基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白數(shù)量的測量和統(tǒng)計［13］。

1.2 基因轉(zhuǎn)錄差異分析

在兩個樣本之間，例如正常樣本和癌癥樣本，為了分析各個基因的轉(zhuǎn)錄是否發(fā)生變化以及變化程度，這里分為基因轉(zhuǎn)錄變化和基因轉(zhuǎn)錄不變化兩類。然后，依據(jù)不同樣本中分布在各基因轉(zhuǎn)錄區(qū)內(nèi)PolⅡ蛋白的個數(shù)及其個數(shù)變化程度來判斷各個基因?qū)儆谀囊活愐约半`屬程度。為此，提出基因轉(zhuǎn)錄模型及其轉(zhuǎn)錄差異模型并對此進(jìn)行分析。

1.2.1 基因轉(zhuǎn)錄模型

由于細(xì)胞中對全部n個基因的調(diào)控強(qiáng)度不同，因此一個PolⅡ蛋白結(jié)合在不同基因轉(zhuǎn)錄區(qū)內(nèi)以對其進(jìn)行轉(zhuǎn)錄的概率是不同的，設(shè)該概率為Pi(i=1，…，n)。假設(shè)細(xì)胞中PolⅡ蛋白的個數(shù)處于較低的水平，設(shè)其個數(shù)為m，則在第i個基因的轉(zhuǎn)錄區(qū)內(nèi)結(jié)合k個PolⅡ蛋白的概率可用二項分布表示:

按照這一概率分布，隨著細(xì)胞中PolⅡ蛋白的個數(shù)增多，結(jié)合在各基因轉(zhuǎn)錄區(qū)內(nèi)的PolⅡ蛋白個數(shù)將隨之增加。但隨著細(xì)胞中PolⅡ蛋白的個數(shù)繼續(xù)增多，由于基因轉(zhuǎn)錄區(qū)上可容納的PolⅡ蛋白的個數(shù)有限，新的PolⅡ蛋白結(jié)合在基因轉(zhuǎn)錄區(qū)內(nèi)的概率將隨之下降，最終，基因轉(zhuǎn)錄區(qū)內(nèi)結(jié)合的PolⅡ蛋白個數(shù)將出現(xiàn)飽和，即mPi(i=1，…，n)將保持不變。事實(shí)上，為了保證正常生理功能，細(xì)胞中PolⅡ蛋白的個數(shù)多處于飽和狀態(tài)。此時，各基因轉(zhuǎn)錄區(qū)內(nèi)結(jié)合PolⅡ蛋白的個數(shù)應(yīng)遵循泊松分布:

式中，yi，j是第 j(j=1，2)個樣本中第 i個基因轉(zhuǎn)錄區(qū)內(nèi)結(jié)合 Pol Ⅱ蛋白的個數(shù)，而λi，j=mjPi，j為分布常數(shù)。按照泊松分布特性，實(shí)際測量的PolⅡ蛋白個數(shù) yi，j是基因 i在樣本 j中實(shí)際轉(zhuǎn)錄程度 λi，j的表象，但如何通過 yi，j準(zhǔn)確 計算出 λi，j，僅依靠yi，j是無法做到的。此外，即使計算出 λi，j，如何評估其樣本間變化程度也是很困難的。一個基因轉(zhuǎn)錄從100變?yōu)?00和另一個基因轉(zhuǎn)錄從1 000變?yōu)? 000，其表示的轉(zhuǎn)錄差異程度是不同的，前者更劇烈一些。為此，結(jié)合樣本中所有基因的表達(dá)分布特性，設(shè)計了基因轉(zhuǎn)錄差異模型分析樣本間基因的轉(zhuǎn)錄差異及其差異程度。

1.2.2 基因轉(zhuǎn)錄差異模型

據(jù)Newton實(shí)驗分析，在同一樣本中，所有基因的表達(dá)水平將遵循伽瑪分布［14～15］。基于此，可認(rèn)為:同一樣本內(nèi)，基因轉(zhuǎn)錄實(shí)際值 λi(i=1，…，n)應(yīng)遵循伽瑪分布:

其中，α，β為模型常數(shù)。

兩個樣本之間，如果第i個基因的轉(zhuǎn)錄未發(fā)生變化，則在兩個樣本PolⅡ總數(shù)歸一化前提下，其在兩個樣本內(nèi)均遵循同一伽瑪分布，即:

否則:

為了計算基因轉(zhuǎn)錄差異程度，本研究定義變量Zi來描述兩個樣本之間第i個基因的轉(zhuǎn)錄是否發(fā)生變化，如果完全變化則Zi=1;完全不變化則Zi=0。

基于以上分析，本研究最終構(gòu)建了一個概率模型PA來綜合描述細(xì)胞中所有基因出現(xiàn)當(dāng)前轉(zhuǎn)錄變化情形的概率:

其中，yi1，yi2為基因 i在兩個樣本中測得的轉(zhuǎn)錄值。λi1，λi2為基因 i在兩個樣本中的實(shí)際轉(zhuǎn)錄值。α1，β1，α2，β2，Zi為模型參數(shù)。P 為任一個基因在兩個樣本之間發(fā)生轉(zhuǎn)錄變化概率。

如果概率模型PA的模型參數(shù)得到準(zhǔn)確求取，則通過Zi值就可以分析出第i個基因在不同樣本間，其轉(zhuǎn)錄是否發(fā)生變化以及變化程度。Zi越大則變化程度越大，反之亦然。

在上述概率模型中，如果Zi已知，即第i個基因在兩個樣本間轉(zhuǎn)錄變化程度是已知的，則可以通過調(diào)整 α1，β1，α2，β2的取值，直接計算出 λi，λi1，λi2值并使得PA最大化，從而實(shí)現(xiàn)對兩個樣本中基因轉(zhuǎn)錄情形的透徹分析。但由于Zi的取值是未知的，無法一次對所有模型參數(shù)在PA最大化前提下進(jìn)行求取，因此，采用統(tǒng)計學(xué)中預(yù)測概率最大化迭代算法來對 α1，β1，α2，β2進(jìn)行 PA概率最大化條件下的求取，并最終確定 λi，λi1，λi2及其 Zi值，以確定每個基因是否在不同樣本間發(fā)生轉(zhuǎn)錄變化及其變化程度。

這里采用的預(yù)測概率最大化算法分為反復(fù)迭代的兩步。

第一步是概率預(yù)測，即假設(shè) α1，β1，α2，β2已知的情況下求取Zi值。

首先對 λi，λi1，λi2進(jìn)行估算，即在 yi1，yi2，α1，β1，α2，β2已知的情況下尋找 λi，λi1，λi2值使 PA最大。

若Zi=0:

若Zi=1:

隨后，Zi可估算為:第二步是概率最大化，即在計算出的情況下尋找最優(yōu)的 α1，β1，α2，β2使 PA最大。通過對 PA表達(dá)式(6)的分析發(fā)現(xiàn)，這相當(dāng)于通過尋找 α1，β1，α2，β2值使下兩式最大

最后，對基因在兩個樣本之間發(fā)生轉(zhuǎn)錄變化的概率P值進(jìn)行更新:

以上過程迭代進(jìn)行，直至收斂。

1.3 粒子群優(yōu)化

在應(yīng)用預(yù)測概率最大化算法進(jìn)行概率模型PA的模型參數(shù)求取過程中發(fā)現(xiàn)，針對第一步，即概率預(yù)測，其計算比較容易實(shí)現(xiàn)。但對于第二步，即概率最大化，如果進(jìn)行精確的解析求解將很困難且計算量相當(dāng)大。與之相對應(yīng)，在復(fù)雜解空間內(nèi)搜尋最優(yōu)解恰好是智能尋優(yōu)算法的優(yōu)勢［16］，因此，這里采用智能尋優(yōu)算法中的粒子群算法，在概率最大化步驟中對模型參數(shù)直接進(jìn)行數(shù)值解尋優(yōu)求取。通過初步計算，發(fā)現(xiàn)模型參數(shù) α1，β1和 α2，β2的解空間比較簡單，因此采用基本粒子群算法進(jìn)行 α1，β1和α2，β2的數(shù)值解求取。另外，在式(10)中，由于 α1，β1和 α2，β2不具有相關(guān)性，因此這里可對 α1，β1和α2，β2分別進(jìn)行尋優(yōu)。

具體粒子群算法設(shè)計如下:

算法初始化時，取h個粒子，其初始值均勻分布在二維解空間內(nèi)。對于每個粒子，其解空間內(nèi)當(dāng)前位置矢量的評估函數(shù)采用式(10)。對于第k(k=1，…，h)個粒子，為了在解空間內(nèi)搜索更優(yōu)的位置矢量，其移動速度V′k采用下式進(jìn)行更新

式中，Pk為第k個粒子在解空間內(nèi)當(dāng)前位置矢量。Pglobal為所有粒子中解空間內(nèi)當(dāng)前最優(yōu)的位置矢量。Pk-local為第 k個粒子解空間內(nèi)當(dāng)前最優(yōu)的位置矢量。Vk為第 k個粒子當(dāng)前移動速度。C0，C1，C2為各部分權(quán)值以控制第k個粒子在解空間內(nèi)的移動方向和速度。

最后，第k個粒子在解空間內(nèi)的新位置矢量按下式計算:

以上過程重復(fù)進(jìn)行，直至無法在解空間內(nèi)找到更優(yōu)的全局位置矢量。

2 實(shí)驗測試結(jié)果和討論

2.1 實(shí)驗設(shè)計

為了對本研究所提出基因轉(zhuǎn)錄差異分析模型的有效性進(jìn)行測試，這里將選取美國俄亥俄州立大學(xué)公開提供的兩種MCF-7乳腺癌樣本PolⅡ測量數(shù)據(jù)。一種樣本是普通的MCF-7乳腺癌樣本，另一種是具有抗藥性的 MCF-7乳腺癌樣本。該數(shù)據(jù)共分為四組，分別是加藥前普通乳腺癌、加藥后普通乳腺癌、加藥前抗藥乳腺癌和加藥后抗藥乳腺癌樣本中利用DNA高通測序技術(shù)測得的25470個基因轉(zhuǎn)錄區(qū)內(nèi)結(jié)合的PolⅡ蛋白數(shù)量。

按預(yù)計，在加藥的情況下，普通乳腺癌樣本內(nèi)基因的轉(zhuǎn)錄變化將比抗藥乳腺癌樣本更劇烈，這構(gòu)造出測試樣本分析模型有效性的良好環(huán)境。

對這四組實(shí)驗數(shù)據(jù)在進(jìn)行PolⅡ總數(shù)歸一化后的統(tǒng)計分析結(jié)果如表1所示。

表1 測試數(shù)據(jù)統(tǒng)計分析(單位:個)Tab.1 Statistical analysis of experiment data(Unit:Number )

由表1可以分析出，在4組測試數(shù)據(jù)歸一化后，不同樣本間的測試數(shù)據(jù)在均值一致的前提下方差變化有著明顯的差異。其中，普通乳腺癌樣本在加藥后基因整體的轉(zhuǎn)錄差異明顯降低，這可能是由于部分轉(zhuǎn)錄水平高的基因加藥后轉(zhuǎn)錄水平得到抑制導(dǎo)致的;而抗藥乳腺癌樣本加藥后基因的整體轉(zhuǎn)錄水平變化則不明顯，這說明其對藥物不敏感。這一現(xiàn)象是符合兩種樣本的自然特性的，因此，可證明實(shí)驗測得數(shù)據(jù)是有意義的。

在獲取所有基因整體的轉(zhuǎn)錄差異變化后，具體到每一個基因，如何評估其轉(zhuǎn)錄變化則將用本研究所提出分析模型進(jìn)行分析。

2.2 基因轉(zhuǎn)錄差異分析

首先，如果某個基因所在DNA區(qū)域內(nèi)組蛋白修飾或DNA甲基化等原因?qū)е翽olⅡ蛋白無法在其轉(zhuǎn)錄區(qū)內(nèi)結(jié)合，則該基因轉(zhuǎn)錄將不遵循泊松分布。因此，為保證計算精度，在測試數(shù)據(jù)中將去除加藥前后轉(zhuǎn)錄區(qū)內(nèi)均沒有結(jié)合PolⅡ蛋白的基因。這樣處理后，在全部25 470個基因中，普通乳腺癌樣本測試數(shù)據(jù)中將有24 020個基因參與后續(xù)分析，而抗藥乳腺癌樣本測試數(shù)據(jù)中有24 230個基因參與后續(xù)分析。

隨后，將4組轉(zhuǎn)錄數(shù)據(jù)分成兩個大組。每個大組含有一個乳腺癌樣本加藥前后的兩組基因轉(zhuǎn)錄數(shù)據(jù)。然后分別采用本文所提出的基因轉(zhuǎn)錄差異分析模型對其進(jìn)行個體基因轉(zhuǎn)錄差異分析。

在預(yù)測概率最大化算法的迭代參數(shù)優(yōu)化過程中，在概率預(yù)測步驟中，由于針對n個個體基因構(gòu)造的泊松模型參數(shù) λi，λi1，λi2(i=1，…，n)的最優(yōu)參數(shù)計算式(7)和式(8)也沒有解析解，依據(jù)其所遵循的伽瑪函數(shù)凸特性設(shè)計參數(shù)尋優(yōu)算法如下:λi，λi1，λi2分別由1開始以1為步長進(jìn)行增加并計算其對應(yīng)的式(7)和式(8)中伽瑪函數(shù)值，當(dāng)函數(shù)值由小變大過程中首次出現(xiàn)下降時，在其對應(yīng)的 λi，λi1，λi2值附近進(jìn)行尋優(yōu)，并找到 λi，λi1，λi2的最優(yōu)值。在概率最大化步驟中，將采用粒子群優(yōu)化算法。粒子群優(yōu)化算法的參數(shù)設(shè)置為:粒子個數(shù)為100，C0=1，C1=2，C2=2。

對于普通乳腺癌和抗藥乳腺癌兩大組基因轉(zhuǎn)錄數(shù)據(jù)，預(yù)測概率最大化算法的參數(shù)尋優(yōu)均在迭代4～5次后收斂，收斂后的 Zi(i=1，2，…，n)值將反映第i個基因轉(zhuǎn)錄變化的程度。其最終分析結(jié)果見表2和表3。

表2 普通乳腺癌分析結(jié)果Tab.2 Analytical results of normal breast cancer

表3 抗藥乳腺癌分析結(jié)果Tab.3 Analytical results of drug resist brea st cancer

結(jié)合表1和表2，普通乳腺癌樣本加藥后，盡管基因間整體轉(zhuǎn)錄差異有所降低，但就個體基因而言，在所考慮的24 020個基因中，仍有6505個基因的轉(zhuǎn)錄發(fā)生了明顯變化，占27.1%。其中，轉(zhuǎn)錄提高的基因共有4 346(18.1%)個，提高幅度大的(Z≥0.9)有2 474(10.3%)個，提高幅度小的(0.1≤Z＜0.9)有1 872(7.8%)個。而轉(zhuǎn)錄降低的基因僅為轉(zhuǎn)錄提高基因數(shù)目的一半，有2 159(9.0%)個，其中降低大的為1 373(5.7%)，降低小的為786(3.3%)?？梢娂铀幒?，普通乳腺癌樣本的反應(yīng)還是較劇烈的，更多的反應(yīng)是基因的轉(zhuǎn)錄水平發(fā)生了提高。

結(jié)合表1和表3，抗藥乳腺癌樣本加藥后，盡管基因間整體轉(zhuǎn)錄差異變化不大，但就個體基因而言，仍有6 505(8.4%)個基因的轉(zhuǎn)錄發(fā)生了明顯變化。其中，轉(zhuǎn)錄提高的基因較少，僅有489(2.0%)個，提高幅度大的有123(0.5%)個，提高幅度小的有366(1.5%)個。而轉(zhuǎn)錄降低的基因有1 543(6.4%)個，其中降低大的為1 273(5.3%)，降低小的為270(1.1%)?？梢娂铀幒?，抗藥乳腺癌樣本的反應(yīng)相比普通乳腺癌樣本明顯減弱。而這種減弱主要體現(xiàn)在加藥后基因的轉(zhuǎn)錄反而受到了抑制，這從基因轉(zhuǎn)錄降低的數(shù)目明顯多于基因轉(zhuǎn)錄提高數(shù)目中可以看出。

從以上分析可以看出，加藥對普通乳腺癌樣本和抗藥乳腺癌樣本的影響有著明顯的不同，相比較而言，普通乳腺癌樣本比抗藥乳腺癌樣本對加藥的反應(yīng)要強(qiáng)烈得多。而這種差異與實(shí)驗設(shè)計的預(yù)期結(jié)果是相吻合的。這充分反映出本研究所提出基因轉(zhuǎn)錄差異分析模型的有效性。

3 結(jié)論

利用可直接反映基因轉(zhuǎn)錄水平的基因轉(zhuǎn)錄區(qū)內(nèi)所結(jié)合PolⅡ 蛋白的個數(shù)，提出一個分析模型，分析每一個基因在兩個相關(guān)樣本間轉(zhuǎn)錄水平的變化程度，通過實(shí)際數(shù)據(jù)測試證明了該模型的有效性。

由于不同樣本不同實(shí)驗下所測得PolⅡ 蛋白整體個數(shù)差異較大，因此，即使通過歸一化進(jìn)行了校正，仍然會對本研究所提出模型的分析精度造成影響。目前，隨著DNA測序技術(shù)的發(fā)展，一次實(shí)驗所測得的數(shù)據(jù)越來越大。據(jù)估計，2009年年底該數(shù)字將達(dá)到1000億個堿基和30億條讀數(shù)。隨著測試數(shù)據(jù)量的提高，不同實(shí)驗下所測得PolⅡ 蛋白整體個數(shù)差異對本模型分析精度的影響將會越來越小。

另外，為了使分析結(jié)果有明顯的生物含義，所提出的模型是基于這樣一個前提設(shè)計的:即如果兩個樣本間基因轉(zhuǎn)錄水平未發(fā)生明顯變化，則其轉(zhuǎn)錄分布特性是基本不變的。這樣處理的好處是如果找到兩個樣本間轉(zhuǎn)錄分布特性發(fā)生明顯變化的基因，該基因就很可能與兩個樣本的自然差異有關(guān)(如加藥反應(yīng)，癌病變等)。因此，只要遵從這一前提，本研究所提出模型即可用于對兩個樣本間基因轉(zhuǎn)錄水平的差異進(jìn)行分析，也可以用于多樣本分析。具體方法可以以一個樣本為基準(zhǔn)，其它樣本對其分別利用本研究所提出模型進(jìn)行計算，并將所有計算結(jié)果集中在一起進(jìn)行綜合分析即可。

此外，由于PolⅡ 蛋白的個數(shù)只能反映基因轉(zhuǎn)錄水平，不能直接反映更有生物意義的基因表達(dá)水平，因此本模型還有待進(jìn)一步發(fā)展，以最終用于基因表達(dá)差異分析。目前，最新發(fā)展的生物技術(shù)已使RNA高通測序成為可能。相比于DNA高通測序僅在DNA層次反映樣本組織的基因轉(zhuǎn)錄模式，RNA高通測序則可直接測量樣本組織中所有基因是否表達(dá)以及表達(dá)的程度。這就為直接對基因表達(dá)差異分析提供了可能。針對這一新技術(shù)，通過將分析對象由基因轉(zhuǎn)錄區(qū)內(nèi)所結(jié)合 PolⅡ 蛋白個數(shù)替換為反映該基因表達(dá)水平的RNA高通測序序列讀數(shù)，本模型即可用于分析兩個相關(guān)樣本間基因表達(dá)差異分析，但具體實(shí)現(xiàn)方式還有待針對RNA高通測序技術(shù)進(jìn)行進(jìn)一步研究。

［1］張曉龍，楊艷霞.機(jī)器學(xué)習(xí)在生物信息學(xué)中的應(yīng)用［J］.武漢科技大學(xué)學(xué)報，2005，28(2):201-204

［2］Kerr MK，Martin M，Churchill GA.Analysis of variance for gene expression microarray data［J］.Journal of Computational Biology，2000，7:819-837.

［3］Wolfinger RD，Gibson G，Wolfinger ED.et al.Assessing Gene Significance from cDNA Microarray Expression Data via Mixed Models［J］.Journal of Computational Biology，2001，8:625-637.

［4］Dudoit S，Yang YH，Callow MJ，et al.Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments［J］.Statistica Sinica，2002，12:111-139.

［5］Schmidt D，Wilson MD，Spyrou C，et al.ChIP-seq:Using high-throughput sequencing to discover protein-DNA interactions［J］.Methods，2009，48(3):240-247.

［6］Johnson DS，Mortazavi A，Myers RM，et al.Genome-wide mapping of in vivo protein-DNA interactions［J］.Science，2007，316:1497-1499.

［7］Ansorge WJ.Next-generation DNA sequencing techniques［J］.N Biotechnol，2009，25(4):195-203.

［8］Fullwood MJ，Wei CL，Liu ETB，et al.Next-generation DNA sequencing of paired-end tags(PET)for transcriptome and genome analyses［J］.Genome Res，2009，19(4):521-32.

［9］Cole T，Lior P，Steven LS.TopHat:discovering splice junctions with RNA-Seq［J］.Bioinformatics，2009，25(9):1105-1111.

［10］Daniel RZ，Ewan B.Velvet:Algorithms for de novo short read assembly using de Bruijn graphs［J］.Genome Res，2008，18(5):821-829.

［11］David JS，William GR.Transcriptome sequencing of malignant pleural mesothelioma tumors［J］.PNAS，2008，105(9):3521-3526.

［12］特納PC著.分子生物學(xué)［M］.(第2版).北京:科學(xué)出版社，2001.

［13］Raja J，Suresh C.Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data［J］.Nucleic Acids Research，2008，36(16):5221-5231.

［14］Newton MA，KendziorskiCM，RichmondCS，etal.On differential variability of expression ratios:Improving statistical inference about gene expression changes from microarray data［J］.Journal of Computational Biology，2001，8:37-52.

［15］Kendziorski C，Newton M，Lan H，et al.On parametric empirical bayes methods for comparing multiple groups using replicated gene expression profiles［J］.Stat Med，2003，22(24):3899-3914.

［16］Parsopoulos KE，Vrahatis MN.Recent approaches to global optimization problems through Particle Swarm Optimization［J］.Natural Computing，2002，1:235-306.

A Model for Analyzing Gene Transcription Difference Using Particle Swarm Optimization

FENG Wei-Xing*WANG Ke-Jun
(Pattern Recognition and Intelligent System Institute of Harbin Engineering University，Harbin，Heilongjiang 150001，China)

Q332

A

0258-8021(2010)02-0229-06

10.3969/j.issn.0258-8021.2010.02.012

2009-07-10，

2009-11-04

國家高技術(shù)研究發(fā)展(863)計劃(2008AA01Z148)

*通訊作者。 E-mail:fengweixing@hrbeu.edu.cn