拷貝數變異檢測算法優化研究

林勇

摘? 要： 拷貝數變異與多種復雜疾病密切相關，具有重要的研究意義。本文利用基于測序數據的拷貝數變異檢測過程中丟棄的不匹配讀數據，采用裂讀法和單端匹配法對已有檢測算法的結果進行過濾優化。模擬和實驗數據檢測結果表明，本文方法優化后能得到了更高的檢測性能。

關鍵詞： 拷貝數變異檢測;算法優化;裂讀法; 配對末端讀數

【Abstract】： Copy number variation is closely related to many complex diseases and has important research significance. In this paper， the mismatched read data discarded in the process of copy number variation detection based on sequencing data were used to filter and optimize the results of existing detection algorithms by split reading method and one end matching method. The detection results of simulation and experimental data showed that the proposed method can achieve higher detection performance after optimization.

【Key words】： Copy number variant detection; Algorithm optimization; Split read; Paired-end read

0? 引言

拷貝數變異（Copy Number Variation， CNV）是大小超過1 kb的亞顯微突變，表現為DNA片段缺失或重復，重復又分為散落重復和串聯重復?？截悢底儺悘V泛地存在于人類基因組上[1]，雖然拷貝數變異的頻率較低，但累積的堿基數量卻大大超過單核苷酸多態[2-4]。許多研究證明了拷貝數變異與包括乳腺癌、孤獨癥、肥胖和骨質疏松癥等多種疾病相關[5-8]。隨著下一代測序技術（Next-Generation Sequencing， NGS）的出現和發展，拷貝數變異檢測也出現了全新的方法，下一代測序技術能夠檢測DNA序列上堿基量級的序列信息，通過對測序數據進行分析和檢測，拷貝數變異長度、位置等信息就的可能檢測獲得?，F有的基于測序數據的拷貝數變異檢測方法主要是基于讀深度法（Read-Depth Method），它的基本原理是基于拷貝數的變化將引起測序數據匹配區域的讀數據累計量會出現顯著地增加或減少。如圖1所示。目前，研究人員已經開發了許多基于度深度法的檢測工具，比如：CNV-seq[9]，ReadDepth[10]，CONTRA[11]和CNVnator[12]。

目前常規方法檢測過程中，首先將讀數據與 參考序列進行比對，由于讀數據存在少量的測序錯誤以及序列上的SNP問題，比對時允許少量堿基的錯配，當大于錯配閾值時，讀數據將被丟棄。實際上這些被丟棄的讀數據并不都是“垃圾數據”，本文對讀數據匹配不成功的成因進行過分析，加以分類，記錄裂讀匹配（Split Read Match，SRM）和單端匹配（One End Match，OEM）的數據，對匹配位置進行聚類后對拷貝數檢測的初步結果進行? 篩選，從而得到更精確的檢測結果，從而提高檢測性能。

1? 基于OEM和SRM的拷貝數變異檢測優化

本文檢測拷貝數變異的整體流程如圖2所示，該方法僅針對末端配對（Paired-End）的Read數據進行CNV檢測，目前通用的測序儀產生的Read以Paired-End數據為主，因此本文算法適用范圍較廣。首先對原始Read數據進行質控，過濾掉Q值低于20的Read，然后將質控結束后得到的read數據與參考序列進行比對，本文采用的比對工具為BWA[13]，比對得到的數據通過samtools轉換成SAM文件，便于分析比對結果。比對得到的結果分為兩部分：

（1）匹配成功的Read數據用于初步的拷貝數變異檢測，本文采用的拷貝數檢測工具為CNVnator[14]，CNVnator通過檢測序列的深度分布情況，結合GC校正、均值漂移法、多帶寬分割等方法進行拷貝數變異的檢測，是目前較常用的一種基于讀深度法檢測拷貝數變異的工具。在檢測過程中，我們將信號強度閾值降低，這有利于檢測出更多可能的CNV。檢測獲得的結果作為候選結果用于后期篩選;

（2）匹配成功的Read數據，又分為兩種，一種是OEM Read，是指Paired-End Read的兩個單端皆能與參考序列匹配但匹配位置距離與Paired-End的插入距離（insert length）差異很大，因此在序列比對的時候被丟棄的Paired-End Read;另一種是非OEM Read，對于這類Read我們進行裂讀處理，然后將裂讀后的數據進行再比對，獲得裂讀匹配位置。

接著我們對上述獲得的兩種匹配位置進行聚類，最后使用聚類結果完成對候選CNV的篩選得到最后的CNV檢測結果。下面我們對OEM、SRM和聚類分析的工作原理和技術細節進行詳細的描述。

1.1? OEM檢測

單端匹配（OEM）要求Paired-End Read的兩端都能匹配到參考序列上，由圖4可知，在檢測拷貝數變異時根據拷貝數增加和拷貝數減少時情況有所不同。為了方便討論，這里做如下定義：Paired-End Read雙端間的距離稱為插入長度（insert length），所有的Read的插入長度滿足正態分布，令均值為IL，標準差為σ。進行匹配時，令前端匹配至參考序列上的位置為P1，后端匹配至參考序列的位置為P2，拷貝數區域的估計長度為CNV_Len。

（1）拷貝數增加時，若為散落重復，跨越拷貝數邊緣的Paired-End Read匹配到參考序列上，一端會匹配至拷貝數序列區域內，另一端則可能會匹配至正常區域。這里判斷該情況的標準為：

1.2? SRM檢測

當被檢測序列中存在結構變異時，位于結構變異邊緣的讀數據與參考序列無法整體匹配，但該讀數據的前綴或者后綴則可以。而CNV屬于結構變異，拷貝數增加本質上是一種插入變異，而拷貝數減少則是刪除變異。SRM的核心思想是通過將原有的讀數據分裂，將前綴和后綴匹配至參考序列，根據這些前后綴的匹配信息完成拷貝數變異檢測的優化[15]，圖3是裂讀法檢測結構變異中插入和刪除的示意圖。本算法通過讀數據分裂，根據前后綴的匹配信息驗證現有工具檢測得到的CNV，對于信號強度較弱的讀計數異常進行過濾，從而提高檢測精度。

為了方便描述，這里定義讀數據Read的讀長為RL，前綴為Readpre，長度為PL，PL=（RL-IL）/2，IL是可變參數，為Read中間部分忽略子串的長度，本方法中IL=RL/5，讀數據的后綴為Readsuf，長度為SL，SL=PL。對于拷貝數增加，CNV邊界將出現Readpre能匹配至參考序列非CNV所在區域而Readsuf能匹配至CNV區域的情況，而對于拷貝數的減少，則出現前后綴都匹配至非CNV所在區域的情況?；谝陨戏治?，本文將非OEM的丟棄讀數據分別取其PL長度的前后綴與參考序列進行比對，獲得匹配坐標，對于多點匹配的情況采用加權隨機分配方法，然后記錄匹配坐標。

1.3? 聚類分析與候選CNV的篩選

為了實現OEM和SRM檢測獲得的匹配坐標對候選基因的篩選，首先將上面記錄的坐標信息進行聚類，本文采用的聚類方法為k均值聚類，由于OEM和SRM都是基于CNV邊界附近的read信息來完成的，故k均值聚類的聚類中心數量k設置為候選CNV數量的兩倍，分別對應于候選CNV的兩個邊界，且初始聚類中心指定為相應的邊界坐標。OEM記錄的兩個匹配坐標僅使用其中靠近聚類中心的一個實施聚類。聚類分析結束后對每個聚類包含的點進行計數，表示為該聚類的置信度。

本文選用CNVnator作為產生候選CNV的工具，為了提高CNV檢測的靈敏度，在候選CNV檢測階段，降低讀深信號（RD signal）的閾值以檢測出更多的候選CNV。然后根據候選CNV的質量值e-val結合聚類置信度進行最終CNV的篩選，篩選原則是：對于質量值較高的CNV直接保留，而質量值較低的CNV則進一步考慮相應的聚類置信度，置信度較高的也保留作為最終CNV，否則丟棄，不納入最終結果。

2? 實驗驗證與分析

為了測試本算法的有效性，本文同時采用模擬數據和真實的測序數據進行驗證。采用模擬數據可以定制在序列中注入的拷貝數變異所在的位置以及拷貝數量，能夠對算法的檢測性能進行量化分析和比較;采用真實測序數據能夠有效地反映算法應用于實際數據的檢測效果，更接近于實際應用。作為比較，使用CNV檢測工具CNV-Seq和CNVnator與本文提出的優化方法進行比較，檢測性能利用精確性（Precision）、靈敏性（Sensitivity）和F1分數這三個參數進行評估。精確性表示檢測結果被判定為正確的數據占全部檢測結果的百分比;靈敏性表示的是判定為正確的檢測結果占對應實際全部的CNV的百分比;F1分數是精確性和靈敏性的綜合指標。令TP表示真陽性（True Position），表示檢測結果和實際全部的CNV數據都是正類。FP表示假陽性（False Position），表示工具的檢測結果是正類但實際是負類。FN表示假陰性（False Negative），表示工具的檢測結果是負類但實際結果是正類。三個指標的公式分別為：精確性Pre=TP/（TP+FP），靈敏性Sen=TP/（TP+FN），F1=2*Pre*Sen/（Pre+Sen）。

2.1? 模擬數據驗證

模擬的CNV數據根據真實測序數據中拷貝數變異數據的不同類型和不同長度的特點進行設計，一共設計了兩種拷貝數變異子類型：INS，DEL，INS對應于拷貝數的增加而DEL對應的是拷貝數的減少。同時設計了2個不同的長度區間：500—1000bp，1000—10000bp。每種變異長度都設計100變異。參考序列是hg19的11號染色體，模擬CNV數據使用的工具是Svsim （https：//github.com/GregoryFaust/ SVsim），對該工具進行改寫，生成包含準確變異信息的序列文件，以fasta格式進行存儲。由變異序列文件再利用ART[16]生成Paired-End讀文件，讀長為100bp，以FASTQ格式進行存儲。本文共生成三種測序深度的讀數據，分別為30X，50X和70X。

本文選用的用于實驗驗證的工具為CNV-Seq和CNVnator與本文提出的優化方法進行比較，由于本文方法基于CNVnator，簡稱為CNVNOP，實驗結果如表1所示。

由表1實驗結果可以看出，無論是拷貝數增加和拷貝數減少的情況，CNVNOP方法比CNVnator和CNVSeq在精確性、靈敏度和F1值都有較好的表現，而當測序深度增大時，相應工具的檢測性能都有所提升，相對而言50X的測序深度較30X的檢測精度提高更明顯一些。本文方法首先降低了CNVnator的信號篩選閾值已獲得更多的候選結果，然后通過ORM和SRM進行篩選，實驗結果表明該方法能夠得到更高的檢測性能。

2.2? 真實數據驗證

真實數據采用千人基因組項目中具有較高覆蓋度的樣本NA19240，這里僅取第1、2、10、11、20和22號染色體進行檢測，從DVG數據庫中查詢可知這6條染色體的INS共1330個，DEL共777個。比較三種工具后的實驗結果表2。由表2可以得到與模擬實驗相同的的結論，優化后的方法具有較好的CNV檢測性能。

3? 總結與展望

本文提出了一種利用比對過程中丟棄的讀數據對拷貝數檢測進行優化的方法，以CNVnator為基礎產生候選拷貝數，基于SRM和OEM的聚類結果進行篩選，有效地提高了檢測精度和靈敏性，模擬數據和真實數據實驗的結果驗證了該方法的有效性。本文方法還具有兩個重要的潛在優點：第一，它具有較好的通用性，本文方法中使用的CNVnator可以使用其它方法替代，例如：Speedseq，ReadDepth等。隨著技術的進步，可能會有更高檢測性能的方法，同樣可以進行候選CNV檢測方法的替換然后采用本文的優化方法;第二，通過本文方法可能計算出拷貝數變異的軟切位點的位置，由OEM和SRM技術原理可知，其匹配位置通常會有兩個，其中一個位于拷貝數變異區域，另一個則是增加的拷貝數區域位置，這對下游的功能分析能提供較好的幫助。本文方法也存在著一些不足之處，由于裂讀匹配時需要將單個Read進行分割，因此要求的讀長不能太小，否則匹配時非常容易產生多位置匹配，降低算法的檢測精度;另外由于OEM利用了配對末端的insert size的信息，對于一些測序儀產生的非paired-end Read數據，本文方法也無法處理。

本文的優化方法能夠有效提高檢測精度，但還有值得完善的地方，首先本文采用的候選拷貝數的工具只采用一種工具進行檢測，實際上這里可以多使用幾種方法同時進行，利用群體優勢來獲得高可靠性的候選變異;其次在額外信息的使用中，除了OEM和SRM外，還有一些技術也可以被引入進行優化，例如：local assembly方法，在拷貝數變異區域附近進行局部拼接有助于獲得進一步的變異細節，也能夠提高檢測的精度。這些方法將在納入將來的研究中，進一步提高檢測性能。

參考文獻

[1] McCarroll， S. A.， Extending genome-wide association studies to copy-number variation[J]. Hum Mol Genet， 2008. 17（R2）： p. R135-42.

[2] Hinds， D. A.， et al.， Common deletions and SNPs are in linkage disequilibrium in the human genome[J]. Nat Genet， 2006. 38（1）： p82-5.

[3] Redon， R.， et al.， Global variation in copy number in the human genome[J]. Nature， 2006. 444（7118）： p. 444-54.

[4] Wong， K. K.， et al.， A comprehensive analysis of common copy-number variations in the human genome[J]. Am J Hum Genet， 2007. 80（1）： p. 91-104.

[5] Bochukova， E. G.， et al.， Large， rare chromosomal deletions associated with severe early-onset obesity[J]. Nature， 2010. 463（7281）： p. 666-70.

[6] Diskin， S. J.， et al.， Copy number variation at 1q21. 1 associated with neuroblastoma[J]. Nature， 2009. 459（7249）： p. 987-91.

[7] Fanciulli， M.， et al.， FCGR3B copy number variation is associated with susceptibility to systemic， but not organ- specific， autoimmunity[J]. Nat Genet， 2007. 39（6）： p. 721-3.

[8] Stefansson， H.， et al.， Large recurrent microdeletions associated with schizophrenia[J]. Nature， 2008. 455（7210）： p. 232-6.

[9] McKernan， K. J.， et al.， Sequence and structural variation in a human genome uncovered by short-read， massively parallel ligation sequencing using two-base encoding[J]. Genome Res， 2009. 19（9）： p. 1527-41.

[10] Miller， C. A.， et al.， ReadDepth： a parallel R package for detecting copy number alterations from short sequencing reads[J]. PLoS One， 2011. 6（1）： p. e16327.

[11] Li， J.， et al.， CONTRA： copy number analysis for targeted resequencing[J]. Bioinformatics， 2012. 28（10）： p. 1307-13.

[12] Abyzov， A.， et al.， CNVnator： an approach to discover， genotype， and characterize typical and atypical CNVs from family and population genome sequencing[J]. Genome Res， 2011. 21（6）： p. 974-84.

[13] H.， L.， Aligning sequence reads， clone sequences and assembly contigs with BWA-MEM[J]. eprint arXiv： 1303. 3997， 2013.

[14] Abyzov， A.， et al.， CNVnator： an approach to discover， genotype， and characterize typical and atypical CNVs from family and population genome sequencing[J]. Genome Res， 2011. 21（6）： p. 974-84.

[15] Wang， J.， et al.， CREST maps somatic structural variation in cancer genomes with base-pair resolution[J]. Nat Methods， 2011. 8（8）： p. 652-4.

[16] Huang， W.， et al.， ART： a next-generation sequencing read simulator[J]. Bioinformatics， 2012. 28（4）： p. 593-4.