崔競(jìng)松，薛 慧，王蘭蘭，郭 遲

1.空天信息安全與可信計(jì)算教育部重點(diǎn)實(shí)驗(yàn)室，武漢大學(xué)國(guó)家網(wǎng)絡(luò)安全學(xué)院，湖北武漢 430072；

2.河海大學(xué) 理學(xué)院，江蘇 南京 211100；

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

0 引言

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

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

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

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

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

1 Levenshtein距離

1.1 Levenshtein距離定義

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

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

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

1.2 Levenshtein距離求解原理

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

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

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

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

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

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

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

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

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

2 LEDA算法

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

2.1 算法原理

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

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

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

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

2.2 算法流程

設(shè)

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

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

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

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

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

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

2）矩陣的第i行第0列，代表由overlap1[1～i]變換到‘’，需要經(jīng)過i次刪除堿基的操作，共需要?jiǎng)h除i個(gè)堿基，即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三者中的最小值。其中，當(dāng)且僅當(dāng)overlap1i和overlap2j相等時(shí)，flag為1；否則，flag為0。

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

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

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

矩陣op中的每一個(gè)元素opi，j為一個(gè)集合，代表從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其中每一個(gè)元素error_correcti代表一項(xiàng)具體的操作，可能是無錯(cuò)/替換/刪除/插入，即

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

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

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

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

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

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

②若overlap1i≠overlap2j，則將opi-1，j-1里包含的所有的集合均添加一個(gè)替換操作（操作‘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里包含的所有的集合均添加一個(gè)刪除操作（操作‘2’）后，添加到opi，j中；

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

2.2.3 糾錯(cuò)恢復(fù)正確的overlap序列t

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

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

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

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

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

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

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

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

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

上述操作如圖2所示。

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

3）當(dāng)遍歷完opm，m的所有元素時(shí)，即可得到所有可能正確的原始的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í)驗(yàn)部分

為了驗(yàn)證本文的研究成果，我們做了仿真數(shù)據(jù)實(shí)驗(yàn)與真實(shí)數(shù)據(jù)實(shí)驗(yàn)，在仿真數(shù)據(jù)實(shí)驗(yàn)中從改變錯(cuò)誤類型與改變錯(cuò)誤數(shù)量?jī)蓚€(gè)維度進(jìn)行實(shí)驗(yàn)，以評(píng)估本文所提出的拼接算法的糾錯(cuò)能力。在真實(shí)數(shù)據(jù)實(shí)驗(yàn)中，主要通過拼接正確率和時(shí)間復(fù)雜度評(píng)估本算法的可行性。

3.1 仿真數(shù)據(jù)實(shí)驗(yàn)

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

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

1）錯(cuò)誤類型的影響

對(duì)于Read1和Read2的overlap部分的序列分別進(jìn)行1個(gè)錯(cuò)的隨機(jī)注錯(cuò)測(cè)試，分6種不同錯(cuò)誤類型討論，見表1。

由于DNA雙鏈的互補(bǔ)與對(duì)稱性，6種分類對(duì)于Read1與Read2的順序沒有要求，因此無需交換Read1與Read2的順序再次實(shí)驗(yàn)。

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

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

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

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

對(duì)于Read1和Read2的overlap部分的序列分別進(jìn)行錯(cuò)誤數(shù)量為1、2、3、4個(gè)錯(cuò)誤的隨機(jī)注錯(cuò)測(cè)試。

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

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

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

3.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

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

3.2.1 實(shí)驗(yàn)過程

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

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

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

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

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

2）程序運(yùn)行

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

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

③將運(yùn)行得到的所有可能正確的DNA序列與測(cè)序文庫(kù)比較，檢驗(yàn)算法拼接正確率。

3.2.2 實(shí)驗(yàn)結(jié)果

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

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

4 討論

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

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

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

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

5 結(jié)語(yǔ)

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

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