拷貝數(shù)變異檢測算法優(yōu)化研究

林勇

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

關(guān)鍵詞： 拷貝數(shù)變異檢測;算法優(yōu)化;裂讀法; 配對末端讀數(shù)

【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? 引言

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

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

1? 基于OEM和SRM的拷貝數(shù)變異檢測優(yōu)化

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

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

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

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

1.1? OEM檢測

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

（1）拷貝數(shù)增加時(shí)，若為散落重復(fù)，跨越拷貝數(shù)邊緣的Paired-End Read匹配到參考序列上，一端會匹配至拷貝數(shù)序列區(qū)域內(nèi)，另一端則可能會匹配至正常區(qū)域。這里判斷該情況的標(biāo)準(zhǔn)為：

1.2? SRM檢測

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

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

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

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

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

2? 實(shí)驗(yàn)驗(yàn)證與分析

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

2.1? 模擬數(shù)據(jù)驗(yàn)證

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

本文選用的用于實(shí)驗(yàn)驗(yàn)證的工具為CNV-Seq和CNVnator與本文提出的優(yōu)化方法進(jìn)行比較，由于本文方法基于CNVnator，簡稱為CNVNOP，實(shí)驗(yàn)結(jié)果如表1所示。

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

2.2? 真實(shí)數(shù)據(jù)驗(yàn)證

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

3? 總結(jié)與展望

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

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

參考文獻(xiàn)

[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.