對基因組內(nèi)不同序列分布差異進行量化的探討△

李 治 崔躍華,2 張巖波，△

基因組常被比喻為一本“生命天書”，如果有幾個關(guān)鍵詞在書中出現(xiàn)的位置比較接近，則其關(guān)系可能比較密切。反之，若其出現(xiàn)的位置具有明顯差異，則其可能關(guān)系較遠。如果把序列看作是基因組內(nèi)的關(guān)鍵詞，則序列在基因組上出現(xiàn)的位置差異則可能說明這些序列之間的功能差異。如果兩個重復(fù)序列在同一基因組上的分布相同，則說明這些序列之間可能有非常密切的關(guān)系。如果不同，則可以通過分布一致性檢驗來判斷其差異是否具有統(tǒng)計學(xué)意義。當其差異具有統(tǒng)計學(xué)意義時，我們一般很想知道其差別到底有多大，以此來推斷兩者之間功能的差異。由于P值受樣本含量的影響，其大小難以反應(yīng)分布之間的差異，因此需要尋找一個合適的指標對分布差異進行量化。相對熵(relative entropy)，又稱KL散度(Kullback-Leibler divergence)，是衡量不同分布之間差異的常用方法。但其有兩個缺陷：①當計算的順序不一樣時，其結(jié)果不同；②對于定量數(shù)據(jù)，一般需將抽樣數(shù)據(jù)進行適當分組再計算，而分組會損失一部分樣本信息。

Kolmogorov-Smirnov檢驗(KS檢驗)一般用于兩種分布之間是否有差異的假設(shè)檢驗[1]。該方法完全避免了相對熵計算所存在的兩個缺陷。那么是否可以利用KS檢驗的統(tǒng)計量對分布之間差異進行量化呢？另外，不同的分布具有不同累積概率曲線，而圖心(centroid)可以視為一個圖形的中心，那么累積概率曲線下圖形的圖心差異也有可能用于衡量分布之間的差異。本文就對這兩個指標進行了一些探討。

研究方法

1.賦值方法

如果將基因組看作[0,1]，則基因組上序列出現(xiàn)的位置即可表示為[0,1]內(nèi)的數(shù)字，如：

基因組序列:CTTTACGCCT堿基T位置:23410數(shù)字轉(zhuǎn)換:0.20.30.41.0

其他堿基以此類推。當基因組很長時，若干堿基組成的序列就可能大量重復(fù)出現(xiàn)在基因組的不同位置上，這些位置也都可以轉(zhuǎn)為[0,1]上的數(shù)字，這樣不同序列的分布就可轉(zhuǎn)化為[0,1]上的各種分布。基因組越長，所需要的精度越高，如基因組長度為106的，則精度為10-6，這樣就可對不同序列的分布進行比較。

2.KS檢驗原理

KS檢驗是比較樣本與理論分布之間，或兩樣本之間，累積概率的最大差異。在R語言[2-3]中KS檢驗的統(tǒng)計量表示為D。該統(tǒng)計量D即累積概率的最大差異，由D值計算兩樣本相同的概率。

3.圖心的計算

圖心橫坐標的計算：作一條垂直于橫軸的線，若該線可以將圖形面積二等分，則該線與橫軸的交點為圖心的橫坐標。圖心縱坐標的計算：作一條垂直于縱軸的線，若該線可以將圖形面積二等分，則該線與縱軸的交點為圖心的縱坐標。如圖1所示，圖中曲線為β(3,3)的累積概率曲線，其曲線下面積分別被兩條平行于x軸和y軸的虛線平分，這兩條虛線的交點即該圖形的圖心。對于理論分布，采用迭代的方法查找圖心坐標，最終坐標的誤差小于10-10。對于抽樣數(shù)據(jù)，則直接計算二分面積的位置。

圖1 圖心示意圖

數(shù)值模擬結(jié)果分析

采用數(shù)值模擬的方法對[0,1]上的若干分布進行了分析。所有隨機數(shù)字的獲取、統(tǒng)計分析和作圖都使用R語言完成。數(shù)據(jù)模擬形式如上述。

1.樣本含量對兩個指標的影響

(1)樣本含量對統(tǒng)計量D的影響

數(shù)值模擬：分別根據(jù)5個β分布生成隨機數(shù)字，然后利用KS檢驗進行均勻分布一致性檢驗。5個分布分別為β(1,1)、β(3,3)、β(1.5,0.8)、β(2,0.5)、β(0.5,2)，分別代表5條不同的序列。β(1,1)即均勻分布。5個分布的累積概率曲線如圖2所示。

圖2 五種分布的累積概率曲線

抽樣的樣本含量從100開始，逐次加100，直至10000。每個樣本含量下進行1000次抽樣。記錄每組抽樣數(shù)據(jù)均勻分布一致性檢驗所得的D值。每個樣本含量根據(jù)所獲得的1000個D值計算5個百分位數(shù)，即2.5%，25%，50%，75%，97.5%。計算結(jié)果如圖3所示。

圖3 樣本含量與D統(tǒng)計量的關(guān)系(100～10000)

在圖3中，每個分布由5條曲線構(gòu)成，從下到上，分別表示D值的5個百分位數(shù)。從圖中可見，隨著樣本含量的增加，D值的離散趨勢逐漸減小，但其集中趨勢并沒有受到明顯影響，且不同的分布集中于不同的位置。表明D值可以看做分布差異的一種量化指標。另外，從圖中可以看到β(2,0.5)和β(0.5,2)兩個分布的數(shù)據(jù)基本重合。

(2)樣本含量對圖心差異的影響

圖心差異采用圖心之間的距離表示，數(shù)值模擬過程同前。計算結(jié)果如圖4所示。

圖4 圖心差異與樣本含量的關(guān)系(100～10000)

從圖中可見，樣本含量對圖心之間距離的影響和對統(tǒng)計量D的影響相似。即隨著樣本含量的增加，圖心距離的離散趨勢逐漸減小，但集中趨勢并沒有受到明顯影響，且不同的分布集中于不同的位置。表明圖心距離也可以看做分布差異的一種量化指標。

2.不同樣本含量下兩指標能判別的最小差異

(1)KS檢驗的統(tǒng)計量D

從表1可見，當樣本含量大于100時，所能判別的最小D值差異已小于0.1。雖然隨著樣本含量的增加，所能判別的最小D值差異不斷減小，但其差異是否具有實際意義可能值得商榷。

(2)圖心距離

計算方法同前，計算結(jié)果見表1。當樣本含量為100時，其最小判別差異在0.02附近。雖然在絕對值上比統(tǒng)計量D小，但需注意理論上此時圖心距離的最大差異小于0.5，而統(tǒng)計量D的最大差異是小于1。另外，從變異程度上看，圖心距離的變異程度約為統(tǒng)計量D的兩倍。

3.在同一坐標系內(nèi)標記不同分布

(1) KS檢驗的統(tǒng)計量D

從前一部分結(jié)果可以看到，雖然統(tǒng)計量D可以用于區(qū)分兩個分布的差異，但正如β(2,0.5)和β(0.5,2)兩個分布和β(1,1)的關(guān)系一樣，當使用一個分布作為基準，存在兩個不同分布的D值相同的情況。為能夠區(qū)分各種不同的分布，可以選擇兩個分布作為基準分布。本實驗選擇β(2,1)和β(1,2)為基準分布。

數(shù)值模擬：分別根據(jù)5個分布生成隨機數(shù)字，然后針對兩個基準分布分別進行分布一致性檢驗。5個分布分別為β(1,1)、β(3,3)、β(1.5,0.8)、β(2,0.5)、β(0.5,2)。抽樣的樣本含量為1000，每個分布進行1000次抽樣，記錄所得D值的平均值。計算結(jié)果如圖5所示。

在圖2中，由于β(1,1)累積概率曲線是β(2,0.5)和β(0.5,2)兩個累積概率曲線的對稱軸，因此使用β(1,1)作為基準就很難區(qū)分二者的差別。從圖5可見，當采用β(2,1)和β(1,2)兩個基準分布時，即可將5個分布分別標記在不同的位置。

(2)圖心

5個分布的累積概率曲線下圖形的圖心采用迭代法計算，計算結(jié)果如圖6所示。由于不同分布具有不同的累積概率曲線，且累積概率曲線只能是單調(diào)遞增曲線。因此不同的分布會對應(yīng)不同的圖心。

實例分析

選取大腸桿菌K-12(MG1655)參考基因組序列，該序列從NCBI網(wǎng)站獲取，全長4639675bp。待測序列為GTGGGCG、GGCGGGG、CCCTCTT、CCCTGCC、CCGCCCT、GCCACCG。6條序列可分別獲取其在基因組中的位置，其重復(fù)出現(xiàn)的次數(shù)即檢驗所需的樣本含量。根據(jù)這些位置進行均勻分布一致性檢驗，并在轉(zhuǎn)換為[0,1]內(nèi)的數(shù)字后計算圖心。

6條序列的統(tǒng)計結(jié)果見表2。在表2中，P和D分別為均勻分布一致性檢驗的P值和統(tǒng)計量，x和y分別為計算所得圖心的橫坐標和縱坐標。其圖心在直角坐標系中的表示見圖7。從圖中可見6條序列大致分為3群。

表1 不同樣本含量下統(tǒng)計量D和圖心距離能判別的最小差異

圖5 利用兩個基準分布和統(tǒng)計量D 顯示五個不同分布的差異

圖6 利用圖心顯示五個不同分布

表2 6條序列的分布差異

圖7 6條序列的圖心表示

討 論

如前所述，基因組上任一序列的分布都可轉(zhuǎn)為[0,1]分布，這樣就可以比較不同序列在基因組上的分布差異。但在計算分布差異時，若使用相對熵來衡量分布之間的差異，由于其對概率差異進行了加權(quán)平均，導(dǎo)致計算順序不同會得出不同的結(jié)果，而這樣的結(jié)果往往不是我們所預(yù)期的。KS檢驗的統(tǒng)計量D以及累積概率曲線下圖形的圖心差異則與此不同。由前面實驗結(jié)果可知，二者隨著樣本含量的增加，離散趨勢逐漸減小，而集中趨勢沒有明顯影響。因此認為這兩個指標可以用于衡量不同分布的差異。

當考察同一基因組上多個序列之間的關(guān)系時，若能將所有序列的分布畫在一個坐標系內(nèi)，則可以使序列分布之間的差異比較直觀。因此本實驗還探討了在同一坐標系標記任意兩個分布差異的可行性。當使用D值差異時，由于對稱性的存在，采用單個基準分布將有可能出現(xiàn)兩個分布D值差異相同的情況。此時采用兩個基準分布可消除這一缺陷，而且可以將不同分布表示在同一坐標系，但該坐標系內(nèi)分布之間的差異不一定等于D值差異。但此時需注意，若采用不同的基準分布，則相當于換了不同的坐標系，分布之間的關(guān)系也將發(fā)生改變。若采用圖心指標，不同的累積概率曲線會對應(yīng)不同的圖心，而且圖心不會像D值那樣隨基準分布而變。另外由于D值反映的是最大累計概率差異，而圖心則與整個圖形有關(guān)，因此采用圖心指標作圖可能會更好地反映分布之間的差異。

在實例分析中，我們選取了同一基因組上的6條序列進行了分析。正如本文開頭所言，分布之間的關(guān)系可能作為功能聯(lián)系密切程度的一個指標，因此該實例分析結(jié)果顯示圖心有可能作為序列功能聚類的一種指標。

最后，由于本研究采用數(shù)值模擬的方法，并不能完全反應(yīng)各種實際情況，有待在實際運用中探討證實。

參 考 文 獻

1.Frank J,Massey Jr.The Kolmogorov-Smirnov test for goodness of fit.Journal of the American Statistical Association,1951,46(253):68-78.

2.鄧宇,彭斌,田考聰.R統(tǒng)計計算環(huán)境簡介.中國衛(wèi)生統(tǒng)計,2008,25(4): 421-422.

3.R Core Team (2013).R: A language and environment for statistical computing.R Foundation for Statistical Computing,Vienna,Austria.ISBN 3-900051-07-0,URL http://www.R-project.org/.

4.Pham-Gia T,Hung TL.The mean and median absolute deviations.Mathematical and Computer Modelling,2001,34(7): 921-936.