摘要
为了适配非天然核酸（dP、dZ、dS、dB）和天然核酸（A、T、C、G）的存储，进一步提高DNA存储密度，本文提出了一种基于非天然核酸的高密度DNA存储编码方法。该方法首先设计了控制GC含量和均聚物长度的动态轮转映射算法，然后根据动态轮转映射算法，提出了适配八核酸的七进制哈夫曼算法，提升了DNA存储密度，最后利用七进制汉明码实现了数据纠错。对文本、图片、音频不同格式的文件利用该方法编码后，该方法平均编码密度可达3.38 bits/nt，各碱基对的含量可稳定保持在23%~26.5%，且不存在单碱基重复的情况。仿真实验结果表明，所提出的方法提高了DNA存储密度，约束了GC含量和均聚物长度且具有良好的纠错能力。
Abstract
In order to adapt the storage of non-natural nucleic acids (dP, dZ, dS, dB) and natural nucleic acids (A, T, C, G), and to further improve the DNA storage density, this paper proposes a high-density DNA storage coding method based on non-natural nucleic acids. The method firstly designs a dynamic rotation mapping algorithm to control the GC content and homopolymer length, then based on the dynamic rotation mapping algorithm, a heptadecimal Huffman algorithm adapted to octanucleic acid is proposed to enhance the DNA storage density, and finally, the data error correction is realized by utilizing the heptadecimal Hamming code. After utilizing this method to encode files of different formats of text, audio and pictures, the average coding density of this method can reach 3.38 bits/nt, and the content of each base pair can be maintained stably between 23% and 26.5%, and there is no single base repetition. Simulation experiments show that the proposed method improves the DNA storage density, constrains the GC content and homopolymer length and has good error correction capability.
Keywords
随着科学技术和信息技术的飞速发展，全球的数据信息总量正在以惊人的速度增长。据估计，到2040年全球信息总量将达到3×1024 bits[1]。目前，主流存储介质（磁带、硬盘、闪存等）已无法满足数据的存储需求[2]，导致一些重要的数据因无法存储而被丢弃[3]。因此，寻找可以解决未来数据存储的介质迫在眉睫。由于DNA具有容量大、密度高和寿命长等优点[4-7]，DNA被认为有望成为解决未来数据存储危机的存储介质[8]。
DNA是生物体携带遗传信息的物质，由腺嘌呤（A）、鸟嘌呤（G）、胸腺嘧啶（T）和胞嘧啶（C）四种碱基构成[9-10]，而DNA存储是以DNA等生物分子作为存储介质的一种新型技术[9，11]。从2012年Church团队成功在DNA存储了650 KB数据后[12]，DNA存储技术迎来了飞速发展。2019年，哥伦比亚大学联合苏黎世联邦工学院利用DNA存储3D打印出了斯坦福兔子[13]，2021年北京大学利用DNA存储构建了一种DNA“纳米弹弓”[14]，同时，DNA存储作为一项数据存储和读取的技术，随机读取也是近年来研究的热点[4，6，15]。
DNA存储数据的质量主要受DNA序列长度、DNA的生化性质和存储过程中的生化操作（DNA 分子的合成、保存和测序等）等影响[16]。DNA的生化性质（GC含量、均聚物长度）会影响生化操作，碱基序列中GC含量异常（高于60%或低于40%）和均聚物会增加DNA合成和测序过程中的错误率[17-19]。而DNA编码算法除了将数字信息映射到 DNA 碱基序列（编码）[3]，还要保证DNA序列中的GC含量、均聚物长度[20-21]。同时由于DNA分子合成、保存和测序过程中，可能会发生碱基删除、替换、插入的错误[18]，DNA的编码算法还要引入纠错方案纠正DNA序列中发生的错误[22-24]。DNA解码是将DNA无错误的恢复成数字信息，因此DNA编解码是实现高效DNA存储的关键步骤。
DNA编解码的研究主要聚焦于如何提升存储密度、控制生化约束条件和保证DNA存储数据的准确性。在提升存储密度方面，部分研究采用了压缩编码算法来减少数据存储占用的空间，提高存储密度[25-26]。研究者以四种天然核苷酸按比例混合成复合DNA字母，增加编码碱基数来提高存储密度[27-28]。在生化约束方面，Church算法[12]、Goldman算法[29]、Grass算法[30]和Blawat算法[31]通过设计映射规则控制均聚物长度和GC含量，但是存在着GC含量低于40%或高于60%的情况。DNA-Fountain算法[32]和 Yin-Yang算法[33]通过筛选将均聚物长度和GC含量控制在指定范围内，消除了GC含量低于40%或高于60%的情况。在保证DNA存储数据的准确性方面，Grass算法、DNA-Fountain算法和Yin-Yang算法均使用里德所罗门（Reed-solomon）纠错码进行纠错，Goldman算法采用物理冗余，分段重复保存来进行纠错。
DNA天然核苷酸的个数限制了每个核苷酸最大编码效率为2比特，为了提高DNA的存储密度，科学家们开始探索将非天然核苷酸引入到DNA存储序列中的方法。非天然核苷酸是人工合成的核苷酸，它们具有和活细胞中天然核苷酸相同的特性[34]，非天然核苷酸不仅可以和天然核苷酸一样进行DNA扩增过程的聚合酶链反应，还具有参与转录的能力[35]。Okamoto等[36]证明了非天然碱基Ds-Px具有稳定性和特异性，在聚合酶链反应扩增中具有高保真度和高效率性，Saito-Tarashima等[37]提出的非天然碱基对Im-Na除了具有很高的保真度和选择性，还拥有形成四个分子间氢键的能力，可以防止核苷酸碱基对错配的情况。Hoshika & Dien等[38-39]发现了四种非天然核苷酸（dP、dZ、dS、dB），dP-dZ和dS-dB非天然碱基对不仅不会破坏天然核苷酸碱基对的结构和稳定性，还具有和天然碱基对相同的行为。
由于非天然核苷酸的映射规则、GC含量范围和纠错码不同于天然核苷酸，天然核苷酸的编码算法无法适配于非天然核苷酸和天然核苷酸DNA存储的情况。因此面对非天然核苷酸的引入，Biswas等[40]于2020年提出了基于天然和非天然遗传核苷酸的编码方案，将二进制数据流成功映射到天然和非天然遗传核苷酸，理论存储密度为3 bits/nt，但是该方案没有考虑碱基序列中均聚物长度、GC含量和碱基发生错误后数据的恢复性。在2022年Biswas等[41]提出了密码约束编码算法，该算法约束了均聚物长度和GC含量，但没有添加纠错码来保证数据的准确性，理论存储密度可达到2.99 bits/nt。同年清华大学刘凯等[16]开发了适用于4、6和8个核苷酸的高可靠性编码系统，该系统控制了均聚物长度和GC含量，同时引入了纠错算法，可以纠正DNA序列丢失、替换错误，在8个碱基的编码中存储密度为3.64 bits/nt。
针对非天然核苷酸编码算法的存储密度还有提高空间，因此本研究提出了一种基于非天然核苷酸的高密度约束性DNA存储编码方法。该方法首先设计出可以控制GC含量和均聚物长度的动态轮转映射算法，然后根据映射算法改进了二进制哈夫曼编码，采用七进制哈夫曼编码对文件进行压缩，提升了存储密度，最后利用七进制汉明码保证了数据的正确性。
1 方法设计及原理
基于非天然核苷酸的高密度DNA存储编码方法是基于天然核苷酸（A、T、C、G）和非天然核苷酸（dP、dZ、dS、dB）的DNA存储方法，具体实现流程如图1所示。首先将需要编码的文件进行读取，接着采用哈夫曼编码对读取的字符进行编码压缩，编码得到七进制数据流。然后将数据流进行分段，为了便于数据流内部排序，为每段数据流添加地址，对加入地址的数据流添加七进制汉明纠错码，以保证数据的正确性。最后对数据流动态轮转映射进行编码，得到满足生化约束的DNA序列。DNA解码是编码的逆过程，将测序后的DNA序列恢复成数字信息。将DNA序列通过动态轮转映射解码获得七进制数据流，然后将七进制数据流去除七进制汉明纠错码和地址后，通过相应的哈夫曼编码进行解码恢复文件内容。
图1基于非天然核苷酸高密度约束性编码方法流程图
Fig.1Flowchart of high-density constrained coding method based on non natural nucleotides
1.1 动态轮转映射算法
在DNA存储编码中，映射算法是将数据流转变为碱基序列的关键步骤，它不仅需要适配文件压缩后数据流，还要控制碱基序列中GC含量和均聚物长度，因此本文第一步设计映射算法。目前基于天然核苷酸（A、T、C、G）编码方案中的碱基序列的GC含量在50%左右，均聚物的长度不大于3[24]。而基于非天然核苷酸和天然核苷酸的编码方案对碱基序列GC含量要求不同于天然核苷酸序列，2022年Biswas等[41]提出的密码约束编码算法中，每条碱基序列中单个碱基的含量控制在10.7%~14.2%，均聚物长度不大于7；同年清华大学刘凯团队[16]设计的编码系统中，每条碱基序列中单个碱基的含量约束在11.1%~13.6%，单碱基重复不超过2。在本文中，采用P、Z、S、B和A、T、C、G进行编码（P、Z、S、B为dP、dZ、dS、dB），因为碱基对A-T、C-G、P-Z、S-B是互补碱基对，所以只需控制每条碱基序列中A-T、C-G、P-Z、S-B的含量在25%左右。
Goldman编码[29]采用三进制数据流和前一位碱基映射编码完全消除了单碱基重复的情况，但固定的映射规则无法实现对GC含量调控。在Goldman编码基础上，本文提出了动态轮转映射算法，该算法采用七进制数据流和前三位碱基进行编码，消除了均聚物；同时设计了首编码碱基组轮换和筛选操作，通过动态改变首编码碱基组进行轮转编码，增加编码成功的可能性，并在编码后增加筛选操作，达到对GC含量调控的目的。
动态轮转映射算法的核心思想是建立一个多列映射编码表，增加数据流编码的多样性，通过筛选操作控制GC含量。为了增加GC含量调控成功的可能性，需要尽可能多的编码序列。如果首编码碱基为一位的话，动态改变首编码碱基，编码表序列只有八种可能，如果首编码碱基组为两位，编码选择有56种，在调控GC含量上都存在局限性。因此本映射算法采用三位的首编码碱基组，同时保证不存在均聚物，编码表中的序列有392种可能性。首编码碱基组是在“ATCGPZSB”八个核苷酸上进行的扩展，就“A”而言，在保证不存在单碱基重复的情况，以“A”为最后的一位的碱基组有49个，依次为TCA、TGA、TPA、TZA、TSA、TBA、CTA、CGA、CPA、CZA、CSA、CBA、GTA、GCA、GPA、GZA、GSA、GBA、PTA、PCA、PGA、PZA、PSA、PBA、ZTA、ZCA、ZGA、ZPA、ZSA、ZBA、STA、SCA、SGA、SPA、SZA、SBA、BTA、BCA、BGA、BPA、BZA、BSA、ATA、ACA、AGA、APA、AZA、ASA、ABA。
在本研究中，编码表的每一列由八个碱基字母“ACDITGPN”构成，其中通过去除首编码碱基组中最后一个碱基并对剩余碱基进行特定的顺序调整（即错序操作），最终生成七进制数字0~6的编码序列。具体而言，去除首编码碱基组最后一个碱基后，编码表每列中的碱基序列发生变化，并且这些序列不重复，从而确保每列碱基序列在编码过程中具有唯一性。如图2所示，编码过程中，当前的数字信息和前三位碱基的组合共同决定了所需的编码碱基，每个碱基与相应的七进制数值（0~6）一一对应。
动态轮转映射算法实现过程如算法1所示，获取待编码的数据列表，依次编码numer_list中每段数据。在编码每段数据时，依次选择首编码碱基组，循环编码each numer列表每个元素，计算编码后每段DNA序列中A-T、C-G、P-Z、S-B的含量，直到各类碱基含量符合条件，中断循环首编码碱基组，开始编码numer_list中下一段数据。在该算法中，最坏情况下，时间复杂度为：O（n*l*k），其中：n是输入的编码列表长度，l是首编码碱基列表的长度（392），k是each code列表中数据的个数（186），所以时间复杂度为O（n）。在空间复杂度方面，由于只会生成一个符合条件的DNA序列，因此它的空间复杂度为O（1）。
图2动态轮转编码映射图
Fig.2Dynamic rotation mapping table
1.2 七进制哈夫曼压缩
哈夫曼编码（霍夫曼编码）是一种无损压缩的熵编码算法，和其他熵编码算法一样，该算法在编码过程中利用符号之间出现的概率进行编码[42]。在哈夫曼编码中，静态哈夫曼编码在处理固定数据集或已知统计特征的数据时，能够提供很好的压缩率。但静态哈夫曼编码在压缩过程中需要对数据进行两次扫描，无法对动态数据进行压缩编码，同时静态哈夫曼编码需要保存额外的信息（比如字符的频率表和构建哈夫曼树的关键节点等）用于解码，使得整体的压缩效率降低[43]。自适应哈夫曼编码是基于静态哈夫曼编码的动态编码技术，不需要读取整体数据，可以随着读取的字符动态调整字符编码长度和改变哈夫曼树，可以压缩动态数据[44]。每当读入一个待编码字符后，都需要更新哈夫曼树中节点的权重值，以保证哈夫曼树的兄弟性质，因此编码过程消耗的时间比较长。
针对这两种压缩算法的特点，本文为了适配动态轮转映射算法的编码规则，提出了两种压缩算法，分别是七进制静态哈夫曼编码算法和七进制自适应哈夫曼编码算法。
读取文件是哈夫曼压缩前的关键一步，读取文件的方式直接影响哈夫曼的压缩率和所需时间。在读取文本文件过程中，以Unicode编码整数的形式读取文字，会使得编码字符数量增多，导致压缩率较高，而直接读取字符会使编码字符种类多于常规读取的0-255整数，当字符种类越多时，哈夫曼树的叶子节点越多，层数越大，在更新和遍历哈夫曼树时，需要消耗的时间越长。本文在读取文件时，对文本直接读取字符，称为直接读取；对其他文件读取UTF-8编码的整数，称为常规读取。
七进制静态哈夫曼编码与二进制哈夫曼树的原理一样，首先统计每个字符或者整数出现的频率，并按频率大小将字符（或整数）从小到大排序。然后利用频率构建哈夫曼树，为了保证哈夫曼树为满七叉树，在节点不够的情况下用“0”补上，且作为频率为“0”的节点。最后根据哈夫曼树将数据编码成七进制信息数据流。为了帮助顺利解码，需要在信息数据流前保存频率表（即编码字符、整数和出现的次数）、信息表元素个数等信息。
七进制自适应哈夫曼编码分为初始化树、对输入每个符号进行编码，更新树三部分。初始化树：哈夫曼树仅有一个NYT（Not yet transmitted）节点，表示树的初始状态。对输入每个符号进行编码：当编码端接收到一个待编码字符后，判断该字符是否首次出现，如果是首次出现，查询现有Huffman树是否存在空节点，将该字符保存到编号最大的空节点中；没有空节点，将NYT节点分出七个子节点，最左侧节点为新的NYT节点，最右侧节点保存字符，编码当前状态下字符的位置。更新树：在编码一个字符后，更新一次哈夫曼树状态，保证哈夫曼树中的节点始终满足兄弟性质（权重值越大的节点，对应的节点编号越大；子节点的节点编号始终小于父节点的节点编号）。
在更新哈夫曼树时，从读取字符对应的节点或读取字符将要存储的节点开始，一路向上更新到根节点，如图3所示，具体步骤为：
图3自适应哈夫曼树更新过程
Fig.3Adaptive huffman tree update process
Step1：判断当前节点是否为空，如果当前节点为空，算法结束；否则执行Step2。
Step2：把上一步访问的节点开始作为当前节点，寻找具有权值N的最远的节点（假设当前节点的权值为N）。
Step3：判断具有权值N的最远的节点是否为当前节点或者当前节点的父节点，如果是，执行Step4，如果不是则进行交换。
Step4：当前节点的权值加一，变为具有权值N+1的节点。
Step5：追溯当前节点的父亲节点，执行Step1。
为了缩短哈夫曼树更新所用的时间，提高该算法性能，在执行更新操作时，将更新的节点进行分类优化算法。具体步骤为：判断该节点存储的字符是否是第一次出现，如果是第一次出现，当前节点跳过查找相同权值节点和判断交换操作，直接将该节点的权值加一，然后追溯到当前节点的父亲节点进行更新操作（图4）。
从图4可以看出，采用七进制自适应哈夫曼编码对文件进行压缩编码的步骤为：
Step1：初始化树，生成NYT节点。
Step2：读取字符。
Step3：判断该字符是否首次出现，如果是首次出现，执行Step4，否则执行Step7。
Step4：查询现有Huffman树是否存在空节点。
Step5：如果是，则在所有空节点中选取节点编号最大的节点，将字符保存该节点。
Step6：如果否，将NYT节点分出七个子节点，最左侧节点为新的NYT节点，最右侧节点插入字符。
Step7：编码（编码当前状态下字符的位置）。
Step8：更新Huffman树。
Step9：判断是否存在未编码数据，如果是执行Step2，否则结束。
1.3 纠错码设计
在DNA存储过程中，数据的读写离不开DNA合成、PCR复制和DNA测序等生化过程。这些过程中均可能产生碱基插入、删除、替换等问题[45-47]。DNA编码算法采用物理冗余或逻辑冗余进行数据恢复，保证信息的正确性[5，9]。在逻辑冗余中，经常通过引入RS码[48-49]、喷泉码[25，32]、LDPC码[50]等进行纠错。本文在汉明码的基础上，提出了一种七进制汉明码的纠错机制。七进制汉明码相较于二进制汉明码，可以用于上述的映射方法，保证信息存储的正确性。
汉明码是一类用于纠错的线性分组码，通过添加冗余监督码原来检测和纠正数据传输中的错误[51]。七进制汉明码的信息码元和监督码元满足汉明码的关系，如公式（1）：
(1)
（n，k）汉明码表示在k个信息码元的基础上插入（n-k）个冗余码元，（n-k）为监督码元的长度。
图4自适应哈夫曼编码过程
Fig.4Adaptive huffman encoding process
以（7，4）七进制汉明码为例，按照添加监督位、添加纠错码、检错、纠错四步说明纠错过程。
1）添加监督位（校验位）
在（7，4）七进制汉明码中，码元总长度为7，信息码元长度为4，则监督码元长度为3。监督码元的监督位分别位于第1位、第2位和第4位。每个监督码元所负责监督的传输码元是根据其在汉明码中的位置来确定的。具体而言，监督码元（pi）（从1开始计数）负责监督第i位为1的位置上的传输码元（包括监督码元本身）。例如，第一个监督码元（p1）监督的是第一位为1的传输码元和其本身。以数据“5 236”为例，p1、p2、p3表示为监督码元，a1、a2、a3、a4表示为信息码元，位置如表1所示。p1负责校验p1、a1、a2、a4；p2负责校验p2、a1、a3、a4；p3负责校验p3、a2、a3、a4。
表1信息码元和监督码元的位置表
Table1Location table of information symbols and supervision symbols
2）添加纠错码（即计算监督码元）
由监督码元的校验关系可写出监督方程：
(2)
根据监督方程可以获得矩阵Q，矩阵Q可以计算监督码元和构建生成矩阵G。
(3)
(4)
二进制汉明码生成的监督码元和传输码元都是二进制数字0和1，而七进制汉明码为保证传输码元为七进制数字，规定生成监督码元和传输码元必须经过“模7”和“7减”两步运算（在“取模7”后，如果某位得数为0，则该位不进行“7相减”）。三个监督码元的计算如下：
(5)
(6)
同理，生成矩阵G和信息码元相乘，经过“模7”和“7减”，获得最终的传输码元。经过计算，“5236”插入监督码元，最终生成的传输码元为“1053236”。
3）检错
由监督方程（3）得到监督矩阵H，监督矩阵H用于检验传输码元是否发生错误。
(7)
在汉明码检错中，监督矩阵H与读取码元相乘计算出校正子，S1，S2，S3表示三个校正子，如果监督矩阵H与读取码元相乘得到S1S2S3，经过“模7”运算后结果为零向量，则表示各个传输码元没有错误。如下式所示：
(8)
如果S1S2S3经过“模7”运算后结果不是零向量，则传输码元出现了错误。假如传输码元为“1053236”，第三位发生错误后为“1033236”，则S1S2S3计算如下式所示：
(9)
4）纠错
通过校正子组合（S1S2S3）计算出错误发生的位置，错误位置与S1S2S3关系如表2所示。确定错误位置后，用误码和校正子非零数（X）相减计算正确的码元，相减计算时允许借位。即通过“550”可知a1发生错误，经过计算可得a1正确的码元为5。
表2错误码元位置鉴别表
Table2Error bit identification table
2 实验结果及分析
为了验证本文基于非天然核苷酸的高密度约束性DNA存储编码方法的有效性，以文本、图片和音频等不同格式的文件为测试对象进行了仿真实验。在仿真实验中，选用处理器为Intel（R）Core（TM）i5-6300HQ CPU@2.30GHz、RAM内存8.0 GB、64位windows10计算机，语言环境为python3.9.7，其中测试对象如表3所示。
表3测试文件细节表
Table3Test file details table
2.1 动态轮转映射算法效果分析
在分析该算法GC含量调控成功率时，首先计算出单次编码后CG、AT、PZ、SB含量调控的成功率，进而验证动态轮转映射算法的效果。CG、AT、PZ、SB在动态轮转映射算法中具有以下规律：以GC为例说明，当确定当前数字时，C、G在最后碱基为A、T、P、I、D、N列中的出现概率为2/7，非C、G出现概率为5/7；C、G在最后碱基为C、G的列中出现概率为1/7，非C、G出现的概率为6/7。同时在编码过程中，下一时刻出现是C、G的概率依赖于前一个碱基。比如前一个碱基为A、T、P、I、D、N之一时，则出现C、G的概率为2/7，非C、G的概率为5/7，前一个碱基为G或C时，则出现C或G的概率为1/7，非C、G的概率为6/7，这和马尔可夫链的定义不谋而合。因此可以将GC和非GC两个状态的变换通过有向图（图5）表达出来，
图5GC和非GC状态的有向图
Fig.5Directed graphs of GC and non GC states
并写出概率转移矩阵
(10)
根据首编码碱基可知，初始状态非GC和GC的概率分别为1，0或0，1，记为初始状态量X0=[1，0]或X0=[0，1]。
通过公式Xn=X0Pn计算，当X0=[1，0]，X1=[0.714，0.286]，可以计算出X9、X10···X100=[0.74978128，0.25021872]。当X0=[0，1]，X1=[0.857，0.143]，X2=[0.734449，0.265551]可以计算出X10、X11···X100=[0.74978128，0.25021872]，可知非CG碱基出现的概率收敛于0.74978128，GC碱基出现的概率收敛于0.25021872。
控制碱基序列中AT、CG、PZ、SB的含量在23%~26.5%时，当碱基序列为200 nt，需要保证GC的碱基个数为46~53。我们设定非CG概率为0.75，CG概率为0.25，通过二项分布累积函数计算，
(11)
通过计算可知，在利用一个首编码碱基组编码长度为200 nt的碱基序列中，GC含量控制在23.5%~26.5%的功率为48.579%。当AT、PZ、SB和GC含量都控制在23%~26.5%时，可用公式X计算出概率。
(12)
那么动态轮转映射算法成功控制各碱基对含量在23%~26.5%的成功率为
(13)
上述数学过程验证了通过动态轮转映射算法，每200个碱基AT、CG、PZ、SB的含量为23.5%~26.5%的成功率。接下来通过数学计算给出碱基序列长度和碱基控制成功率的关系。假设碱基序列长度为40、60、80、100、120、140、160、180、200，分别计算AT、CG、PZ、SB的含量为23.5%~26.5%的成功率。如图6所示，各碱基对含量调控的成功率与碱基序列的长度有关，碱基序列长度越长AT、CG、PZ、SB的含量越好控制。在应用动态轮转映射算法时，我们可以综合考虑各碱基对含量的控制范围、DNA合成长度和DNA存储成本等因素。在DNA存储过程中，主要包括DNA合成和测序两个关键步骤。DNA测序的成本受测序深度、技术类型（如Sanger测序或高通量测序）、目标（例如单基因或全基因组）以及样本数量等因素的影响[27]。而DNA合成的成本则与合成序列的长度、复杂性和数量密切相关，较长且复杂的DNA序列通常需要更高的合成费用[12]。因此，DNA合成的长度直接影响DNA存储的总体成本。在碱基对含量可控的范围内，合理选择合适的碱基序列长度，是优化DNA存储成本的关键策略。
图6碱基序列长度与碱基控制成功率的关系图
Fig.6Graph of base sequence length against base control success rate
为了支持上述的数学分析，我们对压缩后Mona Lisa_color.bmp获取的七进制数据流进行仿真，其中碱基序列长度为200 nt，控制每条碱基序列中AT、CG、PZ、SB含量范围为23%~26.5%。结果如图7所示，643条碱基序列的AT、CG、PZ、SB都符合设定范围。其中随机选取10条碱基序列，如图8所示，每条碱基序列中的AT、CG、PZ、SB含量都保持在25%左右。
图7每条碱基的AT、CG、PZ、SB含量
Fig.7AT, CG, PZ, SB content per base
图8十条碱基AT、CG、PZ、SB的含量
Fig.8Content of ten bases AT, CG, PZ, SB
2.2 压缩效果分析
2.2.1 哈夫曼编码压缩效果分析
测试七进制静态哈夫曼和自适应哈夫曼对各种文件的压缩效果。由于压缩后的数据流无法直接与原文件二进制数据流比较大小，将（01）数据流除以3作为原文件数据长度（八个核苷酸的理论存储密度为3 bits/nt）。选用文本、图片、音频分别进行压缩测试，记录每个文件的压缩率（压缩后的七进制数字个数/文件数据长度），在文本压缩时，为了追求低的压缩率，该实验采用直接读取的方式。如图9所示，七进制自适应哈夫曼编码和静态哈夫曼编码对文本、图片和音频都有压缩效果，其中对文本压缩效果最好。通过静态和自适应哈夫曼编码两条压缩率折线，可发现两种压缩编码对相同文件的压缩率相差不到1%。
图9哈夫曼编码压缩效果图
Fig.9Huffman code compression effect diagram
2.2.2 读取方式压缩效果分析
测试直接读取和常规读取两种方式对文本压缩的影响，测试文件选用Pride and Prejudice英文版（节选）和《傲慢与偏见》中文版（节选）。第一基于七进制自适应哈夫曼编码分别采用直接读取和常规读取进行压缩。如图10所示，采用直接读取的方式下，《傲慢与偏见》的压缩率为43.7%，消耗时间为488.34 s，Pride and Prejudice的压缩率为62.6%，消耗时间为246.22 s；而采用常规读取的方式下，《傲慢与偏见》的压缩率为80.6%，消耗时间为277.9 s，Pride and Prejudice的压缩率为63.2%，消耗时间为264.12 s。
图10自适应哈夫曼编码下不同读取方式的压缩效果图
Fig.10Compression effect graph of different reading methods under adaptive Huffman coding
第二基于七进制静态哈夫曼编码分别采用直接读取和常规读取进行压缩。如图11所示，采用直接读取的方式下，《傲慢与偏见》的压缩率为43.6%，消耗时间为0.81 s，Pride and Prejudice的压缩率为62.7%，消耗时间为0.36 s；而采用常规读取的方式下，《傲慢与偏见》的压缩率为80.56%，消耗时间为0.62 s，Pride and Prejudice的压缩率为63.18%，消耗时间为0.76 s。根据第一和第二实验可知，直接读取文本的压缩率低于常规读取的压缩率，在中文文本压缩中最为明显，直接读取中文文本的消耗时间高于常规读取的消耗时间。
图11静态哈夫曼编码下不同读取方式的压缩效果图
Fig.11Compression effect graph of different reading methods under static Huffman coding
综上所述，在压缩大部分文件时，静态哈夫曼编码的压缩率略低于自适应哈夫曼编码，消耗的时间少于自适应哈夫曼编码；但静态哈夫曼编码在不确定数据整体特征时，无法进行数据压缩，当额外信息的大小大于原信息大小时，静态哈夫曼编码的压缩率会高于自适应哈夫曼编码；而自适应哈夫曼编码可以适用于实时编码。因此在压缩动态数据时，选用自适应哈夫曼编码，压缩静态数据时，选用静态哈夫曼编码；在压缩文本文件时，采用直接读取，压缩其他文件时采用常规读取。此外，哈夫曼编码是一种无损压缩算法，在处理压缩过的文件时，哈夫曼编码的压缩效果会大大降低，甚至可能无法进一步压缩，这将导致该算法适合存储无损文件。
2.3 七进制汉明纠错码效果分析
接着对七进制汉明纠错码进行仿真。实验选用77715位数据流，分别采用七进制汉明（7、4）码、七进制汉明（15、11）码、七进制汉明（31、26）码、七进制汉明（63、57）码、七进制汉明（127、120）码进行纠错。将七进制汉明码纠错后的数据流编码成DNA序列，分别引入0.01%、0.02%、0.04%、0.08%、0.1%的出错率，模拟不同信道质量条件下的错码情况。将引入错误率的DNA序列进行解码，解码得到的七进制序列和未加入纠错码前的七进制序列进行对比，统计解码错码的数量。最后将错码数量除以序列总个数，得到错码频率。为防止实验结果具有偶然性，将所有随机错误频率下的错码引入和纠错操作均进行了100次试验，最终获得的错误频率是纠错之后100次试验的平均值，如表4所示。
从表4可以看出，在纠错77715位数据流时，七进制汉明码（7，4）、（15，11）、（31，26）均可以应对替换错误率在0.02%及以下的情况，可以纠正DNA序列中一个替换错误。DNA碱基序列产生错码具有随机性，当DNA序列越长，在每条引入纠错码的DNA序列上发生错误的位置可能越多。因此在DNA存储生物实现时，要根据信道质量（DNA合成、测序出现的错误率）、存储信息的大小（DNA序列的长度），纠错精度（错码频率）和DNA存储成本等综合考虑，选取最合适的七进制汉明码。
2.4 编码方法的综合效果分析
为了验证该方法的可行性和适用性，对整个流程进行了测试。由于在上述哈夫曼编码压缩实验中，已知自适应哈夫曼编码和静态哈夫曼编码的压缩效果和适用情况，因此在整体实验中，仅使用自适应哈夫曼编码进行实验。在实验中，对压缩转换后的七进制数据流进行分段，每150位数据为一段，每段添加6位地址，选用七进制汉明（31，26）码纠错，此时每段七进制数据流长度为186。经过映射编码后，每段DNA序列前添加三位首编码碱基组，如图12所示，每条DNA序列有189个碱基，150-nt存储信息，6-nt存储地址，30-nt存储纠错码，3-nt存储首编码碱基组。
表4七进制汉明码纠错结果
Table4Error Correction Results of Sevenary Hamming Code
图12短链DNA模型图
Fig.12Short-stranded DNA modeling diagram
以编码文件的比特数与编码的碱基个数（不包括引物区）的比值作为存储密度，即：存储密度=编码文件比特数/编码碱基数。本文以编码时间和存储密度作为衡量指标测试了该方法的存储性能，表5为自适应哈夫曼编码作为压缩编码存储不同文件的存储性能参数。
表5存储性能参数表
Table5Storage Performance Parameters Table
为了全方面的评价该算法的性能，将该算法和前人在非天然核苷酸存储编码算法的性能参数进行比较，对比结果如表6所示。Biswas等[40]（2020年）提出的方法没有考虑均聚物和GC含量对DNA序列的影响，也没有引入纠错码保证DNA存储信息的可靠性，理论存储密度为3 bits/nt。在2022年Biswas等[41]提出的算法理论存储密度可达到2.99 bits/nt，同时约束了均聚物长度和GC含量，但均聚物长度为小于等于7，还是会增加DNA测序的错误率。刘凯团队开发的编码系统控制了均聚物长度和GC含量，引入了纠错算法，可以纠正DNA序列丢失、缺失、替换多个错误，在8个碱基的编码中存储密度为3.64 bits/nt[16]。本文在添加纠错、地址和考虑均聚物长度和GC含量的前提下，最高存储密度为5.45 bits/nt，最低存储密度为2.45 bits/nt，平均存储密度为3.38 bits/nt；在生化约束方面，AT、CG、PZ、SB含量稳定保持在23%~26.5%，没有单碱基重复；值得一提的是，由于该算法内置七进制静态哈夫曼和自适应哈夫曼编码，因此可以对静态数据和实时数据进行存储，且在存储中文文本时，具有超越这三种算法的存储密度。但本算法的七进制汉明码只能纠正DNA序列单个替换错误，当DNA序列发生丢失或多个错误时，该算法DNA信息存储的可靠性会降低。
表6不同DNA编码算法比较表
Table6Comparison Table of Different DNA Encoding Algorithms
3 讨论
在基于非天然核酸的高密度DNA存储编码方法中，由于压缩部分使用了哈夫曼编码，因此该方法在存储文件类型上受哈夫曼编码影响，存储大型无损文本、图片、音频文件具有较高的存储密度，而在存储具有较高的时间和空间冗余性的视频数据时，存储密度较低。在映射部分，动态轮转映射算法通过遍历循环使GC含量和均聚物长度符合生化约束，增加了编码时间。在纠错部分，该方法的七进制汉明码只能纠正一位错误，当DNA序列发生缺失、插入或多个错误时，会降低DNA信息存储的可靠性。
综合考虑，本文基于非天然核酸的高密度DNA存储编码方法在存储无损文件具有较好的存储性能；同时在DNA存储方面，非天然核苷酸的出现不仅突破了每个天然核苷酸存储信息量的限制，提高了存储密度，而且有望进一步提高存取和读写效率、增强存储稳定性。鉴于非天然核苷酸的优势，本文采用了静态哈夫曼编码用来存储静态数据，内置了自适应哈夫曼编码应对存储流式数据的场景。然而，在使用自适应哈夫曼编码压缩时，由于需要实时更新哈夫曼树结构，消耗时间较长，后续我们将会对该算法的运行步骤进行并行化处理以及提高硬件环境，来减少运行时间。针对八进制汉明码纠错的局限性，下一步工作，我们将用RS码和LDPC码等纠错码来提高DNA存储的效率和可靠性。未来，随着生物合成学和DNA存储技术不断发展，DNA存储读取效率将逐步提升，本文提出的存储方法为流式数据和实时数据存储提供了一个解决方案，也为DNA信息存储提供了一种新思路。
