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

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

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

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

研究方法

1.賦值方法

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

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

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

2.KS檢驗(yàn)原理

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

3.圖心的計(jì)算

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

圖1 圖心示意圖

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

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

1.樣本含量對(duì)兩個(gè)指標(biāo)的影響

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

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

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

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

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

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

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

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

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

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

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

(1)KS檢驗(yàn)的統(tǒng)計(jì)量D

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

(2)圖心距離

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

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

(1) KS檢驗(yàn)的統(tǒng)計(jì)量D

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

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

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

(2)圖心

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

實(shí)例分析

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

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

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

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

圖6 利用圖心顯示五個(gè)不同分布

表2 6條序列的分布差異

圖7 6條序列的圖心表示

討 論

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

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

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

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

參 考 文 獻(xiàn)

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)計(jì)計(jì)算環(huán)境簡(jiǎn)介.中國(guó)衛(wèi)生統(tǒng)計(jì),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.