一种高效的前向纠错码桶分配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码的最优分配，进一步提高分配的准确率将是未来研究的一个主要的方向。