改進(jìn)的基因拷貝數(shù)變異檢測算法

李 平，楊洪斌，吳 悅

(上海大學(xué)計(jì)算機(jī)工程與科學(xué)學(xué)院，上海 200072)

1 概述

人類基因組在 2004年意外發(fā)現(xiàn)正常個(gè)體間部分基因的拷貝數(shù)存在差異，其中，2個(gè)科研小組率先公布了正常個(gè)體人類全基因組的拷貝數(shù)變異情況?？截悢?shù)變異(Copy Number Variation, CNV)是指較之于參照基因組，DNA片段位點(diǎn)缺失或復(fù)制大于 1 Kb至1 Mb的結(jié)構(gòu)變異[1]。拷貝數(shù)變異形式多樣，包括DNA片段的缺失、嵌入、復(fù)制和復(fù)合多位點(diǎn)變異等。此后，相繼有報(bào)道指出拷貝數(shù)變異對正常人群間的遺傳變異作用。文獻(xiàn)[2]觀察基因中影響和調(diào)節(jié)細(xì)胞生長、新陳代謝的幾種基因都與CNV有一定的聯(lián)系?？截悢?shù)變異檢測的成果不斷豐富，其中，拷貝數(shù)變異圖是一新的人類基因圖譜，它補(bǔ)充了被稱為“生命之書”中殘缺的章節(jié)，目的是使研究人員可以分辨基因的增加、減少和變化，從新的角度解釋基因與疾病的關(guān)聯(lián)[3]。

迅速發(fā)展的全基因組拷貝數(shù)掃描平臺使 CNV用于全基因組關(guān)聯(lián)分析成為可能，產(chǎn)生了大量不同的研究算法。常見的技術(shù)平臺有基于大插入片段的比較基因組雜交(Comparative Genomic Hybridization, CGH)，代表性寡核苷酸微陣列分析(Representational Oligonucleotide Microarray, ROMA)，和單核苷酸多態(tài)性(Single Nucleotide Polymorphisms, SNP)芯片[4]。SNP 芯片是目前通量最高、使用范圍最廣的基因組拷貝數(shù)變異分析平臺。它不僅可以有效地提高拷貝數(shù)變異檢測精度，而且在一定的范圍內(nèi)可以增加樣本，不需要額外的費(fèi)用，其樣本也可反復(fù)使用，增加效率。在高通量的基因檢測中，數(shù)據(jù)量更為龐大，檢測精確度更高[5]。最常見的 SNP陣列數(shù)據(jù)來自商業(yè)公司 affymetirix和illumina，他們所出售的 SNP陣列數(shù)據(jù)在檢測 CNV區(qū)域有著一定的競爭力[6]。

cnvPartition算法是基于 SNP芯片上的常見算法[7]，可以有效地檢測出隱藏性小片段變異 CNV和新的CNV。但是由于SNP數(shù)據(jù)位點(diǎn)較多、數(shù)據(jù)量龐大、運(yùn)行負(fù)擔(dān)較大且速度緩慢，普通地運(yùn)行檢測算法也需要很好的軟硬件支持。因此，改進(jìn) cnvPartition算法，提高其速度又能保持較好的運(yùn)行結(jié)果是十分必要的。cnvPartition算法是基于三元分割的原理，為了找出顯著變異點(diǎn)，需要不斷進(jìn)行測試。二元分割是一種簡單快速的分割方法，可以快速地找出分割斷點(diǎn)，但忽略了隱藏在大片段中的小片段檢測。因此，本文基于分割原理[8]提出一種改進(jìn) cnvPartition算法的基因拷貝數(shù)變異檢測算法，有效地結(jié)合了二元分割和三元分割[9]。

2 相關(guān)方法

2.1 cnvPartition算法

cnvPartition算法是使用 LRR(Log R Ratio))和BAF(B Allele Frequency)數(shù)據(jù)定義出基因組上拷貝數(shù)發(fā)生變異的區(qū)域。cnvPartition算法通過循環(huán)查找基因數(shù)據(jù)斷點(diǎn)的方法去確定拷貝數(shù)變異和其存在的位置，其目的是在染色體內(nèi)得到均值不同于其余區(qū)域的CNV區(qū)域。該算法可以分為：檢測斷點(diǎn)得到變異段和計(jì)算斷點(diǎn)間片段的CNV值。在查找斷點(diǎn)時(shí)，通過循環(huán)比較，確定可疑變異區(qū)域，進(jìn)而查找確定變異斷點(diǎn)[10]。該方法可以有效地檢測出隱藏在大變異片段中的小型片段，但是需要重復(fù)比較，對軟硬件要求較高，對高通量數(shù)據(jù)來說運(yùn)行速度較慢。

2.2 二元分割

在二元分割中，使用Z-test公式檢查出斷點(diǎn)所在的位置，然后確定斷點(diǎn)所做的分割區(qū)域[11]。斷點(diǎn)將染色體片段分成2個(gè)部分，若斷點(diǎn)所分割的2片段的差異度大于片段分割閾值，則肯定有一部分為變異片段。通過循環(huán)查找所有顯著斷點(diǎn)，以至于找到所有CNV變異片段。其基本過程為：在片段中根據(jù)Z-test公式查找其最大斷點(diǎn)，判斷此斷點(diǎn)分成的2片段差異度是否大于片段閾值，以確定是否為顯著斷點(diǎn)。在顯著斷點(diǎn)分割形成的兩側(cè)片段中，再次查找出所有顯著斷點(diǎn)，最后根據(jù)顯著斷點(diǎn)所劃分的變異區(qū)域計(jì)算其CNV值。

由于染色體上基因位點(diǎn)個(gè)數(shù)在數(shù)萬個(gè)之上，因此使用Z-test算法進(jìn)行計(jì)算，易發(fā)生小數(shù)被大數(shù)吞沒而無法檢測出隱藏在大片段中的小變異片段，有失準(zhǔn)確性[12]。

2.3 環(huán)狀二元分割算法

環(huán)狀二元分割(Circular Binary Segmentation,CBS)算法是應(yīng)用在CGH平臺上的拷貝數(shù)變異檢測算法，它是對二元分割方法的修改以提高CNV的檢測效果[12]。為檢測出隱藏在大片段中的小型變異區(qū)域，CBS算法將染色體兩端連接起來形成環(huán)狀分割，通過生成重置排列方法比較非正常數(shù)據(jù)片段。CBS可以有效地檢測出小片異區(qū)域，但其運(yùn)行速度易受到數(shù)據(jù)量大小限制。在 CBS算法中，重置排列計(jì)算時(shí)間與數(shù)組中位點(diǎn)數(shù)目的二次方位增長，對于高通量的 CGH數(shù)據(jù)將是巨大的負(fù)擔(dān)。因此，這種方法在有效地檢測CNV區(qū)域上也需要相應(yīng)改進(jìn)。

定義1(Z-test檢驗(yàn)) Z檢驗(yàn)是一般用于大樣本(即樣本容量大于 30)平均值差異性檢驗(yàn)的方法。它是用標(biāo)準(zhǔn)正態(tài)分布的理論來推斷差異發(fā)生的概率，從而比較2個(gè)樣本的平均數(shù)的差異是否顯著。

如果檢驗(yàn)一個(gè)已知樣本平均數(shù) X與一個(gè)已知的總體平均數(shù)(μ0)的差異是否顯著。Z值計(jì)算公式為：

其中，X是檢驗(yàn)樣本的平均數(shù)；μ0是已知總體的平均數(shù)；S是樣本的標(biāo)準(zhǔn)差；n是樣本容量。

如果檢驗(yàn)2組樣本平均數(shù)的差異性，從而判斷它們各自代表的總體的差異是否顯著。Z值計(jì)算公式為：

其中，X1和 X2是樣本 1和樣本2的平均數(shù)；S1和S2是樣本 1和樣本 2的標(biāo)準(zhǔn)差；n1和 n2是樣本 1和樣本2的容量。

定義2(斷點(diǎn)) 假設(shè)X1, X2,…是一組隨機(jī)變量，若X1, X2,…, Xv有符合分布函數(shù) F1, Xv,…符合另一分布函數(shù) F2，而 F1和 F2不同，則索引位置 v點(diǎn)即為斷點(diǎn)[12]。

定義3(BAF) BAF表示基因位點(diǎn)中攜帶了B等位基因的雜交樣本比例。在一個(gè)正常樣本中，BAF有3種表達(dá)值，即0.0、0.5和 1，這分別表示位點(diǎn)的基因型為AA、AB和BB。當(dāng)計(jì)算出BAF的值與這3個(gè)值有偏差時(shí)，說明該位置可能有拷貝數(shù)變異存在[13]。

定義4(LRR) LRR是基因位點(diǎn)的實(shí)際觀測信號強(qiáng)度與期望信號強(qiáng)度比的lbn值。若LRR不等于0，則說明拷貝數(shù)有變化。這個(gè)函數(shù)可以測試出發(fā)生變異片段的信號強(qiáng)度與正常信號強(qiáng)度的偏異量，估計(jì)出變異強(qiáng)度[14]。

3 本文算法

本文將二元算法的快速性與 cnvPartition算法的準(zhǔn)確性結(jié)合起來，可以有效地減少高通量數(shù)據(jù)中的循環(huán)比較，又保證對大型數(shù)據(jù)的詳細(xì)檢測，既提高速度又可保持良好的檢測結(jié)果。

3.1 本文算法結(jié)構(gòu)

本文算法總體流程如圖1所示。

圖1 本文算法總體流程

本文算法主要分為3個(gè)階段：(1)在片段中使用二元分割方法查找斷點(diǎn)且判斷該斷點(diǎn)是否為顯著斷點(diǎn)，是顯著斷點(diǎn)則轉(zhuǎn)到階段(2)，否則轉(zhuǎn)到階段(3)。(2)根據(jù)顯著斷點(diǎn)劃分的片段中用二元分割方法再次查找斷點(diǎn)，直到無顯著斷點(diǎn)為止。(3)在無顯著斷點(diǎn)的片段中，用 cnvPartition算法查找是否有小變異片段，若有小變異片段則對新劃分的片段轉(zhuǎn)階段(1)再次按二元分割查找數(shù)斷點(diǎn)；否則將最后生成的無變異片段的片段計(jì)算其CNV值。

3.2 顯著斷點(diǎn)查找

在染色體上，每個(gè)位點(diǎn)的信號強(qiáng)度沿著染色體上的位置連續(xù)自然分布，用來檢測拷貝數(shù)變異的陣列數(shù)據(jù)，也是按位點(diǎn)索引排序的 LRR值。若染色體上有多個(gè)斷點(diǎn)，必然相應(yīng)位點(diǎn)的 LRR值也發(fā)生變化。因此，本文算法最主要的內(nèi)容就是根據(jù)染色體上的LRR值分布查找出所有顯著CNV斷點(diǎn)。

根據(jù)Z檢驗(yàn)方法得到斷點(diǎn)，將斷點(diǎn)所劃分的兩區(qū)域比較其數(shù)據(jù)差異度，若差異度大于規(guī)定片段分割閾值 segment-threshold，則此斷點(diǎn)為該片段內(nèi)的顯著斷點(diǎn)。根據(jù)Z檢驗(yàn)方法，此閾值可以應(yīng)用尾部檢測概率進(jìn)行計(jì)算，即P值小于P(0.01)。其基本過程為：

(1)定義染色體片段X1, X2,…, Xm上T值最大的位點(diǎn)索引位置v為此區(qū)域內(nèi)的斷點(diǎn)B1，T的公式為：

T=max1≤i＜m|Ti|

(2)計(jì)算由v所劃分的染色體上2片段的差異度，并與片段分割閾值相比較。根據(jù) Z檢驗(yàn)方法，2片段的差異度計(jì)算方法的公式為：

其中，Yij和 Zij表示 2片段 i到 j區(qū)域及 1到 i及j到m區(qū)域數(shù)據(jù)的平均值；Sij表示對應(yīng)片段數(shù)據(jù)的平均差。若差異度大于 segment-threshold，則此斷點(diǎn) B1為顯著斷點(diǎn)，根據(jù)顯著斷點(diǎn)所分割的區(qū)域轉(zhuǎn)步驟(2)，否則轉(zhuǎn)步驟(3)。

(3)在B1的2側(cè)區(qū)域內(nèi)根據(jù)步驟(1)描述查找斷點(diǎn)B2和B3，根據(jù)步驟(2)描述分別比較B2B1段與該染色體區(qū)域中其余片段的差異度和B1B2段與其余片段的差異度，選擇差異度值最大的區(qū)域端點(diǎn)定義為P1、P2。

(4)再次比較 P1和 P2所劃分的 2相鄰片段的差異度是否大于segment-threshold以確定P1和P2是否是顯著斷點(diǎn)。若P1和P2中至少有一個(gè)顯著斷點(diǎn)，則轉(zhuǎn)步驟(1)，否則轉(zhuǎn)步驟(5)。

(5)該片段內(nèi)沒有顯著斷點(diǎn)，也即此區(qū)域內(nèi)沒有拷貝數(shù)變異片段。計(jì)算此片段的CNV值[12]。

3.3 CNV值計(jì)算

對于無顯著斷點(diǎn)的區(qū)域內(nèi)，計(jì)算其CNV值，檢測出CNV變異區(qū)域。根據(jù)染色體位點(diǎn)檢測的LRR值和BAF值，與14個(gè)標(biāo)準(zhǔn)的正常位點(diǎn)值相比較計(jì)算每位點(diǎn)的初步拷貝值。對每一分割片段區(qū)域內(nèi)，對每一假定的拷貝值(0～4)，計(jì)算所有位點(diǎn)的 lbLk值之和，和值最大的 k值即為這一區(qū)域的拷貝數(shù)值。對于CNV值不為2的區(qū)域，對每個(gè)CNV值計(jì)算一個(gè)置信度。置信度的計(jì)算形式為區(qū)域?yàn)樗衛(wèi)bLK之和減去區(qū)域內(nèi)所有位點(diǎn)的lbL2值[11]。

3.4 時(shí)間復(fù)雜度分析

在 cnvPartition算法中，確定一個(gè)顯著斷點(diǎn)至少需要比較7次。若一個(gè)區(qū)域內(nèi)只有n個(gè)顯著斷點(diǎn)，至少需要比較 7×n次，同時(shí)為檢查 n個(gè)斷點(diǎn)所劃分的n+1個(gè)區(qū)域內(nèi)是否還有斷點(diǎn)，共需要7×n+7×(n+1)次比較。在本次segcnv算法中，n個(gè)顯著斷點(diǎn)的區(qū)域中需要比較至少n+7×(n+1)次。當(dāng)n值不斷變大時(shí)，其算法運(yùn)行時(shí)間提高7×(2n+1)/(8n+7)≈2倍。

4 實(shí)驗(yàn)結(jié)果與分析

本文將描述使用上述算法的運(yùn)算結(jié)果，在Windows系統(tǒng)下使用 C++語言編寫運(yùn)行程序。通過infinium HD芯片得到的白血病樣本數(shù)據(jù)陣列來進(jìn)行模擬本文算法。這里選擇了8份3號染色體下同種樣本的數(shù)據(jù)去考慮在不同數(shù)據(jù)量下 cnvPartition和本文算法在運(yùn)行時(shí)間和運(yùn)行結(jié)果上的比較。

本文實(shí)驗(yàn)從運(yùn)行時(shí)間和運(yùn)行結(jié)果兩方面考查本文算法的效率。圖 2為本文算法與cnvPartition算法的相對運(yùn)行時(shí)間比較。

圖2 2種算法的相對運(yùn)行時(shí)間比較

從圖2可以看出，cnvPartition算法在數(shù)據(jù)量較少時(shí)，可以維持著較短的運(yùn)行時(shí)間。但是隨著數(shù)據(jù)量的增加，運(yùn)行時(shí)間不斷增長，以至于在成千上萬的高精度數(shù)據(jù)中，成為限制本文算法效率的重要原因。由二元分割算法和 cnvPartition算法結(jié)合的本文算法，其運(yùn)行時(shí)間低于 cnvPartition算法。當(dāng)數(shù)據(jù)量不斷增大時(shí)，其運(yùn)行時(shí)間并沒有表現(xiàn)出指數(shù)增長，僅表現(xiàn)出一定的比率增長。相對于原算法，本文算法在高數(shù)據(jù)量運(yùn)行時(shí)可以保持較低的運(yùn)行時(shí)間，在位點(diǎn)數(shù)增加時(shí)，最高可提高將近50%的速度。這對于提高cnvPartition算法的效率是至關(guān)重要的。

2種算法分割片段數(shù)目比較如圖3所示。

圖3 2種算法分割片段數(shù)目比較

由圖3可以看出，由于cnvPartition算法在選擇顯著變異斷點(diǎn)時(shí)，進(jìn)行詳細(xì)檢查，其可以有效地將大片段數(shù)據(jù)根據(jù)其所在的位置進(jìn)行有效的劃分，得到更多小型片段。且本文算法即保留了 cnvPartition良好的仔細(xì)檢測模式，同時(shí)，對于顯著斷點(diǎn)的檢測提高了檢測速度，維持了良好的檢測算法，非常有利于檢測出小片段的CNV變異。

5 結(jié)束語

基因拷貝數(shù)變異作為影響人類疾病的一個(gè)重要因素已經(jīng)得到了廣泛關(guān)注。高通量的檢測平臺和有效的檢測算法已經(jīng)被提出，目前著重于提高算法的效率。研究發(fā)現(xiàn)，CNV的基本類型都為小型變異片段，因此，有效地檢測出小型變異對CNV的研究具有重要意義。本文提出一種改進(jìn)的基因拷貝數(shù)變異檢測算法。通過改進(jìn) cnvPartition算法，減少部分內(nèi)部循環(huán)比較，提高運(yùn)行的速度，使用斷點(diǎn)分割的原理，又重視對于小變異片段的檢測。實(shí)驗(yàn)結(jié)果表明，該算法在高數(shù)據(jù)量運(yùn)行時(shí)可以保持較低的運(yùn)行時(shí)間，保持cnvPartition算法在檢測隱藏性小片段變異上的良好效果。今后將著重研究染色體上位點(diǎn)數(shù)量龐大所造成計(jì)算速度下降的問題。

[1]吳志俊, 金 瑋.拷貝數(shù)變異: 基因組多樣性的新形式[J].遺傳, 2009, 31(4): 339-347.

[2]陳執(zhí)中.人類基因組拷貝數(shù)變異與新藥研究開發(fā)[D].上海: 中國科學(xué)院上海冶金研究所, 2000.

[3]譚 琪, 曾凡一.遺傳變異的又一來源: 拷貝數(shù)變異[J].生物技術(shù)通迅, 2009, 20(3): 396-398.

[4]孫玉琳, 劉 飛, 趙曉航.拷貝數(shù)變異的全基因組關(guān)聯(lián)分析[J].生物化學(xué)與生物物理進(jìn)展, 2009, 36(8): 968-977.

[5]Daniel A P, Jennie M L, Frank J S, et al.High-resolution Genomic Profiling of Chromosomal Aberrations Using Infinium Whole-genome Genotyping[J]. Genome Research, 2006, 16(9): 1136-1148.

[6]Winchester l, Christopher Y, Ragoussis J.Comparing CNV Detection Methods for SNP Arrays[J].Briefings in Functional Genomics and Proteomics, 2009, 8(5):353-366.

[7]Myungjin M, Jaegyoon A, Chihyun Y Y.A Computational Approach to Detect CNVs Using High-throughput Sequencing[C]//Proc.of the 9th IEEE International Conference on Bioinformatics and Bioengineering.[S.l.]:IEEE Press, 2009.

[8]Chihyun P, Youngmi Y, Jaegyoon A, et al.A Novel Approach to Detect Copy Number Variation Using Segmentation and Genetic Algorithm[C]//Proc.of ACM Symposium on Applied Computing.New York, USA:ACM Press, 2009.

[9]Erez B Y, Yonina E.A Fast and Flexible Method for the Segmentation of aCGH Data[J].Bioinformatics, 2008,24(16): 139-145.

[10]Gaellla M, Benjanim R S, Montserrat G C, et al.Assessment of Copy Number Variation Using the Illumina Infinium 1M SNP-array: A Comparison of Methodological Approaches in the Spanish Bladder Cancer/EPICURO Study[J].Human Mutation, 2011, 32(2): 240-248.

[11]Illumina Corporation.DNA Copy Number and Loss of Heterozygosity Analysis Algorithms[EB/OL].(2010-11-21).http://www.illumina.com.

[12]Adam B, Olshen E, Venkatraman S.Circular Binary Segmentation for the Analysis of Array-based DNA Copy Number Data[J].Biostatistics, 2004, 5(4): 557-572.

[13]Venkatraman E S, Adam B O.A Faster Circular Binary Segmentation Algorithm for the Analysis of Array CGH Data[J].Bioinformatics, 2007, 23(6): 657-663.

[14]Illumina Corporation.DNA Copy Number Analysis Algorithms[EB/OL].(2010-08-21).http://www.pasteur.fr/ip/portal/action/WebdriveActionEvent/oid/01s-00003f-00.