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

昝鄉鎮 姚翔宇 許 鵬 陳智華 石曉龍 李樹棟 劉文斌*

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

②(廣州大學網絡空間技術先進研究院 廣州 510006)

1 引言

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

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

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

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

2 方法

2.1 編碼方法

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

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

2.2 解碼方法

2.2.1 分組策略

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

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

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

2.2.2 桶式糾錯策略

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

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

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

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

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

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

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

(3)如果不存在可識別DNA碼，則從當前位置前進一個堿基，重復(1)。

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

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

3 實驗結果

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

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碼，最終導致投票失敗。因此，“索引+CRC+索引”分組策略既保障了精度又適當提升了召回率，為后面桶式糾錯奠定了關鍵的基礎。

3.2 糾錯策略性能分析

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

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

3.3 與其他方法的比較

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

表2 與其他方法的比較

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

4 結束語

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