基于高通量測(cè)序數(shù)據(jù)的結(jié)構(gòu)變異檢測(cè)方法的研究

張雨豪　王亞?wèn)|

摘要：隨著高通量測(cè)序數(shù)據(jù)技術(shù)的發(fā)展，人類全基因組的測(cè)序成本在不斷降低，測(cè)序速度也有了較為顯著地提升。運(yùn)用生物信息學(xué)的手段處理這些海量基因組數(shù)據(jù)的需求也越來(lái)越迫切，而對(duì)于基因組結(jié)構(gòu)變異的檢測(cè)更是這個(gè)領(lǐng)域的核心內(nèi)容。由高通量測(cè)序數(shù)據(jù)特征入手，介紹了當(dāng)前主流的生物信息學(xué)結(jié)構(gòu)變異檢測(cè)方法，并闡述了有關(guān)基因組結(jié)構(gòu)變異檢測(cè)結(jié)果的評(píng)測(cè)指標(biāo)和手段，最后，結(jié)合個(gè)人基因組的發(fā)展，對(duì)于該領(lǐng)域未來(lái)的發(fā)展提出了改進(jìn)建議。

關(guān)鍵詞：高通量測(cè)序； 結(jié)構(gòu)變異檢測(cè)； 生物信息學(xué)

中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：2095-2163（2013）05-0035-04

0引言

隨著人類基因組計(jì)劃的宣告完成，對(duì)于人類基因組海量數(shù)據(jù)的研究工作也逐步拉開了序幕，這給生物信息學(xué)的發(fā)展提供了很好的發(fā)展機(jī)遇，同時(shí)也帶來(lái)了諸多挑戰(zhàn)。之后的千人基因組計(jì)劃更提供了大量第一手的人類基因組數(shù)據(jù)，這些數(shù)據(jù)既可以為生物學(xué)數(shù)據(jù)處理提供原始輸入，又能為處理生物學(xué)數(shù)據(jù)所得的結(jié)果提供了良好的驗(yàn)證。

當(dāng)利用高通量測(cè)序數(shù)據(jù)來(lái)檢測(cè)結(jié)構(gòu)變異時(shí)，主要有以下幾種思路。第一種是單純依靠覆蓋率信息的方法，這種方法是最早提出檢測(cè)結(jié)構(gòu)變異的方法，現(xiàn)在已很少單獨(dú)利用。第二種主要是依靠雙末端測(cè)序數(shù)據(jù)中非一致序列并通過(guò)聚類來(lái)發(fā)現(xiàn)結(jié)構(gòu)變異信息，這種方法很難發(fā)現(xiàn)具體的機(jī)構(gòu)變異位點(diǎn)信息。第三種方法是利用split read來(lái)精確發(fā)現(xiàn)結(jié)構(gòu)變異，這種方法可以精確發(fā)現(xiàn)結(jié)構(gòu)變異信息，但是重復(fù)序列對(duì)其影響很大。現(xiàn)在大多數(shù)結(jié)構(gòu)變異檢測(cè)軟件都會(huì)集成整合上述幾種方法，取長(zhǎng)補(bǔ)短，并會(huì)相應(yīng)地構(gòu)建一套獨(dú)特的數(shù)據(jù)篩選處理流程，運(yùn)用更快捷更有效的算法，由此而不斷提高基因組結(jié)構(gòu)變異檢測(cè)的能力。

1高通量測(cè)序數(shù)據(jù)介紹

1.1高通量測(cè)序技術(shù)的介紹[HTSS]

對(duì)于人類基因組的全測(cè)序技術(shù)是解決基因組生物信息學(xué)的一個(gè)至關(guān)重要的前提。傳統(tǒng)意義上最著名應(yīng)用、最廣泛的測(cè)序方法是Sanger測(cè)序法[1]，這種方法起源于上世紀(jì)70年代，已經(jīng)過(guò)不斷地改進(jìn)而逐步趨于完善。而且，在2001年得到的第一條人類全基因組序列主要采用的就是這種方法，不過(guò)，這一過(guò)程是通過(guò)全球多個(gè)研究機(jī)構(gòu)的共同努力，且耗費(fèi)了數(shù)年時(shí)間花費(fèi)巨資才完成。

隨著對(duì)于更廉價(jià)、更快捷測(cè)序技術(shù)的需求激增，并經(jīng)過(guò)該領(lǐng)域科學(xué)家的通力協(xié)作，高通量測(cè)序技術(shù)應(yīng)運(yùn)而生。高通量測(cè)序技術(shù)的出現(xiàn)極大地降低了全基因組的測(cè)序時(shí)間以及測(cè)序花費(fèi)。

表1中顯示了幾種高通量測(cè)序技術(shù)的花費(fèi)和優(yōu)缺點(diǎn)，最后一列是第一代Sanger測(cè)序技術(shù)。從表中可以發(fā)現(xiàn)，雖然設(shè)備較貴，但是Illumina測(cè)序儀還是有相對(duì)便宜的價(jià)格和時(shí)間開銷，并且由于Illumina測(cè)序儀可以使用戶根據(jù)其需求生成不同的測(cè)序數(shù)據(jù)，因此，在結(jié)構(gòu)變異檢測(cè)中，原始數(shù)據(jù)大多數(shù)是通過(guò)Illumina測(cè)序儀得到的。

1.2雙末端測(cè)序數(shù)據(jù)介紹

在Illumina測(cè)序儀的結(jié)果中主要會(huì)產(chǎn)生兩種數(shù)據(jù)，一種是單末端數(shù)據(jù)（single end），一種是雙末端測(cè)序數(shù)據(jù)（pair end）。這兩種數(shù)據(jù)分別是根據(jù)不同的生物學(xué)手段得到的，其中雙末端測(cè)序數(shù)據(jù)不僅有短序列（read）信息，而且還包含了插入距離信息，這對(duì)于同一組序列的位置關(guān)系提供了新的一種依靠和保證。在此重點(diǎn)介紹有關(guān)雙末端測(cè)序數(shù)據(jù)的相關(guān)信息。

在雙末端測(cè)序數(shù)據(jù)中，主要包含了相對(duì)基因組的上游序列信息、下游序列信息和插入距離信息，而且數(shù)據(jù)總是成對(duì)出現(xiàn)。由于在處理單末端數(shù)據(jù)時(shí)，主要通過(guò)短序列覆蓋率信息和短序列自身信息來(lái)檢測(cè)結(jié)構(gòu)變異，在利用雙末端測(cè)序時(shí)，不僅可以使用單末端數(shù)據(jù)中的信息，更能通過(guò)對(duì)于插入距離的信息來(lái)有效地檢測(cè)結(jié)構(gòu)變異，因此，在檢測(cè)結(jié)構(gòu)變異的時(shí)候大量采用了雙末端測(cè)序數(shù)據(jù)。

2基因組結(jié)構(gòu)變異類型介紹

隨著人類基因組測(cè)序技術(shù)的進(jìn)步，全基因組的數(shù)據(jù)每天都以海量的規(guī)模在增長(zhǎng)。即使是兩個(gè)不同人種的同性個(gè)體，其基因組之間的差別也是相當(dāng)小的，雖然比例非常低，但是由于人類全基因組有30億堿基序列，所以其數(shù)目仍是非?？捎^的，也正是這些差別導(dǎo)致了人類所有個(gè)體之間的萬(wàn)千差別。因此，開展這些差異的研究對(duì)于無(wú)論是疾病、或是醫(yī)學(xué)等其他領(lǐng)域都有著至關(guān)重要的深遠(yuǎn)意義。

將參考基因組作為比對(duì)依據(jù)，由此得到的差異信息主要分為兩類。第一類是SNP（單核苷酸多態(tài)性）；第二類是結(jié)構(gòu)變異，在結(jié)構(gòu)變異中較為常見的則是如圖1所示的片段刪除和片段插入。

一般來(lái)說(shuō)，將某個(gè)體的基因組序列同參考序列進(jìn)行比對(duì)，如果在一段序列區(qū)間內(nèi)僅有一個(gè)位點(diǎn)不同，就將認(rèn)定為SNP信息。如今的主要檢測(cè)方法是基于貝葉斯估計(jì)進(jìn)行分類，這種方法當(dāng)1-5bp的結(jié)構(gòu)變異時(shí)，就會(huì)產(chǎn)生一個(gè)基于統(tǒng)計(jì)學(xué)的較準(zhǔn)確的結(jié)果，不過(guò)對(duì)于長(zhǎng)序列問(wèn)題的復(fù)雜度卻會(huì)迅速增加，分析難度也會(huì)顯著加大，此時(shí)該方法就不再可取。

3主流結(jié)構(gòu)變異檢測(cè)方法

相較于數(shù)量眾多、性能優(yōu)良的DNA序列比對(duì)工具，結(jié)構(gòu)變異的檢測(cè)工具一方面由于其發(fā)展起步較晚的影響，另一方面則由于結(jié)構(gòu)變異事件的情況相對(duì)SNP和DNA比對(duì)更為復(fù)雜，因此，直到雙末端測(cè)序數(shù)據(jù)的大量使用才出現(xiàn)了很多基于雙末端測(cè)序數(shù)據(jù)的方法。

針對(duì)雙末端測(cè)序數(shù)據(jù)，當(dāng)前主流結(jié)構(gòu)變異檢測(cè)工具都已進(jìn)行了較為充分的利用。SVseq2[2]是一款在低覆蓋率情況下，主要通過(guò)對(duì)雙末端測(cè)序數(shù)據(jù)中產(chǎn)生的split read和非一致的序列的分析來(lái)精確檢測(cè)結(jié)構(gòu)變異的工具。Pindel[3]是一款主要根據(jù)雙末端測(cè)序數(shù)據(jù)來(lái)檢測(cè)大長(zhǎng)度的片段刪除和中長(zhǎng)度的片段插入的結(jié)構(gòu)變異的軟件。Delly[4]是一款使用短插入距離雙末端測(cè)序數(shù)據(jù)、破碎序列片段數(shù)據(jù)相結(jié)合，既能檢測(cè)基因組片段刪除、片段倒置、片段連續(xù)重復(fù)，又能在平衡的基因重排數(shù)據(jù)中檢測(cè)倒易移位事件（reciprocal translocations）的軟件。PRISM[5]是一款既采用了雙末端測(cè)序數(shù)據(jù)的聚類分析又使用了破碎片段序列的比對(duì)分析結(jié)果，并使用改良的Needleman-Wunsch算法而以精確到位點(diǎn)的方式來(lái)檢測(cè)片段刪除、片段插入為主的結(jié)構(gòu)變異信息的軟件。

其中，可利用的雙末端測(cè)序數(shù)據(jù)都是經(jīng)過(guò)BWA等[6]軟件比對(duì)之后的SAM格式文件。文中將可利用的數(shù)據(jù)主要分為兩類：非一致短序列對(duì)（discordant pair）和單映射雙末端測(cè)序數(shù)據(jù)（hanging pair）。如果這兩個(gè)序列片段的映射距離被認(rèn)為是在插入距離的可接受范圍內(nèi)，而且兩個(gè)片段的朝向都沒有發(fā)生改變，即可認(rèn)為這種序列對(duì)為一致的序列對(duì)（concordant pair），該種序列在絕大多數(shù)情況下均不會(huì)被認(rèn)為覆蓋了一個(gè)結(jié)構(gòu)變異。除此之外，其他的雙末端測(cè)序數(shù)據(jù)，無(wú)論是序列朝向問(wèn)題、插入距離問(wèn)題或者CIGAR值異常等問(wèn)題發(fā)生時(shí)，均可認(rèn)為產(chǎn)生的是非一致的序列對(duì)（discordant pair）。除此之外，一種特殊情況，就是雙末端測(cè)序數(shù)據(jù)中僅有一個(gè)序列片段比對(duì)到參考序列上，而另一個(gè)卻未能比對(duì)到參考序列上，由此將沒有CIGAR值，這類特殊的序列可稱為單映射雙末端測(cè)序數(shù)據(jù)對(duì)[7]。

在此，以最近發(fā)布的一篇檢測(cè)結(jié)構(gòu)變異信息的軟件PRISM為例，來(lái)講述檢測(cè)結(jié)果變異的大致流程：

（1）數(shù)據(jù)篩選。DNA全序列比對(duì)是檢測(cè)結(jié)構(gòu)變異的前提條件，因此輸入是標(biāo)準(zhǔn)SAM格式文件。而其中的大多數(shù)序列不會(huì)覆蓋結(jié)構(gòu)變異信息，因此對(duì)于discordant pair和hanging pair的篩選工作將是后續(xù)的研究基礎(chǔ)。[JP2]這里主要還是根據(jù)對(duì)SAM格式中的CIGAR值、插入距離、序列對(duì)正反朝向這三個(gè)方面的判斷，其中之一有異常的狀況，就將其歸類并分別輸出。

（2）數(shù)據(jù)聚類。雖然各個(gè)軟件在聚類時(shí)采用的方法不盡相同，比如：PRISM采用了CNVer[8]作為聚類工具，并采用貪心策略將相似映射距離和朝向的序列對(duì)進(jìn)行劃分聚類，但卻基本都采用了一種聚類手段以實(shí)現(xiàn)對(duì)數(shù)據(jù)進(jìn)行更好地劃分。

（3）Split read比對(duì)。通過(guò)上述兩步得到了Split read數(shù)據(jù)，不同軟件采用各自的算法將split read重新比對(duì)到參考序列上，再結(jié)合上一步的聚類結(jié)果，共同判斷是否可以支持一個(gè)結(jié)構(gòu)變異信息，并分別將其完整記錄。

（4）發(fā)現(xiàn)結(jié)構(gòu)變異信息。大多數(shù)軟件會(huì)根據(jù)序列的質(zhì)量數(shù)，并結(jié)合已經(jīng)設(shè)定的可支持序列數(shù)的具體值來(lái)判定一個(gè)結(jié)構(gòu)變異事件。如PRISM一般采用5作為支持一個(gè)結(jié)構(gòu)變異信息的短序列最小值，由此一個(gè)結(jié)構(gòu)變異事件只有被5個(gè)或者更多短序列支持的時(shí)候，才能判斷是發(fā)生了結(jié)構(gòu)變異。但Pinel卻因其設(shè)計(jì)本就是為低覆蓋率數(shù)據(jù)而實(shí)現(xiàn)準(zhǔn)備的，因而通常選用2作為支持結(jié)構(gòu)變異的短序列最小值。這些參數(shù)既跟數(shù)據(jù)的選用有密切關(guān)系，又與具體采用的聚類方法以及split read比對(duì)算法有關(guān)，因此往往需要不斷調(diào)整以獲得最理想數(shù)據(jù)。

4結(jié)構(gòu)變異評(píng)測(cè)手段

在評(píng)測(cè)中，主要關(guān)注在一定覆蓋率情況下該軟件發(fā)現(xiàn)結(jié)構(gòu)變異的精確率和召回率。首先，可認(rèn)為包含實(shí)際存在的結(jié)構(gòu)變異信息的集合為答案集，在模擬數(shù)據(jù)中，答案集可以

通過(guò)人為手段實(shí)現(xiàn)驗(yàn)證，并通過(guò)對(duì)于植入的結(jié)構(gòu)變異信息的監(jiān)控和追蹤而對(duì)結(jié)構(gòu)變異事件實(shí)時(shí)記錄。同樣，也可認(rèn)為由結(jié)構(gòu)變異檢測(cè)工具報(bào)告的結(jié)構(gòu)變異信息集合為結(jié)果集。共同集則表示既認(rèn)定其存在、而又被結(jié)構(gòu)變異檢測(cè)工具報(bào)告的共同結(jié)果，三者之間的關(guān)系如圖2所示。

其中，召回率=共同集/答案集；精確率=共同集/結(jié)果集。

文中以40x覆蓋下chr1染色體上的模擬數(shù)據(jù)為基準(zhǔn)通過(guò)使用PRISM、Pindel兩種軟件主要針對(duì)片段刪除和片段插入進(jìn)行檢測(cè)。首先，人為地將Venter的結(jié)構(gòu)變異信息植入hg18版本的1號(hào)染色體中，通過(guò)生成的序列與參考基因組相比較生成有關(guān)結(jié)構(gòu)變異信息的答案集，再通過(guò)結(jié)果變異檢測(cè)工具生成結(jié)果集，其后根據(jù)序列名稱間的對(duì)應(yīng)關(guān)系，由此而形成評(píng)估結(jié)果。其具體流程如圖3所示。

這里分別使用15x和40x覆蓋率的chr1染色體作為輸入，對(duì)Pindel和PRISM兩款結(jié)構(gòu)變異檢測(cè)軟件的精確率和召回率進(jìn)行比較。

結(jié)合圖4和圖5可以發(fā)現(xiàn)，當(dāng)覆蓋率較高的時(shí)候，主流的結(jié)構(gòu)變異軟件都具有良好的召回率和精確率。并且隨著覆蓋率的降低，精確率有小幅度的增加，召回率卻有較明顯的減少，因?yàn)槠渲邪嗽S多小型結(jié)構(gòu)變異信息。這些信息采用split read比對(duì)的方法進(jìn)行實(shí)驗(yàn)時(shí)會(huì)有一些噪音數(shù)據(jù)出現(xiàn)，所以，當(dāng)剔除小型結(jié)構(gòu)變異信息時(shí)，這些軟件所能檢測(cè)的數(shù)目會(huì)有所降低，但是結(jié)果卻會(huì)有更好的精確率。

5基因組結(jié)構(gòu)變異發(fā)展方向

雖然當(dāng)今大多數(shù)軟件在一定輸入下都能夠取得較滿意的結(jié)果，但是，采用split read比對(duì)結(jié)合聚類的方法仍有一定的局限性。例如：對(duì)于被測(cè)對(duì)象序列中的大于插入距離的片段插入，現(xiàn)有的各類方法都是很難檢測(cè)出來(lái)的。這種現(xiàn)象主要是因?yàn)?，沒有一個(gè)雙末端測(cè)序序列能夠覆蓋這個(gè)片段插入事件。此外，由于聚類方法的采用，將很難精確發(fā)現(xiàn)一個(gè)結(jié)構(gòu)變異事件的具體精確位點(diǎn)。而且，現(xiàn)在的方法主要是針對(duì)片段刪除和片段插入的情況，對(duì)于其他結(jié)構(gòu)變異信息卻還未得到很優(yōu)異的結(jié)果。

6結(jié)束語(yǔ)

隨著高通量測(cè)序技術(shù)的迅猛發(fā)展，對(duì)于人類基因組的研究正在逐步邁入后基因組時(shí)代，對(duì)于測(cè)序數(shù)據(jù)的處理工作正變得尤為重要，利用生物信息學(xué)手段檢測(cè)人類基因組結(jié)構(gòu)變異這一核心領(lǐng)域的研究也將越來(lái)越深入和細(xì)致。本文從測(cè)序數(shù)據(jù)的介紹入手，對(duì)主流結(jié)構(gòu)變異檢測(cè)軟件的流程和方法進(jìn)行了分析，并且根據(jù)PRISM和Pindel在不同覆蓋率情況下的輸出結(jié)果對(duì)檢測(cè)工具的性能評(píng)估進(jìn)行了敘述。