DNA測序技術及其拼接算法綜述

李艷慧，張少強

（天津師范大學計算機與信息工程學院，天津300387）

DNA測序技術用以分析特定DNA片段的堿基序列（腺嘌呤（A）、胸腺嘧啶（T）、胞嘧啶（C）與鳥嘌呤（G））的排列方式.DNA序列的測定是分子生物學研究的一項重要且關鍵的內容，如基因的定位、基因結構和功能的研究、基因表達與調控、基因與疾病的關系等，都需要對DNA結構有詳細的了解.因此，DNA測序技術極大推動了生物學和醫學的研究和發展.1990年代啟動的人類基因組計劃，其主要目的就是測定人類DNA的30億個堿基對序列及識別人類DNA中所有的基因，這直接推動了生命科學和現代醫學的發展，也加快了DNA測序技術的日趨成熟.

自DNA測序技術發明以來，人類對各物種基因和基因組結構的研究和探索從未中斷.從第一代測序技術到第四代測序技術，測序的讀序（reads，即測序儀器產生DNA片段的序列）長度由長到短，然后再由短變長，測序技術取得了長足的發展.目前，第一代測序技術產生的基于熒光標記和Sanger的雙脫氧鏈終止法原理的測序儀慢慢退出了歷史舞臺，第二代測序技術在市場上占有絕對的優勢，但第三代和第四代測序技術也在迅猛發展.第二代測序技術目前是測序的主流，其顯著特征是高通量，一次能對幾十萬條DNA片段進行序列測序.第三代是以單分子測序為特點的測序技術，這種技術無需PCR擴增，更能反應樣本的真實情況，測序通量也更高，操作過程更簡單，成本更低.第四代測序技術又稱納米孔測序技術，它具有高通量、低成本，并能快速進行基因檢測的特點.測序技術的發展歷程如圖1所示.

測序完成后的第一步也是最重要的一步就是根據讀序拼接回貼成完整的序列，其測序與拼接過程如圖2所示.拼接算法用于將測出的讀序拼接成完整的染色體序列（chromosomesequence）、轉錄本序列（transcript sequence）或因為測序不完整而形成的支架序列（scaffolds），最終形成基因組、轉錄組序列或其他功能序列.所有的拼接算法基本上是先將讀序根據重疊關系連接成“疊連群”（contigs），然后再將疊連群根據配對（pair-end）讀序和其他信息形成有序集合來構建更長的支架序列.

圖1 測序技術的發展歷程Fig.1 Development of sequencing technology

圖2 DNA測序及拼接過程示意圖Fig.2 Diagram of DNA sequencing and assembly

1 第一代測序

1.1 第一代測序技術

第一代測序技術起源于1970年代，測序采用均一的單鏈DNA分子.第一代測序技術的主要方法有化學降解法、雙脫氧鏈終止法、熒光自動測序技術和雜交測序技術.

化學降解法的優點是準確性較好，易掌握，缺點是過程較繁瑣.雙脫氧鏈終止法又稱Sanger法，該方法簡便、快捷，經過后續的不斷改良，成為了第一代DNA測序的主流.Sanger法測序的原理是：每個反應含有所有4種脫氧核苷酸三磷酸（dNYTP），使之擴增，并混入限量的一種不同的雙脫氧核苷酸三磷酸（ddNTP）使之終止.由于ddNTP缺乏延伸所需要的3′-OH基團，使延長的寡聚核苷酸選擇性地在G、A、T或C處終止，終止點由反應中相應的雙脫氧核苷酸而定.熒光自動測序技術基于Sanger原理，使用熒光標記代替同位素標記，利用成像系統自動監測，進而提高了測序的準確性和速度.雜交測序技術利用DNA雜交原理，將一系列已知序列的單鏈寡核苷酸片段固定在基片上，將待測DNA樣本片段變性后與其雜交，根據雜交情況排列出樣品的序列信息.

1.2 第一代測序算法

第一代測序技術最早應用的拼接算法是OLC（overlap-layout-consensus）算法，最初由 Staden（1980年）[1]開發，隨后許多科學家對其進行了擴展和闡述.OLC算法是一個直觀的拼接算法，它適用于讀序長度較長的測序數據.該算法一般分為3個步驟：（1）找到所有讀序之間的重疊（O）；（2）對所有的讀序根據信息在圖上進行布局（L）；（3）推斷出共有序列（C）.隨著Sanger測序技術的廣泛應用，基于OLC的拼接算法有Arachne[2]、Celera Assembler[3]、CAP3[4]、PCAP[5]、Phrap[6]、Phusion[7].Arachne拼接算法的關鍵特征包括：尋找重疊讀序的高效和靈敏、糾正拼接前的錯誤實現高度準確性的重疊評分步驟、基于正向-反向拼接的讀序合并、通過正向-反向拼接不一致檢測重復疊連群.CAP3是第三代CAP序列拼接算法，它具有修剪5’端和3’端低質量區域讀序的能力.CAP3和Phrap比較而言，Phrap通常比CAP3產生更長的疊連群，而通常CAP3在共有序列中產生的誤差比Phrap少，并且使用正向-反向約束的低通量數據，CAP3比Phrap更容易構建支架序列.Phusion使用讀序數據和配對讀序數據提高了大型基因組數據拼接結果的高度精確性.

2 第二代測序

2.1 第二代測序技術

近幾年，第二代測序技術在基因組學的研究得到了廣泛的應用.第二代測序技術也被稱為新一代測序技術（next generation sequencing，NGS）.第二代測序技術相比于第一代測序技術，成本大幅度下降.第二代測序技術主要包括羅氏454公司的454級焦磷酸測序技術、Illumina公司的Solexa測序技術和ABI公司的SOLiD測序技術.目前全球使用量最大的第二代測序技術是Illumina公司的Solexa和Hiseq測序技術，這2種測序技術的核心原理是相同的，都采用邊合成邊測序的方法，它的測序過程主要分為4個步驟：（1）DNA待測文庫構建；（2）Flowcell（即將基因組DNA的隨機片段附著到光學透明的玻璃表面）；（3）橋式PCR擴增與變性；（4）測序.

第二代測序技術引入表觀遺傳學后，形成了各種表觀遺傳學測序的方法.其中，第二代測序技術與染色質免疫共沉淀技術相結合的ChIP-seq技術，已成為功能基因組學，尤其是基因表達調控領域研究的關鍵技術.當前，ChIP-seq的應用主要包括：（1）DNA序列上轉錄因子結合位點的識別；（2）表觀遺傳學領域中基因組DNA甲基化、組蛋白修飾、核小體定位等.目前，第二代測序平臺中的Illumina/Solexa測序平臺上的RNA-seq技術應用最廣.該技術利用高通量測序對組織或細胞中所有RNA反轉錄而成的cDNA文庫進行測序，統計相關讀序的數目并計算不同RNA的表達量，進而發現新的轉錄本，而不再需要構建克隆文庫.以亞硫酸氫鹽轉換法與第二代測序技術相結合的全基因組亞硫酸氫鹽測序法是高通量檢測基因組中DNA甲基化的標準方法.DNase-seq和FAIRE-seq是2種在全基因組范圍內鑒定染色質可接近性、檢測與調節活性相關的DNA序列的方法.FAIRE-seq是在DNase-seq的基礎上產生的.將第二代測序技術與染色質構象捕獲技術（chromosome conformation capture，即3C技術）相結合，用高通量的方式檢測染色體的高級結構及基因組位點間的相互作用，產生了如4C-seq、5C-seq、Hi-C 和 ChiA-PET 等技術.

2.2 第二代測序算法

2.2.1 RNA-seq數據轉錄組拼接算法

第二代測序RNA-seq數據轉錄組的拼接算法有2種形式（如圖3所示）：基于參考基因組的拼接和從頭（de novo）拼接（不基于參考基因組的拼接）.基于參考基因組的拼接首先將大量的RNA-seq測序短讀序準確地回貼到基因組的相應位置再進行拼接，而從頭拼接是根據RNA-seq測序片段之間的重疊關系進行拼接.

圖3 RNA-seq數據拼接的2種形式Fig.3 Two categories of RNA-seq data assembly

第二代測序RNA-seq轉錄組拼接算法有3種類型：剪接圖（splicing graph）、重疊圖（overlap graph）和De Bruijn圖.基于參考基因組拼接大致包括回貼、構圖和根據圖結構設計算法重構表達轉錄本.其中，在回貼部分常使用的聯配工具有Hisat[8]、Tophat[9]、Tophat2[10]、SpliceMap[11]、MapSplice[12]等，這些工具可以處理可變剪切.回貼結束后，所有測序片段會由于回貼到基因組相應位置而聚類成簇（代表特定基因相應的表達情況），進而形成有向無圈圖結構（剪接圖或重疊圖）.

基于參考基因組的拼接算法包栝StringTie[13]、Bayesember[14]、Cufflinks[15]、 Scripture[16]、IsoLasso[17]、iReckon[18]、CER[19]、Traph[20]、CIDANE[21]和 Scallop[22]等 .比較典型的算法是Cufflinks和StringTie.Cufflinks的主要思想是根據所有回貼后兼容的測序片段的兼容關系建立重疊圖（有向無圈），圖中的點代表雙端測序片段，點之間的有向邊代表點之間的測序片段是兼容的.通過最小路覆蓋模型[23]，從圖中找到數目最少并且能覆蓋所有點的路，將找到的路作為最終預測的表達轉錄本的集合.StringTie的主要思想是根據所有回貼后兼容的測序片段的兼容關系建立剪接圖，圖中的點代表外顯子，點與點之間的有向邊表示2個外顯子之間發生剪接事件.與Cufflinks的不同之處在于，StringTie設計了一個貪婪算法逐條預測表達的轉錄本，然后采用最大流模型逐條估計轉錄本表達量，最后對剪接圖進行更新，迭代進行.Scallop算法能精確拼接識別多外顯子和低表達水平的轉錄本，與StringTie相比可以得到更多正確的多外顯子轉錄本.

轉錄組的從頭拼接算法包括ABySS[24]、SOAP de novo-Trans[25]、Oases[26]、IDBA-Tran[27]、BinPacker[28]、Bridger[29]、Trinity[30]等.其中 ABySS、SOAP de novo-Trans、Oases和IDBA-Tran本質上是基因組拼接的方法.Trinity是應用最廣、公認度最高的從頭轉錄組拼接軟件，也是第一個專門針對轉錄組拼接開發的軟件.Trinity的主要思想是將測序短序打碎成k-mer（連續k長度的片段），并用哈希表存儲每個k-mer及其出現的頻率，然后選擇頻率最高的k-mer構建De Bruijn圖，最后通過動態規劃算法遍歷De Bruijn圖，根據打分機制篩選出符合條件的路.Bridger將Cufflinks和Trinity兩者優勢相結合，它的主要思想是構建k-mer哈希表后，選取頻率最高的k-mer作為種子直接延伸出整個連通分支，作為一個剪接圖，之后在重復選取k-mer延伸連通分支時，被使用過的k-mer根據一定規則被繼續使用，最后對每個剪接圖應用最小路覆蓋算法找出最小路覆蓋作為輸出.

2.2.2 ChIP-seq數據相關算法

ChIP-seq數據分析流程見圖4，其主要階段包括讀序回貼（read mapping）[31]，峰值檢測（peak calling）[32-35]和模體識別[36].現在比較流行的回貼軟件包括Bowtie[37]、BWA[38]、GSNAP[39]、MAQ[40]等.其中，Bowtie 和 BWA 是基于Burrows-Wheeler的轉換法.峰值檢測包括3個步驟：尋找sharp peak和broad peak；生成富集倍數（fold enrichment）軌道；生成logLR軌道.峰值檢測軟件主要有 MACS[41]、CisGenome[42]、SICER[43]、CCAT[44]、ZINBA[45].其中MACS和CisGenome適用于分析尖銳信號，SICER和CCAT適合分析連綿信號，ZINBA和MACS適合分析混合型峰.模體識別是指在給定的一組峰值區域中識別模體并找出轉錄因子結合位點.其相關算法有DREME[46]、MEME-ChIP[47]、WEEDER[48]等.DREME 專門設計用于發現短核DNA結合模體的真核生物的轉錄因子，它可以在幾分鐘內優化分析ChIP-seq數據集.DREME可作為基于模體（motif）的序列分析工具的MEME[49]組件的一部分.MEME-ChIP用于DNA大數據集模體分析，其網絡服務設計用于分析ChIP-seq數據峰值區域.WeederWeb[48]是一個Web接口，用于在一組相關的調控DNA序列中自動發現保守模體.

圖4 ChIP-seq數據分析流程Fig.4 Analysis process of ChIP-seq data

2.2.3 第二代測序基因組拼接算法

目前基于第二代測序技術產生的短讀序序列進行基因拼接的算法，大致分為OLC圖、Greedy圖和De Bruijin圖3種類型.因為第二代測序技術產生的基因片段較短且高通量測序速度較快，使得De Bruijn圖算法在第二代測序技術中占據主要地位.

OLC圖算法最早提出于第一代測序，之后科學家不斷對其進行闡述和擴展.該算法遵循overlap/layout/consensus框架，其原理是利用重疊圖，用短讀序表示重疊圖的結點，讀序之間的重疊關系表示重疊圖的邊，通過計算短讀序重疊關系進行短讀序基因拼接.其優點是將序列拼接問題轉換為圖論問題，適用于長讀序基因拼接，缺點是算法無法解決重復序列問題，而且對于第二代基因測序數據的特點，算法運算量較大.典型 算 法有 Arachne[2]、CAP3[4]、EDENA[50]、Newbler[51].其中CAP3和Arachne是在第一代測序時提出的，CAP3是經典的OLC算法，該算法能快速找出讀序之間的重疊關系，能有效計算和評價重疊關系.EDENA也是經典的OLC算法，其優點是拼接正確的contig概率更大.

Greedy圖算法是最早提出的面向第二代測序技術的拼接算法之一，它遵循貪婪算法的原理，給定一個讀序或seed，只要它滿足條件就對這個seed不斷地延長直到所有的數據用完或超過設定閾值.其優點是原理簡單，操作方便，缺點是算法容易得到局部最優解，難以得到全局最優解.基于Greedy圖的算法有SSAKE[52]、SHARCGS[53]、VCAKE[54]、QSRA[55].SSAKE 和VCAKE是典型的Greedy算法，選擇最大重疊進行組裝.SHARCGS和QSRA的相似之處在于算法主要有濾除、組裝和融合3個步驟.

De Bruijin圖算法遵循Euler超路徑算法，其原理是將讀序打斷成長度為k的核苷酸片段（k-mer），利用核苷酸片段間的overlap關系構建DeBruijn圖，然后得到疊連群，通過遍歷De Bruijn圖得到支架序列，進一步操作后得到基因組序列.其優點是能有效解決重復序列問題，使用哈希表存儲方便快速查找，缺點是構造De Bruijin圖對內存要求較高，消耗計算資源較大.基于 De Bruijn 圖的算法有 ABySS[24]、SOAP de novo[25]、EULER-SR[56]、ALLPATHS[57]、Velvet[58].Velvet和 ABySS是比較經典的基于De Bruijn圖的拼接算法，其中ABySS最大的特點是采用De Bruijn圖的分布式表示，在多臺計算機網絡中并行計算進行拼接.SOAP de novo是一個短讀序的基因拼接算法，該算法能精確拼接大型基因組數據.

3 第三代測序

3.1 第三代測序技術

第三代測序技術是在第二代測序技術基礎上增加讀序長度，加快測序運行速度，其特點是單分子測序（不經PCR直接進行邊合成邊測序）.第三代測序技術包括Heliscope測序技術、SMRT（single molecule real time，單分子實時測序）、離子半導體測序技術（ion torrent）等.其中SMRT測序技術是較為成熟的.SMRT技術以SMRT芯片（帶有很多ZMW（zero-modewaveguides）孔的厚度為100 nm的金屬片）為測序載體進行測序反應，并采用邊合成邊測序的思想.第三代測序技術原理是：首先將長讀序互相比對并識別長讀序之間的重疊區，如一個較短的讀序完全包含在一個較長的讀序內或者2個長讀序之間有相同的序列區域；然后通過不同的擴展策略對長讀序進行拼接.其優點是更容易匹配讀序序列和降低測序成本，缺點是測序方法不夠健壯，錯誤率高，比對時間較長，不能被廣泛使用.

3.2 第三代測序算法

3.2.1 第三代測序拼接算法

第三代測序拼接算法利用長讀序跨過基因組中大部分重復區，該類方法一般首先將長讀序相互比對并識別長讀序之間的重疊區，然后通過不同的擴展策略對長讀序進行拼接.

目前已有的長讀序拼接工具主要包括HGAP[59]、Dazzler[60]、Sprail[61]、MHAP[62]、Sparc[63]、PBcR[64]等 .HGAP首先建立一個加權的有向無環圖來表示比對在一起的多個長讀序的關聯關系，然后通過找到圖中的最大權重來確定一致序列，最后利用Celera[65]拼接糾錯后的長讀序.自拼接的PBcR方法通過BLASR[66]軟件識別長讀序之間帶有噪聲的重疊區，利用重疊區對長讀序進行糾錯，最后同樣利用Celera對糾錯后的長讀序進行拼接.Dazzler方法基于排序和合并策略，設計了一個高敏感度、高效率的過濾器，用來處理長讀序之間局部比對所包含的種子點，并提出了一個線性時間的啟發式算法用來找到經過種子點的長讀序之間的局部比對.Sprail方法首先選擇一定比例的任意長讀序作為種子讀序，利用比對程序將其他長讀序比對到種子讀序，從而對種子讀序進行糾錯，最后使用Celera拼接糾錯后的長讀序.MHAP方法提出了一個概率算法，可以有效地識別帶有噪聲的長讀序之間的重疊區，具體而言，將長讀序重疊區識別簡化為整數序列比較問題，從而大幅減少長讀序重疊識別時間.

3.2.2 第三代測序糾錯算法

第三代測序技術解決高錯誤率策略主要有2個：（1）利用來自另一個測序平臺的低錯誤率的短讀序，在保證足夠的覆蓋倍數下，糾正長讀序中的錯誤，該類方法稱為混合方法；（2）利用長讀序的錯誤分布比較均勻的事實，在保證長讀序足夠的覆蓋率下，利用長讀序自身的信息進行糾錯，該類方法稱為長讀序自糾錯方法.

目前，針對PacBio長讀序的糾錯算法有：（1）自糾錯方法：DAGCon[67]、PBcR[64]、LoRMA[68]；（2）基于短讀序的糾錯方法：PBcR[64]、LSC[69]、PacBioToCA[70]、proovread[71]、ECTools[72]、Cerulean[73]、LoRDEC[74]、Jabba[75].DAGCon 的主要思想是通過長讀序間的多對多比對進行自身糾錯.PBcR的主要思想是通過長讀序間的局部比對構建多重序列比對，然后生成共有序列來進行糾錯.PBcR的最新版本既支持長讀序的自糾錯，也支持基于短讀序的糾錯.LSC的主要思想是先將長讀序中相同的堿基壓縮，然后再將短讀序比對到長讀序，利用比對結果對長讀序進行糾錯.PacBioToCA直接將短讀序比對到長讀序上，利用比對結果對長讀序進行糾錯.

4 第四代測序技術

第四代測序技術又稱為納米孔測序技術，它不需要對DNA進行生物或化學處理，而是基于物理的原理進行測序.其簡要原理是：將經過提取的DNA的一端與一種特殊的蛋白連接，該蛋白主要起到與測序納米孔的連接及控制序列進入速度的作用，同時將DNA的另一端連接，這樣DNA的正義鏈和反義鏈可以依次通過納米孔.通過納米孔的DNA單鏈可以產生微弱的電流變化，感應器感受電信號并導出，通過對隨時間變化的電信號進行解析，可以識別基因中堿基對的排列順序.其優點是超長讀序、高通量、低成本、更少的測序時間和更為簡單的數據分析，缺點是測序準確率較低，納米技術在制造、測序和集成等方面也存在著諸多挑戰.

目前，針對Nanopore長讀序的糾錯方法包括Nanocorr[76]、NaS[77]、MiRCA[78]混合糾錯算法.Nanocorr的主要思想是將短讀序比對到長讀序上，通過最大遞增子序列篩選出最優的短讀序比對片段，過濾掉噪音信息.NaS的主要思想是將短讀序比對到長讀序上，將完全匹配的子序列作為種子，通過種子招募相似的短讀序，最后對招募到的短讀序進行拼接.MiRCA的主要思想是利用拼接工具將短讀序拼接成長的疊連群，然后將疊連群比對到長讀序上，利用疊連群和長讀序的比對結果，對長讀序進行糾錯.

5 總結與討論

隨著基因測序技術的不斷發展，測序拼接的正確率和速度也在不斷提高.第一代測序技術的測序讀長可達1 000 bp～2 000 bp，準確率高達99.999%，但存在測序成本高，通量低等方面的缺點.隨著基因組研究的不斷深入，第一代測序技術已經無法滿足人們大規模測序工作的需求，這直接推動了第二代測序技術的產生.由于第二代測序技術產生的讀序較短，而且基因組中存在大量復雜的重復區，這些重復區使基于第二代測序技術的短讀序序列組裝結果比較零散，無法獲得完整的基因組序列.第三代測序技術所產生的長讀序長度可以達到數萬堿基，這些長讀序可以跨過基因組中大部分的重復區，進而幫助研究人員獲得完整的基因組序列.但是第三代測序技術的測序錯誤率高達15%左右，這對長讀序組裝結果的質量和效率造成了較大的影響，特別是錯誤富集的區域.第四代測序的納米孔技術，其顯著優點為無標簽、超長讀序（104～106個堿基）、高通量和低材料需求.這些優點使實驗過程大大簡化，可以很容易地用于DNA測序應用.四代測序技術的比較見表1.

目前第四代技術尚未大面積推廣，最新的拼接算法研究思路是將第二三代測序的短讀序和長讀序相結合，這類算法主要有LINK[79]、SSPACE-LongRead[80]等，另外很多拼接算法存在占用內存資源較多的缺陷.

目前拼接算法面臨的主要挑戰有：如何有效地分析大型數據集、處理序列錯誤和正確捕獲基因組的重復結構.除此之外，還需要繼續開發新的方法來表示和分析拼接圖，從而正確地拼接重復序列并正確區分共組裝單倍體（co-assembled haplotypes）的基因組變異體.同時，為有效地處理越來越大的數據集，需要對基于De Bruijn圖的拼接算法進行并行化.此外，內存效率也是基因組拼接算法的重點研究內容.為減少內存，相關研究在拼接算法中引入了稀疏圖、數據結構壓縮、布隆過濾器和重疊有效計算的FM索引等.其中，基于布隆過濾器的算法利用了布隆過濾器的空間效率很高的隨機數據結構，布隆過濾器利用位數組簡潔地表示集合，并判斷元素是否屬于集合，這使得基于布隆過濾器算法的空間效率和查詢時間都遠遠優于一般的算法，其缺點是有一定的誤識別率和刪除困難.這些新的高效內存拼接算法允許分析更大的數據集，并且不需要高性能的計算設備.總之，目前拼接算法的應用大部分都處于串行的實現階段，拼接速度不夠高，在序列拼接問題中如何實現并行處理，提高拼接準確率以及減少內存資源占用等方面的研究具有重要意義.

現代測序實驗通常會產生配對數據，這些數據可以限制讀序對之間的相對方向和距離的信息.分析這些信息可以解決序列重復問題，并可以將單個的重疊群連接成長度較長的支架，這是提升拼接算法效率的一個重要領域.同時，在低覆蓋率和重復區域進行局部拼接的算法設計是最困難和亟待解決的問題.

目前，在對拼接問題進行廣泛的數學分析的情況下，依然需要依靠啟發式算法和其他技術，這是因為難以提出實際的拼接數學模型.應用于拼接的大多數數學公式表明，雖然找到最佳的解決方案是可能的，但仍需要大量工作.

表1 測序技術的比較Tab.1 Comparison of sequencing techniques