一種高效的前向糾錯碼桶分配DNA存儲解碼方法

昝鄉(xiāng)鎮(zhèn) 姚翔宇 許 鵬 陳智華 石曉龍 李樹棟 劉文斌*

①(廣州大學(xué)計算科技研究院 廣州 510006)

②(廣州大學(xué)網(wǎng)絡(luò)空間技術(shù)先進(jìn)研究院 廣州 510006)

1 引言

云計算、大數(shù)據(jù)等技術(shù)的發(fā)展，人類存儲數(shù)據(jù)的需求呈現(xiàn)出指數(shù)級增長的趨勢。據(jù)國際數(shù)據(jù)公司預(yù)測[1]，2025年全球數(shù)據(jù)總量預(yù)計達(dá)到175 ZB。傳統(tǒng)存儲介質(zhì)技術(shù)在滿足未來數(shù)據(jù)存儲需求方面逐漸暴露出一系列缺點[2,3]，比如有效存儲時間短、數(shù)據(jù)易損壞以及維護(hù)成本高······與此同時，攜帶有遺傳信息的脫氧核糖核酸(DeoxyriboNucleic Acid,DNA)，因其具有超高的存儲密度、低維護(hù)成本以及數(shù)據(jù)保存持久等特點[4—6]，有望是一種極具潛力的存儲介質(zhì)，解決海量數(shù)據(jù)存儲面臨的困境。

DNA存儲主要涉及合成、聚合酶鏈反應(yīng) (Polymerase Chain Reaction, PCR)擴(kuò)增、測序等生物過程。由于技術(shù)局限，這些過程會導(dǎo)致DNA存儲發(fā)生一系列復(fù)雜組合錯誤，從而給數(shù)據(jù)的可靠恢復(fù)帶來了挑戰(zhàn)[7,8]。據(jù)估計，二代測序技術(shù)和陣列合成技術(shù)的錯誤發(fā)生概率為1%～2%[9]，三代測序技術(shù)將達(dá)10%～15%[10]。與傳統(tǒng)數(shù)字存儲相比，DNA存儲序列的堿基錯誤不僅有替換錯誤，還包括插入錯誤。這些錯誤的復(fù)雜交織遠(yuǎn)遠(yuǎn)超出了傳統(tǒng)糾錯碼(Error Correcting Codes, ECC)[11—14]的糾錯能力。同時，堿基插入還會導(dǎo)致DNA序列的長度與標(biāo)準(zhǔn)長度不一致。有研究表明，三代納米孔測序技術(shù)會產(chǎn)生大約88%的非標(biāo)準(zhǔn)長度序列[15]。以往利用RS糾錯碼[16,17]或噴泉碼[18,19]的DNA存儲方法中，通常會舍棄掉這些測序讀段，從而導(dǎo)致大量浪費，并增加了測序過程的時間和費用。因此，研究面向DNA堿基插入、刪除和替換組合錯誤的高效信息恢復(fù)方法是未來DNA存儲亟待解決的一個重要問題。

Blawat等人[13]對每個字節(jié)設(shè)計了兩套DNA編碼，然后交替使用第1類和第2類DNA編碼來編碼原始信息。當(dāng)插入錯誤發(fā)生時，將會打破兩類DNA編碼交替出現(xiàn)的規(guī)律，由此可以定位發(fā)生插入/刪除的位置。Press等人[20]通過哈希操作將二進(jìn)制序列的每一比特，都與其附近比特、序列索引產(chǎn)生強(qiáng)關(guān)聯(lián)后，再進(jìn)行編碼。在解碼的時候，通過貪心窮舉搜索策略評估每一比特滿足強(qiáng)關(guān)聯(lián)的程度，進(jìn)而完成堿基插入、刪除和替換等錯誤的糾正。但是該方案需要較高的冗余度，且解碼過程復(fù)雜。Xue等人[21]通過添加一些冗余位，將二進(jìn)制序列拆分成兩個子串，使其滿足萊文斯坦碼(Levenshtein code)的數(shù)據(jù)形式，然后利用萊文斯坦碼可以糾正1 bit插入/替換的性質(zhì)，完成DNA序列中1位堿基錯誤的糾正。天津大學(xué)Song等人[22]通過構(gòu)建多拷貝讀段的德布萊英圖(de Bruijn graph)，將一致性序列的確定問題轉(zhuǎn)換為圖中的最大權(quán)路徑搜索問題，從而過濾掉插入、刪除、替換導(dǎo)致的低頻子串路徑。本研究團(tuán)隊[23]最近提出一種基于前向糾錯碼的英文文本3層糾錯方法，基本思想是通過前向糾錯碼對DNA讀長進(jìn)行初步糾錯，然后對轉(zhuǎn)化的字符序列進(jìn)行多序列比對糾錯，最后通過單詞拼寫進(jìn)一步糾正錯誤。該方法在錯誤率0.05情況下，20X測序深度糾錯準(zhǔn)確率達(dá)90%以上。但在錯誤率0.10情況下，60X測序深度僅達(dá)到64%。因此，難以適應(yīng)高錯誤率的情況。

本文在前向糾錯碼的基礎(chǔ)上，通過在序列編碼中采用“索引+CRC哈希+索引”模式，提高測序讀長的聚類精度；然后，提出了基于可識別DNA碼的桶式分配策略的糾錯算法。仿真結(jié)果表明在0.1和0.05錯誤率條件下，平均解碼準(zhǔn)確率在20X測序深度時可達(dá)94%以上；在0.15錯誤率條件下，平均解碼準(zhǔn)確率在60X測序深度時可達(dá)90%以上。

2 方法

2.1 編碼方法

表1為本文使用的英文文本前向糾錯編碼表[23]。該編碼表包含26個常規(guī)DNA編碼序列(表1中白色底紋部分)以及4個特殊DNA編碼序列(表1中陰影部分)。常規(guī)DNA編碼序列用于編碼英文字母、標(biāo)點符號和數(shù)字等字符。特殊DNA編碼序列包括大寫鍵、標(biāo)點符號鍵、數(shù)字鍵和空格鍵。除空格鍵外，其他特殊鍵用于標(biāo)記位于其后的下一個DNA編碼序列編碼何種字符(大寫字母、數(shù)字字符或標(biāo)點符號字符)，例如編碼序列“CTTGTC ACACAC”表示編碼的字符為數(shù)字字符“6”。需要說明的是，前一個編碼序列不是特殊鍵的DNA編碼序列編碼的為小寫英文字母。

該編碼表具有如下特點。首先，編碼表的設(shè)計遵循了生物序列的約束，比如任意兩個DNA編碼序列拼接不會產(chǎn)生長度大于2的均聚物、鳥嘌呤和胞嘧啶 (Guanine and Cytosine, GC)含量保持平衡且分布均勻，有利于減少DNA分子在DNA存儲過程中的錯誤。其次，該編碼表任意兩個DNA編碼序列的漢明距離都至少為3，因此具有1位堿基替換糾錯的能力。此外，該編碼表435對序列的平均編輯(插入、刪除、替換)距離為3.85， 僅有12對編碼間的編輯距離為2。因此，可以近似認(rèn)為該編碼表具有一位的糾錯能力。

2.2 解碼方法

2.2.1 分組策略

DNA存儲解碼的第1步是對測序讀長(reads)分組或聚類，即將屬于同一編碼序列的測序讀長劃分為一組，為后續(xù)基于多序列的一致性序列推斷奠定基礎(chǔ)。由于測序讀長中的各種錯誤，分組的一個目的是精度高，盡可能減少將其他序列的測序讀長錯誤加入分組；另一個是召回率高，即盡可能將屬于同一個序列的測序讀長判別出來并分為一組。在機(jī)器學(xué)習(xí)領(lǐng)域，這兩個指標(biāo)往往難以同時滿足，需要進(jìn)行一定的折中。

本文采用圖1所示的序列設(shè)計，在存儲序列前后各加一個該序列的索引，再加一個索引的循環(huán)冗余校驗碼 (Cyclic Redundancy Check, CRC)。圖1中索引值和CRC校驗值的編碼按照兩個連續(xù)比特編碼成一個DNA堿基[24]。分組原則是：(1)如果索引1與索引2相同，則按索引1分組；(2)如果索引1與索引2不同，則選擇與CRC校驗值一致的索引分組。以上兩個條件不滿足就丟棄該序列。

本質(zhì)上這一方法屬于基于測序讀長索引值直接分組的方法，其時間復(fù)雜度為O(N)(N為測序讀長的總數(shù))。另一種分組方法是直接基于序列比對的相似性分組方法，其時間復(fù)雜度為O(N2n2)(n為測序長度)。相比于前者，后者的召回率相對較高，但時間復(fù)雜度大。特別是DNA存儲中N為百萬級別數(shù)量時，所需時間將難以想象。

2.2.2 桶式糾錯策略

圖2(a)給出了一條測序讀長可能發(fā)生錯誤情況的示意圖。從編碼單元的角度看，按照其受影響的程度可以分為3類：(1)第1類、完全正確，沒有受到錯誤影響(深綠色)；(2)僅受到1位插入/刪除/替換的影響(淺綠色)；(3)大于1位插入/刪除/替換的影響(白色)。由于本文所用的30個前向糾錯編碼具有1位糾錯能力，本文對一個測序讀長觀測到的長度為5,6或7的子串有如下兩個假設(shè)：

(1)如果一個6堿基子串是合法編碼，它以很高的概率屬于第1類編碼；

(2)如果一個子串與合法編碼的編輯距離為1，它以較高的概率屬于第2類編碼。

本文將在一個測序讀長中出現(xiàn)的上述兩種字串稱為可識別DNA碼?；谏鲜稣J(rèn)識，本文提出一個基于可識別DNA碼桶分配策略的糾錯算法，基本思想是搜索一個分組中的所有測序讀長的可識別DNA碼，根據(jù)其在測序讀長的位置分配合適的編碼位置(即桶)，最后根據(jù)每個桶中的編碼投票確定最終的編碼。

可識別DNA碼的搜索方法如下：

(1)搜索長度為5, 6和7的DNA子串并計算其與編碼表的最小編輯距離；

(2)如果存在可識別DNA碼，則確定第1個堿基位置，從該可識別DNA碼后面的堿基重復(fù)(1)；

(3)如果不存在可識別DNA碼，則從當(dāng)前位置前進(jìn)一個堿基，重復(fù)(1)。

重復(fù)上述步驟直到掃描完整個測序讀長，即可得到其中的所有可識別DNA碼。需要說明的是，上述過程每個可識別DNA碼將根據(jù)其最小編輯距離賦值不同的權(quán)重。當(dāng)最小編輯距離為0，則權(quán)重設(shè)置為1，否則設(shè)置為0.1。為了提高計算效率，本文提前編制好一個長度為5, 6, 7的所有DNA串與30個合法編碼的最小編輯距離表，上述搜索過程的最小編輯距離直接查表即可得到，避免了重復(fù)計算。

如果本文將每個合法編碼位置當(dāng)作一個桶，對一個序列分組中所有DNA測序讀長進(jìn)行解碼的過程可以描述為：將測序讀長中的可識別DNA碼按照其對應(yīng)的合法編碼放入對應(yīng)的桶，最后根據(jù)每個桶中編碼權(quán)重進(jìn)行投票，即可確定該序列的可能編碼。

3 實驗結(jié)果

本文仿真實驗采用《老人與?！泛汀读_伯特·路易斯·斯蒂文森評傳》的英文文本，總文件大小為324 kB。編碼英文文本的DNA存儲行的長度為208堿基，其中數(shù)據(jù)域長度為180個堿基，索引1和索引2均為10堿基，CRC校驗值為8堿基。最終形成11637條DNA序列。仿真實驗的錯誤率分別為0.05, 0.1和0.15。每種錯誤率下，插入:替換:刪除的比例設(shè)置為1 : 2 : 2, 1 : 1 : 1以及2 : 2 : 1 3種情況。測序深度分別為20, 30, 40, 50及60。每組參數(shù)重復(fù)1000次。仿真實驗的配置為Intel(R) Xeon(R)Silver 4210 CPU @ 2.20 GHz處理器、30 GB內(nèi)存的服務(wù)器，軟件環(huán)境為CentOS Linux release 7.6系統(tǒng)。

3.1 分組策略性能分析

圖3分別給出了“索引”、“索引+CRC”、“索引+CRC+索引”3種分組策略的平均精度和平均召回率。從圖3(a)可以看出，后兩種分組的平均精度基本接近100%，且明顯高于簡單索引分組策略。這主要是因為CRC或索引的校驗作用明顯提高了索引分組的精度。此外，簡單索引分組的精度隨錯誤率增加而降低。

圖3(b)的平均召回率表明“索引+CRC+索引”分組的召回率明顯高于“索引+CRC”。這主要是因為只要有一個索引通過CRC檢驗，即可以以很高的概率保證分組的正確性，因而提高了召回率。和簡單索引相比，“索引+CRC+索引”的召回率隨錯誤率增加而逐漸降低。

以錯誤率0.15為例，“索引+CRC+索引”分組的平均精度約為100%，平均召回率約為11%；簡單索引的平均精度為33%，平均召回率約為20%。前者的召回率雖然約為后者的1/2，但是基本都是正確分組。而后者約有67%來自其他分組的錯誤測序讀長，這將造成后面一個桶中會有大量不正確可識別DNA碼，最終導(dǎo)致投票失敗。因此，“索引+CRC+索引”分組策略既保障了精度又適當(dāng)提升了召回率，為后面桶式糾錯奠定了關(guān)鍵的基礎(chǔ)。

3.2 糾錯策略性能分析

圖4(a)是不同測序深度情況下的解碼平均準(zhǔn)確率?？梢钥闯觯?1)當(dāng)錯誤率為0.05和0.10，本文方法在測序深度20時平均準(zhǔn)確率就達(dá)到94%以上。但是隨著測序深度增加，平均準(zhǔn)確率基本不變。這可能與DNA存儲測序的分布不均性有關(guān)。這里存儲測序的分布不均勻主要包括兩個方面：一是序列中堿基錯誤分布的不均勻(仿真數(shù)據(jù)里表現(xiàn)為序列中堿基錯誤隨機(jī)發(fā)生的不均勻性)；二是測序序列拷貝數(shù)分布的不均勻(仿真數(shù)據(jù)里表現(xiàn)為編碼序列的抽樣分布不均勻)。DNA存儲測序的不均性，導(dǎo)致了無論測序深度多高，總是有些序列因為拷貝數(shù)太少以及序列隨機(jī)錯誤發(fā)生的不均勻性，導(dǎo)致不能準(zhǔn)確解碼，進(jìn)而影響了平均準(zhǔn)確率。當(dāng)錯誤率增加到0.15，平均準(zhǔn)確率極具下降，測序深度20時的平均準(zhǔn)確率約為70%。隨著測序深度的增加，在測序深度60時就達(dá)到90%。(2)插入刪除比例對糾錯性能有一定影響，當(dāng)錯誤率為0.05和0.10，刪除比例較大時的平均精度較高。這可能是低錯誤率刪除錯誤對合法編碼的影響小于插入的破壞程度。當(dāng)錯誤率為0.15時，刪除比例較大時的平均精度較低。這說明高錯誤率刪除錯誤對合法編碼的影響大于插入的破壞程度。例如，合法編碼“ATGAGC”，兩位堿基缺失后的編碼可能為“ATGA”, “ATGC”,“AGAC”, ···，任意兩個位置的堿基缺失，都造成了合法編碼本身信息的破壞，每種破壞均有可能導(dǎo)致糾正錯誤(注：高錯誤率下插入錯誤或刪除錯誤導(dǎo)致發(fā)生錯誤的合法DNA碼普遍出錯的堿基數(shù)量大于等于2)。而合法編碼兩位堿基插入錯誤的引入，比如“ATGAGCXX”, “XATGAGCX”,“XXATGAGC”, “ATXGXAGC”, ··· (X為插入堿基)，并不會破壞合法編碼本身的信息。此外在合法編碼所有兩位插入錯誤的種類中，存在少量種類是可以正確識別并糾正的，比如“A TGAGCXX”, “XATGAGCX”和“XXATGAGC”。這表明，缺失兩個堿基的合法DNA編碼序列糾正失敗的概率，要大于插入兩個堿基的合法DNA編碼糾正失敗的概率。這就導(dǎo)致給定一高錯誤率，刪除錯誤比例較大的情況下可識別DNA碼的識別準(zhǔn)確率低于插入比例較大情況下的可識別DNA碼的識別準(zhǔn)確率。

圖4(b)是不同測序深度的情況下的解碼平均運行時間?？梢钥闯觯?1)隨著測序深度的增加，算法運行時間近似線性增加；(2)在給定測序深度下，算法運行時間隨錯誤率增加而減少。這主要是由于錯誤率對分組測序讀長數(shù)量的影響導(dǎo)致。圖3(b)顯示在0.05, 0.10和0.15時，分組讀長數(shù)量占測序深度的百分比大致分別為0.60, 0.27和0.11，基本與相應(yīng)錯誤率下的時間比一致。

3.3 與其他方法的比較

表2給出了不同方法的糾錯策略、插入/刪除糾錯、覆蓋率、錯誤率、準(zhǔn)確率和存儲模型。與文獻(xiàn)[25,26]的文本存儲工作相比，本文提出的方法可以對包括插入、刪除和替換在內(nèi)的3種錯誤進(jìn)行糾正，而其他兩種方法則沒有這種能力，這限制了它們只能在噪聲非常低(≤2%)的情況下應(yīng)用且需要較高的測序深度。本文方法在20X測序深度下，在錯誤率10%的條件下能恢復(fù)94%以上的數(shù)據(jù)。而文獻(xiàn)[25]恢復(fù)錯誤率約為2%的數(shù)據(jù)所需要的測序深度約為2700X，這對于未來的大數(shù)據(jù)存儲是不可行的。

表2 與其他方法的比較

與文獻(xiàn)[20—23]中的4種一般DNA存儲方法相比，本文方法比文獻(xiàn)[21]的方法更強(qiáng)大，文獻(xiàn)[21]只能糾正1位堿基的插入。在錯誤率為3%和50X測序深度的情況下，本文方法和文獻(xiàn)[20]中的方法幾乎相同，可以達(dá)到99%以上的平均準(zhǔn)確率。和文獻(xiàn)[22]相比，在錯誤率為10%的情況下，60X測序深度下本文方法的平均準(zhǔn)確率為94%以上，高于文獻(xiàn)[22]92%的平均準(zhǔn)確率。此外，在錯誤率為10%和20X測序深度下，本文方法的平均準(zhǔn)確率依然達(dá)到94%以上，遠(yuǎn)遠(yuǎn)高于文獻(xiàn)[22]50%的平均準(zhǔn)確率。和文獻(xiàn)[23]相比，在錯誤率為5%和20X測序深度的情況下，本文方法的平均準(zhǔn)確率可以達(dá)到97%以上，高于文獻(xiàn)[23]的平均準(zhǔn)確率。此外，本文方法在錯誤率15%和60X測序深度的情況下，平均準(zhǔn)確率可以達(dá)到90%以上，遠(yuǎn)高于文獻(xiàn)[23]64%的平均準(zhǔn)確率。

4 結(jié)束語

如何解決序列中的組合插入、刪除和替換錯誤，是DNA存儲信息可靠性的基礎(chǔ)。DNA存儲的解碼過程主要包括兩個方面：測序讀長的分組和基于組內(nèi)測序讀長的一致性序列的恢復(fù)。為了提高DNA存儲信息的恢復(fù)精度，本文主要在以上兩個方面進(jìn)行了如下的研究：(1)提出了“索引+CRC哈希+索引”的序列索引編碼方法，仿真結(jié)果表明該索引編碼的分組精度可以達(dá)到99%以上，并保證較高的召回率。(2)在文本字符前向糾錯編碼的基礎(chǔ)上，提出一種基于可識別DNA碼的桶分配糾錯算法。影響本文糾錯方法精度的因素主要有3個：一是可識別DNA碼的檢索與糾錯；二是可識別DNA碼桶的分配；三是基于可識別DNA碼權(quán)重分配的多數(shù)投票。仿真結(jié)果表明在0.10和0.05錯誤率條件下，平均解碼準(zhǔn)確率在20X測序深度時可達(dá)94%以上；在0.15錯誤率條件下，平均解碼準(zhǔn)確率在60X測序深度時可達(dá)90%以上。此外，在給定錯誤率的情況下，本文提出的解碼算法為線性時間復(fù)雜度O(N)。因此，適合于未來面向大數(shù)據(jù)的DNA存儲應(yīng)用。最后，如何解決可識別DNA碼的最優(yōu)分配，進(jìn)一步提高分配的準(zhǔn)確率將是未來研究的一個主要的方向。