改進的基因拷貝數變異檢測算法

李 平，楊洪斌，吳 悅

(上海大學計算機工程與科學學院，上海 200072)

1 概述

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

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

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

2 相關方法

2.1 cnvPartition算法

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

2.2 二元分割

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

由于染色體上基因位點個數在數萬個之上，因此使用Z-test算法進行計算，易發生小數被大數吞沒而無法檢測出隱藏在大片段中的小變異片段，有失準確性[12]。

2.3 環狀二元分割算法

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

定義1(Z-test檢驗) Z檢驗是一般用于大樣本(即樣本容量大于 30)平均值差異性檢驗的方法。它是用標準正態分布的理論來推斷差異發生的概率，從而比較2個樣本的平均數的差異是否顯著。

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

其中，X是檢驗樣本的平均數；μ0是已知總體的平均數；S是樣本的標準差；n是樣本容量。

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

其中，X1和 X2是樣本 1和樣本2的平均數；S1和S2是樣本 1和樣本 2的標準差；n1和 n2是樣本 1和樣本2的容量。

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

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

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

3 本文算法

本文將二元算法的快速性與 cnvPartition算法的準確性結合起來，可以有效地減少高通量數據中的循環比較，又保證對大型數據的詳細檢測，既提高速度又可保持良好的檢測結果。

3.1 本文算法結構

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

圖1 本文算法總體流程

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

3.2 顯著斷點查找

在染色體上，每個位點的信號強度沿著染色體上的位置連續自然分布，用來檢測拷貝數變異的陣列數據，也是按位點索引排序的 LRR值。若染色體上有多個斷點，必然相應位點的 LRR值也發生變化。因此，本文算法最主要的內容就是根據染色體上的LRR值分布查找出所有顯著CNV斷點。

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

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

T=max1≤i＜m|Ti|

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

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

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

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

(5)該片段內沒有顯著斷點，也即此區域內沒有拷貝數變異片段。計算此片段的CNV值[12]。

3.3 CNV值計算

對于無顯著斷點的區域內，計算其CNV值，檢測出CNV變異區域。根據染色體位點檢測的LRR值和BAF值，與14個標準的正常位點值相比較計算每位點的初步拷貝值。對每一分割片段區域內，對每一假定的拷貝值(0～4)，計算所有位點的 lbLk值之和，和值最大的 k值即為這一區域的拷貝數值。對于CNV值不為2的區域，對每個CNV值計算一個置信度。置信度的計算形式為區域為所有lbLK之和減去區域內所有位點的lbL2值[11]。

3.4 時間復雜度分析

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

4 實驗結果與分析

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

本文實驗從運行時間和運行結果兩方面考查本文算法的效率。圖 2為本文算法與cnvPartition算法的相對運行時間比較。

圖2 2種算法的相對運行時間比較

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

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

圖3 2種算法分割片段數目比較

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

5 結束語

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

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

[2]陳執中.人類基因組拷貝數變異與新藥研究開發[D].上海: 中國科學院上海冶金研究所, 2000.

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

[4]孫玉琳, 劉 飛, 趙曉航.拷貝數變異的全基因組關聯分析[J].生物化學與生物物理進展, 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.