HMDA：基于HMAC的DNA组装算法

HMDA： HMAC-Based DNA Assembly Algorithm

摘要

DNA组装是基因组学研究中至关重要的步骤。传统的基于第二代测序技术的基因组组装算法存在以下问题：1） 无法保证数据的完整性和准确性；2） 无法处理混合物种基因组序列组装问题；3） 消耗的内存空间大；4） 缺乏安全性保障。针对这些问题，提出了一种基于HMAC的DNA组装算法——HMDA。在DNA组装过程中，运用HMAC技术对组装信息进行编码，并利用密码子的性质，在基因组内嵌入经HMAC加密后的信息内容。当执行组装操作时，这些信息会被提取出来，作为筛选正确序列的关键依据。此外，HMDA将不同密钥与不同物种或者用户一一映射，实现了特定物种或用户的身份验证。实验结果表明，HMDA算法相对于其他组装算法有以下优势：1） 针对理论上能够恢复基因组的测序文库，能够得到准确且完整的组装结果；2） 能从混合物种基因组数据中组装出不同物种；3） 空间消耗至多为其他算法的33.4%；4） 安全性更强，授权用户和普通用户组装结果不同，且该算法能检测数据篡改并返回错误状态码。

Abstract

DNA assembly is a crucial step in genomic research. Traditional genome assembly algorithms based on second-generation sequencing technology have the following problems: 1) they cannot guarantee the integrity and accuracy of the data; 2) they cannot handle the problem of assembling mixed species genome sequences; 3) they consume a large amount of memory space; and 4) they lack security guarantees. A DNA assembly algorithm based on a Hash-based message authentication code(HMAC), HMDA, is proposed to address these issues. In the DNA assembly process, HMAC technology is used to encode assembly information, and the properties of codons are used to embed HMAC-encrypted information content into the genome. This information is extracted as a key basis for screening the correct sequence when performing the assembly operations. In addition, HMDA maps different keys to different species or users individually, enabling species-specific or user-specific authentication. Experimental results show that this algorithm has the following advantages over other assembly algorithms:1) it can generate accurate and complete assembly results for sequencing libraries that theoretically recover genomes;2) it can assemble different species from mixed-species genomic data;3) it utilizes at most 33.4% of the space resources required by other algorithms; 4) it enhances security by producing different assembly results for authorized and ordinary users. HMDA can detect data tampering and return error status codes.

Graphical abstract

关键词

DNA组装 / 基因组学 / 第二代测序技术 / HMAC

Key words

DNA assembly / genomics / second-generation sequencing technology / Hash-based message authentication code (HMAC)

0　 引 言

DNA组装是基因组学研究中至关重要的步骤，其目标是将碎片化的DNA测序片段重新组合成连续的基因组序列。DNA组装在许多领域具有广泛的应用，包括基因组结构研究、活体存储^[1-3]和物种鉴定等。如何准确、完整地组装出DNA序列一直是学界关心的问题。

高通量测序技术的迅猛发展为序列组装带来了巨大的潜力，然而，这一发展也伴随着前所未有的计算挑战^[4-6]。当前的高通量测序技术仅能检测较短的DNA片段，而基因组的实际长度远远超出了这些片段的长度。因此，在测序时需要将基因组随机打碎成小片段，然后进行测序，最后得到这些小片段构成的读段（reads）文库^[7]。序列组装的核心任务是通过reads文库中的重叠信息，将这些reads组装成原始的基因组序列。这种重叠信息表现为不同reads的末端与起始端之间存在一定长度的重叠。因此，如何有效利用高通量测序技术生成的reads文库还原原始序列，成为DNA组装问题的关键挑战。

在基因组序列组装领域，常用的基于第二代测序技术的基因组组装算法包括Velvet^[8-9]、PHRAP^[10-11]和SPAdes^[12-13]。这三种算法均适用于对待测的野生与人工DNA序列进行分析处理。但是，这些方法本质上都是将DNA组装问题转化为复杂的NP问题进行求解，由此得到的仅是近似最优解，无法保证组装结果的准确性和完整性。它们在设计上普遍忽视了待组装reads数据所属物种的来源，同时也未考量组装过程中的用户身份验证问题，从而潜藏了一定的安全隐患。另外，这些算法在运行过程中通常需要维护庞大的图结构数据，在剔除错误组装路径和过滤无效信息方面效能不高，故而容易造成大量空间资源的浪费。

综合来看，上述算法在实践中暴露出四个主要问题：1） 无法从包含整个基因组信息的文库中准确、完整地组装出原始基因组序列；2） 难以处理混合物种基因组组装问题；3） 高内存消耗；4） 缺乏安全性保障。因此，本文在第二代测序技术的背景下，提出了一种基于HMAC（Hash-based Message Authentication Code）的DNA组装算法，简称HMDA算法。通过在DNA组装问题中引入HMAC技术，增强组装的准确性、完整性、可扩展性、空间节约性和安全性。

1　 基于哈希的消息认证码（HMAC）

HMAC是一种基于哈希（Hash）函数和密钥进行消息认证的方法，可以与各种迭代散列函数结合使用，包括但不限于SHA-1、SHA-256、MD5等^[14-15]。该方法能够提供消息完整性验证和数据源认证，确保通信的安全性，防止数据在传输过程中遭受篡改或伪造。

HMAC机制的定义如（1）式：

\mathrm{H}\mathrm{M}\mathrm{A}{\mathrm{C}}_{\mathrm{k}\mathrm{e}\mathrm{y}}(M)=

H((\mathrm{k}\mathrm{e}{\mathrm{y}}^{+}\oplus \mathrm{o}\mathrm{p}\mathrm{a}\mathrm{d})||H((\mathrm{k}\mathrm{e}{\mathrm{y}}^{+}\oplus \mathrm{i}\mathrm{p}\mathrm{a}\mathrm{d})||M))（1）

其中，key表示共享密钥；M表示明文；H表示HMAC机制采用的哈希密码函数，假设其块长l为64字节；opad是和哈希密码函数块长等长的固定值，其每个字节都为0x5C；ipad是和哈希密码函数块长等长度的固定值，其每个字节都为0x36；“\oplus”表示异或操作；“||”表示字节连接操作。

HMAC机制执行过程如下：

1) 判断key长度是否为l的倍数，若不是，则在key末尾填充0直至其长度达到l的倍数，得到key⁺。

2) 对key⁺和ipad执行异或运算，得到一个64字节的数据块，并将该数据块和明文M连接，形成数据块(\mathrm{k}\mathrm{e}{\mathrm{y}}^{+}\oplus \mathrm{i}\mathrm{p}\mathrm{a}\mathrm{d})||M。

3) 将在2）中获得的数据块作为H的输入，输出一个64字节的数据块。

4) 将key⁺和opad进行异或运算，得到一个64字节数据块。

5) 将在4）中获得的数据块和在3）中获得的数据块进行连接，并作为H的输入，最终得到64字节的密文数据块。

HMAC生成消息认证码所需要的时间主要取决于对明文进行哈希运算的过程，因其简单高效并且安全性能优秀，得到了广泛的应用^[16-18]。将HMAC应用于DNA组装具有两大优势。一是HMAC实现了特定物种或用户的身份验证。只有拥有正确密钥的个体才能正确解密DNA序列中隐藏的信息，确保DNA所嵌入的数据的访问受到控制，可用于处理混合物种基因组组装问题和确保组装安全。二是HMAC具有出色的抗篡改性质，可检测DNA序列的篡改或未经授权的修改。在DNA片段组装时，HMAC验证组装过程的完整性，确保生成的DNA序列无损害或错误。这有助于过滤掉错误序列，减少无效生长，降低内存消耗。这两个优势保证了HMDA算法能够从包含整个基因组信息的文库中准确、完整地组装原始基因组序列。

2　 HMDA算法

将HMAC技术运用到DNA序列组装中，通过HMAC编码的方式，将组装的信息嵌入到事先选定的用作信息载体的基因中。在组装过程中，解码信息载体中嵌入的信息，以校验生长序列的正确性，去除无效的生长。HMDA算法包括编码阶段和解码阶段。编码阶段分为两个步骤：HMAC编码和信息嵌入。解码阶段包括三个步骤：确定信息基因、基于重叠信息的生长、剪枝过滤。

2.1　编码阶段

2.1.1　HMAC编码

在生物基因组上选取N_{\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}\mathrm{o}\mathrm{r}}个基因作为信息载体，基因可表示为A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}(A_j表示氨基酸，\mathrm{1}\le j\le N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i},N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}表示\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i中氨基酸的个数)。为了方便描述，将作为信息载体的基因命名为信息基因，信息基因的作用是嵌入真实序列片段信息。假定基因\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i为信息基因(\mathrm{1}\le i\le N_{\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}\mathrm{o}\mathrm{r}})，\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i与其要嵌入的真实序列片段的对应关系如图1所示，\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i需要嵌入其左侧真实序列片段\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}信息和右侧真实序列片段\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}信息。\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i中A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_{⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋}( ⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋表示\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}向下取整)嵌入\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}信息；\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i中A_{⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋+\mathrm{1}}{A}_{⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋+\mathrm{2}}\dots A_{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}则嵌入\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}信息。

以A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_{⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋}为例，密钥为key。如果i为偶数，则氨基酸A_i嵌入的\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i是\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}按照规则划分得到的若干序列片段中的某个片段\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{l}\mathrm{e}\mathrm{f}{\mathrm{t}}_i经HMAC编码后的数值，如（2）式：

\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i=\left(\mathrm{H}\mathrm{M}\mathrm{A}{\mathrm{C}}_{\mathrm{k}\mathrm{e}\mathrm{y}}\left(\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\_\mathrm{l}\mathrm{e}\mathrm{f}{\mathrm{t}}_i\right)\right)\mathrm{\%}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}（2）

其中，\mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}{\mathrm{e}}_i\subset \mathrm{S}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}{\mathrm{e}}_{i+\mathrm{1}}，\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}为一个常数，%表示取模运算。

如果i为奇数，且i>N_{\mathrm{A}\mathrm{A}}（N_{\mathrm{A}\mathrm{A}}表示存储的氨基酸信息的个数，即氨基酸A_i记忆左侧氨基酸序列的能力），则氨基酸A_i嵌入的\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i是A_i前N_{\mathrm{A}\mathrm{A}}个氨基酸构成的氨基酸序列\mathrm{A}{\mathrm{A}}_{N_{\mathrm{A}\mathrm{A}}}和A_i在整个基因组中的位置信息\mathrm{p}\mathrm{o}{\mathrm{s}}_i连接成的字符串经过HMAC编码的数值，如（3）式：

\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i=\left(\mathrm{H}\mathrm{M}\mathrm{A}{\mathrm{C}}_{\mathrm{k}\mathrm{e}\mathrm{y}}\left(\mathrm{A}{\mathrm{A}}_{N_{\mathrm{A}\mathrm{A}}}|\mathrm{p}\mathrm{o}{\mathrm{s}}_i\right)\right)\mathrm{\%}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}（3）

2.1.2　信息嵌入

在DNA中，每三个相邻的碱基构成一个密码子（codon），这些密码子排列紧密且互不重叠，每个密码子对应一个特定的氨基酸。在编辑DNA的碱基序列时，为了保持其生化特性和功能，需要确保生成的氨基酸序列不发生变化。此外，自然界中存在64种密码子，它们对应着20种氨基酸。同一种氨基酸通常与多个密码子对应，这导致信息存在冗余性。此外，同一种氨基酸对应的密码子出现频率不均等，即存在密码子偏好性。在利用氨基酸与其对应密码子之间的映射关系中存在的信息冗余特性来嵌入有助于组装的信息时，需要尽量不改变或最小程度的影响DNA原有的密码子偏好性^[19-21]。

信息嵌入的核心是将数字信息转换为氨基酸的密码子信息，通过“替换同义密码子而不改变氨基酸”的方式将HMAC编码后的值\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i嵌入到所选定的基因(\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i)中，在此过程中，需要确保嵌入信息后基因中各种密码子出现的比例在统计学上与自然界中的密码子比例保持一致。为了方便实际操作，设某个氨基酸A_i有K_i个对应的密码子，它们分别称为\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,\mathrm{0}},\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,\mathrm{1}},\dots ,\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,K_i-\mathrm{1}}，各密码子在自然界中出现的频率比如（4）式所示：

d_{i,\mathrm{0}}:d_{i,\mathrm{1}}:\dots :d_{i,K_i-\mathrm{1}}（4）

其中，d_{i,j}为整数()，且d_{i,j}\ge \mathrm{0}。频率之和L_i如（5）式所示。

L_i={\displaystyle \sum _{j=\mathrm{0}}^{K_i-\mathrm{1}}}{d}_{i,j}（5）

密码子信息嵌入是指将一个小于L_i的非负整数\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i通过选择A_i对应的密码子，嵌入A_i。在本次应用中，令L_i=\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}。\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i嵌入A_i时，如（6）至（8）式所示选用密码子c_i：

\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}(A_i)=\left[\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,\mathrm{0}},\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,\mathrm{1}},\dots ,\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{i,K_i-\mathrm{1}}\right]（6）

D(A_i)=\left[d_{i,\mathrm{0}},d_{i,\mathrm{1}},\dots ,d_{i,K_i-\mathrm{1}}\right]（7）

c_i=\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\left(\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i,\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}(A_i),D(A_i)\right)=

（8）

其中,\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}(A_i),D(A_i)分别表示氨基酸A_i对应的密码子列表及其统计频率。在\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_i足够均匀随机的情况下，这种映射方法可以使得嵌入后的序列中各密码子的出现频率逐渐接近自然状态下各密码子出现的频率，从而使得嵌入数据前后的氨基酸密码子序列在统计学上不可区分，确保HMAC编码后的值嵌入到氨基酸序列中不会产生新的信息序列。

HMDA编码算法描述如算法1。

2.2　解码阶段

2.2.1　确定信息基因

2.2.1.1 确定信息基因的位置

在解码阶段，解码方只能根据第二代测序技术得到的reads文库、与编码方约定的密钥和序列划分规则组装DNA序列。文库中的读段\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i，表示为b_{\mathrm{1}}{b}_{\mathrm{2}}\dots b_z(b_j\in \left\{\mathrm{A},\mathrm{G},\mathrm{C},\mathrm{T}\right\},\mathrm{1}\le j\le z)，根据编码规则，确定\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i中哪些氨基酸存储了位置信息。由于第一个氨基酸的第一个碱基可能是b_{\mathrm{1}}\mathrm{、}{b}_{\mathrm{2}}或者b_{\mathrm{3}}，且第一个氨基酸在编码阶段时下标奇偶性不确定，因此需要考虑这六种情况去解码reads里面存储的位置信息。本节以其中一种情况为例进行解码，以氨基酸为单位，\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i可以表示为A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_{⌊\frac{z}{\mathrm{3}}⌋}，并且A_{\mathrm{1}}在编码时是下标为偶数的氨基酸。

已知原基因组大小为N个碱基，则每个基因的第一个密码子有N个可能位置。针对\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i中的氨基酸A_j，遵循编码阶段时设定的位置处理规则，将N个可能的位置进行相应的HMAC编码和信息嵌入操作，并在此过程中判断A_j的密码子是否发生了变化，如（9）式和（10）式所示计算解码阶段经过HMAC编码后的值\mathrm{l}\mathrm{o}{\mathrm{c}}_j和解码时氨基酸A_j选择的密码子信息\mathrm{A}{\mathrm{A}}_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}。

\mathrm{l}\mathrm{o}{\mathrm{c}}_j=\left(\mathrm{H}\mathrm{M}\mathrm{A}{\mathrm{C}}_{\mathrm{k}\mathrm{e}\mathrm{y}}\left(\mathrm{A}{\mathrm{A}}_{N_{\mathrm{A}\mathrm{A}}}|\mathrm{p}\mathrm{o}{\mathrm{s}}_j\right)\right)\mathrm{\%}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}（9）

\mathrm{A}{\mathrm{A}}_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}=

\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\left(\mathrm{l}\mathrm{o}{\mathrm{c}}_j,\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}\left(A_j\right),D\left(A_j\right)\right)（10）

其中，\mathrm{p}\mathrm{o}{\mathrm{s}}_j\in \{\mathrm{1,2},\dots ,N\}，表示氨基酸A_j可能的位置信息。

本文设计了一套打分规则，旨在确定信息基因组嵌入的位置信息。\mathrm{在处}\mathrm{理}\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i时，创建一个大小为N+1、初始赋值为0的打分数组SA。随后，针对位置\mathrm{p}\mathrm{o}{\mathrm{s}}_j，如果\mathrm{A}{\mathrm{A}}_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}与A_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}是一致的(A_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}代表真实情况中氨基酸A_j选择的密码子信息)，则\mathrm{位置}\mathrm{p}\mathrm{o}{\mathrm{s}}_j-j+\mathrm{1}的分数加1。\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{S}\mathrm{A}\right)记录最高得分，如果\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{S}\mathrm{A}\right)明显高于数组中其他位置的分数，则认为其对应的下标是该\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i中A_{\mathrm{1}}的位置；否则，认为该\mathrm{r}\mathrm{e}\mathrm{a}{\mathrm{d}}_i没有嵌入位置信息。

2.2.1.2 确定信息基因的大小

经过2.2.1.1节操作以后，存储位置信息的reads都被定位到整个基因组的相应位置。为确定信息基因的大小，需找到可能的起始密码子和终止密码子，统称为边界标志。在边界标志前后，选择与reads大小相等的序列片段，然后执行打分操作，以判断是否为真正的边界标志。真正的边界标志前氨基酸嵌入了位置信息，打分是有规律的；真正的边界标志后氨基酸没有嵌入位置信息，打分是随机的。根据此特性，可以确定真正的边界标志，从而确定出基因的范围。

2.2.2　基于重叠信息的生长

利用布谷哈希(cuckcoo Hash)^[22-23]的思想，对reads文库中的所有reads建立哈希表。哈希表的键是每个read尾部overlap长度的序列构成的字符串，经过哈希映射后得到的数值，相应的值则是该read中剩余碱基序列构成的字符串。其中，overlap是一个数值参数，代表reads之间重叠序列长度。

在组装过程时，首先以信息基因尾部overlap长度序列哈希后的值作为初始键，从哈希表中检索匹配的reads进行组装。然后，将组装得到的所有分支序列尾部overlap长度的序列哈希后的值作为新的键，继续搜索匹配的reads并组装。持续迭代这个过程，直到查询结果为空，或者达到与相邻信息基因的重叠长度等于overlap时，组装过程终止。

2.2.3　剪枝过滤

在HMDA算法中，当序列生长到信息基因中某个氨基酸A_j所嵌入的序列信息的范围时，按照氨基酸中嵌入的真实序列信息，对组装得到的所有分支\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{s}执行剪枝过滤操作。

在处理某个分支\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}{\mathrm{h}}_i\in \mathrm{b}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{s}时，进行如（11）式和（12）式所示的操作：

\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_j=\left(\mathrm{H}\mathrm{M}\mathrm{A}{\mathrm{C}}_{\mathrm{k}\mathrm{e}\mathrm{y}}\left(\mathrm{b}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}{\mathrm{h}}_i\right)\right)\mathrm{\%}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}\mathrm{S}\mathrm{u}\mathrm{m}（11）

\mathrm{A}{\mathrm{A}}_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}=

\mathrm{C}\mathrm{h}\mathrm{o}\mathrm{o}\mathrm{s}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\left(\mathrm{D}\mathrm{a}\mathrm{t}{\mathrm{a}}_j,\mathrm{C}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}\mathrm{L}\mathrm{i}\mathrm{s}\mathrm{t}\left(A_j\right),D\left(A_j\right)\right)（12）

若\mathrm{A}{\mathrm{A}}_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}与A_j^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}\mathrm{n}}不一致，则该分支被视为错误序列并予以去除；反之，该分支会被保留并继续生长。需要注意的是，经过剪枝过滤操作后留存的分支尚不能断言为正确，但被丢弃的分支可以确定为错误。被保留的分支会继续生长，直到达到下一个氨基酸A_{j+\mathrm{1}}所嵌入的序列信息的区域，再次进行剪枝过滤操作。若在文库中无法找到与当前分支匹配的reads导致无法继续生长，那么这条分支也将被直接去除。

HMDA解码算法描述如算法2。

3　 算法分析

3.1　打分规则的性能分析

3.1.1　理论分析

读段read表示为A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_L(其中L=⌊\frac{\mathrm{l}\mathrm{e}\mathrm{n}\left(\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}\right)}{\mathrm{3}}⌋，\mathrm{l}\mathrm{e}\mathrm{n}(\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d})表示\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}的碱基个数。假设read中连续L_{\mathrm{1}}个氨基酸来自信息基因区域，连续L_{\mathrm{2}}个氨基酸来自非信息基因区域，L_{\mathrm{1}}+L_{\mathrm{2}}=L。本节将对下述四种情形进行分析，评估正确位置和非正确位置的得分情况。

情形1：若L_{\mathrm{1}}=\mathrm{0},打分阶段,正确位置得分\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}和其他位置的平均得分\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}都如（13）式：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{j=N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}\right)}{N}（13）

情形2：若L_{\mathrm{2}}=\mathrm{0}，打分阶段，正确位置得分\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}如（14）式：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=L-N_{AA}（14）

其他位置的平均得分\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}如（15）式所示：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{j=N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}-\mathrm{1}\right)}{N-\mathrm{1}}（15）

两者得分相差为（16）式：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}-\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{N}{N-\mathrm{1}}\times ({\displaystyle \sum _{j=N_{AA}+\mathrm{1}}^{j=L}}(\mathrm{1}-P_{A_{{}_j}}^{codo{n}_j}))（16）

情形3：若L_{\mathrm{1}}\ne \mathrm{0},L_{\mathrm{2}}\ne \mathrm{0}，且read\mathrm{中}{L}_{\mathrm{1}}个氨基酸在前，L_{\mathrm{2}}个氨基酸在后。打分阶段，\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}得分为（17）式所示：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=L_{\mathrm{1}}-N_{\mathrm{A}\mathrm{A}}+\frac{{\displaystyle \sum _{j=L_{\mathrm{1}}+\mathrm{1}}^{j=L_{\mathrm{1}}+L_{\mathrm{2}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}\right)}{N}（17）

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}得分如（18）式所示：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{j=N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L_{\mathrm{1}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}-\mathrm{1}\right)}{N-\mathrm{1}}+

\frac{{\displaystyle \sum _{j=L_{\mathrm{1}}+\mathrm{1}}^{j=L_{\mathrm{1}}+L_{\mathrm{2}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}\right)}{N}（18）

两者得分相差为（19）式：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}-\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{N}{N-\mathrm{1}}\times ({\displaystyle \sum _{j=N_{AA}+\mathrm{1}}^{j=L_{\mathrm{1}}}}(\mathrm{1}-P_{A_{{}_j}}^{codo{n}_j}))（19）

情形4：若L_{\mathrm{1}}\ne \mathrm{0},L_{\mathrm{2}}\ne \mathrm{0}，且read中L_{\mathrm{2}}个氨基酸在前，L_{\mathrm{1}}个氨基酸在后。打分阶段，\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}得分为（20）式所示：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}=L_{\mathrm{1}}-N_{\mathrm{A}\mathrm{A}}+\frac{{\displaystyle \sum _{j=N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L_{\mathrm{2}}+N_{\mathrm{A}\mathrm{A}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}\right)}{N}（20）

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}得分如（21）式所示：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{j=N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L_{\mathrm{2}}+N_{\mathrm{A}\mathrm{A}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}\right)}{N}+

\frac{{\displaystyle \sum _{j=L_{\mathrm{2}}+N_{\mathrm{A}\mathrm{A}}+\mathrm{1}}^{j=L_{\mathrm{1}}+L_{\mathrm{2}}}}\left(N\times P_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}-\mathrm{1}\right)}{N-\mathrm{1}}（21）

两者得分相差如（22）式：

\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}-\mathrm{P}\mathrm{o}{\mathrm{s}}_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{N}{N-\mathrm{1}}\times ({\displaystyle \sum _{j=L_{\mathrm{2}}+N_{AA}+\mathrm{1}}^{j=L_{\mathrm{1}}+L_{\mathrm{2}}}}(\mathrm{1}-P_{A_j}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}))（22）

其中，P_{A_j}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}表示氨基酸A_j的密码子为\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j的概率值。在情形1中，所有位置均未显示出明显的得分峰值。而在其他情形下，正确位置和非正确位置的得分差异确然存在，如（16）、（19）和（22）式所示，L_{\mathrm{1}}的大小将会决定这种得分差异的显著性，进而影响算法抗噪声干扰的能力。

3.1.2　实验分析

在实验中，信息基因潜在的起始位置共计4 521 562个可能，其正确的起始位置设定在20。实验选用长度为150个碱基对的reads进行研究，并设定N_{\mathrm{A}\mathrm{A}}为5。为了研究嵌入位置信息的氨基酸在整个read中的占比对最终结果的影响，我们设置了四种不同的L_{\mathrm{1}}占比级别：0%（没有任何氨基酸携带位置信息）、32%、68%和100%（所有氨基酸均携带位置信息）。

图2至图5揭示了在不同嵌入比例下，正确位置与非正确位置得分之间的显著差异。通过实验观察表明，随着L_{\mathrm{1}}占比增大，得分分布图呈现出明显的峰值特征，而该峰值位置能够准确指示信息基因的实际起始点。

3.2　剪枝策略的性能分析

同一氨基酸对应的密码子出现概率各不相同，这意味着即使分支序列组装错误，氨基酸A_j的密码子\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j也有P_{A_j}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}大小的概率选择保持不变。因此，错误的分支序列有约P_{A_j}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}的概率被保留。然而，信息基因中每个氨基酸嵌入的序列信息包含前面氨基酸嵌入的序列信息，每次误判留下的分支序列都需要通过后续氨基酸的又一次筛选操作，导致误判的概率呈指数级下降，从而提高了留下序列的准确性。

假设测序文库的平均测序深度为depth，测序过程中碱基替换错误发生的概率为P_{\mathrm{s}\mathrm{u}\mathrm{b}}，碱基增删错误的概率则为P_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{l}}。同时，考虑到基因组内部两个不同区域发生重复碱基段的概率为P_{\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}}。当一段生长的碱基数目为grow_num的序列时，涉及的信息基因被标记为\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i，其氨基酸序列为A_{\mathrm{1}}{A}_{\mathrm{2}}\dots A_{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}，对应的密码子序列则为\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{\mathrm{2}}\dots \mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}。最后，基于这些参数，通过组装所得到的基因序列数量，即contig_num，可通过（23）式得出上限。

\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{g}\_\mathrm{n}\mathrm{u}\mathrm{m}\le

(\mathrm{d}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{h}\times (P_{\mathrm{s}\mathrm{u}\mathrm{b}}+P_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{l}})+P_{\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}}{\times \mathrm{1}+\mathrm{1})}^{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\_\mathrm{n}\mathrm{u}\mathrm{m}}\times

{\displaystyle \prod _{j=\mathrm{1}}^{{}^{⌊\frac{⌊\frac{N_{\mathrm{G}\mathrm{e}\mathrm{n}{\mathrm{e}}_i}}{\mathrm{2}}⌋}{\mathrm{2}}⌋}}}{P}_{A_{{}_j}}^{\mathrm{c}\mathrm{o}\mathrm{d}\mathrm{o}{\mathrm{n}}_j}（23）

在实际的基因组组装过程中，P_{\mathrm{s}\mathrm{u}\mathrm{b}}通常处于{\mathrm{10}}^{-\mathrm{3}}的数量级，而P_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{l}}更低，大约为{\mathrm{10}}^{-\mathrm{6}}的数量级。尽管基因组内有一些区域存在重复序列，可能在组装阶段产生较多的候选碱基序列，但大多数情况下，P_{\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}}相对较小，理论上可近似表达为{(\frac{\mathrm{1}}{\mathrm{4}})}^{\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}}。对于生长碱基的数量，一般情况下可达40万碱基数量级，而单个基因对应的氨基酸序列平均长度约为500个氨基酸。鉴于以上参数，经过综合考量，在平均意义下，组装得到的基因序列数量contig_num可以有效地控制在个位数范围内，这将原本可能呈指数增长的复杂度问题大幅度简化，从而显著减少了对时间和空间资源的需求。

3.3　安全性分析

3.3.1　HMDA的错误检测性

已经编码的数据遭到恶意篡改后，不同于其他算法会返回错误的组装结果，也不会告知组装方结果有错误，HMDA算法会返回组装错误状态码，并直接宣告组装失败。这是因为HMAC具有防篡改性质，一旦数据被改动，其HMAC值也会相应改变，导致信息基因的位置定位出错以及生长序列的剪枝操作出错，进而造成组装失败。

3.3.2　HMDA的密钥相关性

由于HMAC是密钥相关的，只有拥有正确密钥的用户才能解密并利用信息基因中嵌入的组装信息去除错误的生长序列，最终得到个位数的组装结果，其中包含正确的结果。相反，普通用户则只能使用现有的组装算法进行组装，得到的是若干条无序的短碱基序列，这些序列尚未确定其在基因组中的确切位置。

4　 实验与分析

4.1　序列组装问题

4.1.1　实验设计

选择大肠杆菌基因组作为第一个研究对象，在NCBI官方网站获取其基因组信息，从中挑选了yhjK、ftsP、ythD、yeeO、ydCR、asnS、fdrA、yaaJ和thiC共9个基因作为信息基因。使用密钥“bacterium”对信息基因进行HMAC编码。编码后的结果被用作新的基因组，并通过wgsim^[24]工具生成了一个包含70 000 000条read的reads文库，其中包含了大肠杆菌完整的基因组信息。最后对reads文库进行数据清洗，去除一些不完整、字母种类不正确的数据。每个read的长度为150碱基，N_{\mathrm{A}\mathrm{A}}取值为6。

4.1.2　实验执行与结果

组装结果如表1所示。出发基因和终止基因是信息基因，总生长个数代表相邻信息基因间隔的碱基数量，剪枝次数代表过滤掉的错误分支的数量，结果序列个数代表在这一段生长过程中最后剩下的序列个数，是否是原序列代表是否组装得到原始序列。数据显示最终得到两条序列，其中一条是正确且完整的大肠杆菌基因组序列。表1中的剪枝次数表明，在组装过程中出现大量错误分支，即无效地生长，消耗大量资源。剪枝后的结果序列数量显示，剪枝操作成功去除了绝大多数无效的生长，仅保留了至多两条可能正确的生长序列，从而显著减少了资源的消耗，并提高了组装的准确性。

4.2　混合物种基因组组装

4.2.1　实验设计

选择枯草芽孢杆菌基因组作为第二个研究对象，在NCBI官方网站获取其详细的基因组信息，并从中挑选了ybeC、yeeF、carB、pksS、ilvD、levR、pbpD、hpxW和bglP共9个基因作为信息基因，使用密钥“subtilis”对信息基因进行HMAC编码。编码后的结果被用作新的基因组，并通过wgsim工具生成了一个包含70 000 000条read的reads文库，其中包含了整个枯草芽孢杆菌基因组信息。将4.1节生成的reads文库与本节对枯草芽孢杆菌基因组测序生成的reads文库合并成一个新的reads文库。最后对reads文库进行数据清洗，去除一些不完整、字母种类不正确的数据。每个read的长度为150碱基，N_{\mathrm{A}\mathrm{A}}取值为6。

4.2.2　实验执行与结果

在对4.2.1节生成的新reads文库执行解码过程时，为确定混合物种各自的信息基因，需要尝试密钥“bacterium”和“subtilis”。

根据表2和表3显示的组装结果，该算法在混合物种基因组下成功组装出大肠杆菌和枯草芽孢杆菌的基因组，实现了混合物种基因组组装的目标。

4.3　算法比较

为了更全面地分析HMDA的组装效果，本文基于4.1与4.2节中的数据，将HMDA算法与常用的基于第二代测序技术的基因组组装算法（Velvet、PHRAP和SPAdes）的组装效果进行了比较，四种算法的比较结果如表4所示。

通过比较发现，Velvet、PHRAP和SPAdes适用于人工序列和野生序列的组装，而HMDA仅仅适用于人工序列的组装，因为HMDA算法需要在序列中嵌入信息。但是相较于其他三种算法，本文提出的HMDA算法能够从包含整个基因组信息的文库中准确、完整地组装原始基因组，可以处理混合物种基因组组装问题，空间资源消耗至多为其他算法的33.4%，同时提供了安全性保障。因此，HMDA算法在某些应用场景下具有更优的性能。

5　 结 语

本文提出了一种能够保障数据完整性、处理混合物种基因组组装问题、降低空间资源消耗且提供安全性保障的DNA组装算法，是HMAC思想在DNA组装的创新性应用。与其他经典的基于第二代测序技术的基因组装算法相比，这一算法提供了更准确、完整、可扩展、空间节约和安全的DNA组装方法，有助于DNA组装的改进和发展。

在未来的工作中，本课题组计划提高算法的容错性能，允许保留与真实序列在可接受编辑距离内的组装结果，同时丢弃编辑距离超出可接受范围的组装序列。这一改进将有助于提高算法的鲁棒性和结果的准确性。

参考文献

[1] DONG Y M， SUN F J， PING Z， et al. DNA storage： Research landscape and future prospects［J］. National Science Review， 2020， 7（6）： 1092-1107. DOI： 10.1093/nsr/nwaa007 .

[2] RUAN C H， HAN R D， LI Y X， et al. Efficient DNA-based image coding and storage［C］//2023 IEEE International Symposium on Circuits and Systems （ISCAS）. New York： IEEE Press， 2023： 1-5. DOI： 10.1109/ISCAS46773.2023.10182162 .

[3] 许鹏， 方刚， 石晓龙， 等. DNA存储及其研究进展［J］. 电子与信息学报， 2020， 42（6）： 1326-1331. DOI： 10.11999/JEIT190-863 .

[4] XU P， FANG G， SHI X L， et al. DNA storage and its research progress［J］. Journal of Electronics & Information Technology， 2020， 42（6）： 1326-1331. DOI： 10.11999/JEIT190863（Ch ）.

[5] GIANI A M， GALLO G R， GIANFRANCESCHI L， et al. Long walk to genomics： History and current approaches to genome sequencing and assembly［J］. Computational and Structural Biotechnology Journal， 2019， 18： 9-19. DOI： 10.1016/j.csbj.2019.11.002 .

[6] HU T S， CHITNIS N， MONOS D， et al. Next-generation sequencing technologies： An overview［J］. Human Immunology， 2021， 82（11）： 801-811. DOI： 10.1016/j.humimm.2021.02.012 .

[7] KUMAR K R， COWLEY M J， DAVIS R L. Next-generation sequencing and emerging technologies［J］. Seminars in Thrombosis and Hemostasis， 2019， 45（7）： 661-673. DOI： 10.1055/s-0039-1688446 .

[8] 崔竞松， 薛慧， 王兰兰， 等. LEDA：一种基于Levenshtein距离的DNA序列拼接算法［J］. 武汉大学学报（理学版）， 2022， 68（3）： 271-278. DOI： 10.14188/j.1671-8836.2021.0079 .

[9] CUI J S， XUE H， WANG L L， et al. LEDA： A DNA sequence assembly algorithm based on levenshtein distance［J］. Journal of Wuhan University （Natural Science Edition）， 2022， 68（3）： 271-278. DOI： 10.14188/j.1671-8836.2021.0079（Ch ）.

[10] ZERBINO D R， BIRNEY E. Velvet： Algorithms for de novo short read assembly using de Bruijn graphs［J］. Genome Research， 2008， 18（5）： 821-829. DOI： 10.1101/gr.074492.107 .

[11] SHI H H， WU G. Gene sequence assembly algorithm model based on the DBG strategy and its application［J］. Journal of Healthcare Engineering， 2021， 2021： 6676194. DOI： 10.1155/2021/6676194 .

[12] DE LA BASTIDE M， MCCOMBIE W R. Assembling genomic DNA sequences with PHRAP［J］. Current Protocols in Bioinformatics， 2007， 11： 11.4.1-11.4.15. DOI： 10.1002/0471250953.bi1104s17 .

[13] RANADEV P， K.H KURMINDLA. Genome assembly：The art of creating a big thing from millions of small things［J］. Microbiology，2021，（05）：39.

[14] BANKEVICH A， NURK S， ANTIPOV D， et al. SPAdes： A new genome assembly algorithm and its applications to single-cell sequencing［J］. Journal of Computational Biology： A Journal of Computational Molecular Cell Biology， 2012， 19（5）： 455-477. DOI： 10.1089/cmb.2012.0021 .

[15] GUPTA S K， RAZA S， UNNO T. Comparison of de-novo assembly tools for plasmid metagenome analysis［J］. Genes & Genomics， 2019， 41（9）： 1077-1083. DOI： 10.1007/s13258-019-00839-1 .

[16] 谢俊. HMAC-SM3的差分功耗分析方法研究［D］. 上海： 上海交通大学， 2016. DOI： 10.27307/d.cnki.gsjtu.2016.001955 .

[17] XIE J. Research on differential power analysis method of HMAC-SM3［D］.Shanghai： Shanghai Jiao Tong University， 2016. DOI： 10.27307/d.cnki.gsjtu.2016.001955（Ch ）.

[18] 张诚， 范恒英， 刘维嘉. HMAC在IPSec和SSL中的应用［J］. 中国新通信， 2008， 10（11）： 22-25. DOI： 10.3969/j.issn.1673-4866.2008.11.005 .

[19] ZHANG C， FAN H Y， LIU W J. Application of HMAC in the IPSec and SSL［J］. China New Telecommunications， 2008， 10（11）： 22-25. DOI： 10.3969/j.issn.1673-4866.2008.11.005（Ch ）.

[20] 李丹枫， 王飞， 赵国鸿. 一种大流量报文HMAC-SM3认证实时加速引擎［J］. 计算机工程与科学， 2021， 43（1）： 82-88. DOI： 10.3969/j.issn.1007-130X.2021.01.010 .

[21] LI D F， WANG F， ZHAO G H. A real-time HMAC-SM3 acceleration engine for large network traffic［J］. Computer Engineering & Science， 2021， 43（1）： 82-88. DOI： 10.3969/j.issn.1007-130X.2021.01.010（Ch ）.

[22] SRINIVASAN S， SHIVAKUMAR K B， MUAZZAM M. HMAC-RSA： A security mechanism in cognitive radio for enhancing the security in a radio cognitive system［J］. Journal of Intelligent & Fuzzy Systems， 2019， 36（5）： 4449-4459. DOI： 10.3233/jifs-169999 .

[23] 胡乔木. 基于压缩HMAC算法的传感网隐私保护范围查询协议研究与设计［D］. 桂林： 桂林理工大学， 2022. DOI： 10.27050/d.cnki.gglgc.2022.000405 .

[24] HU Q M. Research and design of privacy protection range query protocol for sensor networks based on compressed HMAC algorithm［D］.Guilin： Guilin University of Technology， 2022. DOI： 10.27050/d.cnki.gglgc.2022.000405（Ch ）.

[25] ZHIRNOV V， ZADEGAN R M， SANDHU G S， et al. Nucleic acid memory［J］. Nature Materials， 2016， 15： 366-370. DOI： 10.1038/nmat4594 .

[26] BONNET J， COLOTTE M， COUDY D， et al. Chain and conformation stability of solid-state DNA： Implications for room temperature storage［J］. Nucleic Acids Research， 2010， 38（5）： 1531-1546. DOI： 10.1093/nar/gkp1060 .

[27] GIBSON D G. Synthesis of DNA fragments in yeast by one-step assembly of overlapping oligonucleotides［J］. Nucleic Acids Research， 2009， 37（20）： 6984-6990. DOI： 10.1093/nar/gkp687 .

[28] PAGH R， RODLER F F. Cuckoo hashing［J］. Journal of Algorithms， 2004， 51（2）： 122-144. DOI： 10.1016/j.jalgor.2003.12.002 .

[29] DEVROYE L， MORIN P. Cuckoo hashing： Further analysis［J］. Information Processing Letters， 2003，86（4）： 215-219. DOI： 10.1016/S0020-0190（02）00500-8 .

[30] HENG L. Wgsim［EB/OL］.［2011-10-18］.

基金资助

国家重点研发计划(2022YFB3903801)

湖北省重大科技专项项目(2022AAA009)