姜玥 王亞東

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

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

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

0引言

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

1基本概念

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

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

1.2Read與split read

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

1.3雙末端測序

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

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

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

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

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

其中：

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

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

4實驗結(jié)果與分析

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

5結(jié)束語

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

參考文獻(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.