高 靜 焦 雅 張文廣

摘要：高通量測(cè)序技術(shù)的飛速發(fā)展，給生物信息學(xué)帶來了新的機(jī)遇和挑戰(zhàn)，第二代測(cè)序序列數(shù)量多、長度短使得原來的序列分析手段不再適用。近幾年來，針對(duì)高通量測(cè)序的序列分析算法和軟件日益增多，目前已有上百種，導(dǎo)致選擇合適的軟件成為一個(gè)難題。對(duì)第二代測(cè)序的測(cè)序類型、序列類型以及分析算法進(jìn)行了總結(jié)和歸納，對(duì)現(xiàn)今常用的分析軟件的序列的類型、長度以及軟件應(yīng)用算法、輸入/輸出格式、特點(diǎn)和功能等方面做了詳細(xì)分析和比較并給出建議。分析了現(xiàn)今測(cè)序技術(shù)和序列分析存在的問題，預(yù)測(cè)了今后的發(fā)展方向。

關(guān)鍵詞：高通量測(cè)序：序列比對(duì)；序列作圖；序列比對(duì)工具

中圖分類號(hào)：Q-31

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1007-7847（2014）05-0458-07

以Roche/454焦磷酸測(cè)序（2005年）、lllumina/Solexa聚合酶合成測(cè)序（2006年）和ABI/SOLiD連接酶測(cè)序（2007年）技術(shù)為代表的第二代測(cè)序技術(shù)與Sanger測(cè)序相比，共有的突出特征是單次運(yùn)行產(chǎn)出序列數(shù)據(jù)量大，故而又被通稱為高通量測(cè)序技術(shù)。高通量測(cè)序技術(shù)大大降低了測(cè)序的時(shí)問和成本，因而得到了廣泛應(yīng)用。但隨之而來的短序列（short read）為基因數(shù)據(jù)分析帶來了新挑戰(zhàn)。目前對(duì)短序列一個(gè)常用的分析方法則是將已有基因組序列作為參考基因序列（reference），將短序列與參考基因序列進(jìn)行序列比對(duì)，并在參考基因序列上進(jìn)行定位，這個(gè)過程稱為mapping。序列比對(duì)是基岡數(shù)據(jù)分析最基本的手段，對(duì)其進(jìn)行研究具有重要意義。

1 高通量測(cè)序概述

1.1 測(cè)序分類

大規(guī)模、低成本、快速的高通量測(cè)序技術(shù)被廣泛應(yīng)用于生物研究的各個(gè)方面，當(dāng)前主流的測(cè)序分為：1）基因組學(xué)的全基因組de novo測(cè)序、全基因組重測(cè)序、外顯子&目標(biāo)區(qū)域測(cè)序、簡化基因組測(cè)序；2）轉(zhuǎn)錄組學(xué)的轉(zhuǎn)錄組測(cè)序、數(shù)字基因表達(dá)譜、小RNAs測(cè)序、降解組測(cè)序和長鏈非編碼BNA洲序；3）表觀基囚組學(xué)的全基因組Bisulfile甲基化測(cè)序、RRBS、MeDIP測(cè)序等。

1.2 測(cè)序數(shù)據(jù)類型

現(xiàn)在比較常見的測(cè)序方式有：single -read、paired-end、Mate-pair、color-space。1）Single-read即為單末端測(cè)序，在測(cè)序前先將DNA樣本進(jìn)行片段化處理形成200～500 hp的片段，且將引物序列連接到DNA片段的一端，然后末端加上接頭。這種方法較簡單，但是它只有一端有拼接信息，不利于拼裝；2）Paired-end稱為雙末端測(cè)序或配對(duì)末端測(cè)序，它是指在構(gòu)建待測(cè)DNA文庫時(shí)，在基因片段的兩端都加上測(cè)序引物結(jié)合位點(diǎn)再進(jìn)行測(cè)序。雙末端測(cè)序是在片段兩端都加上接頭，在序列拼裝的時(shí)候更容易定位片段；3）Male-pair也屬于雙末端測(cè)序，與paired-end不同的是它將基因隨機(jī)打成特定為2～10 kb的片段，然后經(jīng)末端修復(fù)，生物素標(biāo)記和環(huán)化等實(shí)驗(yàn)步驟后.再把環(huán)化后的片段打斷成400～600 bp的片段并標(biāo)記測(cè)序。因?yàn)榻?jīng)過了環(huán)化步驟，Mate-pair比paired-end方式測(cè)得的序列片段更長，從而可以將基因序列中大量的repeat包含在內(nèi)，減少拼接難度；4）Color-space read是ABI公司的SOLiD測(cè)序儀檢測(cè)出的read，SOLiD可以同時(shí)檢測(cè)兩個(gè)相鄰的堿基，并利用顏色空間的4種不同的顏色對(duì)兩個(gè)堿基編碼。對(duì)于這種編碼方式，只要知道顏色編碼對(duì)應(yīng)的序列上任何一個(gè)位置的堿基類型，就可以將顏色編碼解碼為原來的堿基序列，圖l為顏色編碼過程。

2 序列比對(duì)方法

對(duì)于百萬條甚至上億條短序列在reference上的定位，傳統(tǒng)的動(dòng)態(tài)規(guī)劃算法不能滿足我們的要求。為了加快序列比對(duì)的速度，應(yīng)用啟發(fā)式算法是必然趨勢(shì)。同前，絕大多數(shù)的啟發(fā)式算法都是通過建立索引加快比對(duì)速度。建立索引就是用一個(gè)輔助的數(shù)據(jù)結(jié)構(gòu)存儲(chǔ)待比對(duì)序列，這種數(shù)據(jù)結(jié)構(gòu)能夠?qū)⑿蛄兄蟹弦欢ㄒ?guī)則的子序列凸顯出來。常用的數(shù)據(jù)結(jié)構(gòu)有哈希表和后綴樹兩種。

2.1 基于種子的哈希表索引方法

這種方法多用于數(shù)據(jù)庫搜索或read在參考序列上的mapplng，以哈希表形式建立序列的索引，這類算法的代表軟件是SSAHA，基于種子的哈希表方法具體步驟：

第一步，建立哈希表。

DNA序列由A、C、T、G四種脫氧核苷酸組成，長度為k的連續(xù)片段也就是“種子”有4K種可能將4K種子存放在哈希表中，我們用一個(gè)公式將所有種子唯一標(biāo)識(shí)，將4個(gè)堿基的值f（x）轉(zhuǎn)換為二進(jìn)制數(shù)據(jù)表示，只占用2 Bit的內(nèi)存。

第二步，將數(shù)據(jù)庫中的序列與哈希表關(guān)聯(lián)

設(shè)一個(gè)含有n條DNA序列的數(shù)據(jù)庫D={S1，S2，…，Sn}，將數(shù)據(jù)庫中的每條序列分解連續(xù)種子w，其長度k=5，從頭開始每次偏移一位直到序列結(jié)束，那么長度為l的序列有（l-k+1）個(gè)種子。并將每一個(gè)“種子”所在的序列，N（N=l，2，…，n）和在序列中的順序號(hào)L （L=1，2，…，l）在對(duì)應(yīng)的哈希表中的“種子”以二維數(shù)組方式記錄下來（N，L）。將S1

第三步，查詢序列Q也分解成長度為k的種子并計(jì)算其標(biāo)識(shí)號(hào)，在哈希表中找到數(shù)據(jù)庫中與之匹配的所有位置。

為了讓算法更高效，算法在各個(gè)方面進(jìn)行了改進(jìn)：1）“種子”長度，允許錯(cuò)配和indel的數(shù)量等因素直接影響比對(duì)結(jié)果的準(zhǔn)確性、敏感度、時(shí)間和空間花費(fèi)等指標(biāo)，所以研究者們一直致力于對(duì)“種子”形式和選擇的改進(jìn)中。Blastn、Blat、SOAP、SeqMap、MAQ、SHRiMP使用不同的種子模型，如空位種子、空位種子組（同時(shí)使用多個(gè)種子模型）等；2）-些軟件如MAQ、RMAP、SeqMap、ZOOM等將reads序列建立索引，而一些軟件將參考序列數(shù)據(jù)庫建立索引；3）設(shè)置閾值F將哈希表中出現(xiàn)頻率低于F的“種子”刪去。

2.2 基于前/后綴樹的索引方法

后綴樹是一種重要的數(shù)據(jù)結(jié)構(gòu)，AVID、MUMmer、MUMmer等序列比對(duì)算法都是基于后綴樹的這個(gè)結(jié)構(gòu)，但與兩序列相似性比對(duì)時(shí)用于尋找最大公共子序列不同，read在參考序列中mapping時(shí)，是將read或其子序列在按序排列的所有后綴中一一比對(duì)找到匹配位置。算法VCAKE、SSAKE、SHARCGS、BWA-SW應(yīng)用前綴樹數(shù)據(jù)結(jié)構(gòu)，它將序列所有的前綴用樹的結(jié)構(gòu)表示，原理與后綴樹相同。如圖3所示為序列S=GATGAC的前綴樹和前綴DAWG （DirectedAcyclic Word Graph，有向無環(huán)詞圖）。

由于后綴樹結(jié)構(gòu)在比對(duì)大型序列時(shí)，空間占用比較大，Abouelhoda等改進(jìn)了Manber和My-ers提出的后綴樹存儲(chǔ)方法，稱為增強(qiáng)后綴數(shù)組來提高空間利用率。Ferragina和Manzini提出的FM （full-text minute-space）-indexc是一種基于BWT （burrows-wheeler transform）的全文本壓縮索引結(jié)構(gòu)，利用BWT矩陣與后綴樹之間的關(guān)系壓縮后綴數(shù)組和索引，減少內(nèi)存占用率。

第一步，BWT矩陣T設(shè)序列S以“$”為結(jié)束符號(hào)，輪換列出序列S所有的后綴記為數(shù)組M，然后將數(shù)組M的每行按照字母順序排序得到BWT矩陣T；

第二步，矩陣T中第i行在矩陣M中的行號(hào)（也是矩陣T第一列字母在序列S中的序號(hào)）為數(shù)組S[i]矩陣T中最后一列的結(jié)果為BWrr數(shù)組B[i]。我們可以根據(jù)這兩個(gè)數(shù)組將序列S還原，這個(gè)算法稱為UNPERMUTE算法；

第三步，創(chuàng)建數(shù)組Oc （A）表示矩陣T中第一列的第一個(gè)A在第幾行，Occ（2，G）表示矩陣T的最后一列中第2行的G是此列的第幾個(gè)G。以這兩個(gè)數(shù)組為索引可以很快找到匹配點(diǎn)；

用這種方法比對(duì)的軟件有：Bowite、BWT-SW、BWA-SW。我們知道人類的基因組為3 Gb，那么數(shù)組Occ和S的每個(gè)值都要用4 Byte表示，3 Gb就要分別用12 GB內(nèi)存存儲(chǔ)。在Bowtie中應(yīng)用了FM-index，每隔數(shù)行建立一個(gè)索引，使得人類基因組的內(nèi)存占用約為1.3 GB，完全可以在普通機(jī)器上運(yùn)行。

基于這兩種數(shù)據(jù)結(jié)構(gòu)的算法都有利于提高比對(duì)效率，兩者的區(qū)別在于后綴樹可以有效減少不精確匹配，并可避免比對(duì)過程中做的無用功，這個(gè)特點(diǎn)適用于相同物種之間相似性高的序列比對(duì)和尋找保守區(qū)。而應(yīng)用哈希表數(shù)據(jù)結(jié)構(gòu)的算法具有較高的敏感性，有利于發(fā)現(xiàn)SNPs和突變?？捎糜诰植科ヅ浠驈拇罅繑?shù)據(jù)巾搜索匹配點(diǎn)以及跨物種序列間的比對(duì)。

3 序列分析軟件比較

表1中列出民目前常用的第二代測(cè)序序列的分析軟件，對(duì)軟件分析的序列類型、序列長度、軟件特點(diǎn)、功能以及可輸入、輸出格式等方面進(jìn)行了詳細(xì)析和比較。

選擇合適的軟件要根據(jù)軟件適用的序列類型、長度以及運(yùn)用的算法、性能特點(diǎn)、輸入/輸出格式等項(xiàng)全而考慮，如表1中列出的常用軟件的各項(xiàng)屬性進(jìn)行篩選：1）根據(jù)第2列中序列類型，如果是Spliced reads可以選擇QPALMA和TopHat軟件，（color-spaced reads可以選擇BWT -SW、BFAST；2）根據(jù)第4列中read的長度，SSAHA、BLAT用于長序列的比對(duì)；SOAP.ZOOM、SeqMap適用于short-reads的mapping；3）按照表中5和6列軟件的特點(diǎn)和功能，如Bowtie、BWT -SW、RMAP、MAQ、QPALMA、SHiRMP等在分析序列對(duì)，應(yīng)用了質(zhì)量分?jǐn)?shù)，這在很大程度上提高了比對(duì)的準(zhǔn)確性。用Bowlie和BWA在分析速度內(nèi)存占用方面略優(yōu)于其他軟件：4）第8列為軟件在比對(duì)時(shí)是否允許錯(cuò)配，gaps以及indels，這也影響比對(duì)結(jié)果。例如，MAQ在比對(duì)Illumina paired-end reads時(shí)允許gaps這樣可以加大paired-end reads拼接的成功率。BWA、BFAST、SHiRMP、PASS等在進(jìn)行比對(duì)時(shí)允許一定或不定數(shù)量的gaps可以提高比對(duì)的敏感度，降低假SNP的出現(xiàn)率，有利于發(fā)現(xiàn)真正的突變位點(diǎn)；5）最后，根據(jù)軟件的輸入格式是否符合待比對(duì)序列的格式、輸出格式是否又符合下一步分析要求的格式進(jìn)行軟件選擇。

此外，在分析數(shù)據(jù)的時(shí)候?qū)蓚€(gè)或多個(gè)軟件結(jié)合使用，例如在處理spliced reads時(shí)，可以使用Bowtie或BWT建立索引文件，再將索引文件使用TopHat、SoapSplice等軟件進(jìn)行splice junction檢測(cè)。

4 結(jié)論

本文對(duì)國內(nèi)外第二代測(cè)序技術(shù)以及數(shù)據(jù)的分析算法進(jìn)行了總結(jié)和分析，對(duì)目前流行的序列分析軟件進(jìn)行詳細(xì)歸納和比較并給出了參考建議，有利于從整體上把握測(cè)序技術(shù)和序列分析的發(fā)展?fàn)顩r.為今后序列分析和研究提供參考。

通過上述研究發(fā)現(xiàn)雖然基因分析軟件已有很多且解決了一定的問題，但仍然存在困難：1）由于序列中大量repeat（重復(fù)片段）的存在，很多短reads作為repeat的一部分，不能對(duì)其準(zhǔn)確定位，增大了拼裝難度； 2）當(dāng)前大多硬件環(huán)境較低，計(jì)算時(shí)間以小時(shí)或天為單位，嚴(yán)重限制了研究進(jìn)展；3）比對(duì)結(jié)果的準(zhǔn)確率和下一步分析的可用性方面有待提高。

針對(duì)上述問題，第三代單分子納米孔測(cè)序技術(shù)目的是在第二代測(cè)序技術(shù)的基礎(chǔ)上提高測(cè)序的準(zhǔn)確性、增長read的讀長。雖然其技術(shù)還不成熟，但必將成為未來幾年的主流技術(shù)。序列比對(duì)方面，現(xiàn)有越來越多的軟件提供并行處理方式，所以算法并行化以及云計(jì)算技術(shù)的應(yīng)用將會(huì)成為序列分析的重要發(fā)展方向。