高通量測(cè)序數(shù)據(jù)中Split read映射方法的研究

姜玥 王亞?wèn)|

摘要：高通量測(cè)序技術(shù)的快速發(fā)展與廣泛應(yīng)用為計(jì)算機(jī)科學(xué)帶來(lái)了新的挑戰(zhàn)，read的映射問(wèn)題是其中非常重要的一個(gè)部分。Split read是一類特殊的read，其出現(xiàn)通常是由基因組中的結(jié)構(gòu)變異造成的。這類read在映射中不再保持連續(xù)序列的形式，而是包含了一定長(zhǎng)度的空位，因此具有較高的映射難度。提出一種利用雙末端測(cè)序數(shù)據(jù)的映射結(jié)果來(lái)指導(dǎo)split read映射的方法，這種方法可以使split read的映射難度不再與其所包含的空位數(shù)量相關(guān)，從而降低了映射過(guò)程中的搜索空間，提高映射效率。

關(guān)鍵詞：split read； 映射； 高通量測(cè)序； 生物信息學(xué)

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

0引言

人類基因組計(jì)劃的完成為人類基因組的研究提供了一套參考基因組序列，大大地簡(jiǎn)化了人類個(gè)體基因組的序列研究，因?yàn)椴煌祟悅€(gè)體基因組序列之間有著極高的相似性，現(xiàn)在的研究主要專注于個(gè)體基因組序列與參考基因組序列的差異，這大大地簡(jiǎn)化了研究的過(guò)程。而高通量測(cè)序技術(shù)的不斷發(fā)展，則為人類基因組研究提供了有力數(shù)據(jù)支持。為了利用高通量測(cè)序數(shù)據(jù)，需要將上億的測(cè)序短序列（read）映射到參考基因組序列上，這些read當(dāng)中大部分可以以連續(xù)序列的形式被映射，但是仍有一部分read由于個(gè)體基因組序列與參考基因組序列的差異，會(huì)在映射中包含一段空位，這樣的read稱為split read，其映射相比于第一類read是更為困難的。Split read的映射往往可以顯示個(gè)體基因組中變異區(qū)域的序列信息，對(duì)研究更快速、準(zhǔn)確的split read映射方法有著重要的意義。

1基本概念

1.1高通量測(cè)序數(shù)據(jù)

高通量測(cè)序是一種測(cè)序DNA序列的技術(shù)。在測(cè)序過(guò)程中，將完整的樣本DNA序列打碎，從中篩選出滿足特定長(zhǎng)度（通常為數(shù)百bp）的片段，然后在每個(gè)片段的一端或兩端各讀取一段長(zhǎng)度為數(shù)十至數(shù)百bp的序列。這些讀取出的序列長(zhǎng)度通常遠(yuǎn)遠(yuǎn)小于被測(cè)樣本DNA序列的長(zhǎng)度，但是高通量測(cè)序技術(shù)可以同時(shí)讀取大量這樣的短序列，使得短序列總長(zhǎng)度達(dá)到樣本DNA長(zhǎng)度的數(shù)倍至數(shù)十倍，從而使獲得樣本DNA序列成為可能。

1.2Read與split read

在高通量測(cè)序中，從打碎的DNA片段上讀取出來(lái)的短序列稱為read。Read是被測(cè)DNA序列的一個(gè)短片段，單個(gè)的read序列長(zhǎng)度遠(yuǎn)遠(yuǎn)短于被測(cè)DNA序列的長(zhǎng)度，但是通過(guò)將大量read映射到參考基因組序列的方式，就可以獲得被測(cè)DNA的序列內(nèi)容，如圖1所示。測(cè)序時(shí)所讀取的read是一段連續(xù)的序列，但是由于DNA結(jié)構(gòu)變異的存在，一些read在映射結(jié)果中不再保持連續(xù)的形式，而是包含了空位，這樣的read稱為split read。

1.3雙末端測(cè)序

在高通量測(cè)序過(guò)程中，從打碎的DNA片段的兩端讀取序列的方法稱為雙末端測(cè)序。雙末端測(cè)序中獲得的讀取自同一片段的一對(duì)read稱為一個(gè)read pair。理論上，如果被測(cè)DNA序列與參考基因組序列完全相同，read pair被映射到參考基因組之后，其中的兩個(gè)read之間的距離與被測(cè)時(shí)DNA片段的長(zhǎng)度應(yīng)當(dāng)是相同的。但是由于被測(cè)DNA與參考基因組序列存在差異，特別是由于結(jié)構(gòu)變異的存在，read pair映射后其一對(duì)read之間的距離會(huì)與被測(cè)的DNA片段長(zhǎng)度產(chǎn)生明顯的差異。

2Deletion對(duì)附近read 與read pair映射所造成的影響Deletion是一種常見的結(jié)構(gòu)變異形式，表現(xiàn)為被測(cè)DNA序列相比參考基因組序列缺失了部分序列。由于這種變異的存在，其附近的read與read pair在映射過(guò)程中會(huì)發(fā)生異常，如圖2所示。從圖2中可以看出，由于deletion的存在（黑色短線段），跨過(guò)deletion的read pair（左）在映射后兩個(gè)read之間的距離要長(zhǎng)于被測(cè)時(shí)兩個(gè)read之間的距離，這個(gè)距離的差異恰好是deletion的長(zhǎng)度。而跨過(guò)deletion邊界的read（右）在映射時(shí)則會(huì)包含與deletion長(zhǎng)度相同的一段空位，形成split read。

3利用read pair映射分析指導(dǎo)split read映射的方法目前的read映射方法出于運(yùn)行效率的考慮，都會(huì)限制映射結(jié)果中所允許的空位數(shù)量與長(zhǎng)度[1-3]。有一些利用雙末端測(cè)序數(shù)據(jù)特性而特別為split read映射所設(shè)計(jì)的映射方法，利用read pair中一個(gè)映射較好的read作為基點(diǎn)，在臨近的一段區(qū)間為另一個(gè)映射效果不好或者無(wú)法連續(xù)映射的read進(jìn)行允許較多空位的映射[4]。這樣的方法存在著映射效果與搜索空間相關(guān)，映射難度大，效率低等問(wèn)題，如圖3所示。

為了改進(jìn)這些不足，本文提出一種利用deletion附近的read pair的映射結(jié)果來(lái)指導(dǎo)split read映射的方法。從圖2中可以看出，受到deletion影響的read pair，雖然其一對(duì)read之間的映射距離發(fā)生了異常，但兩個(gè)read的映射位置距離deletion的邊界并不遠(yuǎn)。通過(guò)將這樣存在映射異常的read pair按照映射位置與每對(duì)read之間的距離進(jìn)行聚類，可以大致獲得deletion邊界的位置。由于split read的映射實(shí)際上只需要deletion邊界處的一小段序列，而與deletion序列本身無(wú)關(guān)，因此可以每個(gè)聚類結(jié)果中的兩處deletion邊界位置為基點(diǎn)，各選擇一段固定長(zhǎng)度的序列作為參考序列進(jìn)行split read映射，選擇序列的長(zhǎng)度只要確保可以包含deletion的分界點(diǎn)即可（圖4上半部分）。通過(guò)這樣的方式，split read的映射將不再與deletion本身的長(zhǎng)度相關(guān)，因?yàn)閰⑴csplit read映射的參考序列只是deletion邊界處固定長(zhǎng)度的兩段序列的組合，其選取與deletion本身的長(zhǎng)度無(wú)關(guān)。

接下來(lái)，需要將每個(gè)聚類結(jié)果附近映射效果較差或無(wú)法映射的read提取出來(lái)，這些read可能是受到了每個(gè)聚類結(jié)果所對(duì)應(yīng)的deletion的影響而無(wú)法實(shí)現(xiàn)良好的映射，因其是候選的split read。將這些read向組合的參考序列映射需要一種序列映射算法，本文提出一種Needleman-Wunsh算法[5， 6]的變種算法來(lái)完成split read映射。變種算法同樣是一種動(dòng)態(tài)規(guī)劃算法，其遞歸表達(dá)式為：

其中：

db是由兩段參考基因組序列組成的橫向序列，段序列的長(zhǎng)度分別為m1和m2。qr是由read序列構(gòu)成的縱向序列，長(zhǎng)度為l。M（i，j）是當(dāng)qr[i]和db[j]對(duì)齊時(shí)單元（i，j）的打分；Iqr（i，j）是qr[i]和一個(gè)空位對(duì)齊時(shí)單元（i，j）的打分；Idb（i，j）是db[j]和一個(gè)空位對(duì)齊時(shí)單元（i，j）的打分。gapopen是開始一段新空位的罰分；gapext是擴(kuò)展一個(gè)空位的罰分。w（a，b）是一個(gè)打分函數(shù)，當(dāng)a和b相同時(shí)打正分，反之打負(fù)分。jumpqr是matrix2中額外計(jì)算的罰分，是從matrix2中單元向matrix1中單元進(jìn)行跳躍的罰分。jmax是matrix2中單元跳躍目標(biāo)單元的橫坐標(biāo)，對(duì)于matrix2中的單元（i，j）來(lái)說(shuō)，其跳躍的目標(biāo)單元坐標(biāo)為（i-1，jmax）。

變種算法與原算法的最大區(qū)別在于，序列比對(duì)的打分矩陣被劃分為了兩個(gè)部分，分別對(duì)應(yīng)著deletion兩個(gè)邊界附近所選擇出的參考序列（圖4下半部分中Part 1與Part 2）。在第一部分中，全部的比對(duì)分?jǐn)?shù)計(jì)算與原算法相同，在第二部分中，為每個(gè)單元計(jì)算分值時(shí)會(huì)多考慮一項(xiàng)，即來(lái)源于第一部分矩陣上一行中具有最高分值的單元（圖4下半部分中NW-MAX單元）的打分。這個(gè)分值的計(jì)算相當(dāng)于將第一部分矩陣中的部分序列比對(duì)結(jié)果與第二部分矩陣中的部分序列比對(duì)結(jié)果相連接，相連接的兩個(gè)單元所在的位置就是這個(gè)映射所對(duì)應(yīng)的一段連續(xù)空位的邊界點(diǎn)。變種算法對(duì)于這種連接給出一個(gè)固定的罰分，這個(gè)罰分與兩個(gè)單元的橫向距離無(wú)關(guān)。在原算法中，這樣的單元之間的“跳躍”是不允許的，相同的映射在原算法中需要依靠相鄰單元的連續(xù)計(jì)算來(lái)完成（圖4下半部分中虛線箭頭所示），由于原算法中引入空位 需要罰分，因此split read的映射結(jié)果的最終分值將會(huì)受到引入的空位數(shù)量的影響，引入的空位越多，分值越低。這可能導(dǎo)致split read的映射結(jié)果由于引入的空位過(guò)多而導(dǎo)致分值過(guò)低，最終被舍棄。

4實(shí)驗(yàn)結(jié)果與分析

本文將所提出的算法進(jìn)行程序?qū)崿F(xiàn)，稱為PRISM。通過(guò)將人類基因組中deletion注釋加入到參考基因組1號(hào)染色體序列中的方式構(gòu)造了一條模擬基因組序列，并使用模擬測(cè)序軟件[7]對(duì)該模擬基因組序列進(jìn)行模擬測(cè)序生成一套模擬數(shù)據(jù)集。在該模擬數(shù)據(jù)集上，本文將所提出的split read映射方法與一種已有的方法Pindel進(jìn)行了比較。首先是運(yùn)行速度上的比較，結(jié)果如表1所示。由于在取得候選split read時(shí)的標(biāo)準(zhǔn)不同，兩種方法作為輸入的read數(shù)量不同，但是從結(jié)果上可以看出，PRISM的輸入規(guī)模略高于Pindel，而運(yùn)行時(shí)間卻遠(yuǎn)遠(yuǎn)短于Pindel，這證實(shí)了PRISM利用read pair分析結(jié)果來(lái)指導(dǎo)split read映射的方法可以大幅地提高split read映射的效率。第二項(xiàng)比較是split read映射效果的比較，具體結(jié)果如圖5所示，可以看出PRISM在正確映射split read的能力上也要優(yōu)于Pindel。

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

本文提出了一種新的split read映射方法，這種方法利用split read附近的read pair映射結(jié)果分析來(lái)指導(dǎo)split read的映射，以達(dá)到縮小映射過(guò)程中搜索空間，提高映射效率與準(zhǔn)確性的目的。在模擬數(shù)據(jù)實(shí)驗(yàn)中，通過(guò)與已有的方法進(jìn)行對(duì)比，證實(shí)了本文所提出的方法在運(yùn)行效率、與split read映射結(jié)果上都具有優(yōu)勢(shì)。

參考文獻(xiàn)：

[1]LI H， DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform [J]. Bioinformatics， 2009， 25（14）： 1754-1760.

[2]LANGMEAD B， SALZBERG S L. Fast gapped-read alignment with Bowtie 2 [J]. Nature methods， 2012， 9（4）： 357-359.

[3]LANGMEAD B， TRAPNELL C， POP M， et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome [J]. Genome biology， 2009， 10（3）： R25.

[4]YE K， SCHULZ M H， LONG Q， et al. Pindel： a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads [J]. Bioinformatics， 2009， 25（21）： 2865-2871.

[5]DU Z H， LIN F. Improvement of the needleman-wunsch algorithm [J]. Lect Notes Artif Int， 2004， 3066：792-797.

[6]NEEDLEMAN S B， WUNSCH C D. A general method applicable to the search for similarities in the amino acid sequence of two proteins [J]. Journal of molecular biology， 1970， 48（3）： 443-453.

[7]HUANG W， LI L， MYERS J R， et al. ART： a next-generation sequencing read simulator [J]. Bioinformatics， 2012， 28（4）： 593-594.