LEDA：一種基于Levenshtein距離的DNA序列拼接算法

崔競松，薛 慧，王蘭蘭，郭 遲

1.空天信息安全與可信計算教育部重點實驗室，武漢大學國家網(wǎng)絡(luò)安全學院，湖北武漢 430072；

2.河海大學 理學院，江蘇 南京 211100；

3.武漢大學衛(wèi)星導航定位研究中心，湖北武漢430079

0 引言

DNA測序技術(shù)是指在給定的基因組中確定堿基序列排列方式，堿基序列包括腺嘌呤（A）、胸腺嘧啶（T）、胞嘧啶（C）與鳥嘌呤（G）。兩條核苷酸序列配對構(gòu)成DNA序列，其中每一對堿基A與T、G與C構(gòu)成堿基對（base pair，bp）［1］。DNA測序技術(shù)是生物學家們了解DNA結(jié)構(gòu)的重要技術(shù)手段。

目前的測序技術(shù)大致可以分為三代。第一代測序技術(shù)——傳統(tǒng)的Sanger測序［2］，早在1977年被提出。該測序技術(shù)能夠產(chǎn)生長度約為1 000的讀段，且測序錯誤率極低，范圍在0.001%～0.01%，但是Sanger測序成本高，速度慢，應(yīng)用范圍受到限制。因此，自2005年以來，大量的第二代測序技術(shù)被提出，第二代測序技術(shù)通常也被稱為下一代測序技術(shù)［3］（next-generation sequencing，NGS）。與第一代Sanger測序相比，下一代測序技術(shù)成本較低，可以在較短的時間內(nèi)對大量的DNA序列片段同時測序，而且其測序錯誤率也比較低，為0.5%～1.0%。但是，下一代測序技術(shù)產(chǎn)生的讀段長度較短，一般僅包含75到300個堿基。因此，為了解決下一代測序技術(shù)存在的讀段長度較短的問題，第三代測序技術(shù)被提出，主要包括Single Mplecule Real Time Sequencing［4］和Oxford Nanopore［5］。第 三代測序技術(shù)能夠產(chǎn)生包含超過10 000個堿基的讀段，但是測序錯誤率過高，在10%～15%，遠遠高于第一代測序技術(shù)和下一代測序技術(shù)。而且，第三代測序技術(shù)的測序成本也明顯高于下一代測序技術(shù)。

由于下一代測序技術(shù)的低測序成本、低錯誤率和高通量等特性，目前該測序技術(shù)為市面上應(yīng)用最為廣泛的測序技術(shù)。由于受技術(shù)限制，下一代測序技術(shù)存在的問題是測序長度較短，隨著讀段長度的增長，測序錯誤率會急劇增加。為了解決這個問題，目前主流的下一代測序平臺均采用雙端測序技術(shù)（paired-end sequencing）［3，6］。雙 端測序從DNA片段的兩端分別開始測序，并生成高質(zhì)量、可比對的測序數(shù)據(jù)。對于測序末端，由于測序長度過長而導致錯誤增加，可以根據(jù)另一條測序序列來進行糾正，因為此處位于另一條序列的前端，具有較低的測序錯誤率。相比于單端測序從DNA片段的一端開始測序，產(chǎn)生的每條讀段（read）之間沒有聯(lián)系，雙端測序技術(shù)在測序過程中是從待測序列片段的兩端各測序一次。在測序過程中，測序的準確率會隨著測序的進行而下降，即reads越往后面越不準確。雙端測序產(chǎn)生的數(shù)據(jù)是成對出現(xiàn)的，每條read都有與其相匹配的另外一條read相對應(yīng)。一般而言，雙端測序的兩條相對應(yīng)的reads的尾部會有重疊區(qū)域，定義為overlap，將雙端的reads進行比對可以拼接出更長的序列。對雙端測序的reads進行拼接是對整個基因組的序列的分析的第一步，并且更長的reads會顯著的提高基因組拼接質(zhì)量，因此，它的準確性對所有下游分析都至關(guān)重要。

為了將雙端測序產(chǎn)生的大量的reads拼接成更長的序列，已有一些拼接算法解決此問題，然而由于測序錯誤的存在，拼接仍然是一項有挑戰(zhàn)性的工作。即使雙端測序所得到的reads已經(jīng)足夠短了，比如常用的Illumina平臺為150 bp，但在overlap區(qū)域仍然有一定的測序錯誤率。當overlap區(qū)域出現(xiàn)了測序錯誤時，將兩條read拼接起來，便難以確定正確的overlap的區(qū)域；當read1和read2的overlap部分的堿基不匹配時，便難以確定哪個才是正確的。這給拼接造成了極大的困難。因此需要一種可以進行糾錯的DNA序列拼接算法。常用的基于雙端測序的拼接算法有Shera［7］、FLASH［8］、COPE［9］。Shera拼接時需要依賴fastq質(zhì)量值，且其速度較差；FLASH算法正確率較高、速度比Shera快，但是在使用條件上，它同樣依賴于fastq質(zhì)量值且有最短的overlap長度要求，同時在拼接之前需要借助工具Bowtie［10］對reads進行糾錯才能拼接；COPE需要利用kmer頻率信息和fastq質(zhì)量值去拼接，且具有較高的內(nèi)存需求和相對較長的執(zhí)行時間。這三種算法都是先遍歷雙端reads的每一個位點，尋找到合適長度的overlap，再通過測序質(zhì)量值去處理重疊區(qū)域里不相同的堿基。由于read的長度一般為75～300 bp，所以時間代價往往比較大，且拼接效果依賴于測序質(zhì)量值，而測序質(zhì)量值在理論上來說只是一種概率性的值，往往并沒有那么準確。因此從使用條件上來說，以上算法限制較多［11］。同時這些算法自身不具備糾錯能力，糾錯需要借助額外的工具［12］。

因此，為了解決現(xiàn)有的拼接算法在拼接時具有諸多限制條件，且拼接時間復雜度高的問題，本文在第二代測序技術(shù)的背景下，針對雙端測序技術(shù)所產(chǎn)生的兩條reads，充分利用兩條reads尾部的重疊序列，將確定overlap長度問題與處理不匹配堿基問題相統(tǒng)一，即將拼接與糾錯相結(jié)合，設(shè)計了一種基于Levenshtein距離的DNA序列拼接算法，簡稱LEDA算法。

1 Levenshtein距離

1.1 Levenshtein距離定義

Levenshtein距離即文本編輯距離，這一概念是由俄羅斯科學家Levenshtein于1965年提出的［13］，它用于兩個字符串之間差異程度的量測，量測方式是計算至少需要多少次的操作才能將一個字符串變成另一個字符串。給定兩個字符串A和B，將A轉(zhuǎn)換成B所需要刪除、插入等操作的集合就叫做A到B的編輯路徑。其中最短的編輯路徑就叫做字符串A和B的編輯距離。Levenshtein距離許可的編輯操作有：將一個字符替換（substitute，sub）成另一個字符，插入（insert，ins）一個字符，刪除（delete，del）一個字符。

編輯距離的應(yīng)用十分廣泛。Unix下的diff及patch命令也是利用編輯距離來進行文本編輯對比。DNA可以視為由A、C、G和T組成的字符串，因此，編輯距離可判斷兩個DNA的類似程度。

由于DNA序列測序過程中可能發(fā)生的三種錯誤sub，ins，del和Levenshtein距離的三種編輯操作相同，因此本文選擇Levenshtein距離設(shè)計相關(guān)算法。

1.2 Levenshtein距離求解原理

Levenshtein距離的求解采用動態(tài)規(guī)劃的思想，其基本思想是將原問題轉(zhuǎn)換為求解多個相似子問題，通過子問題的求解計算出原問題的解。在求解編輯距離時需要構(gòu)建編輯距離矩陣，逐行求出矩陣中每個元素的值，最終得到矩陣中最后一個元素的值，該值即為兩個字符串的編輯距離［14］。

設(shè)有兩個字符串A和B：A=a1，a2，…，am，B=b1，b2，…，bn，其長度分別是m和n。首先建立A和B的規(guī)模(m+1)×(n+1)的編輯距離矩陣D。設(shè)字符串A的前i個字符所組成的子字符串為A[1～i]，字符串B的前j個字符所組成的子字符串為B[1～j]。定義矩陣D如下

矩陣的每一個元素di，j表示字符串A的前i個字符組成的字符串A[1～i]和字符串B的前j個字符組成的字符串B[1～j]的編輯距離。我們可以對任意一個字符串A或B進行插入一個字符、刪除一個字符、替換一個字符三種操作。因此對于兩個字符串，共有六種操作方法。但其中，對字符串A刪除一個字符和對字符串B插入一個字符是等價的，比如當A=“doge”，B=“dog”時，既可以刪除A的最后一個字符e，得到相同的“dog”，也可以在B的末尾添加一個字符e，得到相同的“doge”；同理，對字符串B刪除一個字符和對字符串A插入一個字符是等價的；對A替換一個字符和對B替換一個字符是等價的，比如當A=“bat”，B=“cat”時，修改A的第一個字符b為c，和修改B的第一個字符c為b是等價的。因此，本質(zhì)上不同的操作只有在A中插入一個字符、在A中刪除一個字符、替換A的一個字符這三種。

假設(shè)已知di-1，j-1，di-1，j，di，j-1的值，則可以計算出di，j的值。對于di，j的求解，分別考慮以下三種不同的操作：

1）插入：把B的第j個字符插入到A的第i-1個字符后。由于已知經(jīng)過di-1，j次操作即可將A[1～i-1]轉(zhuǎn)換為B[j]，因此把B的第j個字符插入到A串的第i-1個字符后，即可將A[1～i]轉(zhuǎn)換為B[1～j]，此時di，j最小為di，j-1+1。

2）刪除：把A的第i個字符刪除。由于經(jīng)過di，j-1次操作即可將A[1～i]轉(zhuǎn)換為B[1～j-1]。因此把A的第i個字符刪除，即可使A[1～i]和B[1～j]相等，此時di，j最小為di，j-1+1。

3）替換：把A的第i個字符替換為B的第j個字符。由于經(jīng)過di-1，j-1次操作即可使A[1～i-1]和B[1～j-1]相等，此時把A的第i個字符替換為B的第j個字符，則可使A[1～i]和B[1～j]相等，此時di，j最小為di-1，j-1+1。特別地，若A的第i個字符和B的第j個字符原本就相同，則不需要進行替換操作，這種情況下，di，j最小可以為di-1，j-1。

對于邊界情況，即對于矩陣D的第一行和第一列，表示一個空字符串和一個非空字符串相互轉(zhuǎn)換的編輯距離，此時編輯距離為非空字符串的長度。

基于以上分析，得到求解編輯距離的公式如下

2 LEDA算法

將上述Levenshtein距離原理運用到DNA序列拼接中，原DNA序列overlap區(qū)域發(fā)生了插入和刪除錯誤造成了overlap長度的變化，原DNA序列overlap區(qū)域發(fā)生了替換錯誤造成了overlap區(qū)域堿基的變化，即原DNA序列overlap區(qū)域發(fā)生了插入、刪除和替換錯誤產(chǎn)生了overlap1和overlap2，將DNA序列拼接問題轉(zhuǎn)化成原序列片段發(fā)生了插入、刪除和替換三種錯誤生成的兩條序列的字符串比對問題。通過將DNA序列拼接問題與編輯距離相結(jié)合，比對雙端測序所產(chǎn)生的兩條DNA序列片段overlap1和overlap2，得到overlap1與overlap2的編輯距離，在由一個序列轉(zhuǎn)換成另一個序列的過程中，尋找所有可能正確的DNA序列片段，最后拼接成完整的DNA序列。

2.1 算法原理

本文所提出的LEDA算法是基于下一代測序技術(shù)所產(chǎn)生的讀段。目前下一代測序技術(shù)采用的是深度測序（sequencing depth），即對基因組進行多次測序，有時可達數(shù)百次甚至數(shù)千次。下一代測序技術(shù)所測得reads長度一般為幾十到幾百bp。

測序儀完成DNA序列的測序后會產(chǎn)生一系列原始文件，稱為fastq文件，該文件用于下一代測序數(shù)據(jù)的存儲。本文所提出的拼接算法需要將測序后生成的原始的fastq文件中讀入測序片段。通常，雙端測序所產(chǎn)生的一對reads的序列的關(guān)系如圖1所示。

圖1 雙端測序產(chǎn)生的一對reads序列Fig.1 A pair of reads generated by DNA paired-end sequencing

Read1和Read2的尾部會有overlap序列，可截取得到雙端的overlap序列，分別為overlap1和overlap2，其長度相等。如果在DNA序列沒有產(chǎn)生錯誤的情況下，overlap1和overlap2應(yīng)該與5-3正序overlap相同。如果overlap1和overlap2不相同，則說明DNA序列合成或者測序過程中產(chǎn)生了錯誤，可能是overlap1和overlap2都發(fā)生了替換/刪除/插入錯誤或其中一段無錯，另一段發(fā)生了替換/刪除/插入錯誤。

2.2 算法流程

設(shè)

其長度均為m。設(shè)序列overlap1的前i個堿基所組成的子序列為overlap1[1～i]，序列overlap2的前j個堿基所組成的子序列為overlap2[1～j]。

2.2.1 構(gòu)造編輯距離矩陣min_ed

通過動態(tài)規(guī)劃的思想構(gòu)造overlap1和overlap2的階為(m+1)×(m+1)的編輯距離矩陣min_ed。定義矩陣min_ed為：

矩陣中的每一個元素min_edi，j代表從overlap1[1～i]變換至overlap2[1～j]所需要的最少變換次數(shù)。

矩陣min_ed的構(gòu)造方法如下：

1）矩陣的第0行第j列，代表由‘’變換到overlap2[1～j]，需要經(jīng)過j次插入堿基的操作，共需要插入j個堿基。即min_ed0，j=j，‘’代表空堿基序列；

2）矩陣的第i行第0列，代表由overlap1[1～i]變換到‘’，需要經(jīng)過i次刪除堿基的操作，共需要刪除i個堿基，即min_edi，0=i；

3）矩陣的第i行第j列，代表由overlap1[1～i]變換到overlap1[1～j]，需要經(jīng)過的操作次數(shù)為min_edi-1，j+flag、min_edi，j-1+flag、min_edi-1，j-1三者中的最小值。其中，當且僅當overlap1i和overlap2j相等時，flag為1；否則，flag為0。

根據(jù)以上步驟得到矩陣min_ed，若min_edi，j=0，則說明overlap區(qū)域未發(fā)生任何錯誤，因此直接跳轉(zhuǎn)到2.2.4節(jié)的步驟進行序列拼接。

2.2.2 構(gòu)造編輯路徑矩陣op

利用動態(tài)規(guī)劃的思想構(gòu)造overlap1和overlap2的階為(m+1)×(m+1)的編輯路徑矩陣op。矩陣op定義如下

矩陣op中的每一個元素opi，j為一個集合，代表從overlap1[1～i]變換到overlap2[1～j]共有in_op1，in_op2，…，in_opu一共u種不同的編輯路徑，即

其中每一條編輯路徑in_op代表將overlap1[1～i]變換到overlap2[1～j]的v次操作，in_op={error_correct1,error_correct2,…,error_correctv},v≥1其中每一個元素error_correcti代表一項具體的操作，可能是無錯/替換/刪除/插入，即

其中，‘0’代表無錯操作，‘1’代表替換操作，‘2’代表刪除操作，‘3’代表插入操作。

矩陣op的構(gòu)造需要采用動態(tài)規(guī)劃的思想逐行求出矩陣中每個元素的值，op由矩陣min_ed以及overlap1和overlap2共同決定，其具體構(gòu)造方法如下：

1）矩陣op的第0行第j列，表示從空堿基序列轉(zhuǎn)換到非空堿基序列的編輯路徑，此時只有一種編輯路徑，即op0，j={‘{3’，‘3’，…，‘3’}}，一共j個‘3’。

2）矩陣op的第i行第0列，表示將一個非空堿基序列轉(zhuǎn)換為空堿基序列的編輯路徑，此時只有一種編輯路徑，即opi，0={‘{2’‘，2’，…，‘2’}}，一共i個‘2’。

3）矩陣的第i行第j列，此時已知opi-1，j，opi，j-1，opi-1，j-1的值，則可以求出opi，j的值，opi，j的求解需要先考慮overlap1i與overlap2j是否相同：

①若overlap1i=overlap2j，則將opi-1，j-1里包含的所有的集合均添加一個無錯操作（操作‘0’）后，添加到opi，j中；

②若overlap1i≠overlap2j，則將opi-1，j-1里包含的所有的集合均添加一個替換操作（操作‘1’）后，添加到opi，j中；

接著分別考慮min_edi-1，j和min_edi，j-1是否等于min_edi，j+1：

①若min_edi，j=min_edi-1，j+1，則將opi-1，j里包含的所有的集合均添加一個刪除操作（操作‘2’）后，添加到opi，j中；

②若min_edi，j=min_edi，j-1+1，則將opi，j-1里包含的所有的集合均添加一個插入操作（操作‘3’）后，添加到opi，j中。

2.2.3 糾錯恢復正確的overlap序列t

利用矩陣op的最后一個元素opm，m去糾錯恢復正確的overlap序列，t表示原始的正確的overlap序列，t的構(gòu)造方式如下：

已 知opm，m={in_op1，in_op2，…，in_opu}，遍 歷其中的每一個元素in_op，進行如下操作：

1）初始化t為空堿基序列，即t=‘’。

2）已知in_op={error_correct1，error_correct2，…，error_correctv}，v≥1，遍歷其中的每一個元素error_correcti(1≤i≤v)：

①若error_correcti為‘0’，說 明overlap1i和overlap2i相同，無需進行糾錯，將overlap1i添加到堿基序列t中。

②若error_correcti為‘1’，說 明overlap1i或overlap2i在此處產(chǎn)生了sub錯，在此處可能是overlap1i正確或者是overlap2i正確，因此將overlap1i或者是將overlap2i添加到t中。

③若error_correcti為‘2’，說明可能是overlap1i在此處發(fā)生了ins錯添加了一個堿基，或者是overlap2i在此處發(fā)生了del錯刪除了一個堿基。此時將overlap1的長度減1，或?qū)verlap1i添加到t。

④若error_correcti為‘3’，說明可能是overlap1i在此處發(fā)生了del錯刪除了一個堿基，或者是overlap2i在此處發(fā)生了ins錯添加了一個堿基。此時將overlap2i添加到t中，或?qū)verlap2的長度減1。

每一步操作后便進入一次遞歸，直至窮舉完所有可能的情況。當遍歷完in_op的所有元素時，便將此時的t添加到結(jié)果集中。

上述操作如圖2所示。

圖2 基于op m,m尋找所有可能正確序列的流程圖Fig.2 The flow chart based on op m,m to find all possible correct sequences

3）當遍歷完opm，m的所有元素時，即可得到所有可能正確的原始的overlap序列。

2.2.4 拼接

將所有可能正確的DNA序列片段t與app1、app2拼接，其中，app1是Read1剩下的序列片段，app2是Read2剩下的序列片段經(jīng)過圖3所示的操作后得到的序列片段。如圖3所示。

圖3 DNA序列拼接Fig.3 DNA sequence assembling

2.3 算法偽代碼

上述2.2節(jié)LEDA算法流程的偽代碼如算法1所示。

3 實驗部分

為了驗證本文的研究成果，我們做了仿真數(shù)據(jù)實驗與真實數(shù)據(jù)實驗，在仿真數(shù)據(jù)實驗中從改變錯誤類型與改變錯誤數(shù)量兩個維度進行實驗，以評估本文所提出的拼接算法的糾錯能力。在真實數(shù)據(jù)實驗中，主要通過拼接正確率和時間復雜度評估本算法的可行性。

3.1 仿真數(shù)據(jù)實驗

3.1.1 實驗環(huán)境

操作系統(tǒng)Windows10 64位，處理器Intel（R）Core（TM）i5-10400 CPU@2.90 GHz 2.90 GHz，RAM 16.00 GB，編譯環(huán)境Python3.7。

3.1.2 結(jié)果及分析

仿真實驗中，我們從改變錯誤類型以及改變錯誤數(shù)量兩個方面設(shè)計實驗，以分析LEDA算法對于測序出錯的序列的拼接與糾錯能力。本文構(gòu)造了20 000條含有200個堿基對的隨機DNA序列作為仿真實驗的DNA模板鏈，然后模擬下一代測序技術(shù)生成20 000對reads。對正確的reads進行注錯測試，然后用LEDA算法進行拼接，將拼接后的序列與原先正確的DNA模板序列進行比對以檢驗算法的拼接正確率以及糾錯能力。

由于在雙端測序的過程中，越接近測序序列的尾部出錯率越高，且本文所設(shè)計的算法僅對于overlap部分的序列有糾錯能力，所以實驗中對序列注錯的位置均在overlap部分。在實驗中，將無錯操作看作特殊的替換，即替換為自身。

1）錯誤類型的影響

對于Read1和Read2的overlap部分的序列分別進行1個錯的隨機注錯測試，分6種不同錯誤類型討論，見表1。

由于DNA雙鏈的互補與對稱性，6種分類對于Read1與Read2的順序沒有要求，因此無需交換Read1與Read2的順序再次實驗。

本實驗計算了在6種不同錯誤類型、錯誤發(fā)生位置隨機且錯誤數(shù)量不超過兩個的情況下，20 000條含有200個堿基對的隨機DNA序列利用該算法拼接的正確率以及運行時間。測試結(jié)果如表1。

表1 不同錯誤類型DNA序列拼接正確率和時間Table 1 DNA sequence assembly success rate and time with different error types

分析結(jié)果可得，上述條件下的隨機DNA序列拼接平均正確率為0.975 6，平均運行時間為15.148 min，平均每對reads拼接所用時間為45.444 ms。由此可得本文算法可以處理各種類型的錯誤，拼接高效且正確率高。

2）錯誤數(shù)量的影響

對于Read1和Read2的overlap部分的序列分別進行錯誤數(shù)量為1、2、3、4個錯誤的隨機注錯測試。

本實驗計算了錯誤數(shù)量從1增加到4的情況下，20 000條含有200個堿基對的隨機DNA序列利用本文算法拼接的正確率以及運行時間。測試結(jié)果如表2。

表2 不同錯誤個數(shù)DNA序列拼接正確率和時間Table 2 DNA sequence assembly correct rate and time of differ ent number of err ors

分析結(jié)果可得，上述條件下DNA序列拼接平均正確率為0.971 2，平均運行時間為15.043 min，平均每對reads拼接所用時間為45.129 ms。20 000條含有200個堿基對的隨機DNA序列拼接正確率隨著錯誤數(shù)量的增加而下降，但仍在可以接受的范圍，因為考慮到下一代測序錯誤率僅為0.5%～1.0%，對于實驗讀段長度為150，錯誤個數(shù)不超過1.5個，可見LEDA算法完全可以應(yīng)對實際測序中所出現(xiàn)的錯誤，在實際應(yīng)用中性能較好。

3.2 真實數(shù)據(jù)實驗

為了驗證本文所述算法在實際應(yīng)用中的有效性，我們不僅做了仿真實驗，而且也做了真實的測序?qū)嶒灐Ｕ鎸嵉臏y序數(shù)據(jù)由二代測序的Illumina的Hiseq4000測序平臺測序所得，測序文庫是一個有11 520條長度為180 nt（nt即核苷酸）DNA的文庫，其reads長度為150 nt，為雙端讀段，所有的Read1和Read2分別由兩個fastq文件存儲。Read1和Read2分別有5 003 753條，大小共為3.36 GB。實驗環(huán)境與3.1節(jié)仿真實驗環(huán)境相同。

3.2.1 實驗過程

1）數(shù)據(jù)處理

①處理fastq文件，讀取其中的Read1和Read2，并分別存儲在不同的列表Read1和Read2里。

②考慮到第二代測序的錯誤率很低，僅為0.5%～1.0%，對于實驗中讀段長度為150，錯誤個數(shù)大都不超過2個，而且下一代測序采用的是深度測序，對于本實驗其平均測序深度約為434，即平均每個堿基都會被測序434次。因此，實驗中篩選出Read1與Read2尾部重疊部分即overlap1與overlap2的編輯距離小于等于1的序列，得到3 928 926對Read1和Read2序列，占總序列78.52%。

③為了確保實驗的可靠性，我們也從實驗中篩選出Read1與Read2尾部重疊部分即overlap1與overlap2的編輯距離小于等于2的序列，最終得到4 459 805對Read1和Read2序列，占總序列的89.13%。

可見，編輯距離控制在1或2以內(nèi)，便可以涵蓋fastq文件中的大部分的數(shù)據(jù)，且由于二代測序的深度測序特性，即對于同一段序列可能會測出好幾百倍甚至幾千倍的結(jié)果，能夠保證篩選出來的序列可以覆蓋原始序列的所有堿基。

2）程序運行

①輸入FASTQ文件處理后產(chǎn)生的兩個列表Read1和Read2。

②由于有四種不同的引物，分別為模板鏈兩端原本的引物以及兩端原本引物的反向互補序列。引物不同會造成Read1與Read2拼接的方式會有所不同，所以實驗中需要先判斷Read1的引物類型，然后再分類型討論返回結(jié)果。

③將運行得到的所有可能正確的DNA序列與測序文庫比較，檢驗算法拼接正確率。

3.2.2 實驗結(jié)果

測序文庫中模板DNA序列共有11 520條，對于編輯距離小于等于1的reads，其拼接結(jié)果可以成功比對其中11 519條，即幾乎可以拼接得到原來文庫中的所有DNA序列，正確率為0.999 9，其總運行時間為2 393.809 6 min，平均每對reads拼接所用時間為36.56 ms。

為了保證實驗的可靠性，我們也拼接了編輯距離小于等于2的Reads，其拼接結(jié)果也為11 519條，與編輯距離控制在1的拼接結(jié)果是相同的。

4 討論

為了能更好地分析對比LEDA算法的拼接效果，我們將Shera、FLASH、COPE算法在3.2節(jié)中真實數(shù)據(jù)上運行，4種算法結(jié)果如表3所示。

從表3中可以看出，在使用條件上，LEDA算法無任何限制，而Shera算法需要fastq文件中質(zhì)量值的信息，F(xiàn)LASH算法不僅需要fastq文件中質(zhì)量值的信息并且在拼接之前必須要使用額外的糾錯工具Bowite對reads進行錯誤糾正以提高拼接正確率，COPE算法需要fastq文件中質(zhì)量值的信息以及k-mer頻率信息對reads進行錯誤糾正。

表3 不同算法的使用條件、運行時間、時間復雜度以及拼接正確率Table 3 Use conditions、running time、time complexity and assembly correct rate of different algorithms

在正確率相同均為99.99%的情況下，LEDA算法平均每對reads拼接所用的時間比Shera、FLASH、COPE算法都要短。這是由于相比于其他拼接算法的時間復雜度與read的長度呈二次方的時間復雜度，LEDA算法與read長度并無直接關(guān)系，而是與overlap長度和編輯距離有關(guān)，時間復雜度為O（n·2x）。而實際情況中read的長度都大于overlap長度，且編輯距離都比較小，因此本算法的時間復雜度優(yōu)于其他三種算法。

5 結(jié)語

相比于其他的一些DNA序列拼接算法，需要遍歷Read1和Read2的所有可能組合的位點尋找所有可能的overlap，復雜且代價比較大，且使用時有諸多限制條件，本文提出的LEDA算法，將問題轉(zhuǎn)化為原序列的overlap片段發(fā)生了ins、del和sub三種錯誤生成兩條序列的字符串比對，拼接與糾錯相結(jié)合，使用本文算法時不需要依賴其他任何工具做任何預(yù)處理工作，且對于overlap長度沒有限制，不需要額外讀段信息。在實際應(yīng)用中使用簡單且運行高效。在DNA序列拼接方面首次引入了Levenshtein距離的思想，不依賴于額外的信息，大大簡化了已有的拼接算法依靠測序質(zhì)量值等其他測序信息計算概率進行拼接的思路，LEDA算法本身就具備一定的糾錯能力。

目前LEDA算法僅僅只能針對雙端測序，無法適用于單端測序技術(shù)。另外，LEDA算法的高效是建立在需要嚴格控制overlap1和overlap2的編輯距離的基礎(chǔ)上。LEDA算法無法解決整個基因組的拼接問題，只能解決基礎(chǔ)的雙端測序的拼接問題。未來，我們將繼續(xù)探索將LEDA算法擴展，使其能用于整個基因組的拼接之中。