基因組高通量測序數(shù)據(jù)結(jié)構(gòu)變異識(shí)別的方法分析

1、基因組高通量測序數(shù)據(jù)結(jié)構(gòu)變異識(shí)別的方法分析0 引 言遺傳變異是生命的基本特征，遺傳變異與表型差異之間的關(guān)系，是現(xiàn)代生物學(xué)的一個(gè)基本問題。由基因決定生物體的遺傳特征和主要個(gè)體差異的觀念正在逐漸改變，過去幾年的許多研究顯示，基因組中大尺度的結(jié)構(gòu)變異（structural variation， SV）與個(gè)體的表型差異和疾病等有一定的關(guān)聯(lián)【1】。近年的研究更進(jìn)一步，很多疾病和表型差異與基因組的結(jié)構(gòu)變異有關(guān)系，人們已經(jīng)逐漸認(rèn)識(shí)到結(jié)構(gòu)變異的普遍性和重要性。1 結(jié)構(gòu)變異及其重要作用基因組結(jié)構(gòu)變異通常指基因組DNA分子上整片的片段差異，包括保持?jǐn)?shù)量平衡的倒位（inversion）和易位（translocati

2、on）變化以及導(dǎo)致數(shù)量失衡的插入（insertion）、缺失（deletion）和重復(fù)（duplication）變化，其中數(shù)量失衡的變化也稱為拷貝數(shù)變異（copy number variation， CNV）。一般認(rèn)為結(jié)構(gòu)變異的片段，長度約在一千到幾百萬個(gè)堿基之間，一致性在90%以上。2004年，全球范圍內(nèi)數(shù)個(gè)人類基因組計(jì)劃;研究基地意外地發(fā)現(xiàn)，表型正常的人群中，不同個(gè)體間在某些基因的拷貝數(shù)上存在差異，結(jié)構(gòu)變異的發(fā)現(xiàn)及其對(duì)基因雜合性的重要作用，徹底改變了人們對(duì)基因型研究的認(rèn)識(shí)【5】。個(gè)體基因組之間的差異，除少數(shù)SNP和短序列的插入/缺失（indel）變化之外，大部分以結(jié)構(gòu)變異為主【6】。相關(guān)研

3、究也證實(shí)，在玉米自交系中存在異常高的結(jié)構(gòu)變異【7】。雖然普遍認(rèn)為當(dāng)基因組規(guī)模相對(duì)較小時(shí)，發(fā)生結(jié)構(gòu)變異的數(shù)量會(huì)較少，但對(duì)大豆的研究卻特別發(fā)現(xiàn)，相對(duì)于大豆基因組的大小和其中SNP的數(shù)量來說，結(jié)構(gòu)變異的發(fā)生率較高。而且結(jié)構(gòu)變異在大豆的不同染色體間的分布也不均勻，在基因分布較密的區(qū)域發(fā)生較多，并且與植物的抗性表現(xiàn)有一定的關(guān)聯(lián)。2 相關(guān)研究最近幾年隨著下一代測序技術(shù)中高通量測序平臺(tái)在應(yīng)用上的推廣普及，基因組序列數(shù)據(jù)的獲取速度越來越快。這些平臺(tái)在測序中都會(huì)產(chǎn)生大量的相互交疊的對(duì)端讀片段（paired-end read）。被測序的DNA分子將打斷為很多碎片，其插入長度約2k150kbp，設(shè)備僅測每個(gè)碎片兩端

4、長度為251kbp的對(duì)端讀片段。使用這樣的序列數(shù)據(jù)識(shí)別結(jié)構(gòu)變異，在長度和精確度上均會(huì)高于芯片技術(shù)的使用效果，因而適用于各種類型結(jié)構(gòu)變異的識(shí)別以及新類型的發(fā)現(xiàn)。采用測序技術(shù)的結(jié)構(gòu)變異識(shí)別可以使用重頭組裝和重測序兩種方法。其中，重頭組裝方法最為理想，但由于目前測序技術(shù)的限制，完整的全基因組組裝結(jié)果并非最優(yōu)，不能很好地滿足先組裝、再識(shí)別的要求；而重測序方法則更加實(shí)際，可將得到的大量讀片段定位到參考基因組并實(shí)施比較分析。讀片段定位（read mapping）是序列比對(duì)的基本任務(wù)，對(duì)其研究已可堪稱深入，但通用的讀序列定位算法應(yīng)用于結(jié)構(gòu)變異檢測時(shí)，卻沒有充分考慮斷點(diǎn)（breakpoint）處的比對(duì)，從而會(huì)

5、導(dǎo)致斷點(diǎn)預(yù)測不準(zhǔn)確甚至結(jié)構(gòu)變異檢測失敗。采用重測序技術(shù)的結(jié)構(gòu)變異分析方法主要有三種，具體來說，分別是：（1）覆蓋度（read depth或depth of coverage， DoC）分析法，根據(jù)在參考基因組上讀片段的覆蓋度，判斷是否有拷貝數(shù)變化的結(jié)構(gòu)變異，如CNV-seq等，這類方法無法識(shí)別拷貝數(shù)不變的倒位和易位變化；（2）對(duì)端定位（paired-end mapping， PEM）法，根據(jù)對(duì)端讀片段插入距離的一致性，判斷是否發(fā)生結(jié)構(gòu)變異，如VariationHunter等；（3）讀片段分割法（split read），利用讀片段在參考基因組上定位時(shí)的斷點(diǎn)判斷結(jié)構(gòu)變異，如Pindel，但由于NG

6、S測序的讀長較短，不利于發(fā)現(xiàn)較大的結(jié)構(gòu)變異。也有將這些信息結(jié)合的方法，如CNVer等只關(guān)注了CNV一種結(jié)構(gòu)變異；Break- Dancer和GASV使用了高置信度讀片段，忽略了與重復(fù)區(qū)域有關(guān)的低置信度讀片段；inGAP-sv和Gnome Strip在精確性和特異性兩方面都較好，但卻都基于二倍體人類基因組設(shè)計(jì)和實(shí)現(xiàn)。因此，目前采用NGS數(shù)據(jù)進(jìn)行基因組結(jié)構(gòu)變異的檢測與分析方法還不夠成熟，缺少完整的分析、發(fā)現(xiàn)和評(píng)價(jià)體系。3 結(jié)構(gòu)變異識(shí)別方法本文的基本的結(jié)構(gòu)變異識(shí)別方法如圖1所示，其原理可概括為：以目標(biāo)物種基因組序列為參考，使用多個(gè)供體全基因組高通量末端配對(duì)測序讀片段數(shù)據(jù)，采用PEM的方法識(shí)別結(jié)構(gòu)變異

7、，并利用對(duì)端信息進(jìn)行讀片段分隔方式的斷點(diǎn)預(yù)測。本文設(shè)計(jì)基于單個(gè)參考基因組和多供體基因組重測序數(shù)據(jù)的結(jié)構(gòu)變異識(shí)別方法。從識(shí)別方法的角度，本文可以分為三個(gè)主要方面：數(shù)據(jù)預(yù)處理，結(jié)構(gòu)變異識(shí)別和精確的斷點(diǎn)預(yù)測。下面將給出其具體實(shí)現(xiàn)。3.1 數(shù)據(jù)預(yù)處理數(shù)據(jù)預(yù)處理的主要任務(wù)是讀片段定位，不一致和歧義讀對(duì)（discordant and ambiguous pairs）的過濾，產(chǎn)生的中間結(jié)果為不同類型的候選讀片段集，用于支持可能的多種結(jié)構(gòu)變異的識(shí)別和斷點(diǎn)精確預(yù)測。在此，對(duì)其過程解析可做如下展現(xiàn)。（1）讀片段定位在參考基因組上定位讀片段，對(duì)每個(gè)供體的讀片段分別向參考基因組上定位，相應(yīng)得到各自的定位結(jié)果。定位時(shí)除

8、需考慮一般定位算法涉及的測序錯(cuò)誤、soft clipping等問題之外，還需要考慮兩方面問題。具體描述為： 多處匹配：若同一個(gè)讀片段可能在多處得到較好的匹配，則可能因?yàn)橛锌截悢?shù)變化，需采用基于種子及擴(kuò)展的比對(duì)算法，保存并返回多個(gè)高分比對(duì)位置； 對(duì)讀片段分割的懲罰：若某條讀片段剛好在斷點(diǎn)處，由于對(duì)空位的罰分，可能造成錯(cuò)誤定位或無法定位，因此要調(diào)整比對(duì)算法的罰分策略，允許一處不罰分的長空格。（2）不一致讀片段對(duì)的過濾由測序技術(shù)可知，同一碎片雙端測序的讀片段對(duì)，在參考基因組上的定位如若符合一定的方向和距離要求即稱為一致的（concordant），否則稱為不一致的（discordant），如圖2所示。

9、不一致的原因之一就是由于參考基因組和供體基因組之間的結(jié)構(gòu)差異，主要分為三類，將其綜合概述如下：方向異常的，可能因?yàn)榈刮唬枚ㄎ坏姆较蚩珊Y選出； 距離異常，可能對(duì)應(yīng)插入、缺失或易位，根據(jù)測序參數(shù)中平均插入長度L和標(biāo)準(zhǔn)差σ來處理，當(dāng)距離超過L+3σ錯(cuò)誤！未指定書簽。和小于L-3σ錯(cuò)誤！未指定書簽。的都視為異常并篩選出； 懸掛對(duì)（hanging pairs），即對(duì)端讀片段中有一條未能成功定位的，可能對(duì)應(yīng)較大的結(jié)構(gòu)變異，由定位結(jié)果可直接篩選出來。（3）歧義讀片段對(duì)過濾若一條讀片段在多處均有較好定位，又可將其稱為歧義定位對(duì)。可能由于潛在的結(jié)構(gòu)變異造成，同樣按照類似不

10、一致讀對(duì)過濾的方法進(jìn)行篩選。只是要依據(jù)在第（1）步中讀片段定位算法所提供的結(jié)果，即將歧義讀片段對(duì)的多處匹配結(jié)果統(tǒng)一、完整地實(shí)現(xiàn)返回。（4）合并多供體的候選讀片段由于文庫制備和測序?qū)嶒?yàn)的不同，每個(gè)供體測序的基本參數(shù)也是不同的，如插入長度、讀片段長度等，所以讀片段的定位和過濾即需分別設(shè)計(jì)運(yùn)行。而讀片段的定位結(jié)果參考同一基因組，因此可將過濾得到的候選讀片段按類型合并，并保存供體來源信息，并在發(fā)生分歧時(shí)提供分類依據(jù)。而相應(yīng)于懸掛對(duì)，則由于只有一條讀片段能夠成功定位，需要對(duì)未成功定位的讀片段連同歧義讀片段進(jìn)行聚類，以發(fā)現(xiàn)較大的結(jié)構(gòu)變異。3.2 結(jié)構(gòu)變異識(shí)別（1）基本結(jié)構(gòu)變異識(shí)別由結(jié)構(gòu)變異導(dǎo)致的不一致讀對(duì)

11、如圖3所示。對(duì)于基本的插入、缺失和倒位等，需根據(jù)定位的座位進(jìn)行分類并計(jì)算讀對(duì)共同的比對(duì)得分。對(duì)每個(gè)座位，若有3組以上來自同一供體的同類讀對(duì)，則認(rèn)為是合法的結(jié)構(gòu)變異；或者，另有多個(gè)供體且每個(gè)有2組以上的同類讀對(duì)，也將認(rèn)定為是合法的結(jié)構(gòu)變異。（2）讀片段覆蓋度分析利用多供體數(shù)據(jù)，在提高等效的覆蓋度方面具有顯著優(yōu)勢。覆蓋度分析的原理如圖4所示，當(dāng)覆蓋度發(fā)生異常時(shí)，可能因?yàn)榛蚪M的重復(fù)片段或存在較大的結(jié)構(gòu)變異。分析定位到參考基因組上讀片段的數(shù)量，就需要估計(jì)全部成功定位的讀片段在參考基因組上的分布和特征，再根據(jù)總體及局部的分布以判斷是否呈現(xiàn)有明顯的變化，并在一定大小的窗口內(nèi)計(jì)算讀片段覆蓋度，而當(dāng)信號(hào)發(fā)生

12、相對(duì)顯著的增加和減少時(shí)，即可判斷有結(jié)構(gòu)變異發(fā)生。覆蓋度與對(duì)端讀片段信息的結(jié)合分析，可發(fā)現(xiàn)較大結(jié)構(gòu)變異的原理如圖5所示。具體地，由于幾個(gè)對(duì)端讀片段的距離較大并且同時(shí)此處的覆蓋度超過正常值。在覆蓋度異常的區(qū)域，尋找方向不一致的對(duì)端讀片段，并試圖利用復(fù)制該區(qū)域的方法，將不一致的對(duì)端讀片段一致排列，以發(fā)現(xiàn)結(jié)構(gòu)變異。（3） 復(fù)雜結(jié)構(gòu)變異識(shí)別與重復(fù)片段的分析對(duì)于超過插入長度的結(jié)構(gòu)變異，會(huì)導(dǎo)致懸掛對(duì)和歧義定位。因此對(duì)這些讀片段，使用貪心算法將距離和方向類似的片段實(shí)現(xiàn)聚類后，即需進(jìn)行基于de Bruijn圖的拼接。在此圖中對(duì)數(shù)據(jù)錯(cuò)誤和重復(fù)序列進(jìn)行處理，并根據(jù)路徑與參考基因組實(shí)行比較，從而分析可能存在的結(jié)構(gòu)變異

13、。如圖6所示，利用多供體信息，可以相互驗(yàn)證路徑中的環(huán)路，是否是因?yàn)橹貜?fù)片段。若構(gòu)成路徑的片段相對(duì)供體來源是均勻的，可以認(rèn)為是在供體中普遍存在的重復(fù)序列；否則即認(rèn)為是由于結(jié)構(gòu)變異，再與參考基因組相比較，并確認(rèn)結(jié)構(gòu)變異的類型和位置。3.3 結(jié)構(gòu)變異的斷點(diǎn)預(yù)測斷點(diǎn)預(yù)測需要讀片段分割分析和對(duì)端讀片段分析相結(jié)合，來自多供體的對(duì)端讀片段構(gòu)成差異較大，更利于斷點(diǎn)位置的精確預(yù)測。如果有多個(gè)讀片段在斷點(diǎn)處的一側(cè)發(fā)生完美匹配，則可以通過該組讀片段直接預(yù)測斷點(diǎn)位置，誤差在讀片段間交疊的長度范圍內(nèi)；而如果在讀片段覆蓋范圍之外，則需要使用對(duì)端讀片段的插入長度進(jìn)行分析。因?yàn)椴迦腴L度L有一定誤差范圍，假設(shè)在Lmin和Lma

14、x之間且對(duì)端讀片段的左右斷點(diǎn)分別是（x， y），那么可以估計(jì)斷點(diǎn)發(fā)生的位置（a， b）滿足 Lmin ≤ （a-x） + （y-b） ≤ Lmax，即（a， b）落在x， y， Lmin， Lmax構(gòu)成的梯形區(qū)域當(dāng)中。若有多個(gè)對(duì)端讀片段同時(shí)支持該區(qū)域的結(jié)構(gòu)變異，則該斷點(diǎn)所在位置需要在多個(gè)梯形的交疊區(qū)域中，如圖7中分別為形成插入和缺失時(shí)的情況。此外，為了能夠處理超過插入長度結(jié)構(gòu)變異的斷點(diǎn)，利用斷點(diǎn)兩翼序列的特點(diǎn)，根據(jù)斷點(diǎn)處的微同源性，分析已知結(jié)構(gòu)變異斷點(diǎn)處的微同源特征和兩翼序列的特征，再進(jìn)行斷點(diǎn)預(yù)測的。4 實(shí)驗(yàn)結(jié)果及分析本文實(shí)驗(yàn)利用模擬數(shù)據(jù)生成供體的測序片段，根據(jù)dbVar對(duì)已知結(jié)構(gòu)

15、變異大小和分布，在大豆1號(hào)染色體上構(gòu)造了插入和缺失類型的結(jié)構(gòu)變異，如表1所示。由于易位、倒位和重復(fù)事件現(xiàn)有數(shù)據(jù)中也很少，而且通常各種軟件對(duì)這些事件識(shí)別的準(zhǔn)確率都不高，所以沒有進(jìn)行模擬。本文即比較選擇了主流的幾個(gè)結(jié)構(gòu)變異識(shí)別工具，分別是BreakDancer、Pindel和CREST。對(duì)缺失事件的識(shí)別結(jié)果如表2所示。綜合該結(jié)果可知，總體而言，Pindel、BreakDancer和本文方法的識(shí)別結(jié)果各有優(yōu)勢，不可概述；分別地，Pindel和本文方法在小型的缺失事件上結(jié)果較好，而在100bp以上的范圍SVdetect、BreakDancer和本文方法的結(jié)果則更佳。其中1001000bp的兩組結(jié)果中，

16、本文方法識(shí)別效果較好，而在100500bp和 1 000b兩組結(jié)果中當(dāng)屬本文方法的效能最優(yōu)。主要原因是對(duì)于利用split-read方式的識(shí)別工具，其大小產(chǎn)生split-read的數(shù)量較少，而對(duì)于不一致讀片段對(duì)的識(shí)別工具其長度卻恰好落入插入長度的誤差范圍內(nèi)。綜合分析實(shí)驗(yàn)結(jié)果，本文方法在不同的條件下均可以得出較好的識(shí)別結(jié)果。本文方法能夠準(zhǔn)確識(shí)別出大部分的結(jié)構(gòu)變異事件，因而可以用于NGS數(shù)據(jù)中結(jié)構(gòu)變異的識(shí)別。參考文獻(xiàn)：【1】CRADDOCK N， Wellcome Trust Case Control Consortium， CRADDOCK N， HURLES ME， et al. Genome-

17、wide association study of CNVs in 16，000 cases of eight common diseases and 3，000 shared controls. Nature， 2010， 464（7289）： 713–20.【2】ABECASIS G R， ALTSHULER D， AUTON A， et al. A map of human genome variation from population-scale sequencing. Nature， 2010， 467（7319）：1061–73.【3】周雪崖， 張學(xué)工.

18、基于拷貝數(shù)變異的遺傳關(guān)聯(lián)研究. 科學(xué)通報(bào)， 2011， 56（6）：370-182.【4】TREANGEN T J， SALZBERG S L. Repetitive DNA and next-generation sequencing： computational challenges and solutions. Nature reviews. Genetics， 2012， 13（1）：36–46.【5】FEUK L， CARSON A R， SCHERER S W. Structural variation in the human genome. Nature reviews. Genetics， 2006