High-robustness DNA storage method based on specific base sequence

Abstract

The invention provides a high-robustness DNA storage method based on a specific base sequence. The high-robustness DNA storage method is used for solving the technical problem that local GC stability and global robustness are neglected in an existing storage scheme. The method comprises the following steps: constructing a '0-1' mapping rule and a specific base sequence; performing DNA coding on a binary sequence to be coded by using a 0-1 mapping rule and a specific base sequence to generate a rule index sequence and a data DNA sequence; compressing the rule index sequence by adopting multivariate Huffman coding, and coding the compressed sequence in combination with multivariate rotation coding to generate a DNA coding sequence; the obtained data DNA sequence is sequenced for multiple times and compared with a specific base sequence to detect and correct errors. According to the method, the error correction capability is improved while the local GC content is close to 50%, and the reliability and efficiency of DNA storage are optimized; under the condition that the average storage density is 1.33 bit / nt, the storage robustness is improved, and a new thought is provided for reliable storage of large-scale DNA data.

Description

Technical Field

[0001] The present invention relates to the technical field of DNA sequence coding, and in particular to a highly robust DNA storage method. Background Art

[0002] With the advent of the big data era, the rapid development of the Internet and the Internet of Things has accelerated the generation of large-scale data. 5G communication technology has greatly increased the data transmission rate. According to the International Data Corporation (IDC) forecast in 2021, the global data volume will grow at a compound annual growth rate of 23% from 2020 to 2025, and is expected to reach 180 zettabytes by 2025. However, traditional storage media such as optical disks, magnetic storage, and semiconductor storage are difficult to cope with the challenges of future big data storage due to their low storage density and short life. The DNA chain is composed of four nucleotides: adenine (A), thymine (T), guanine (G), and cytosine (C). With the development of DNA synthesis and sequencing technology, DNA storage provides a new solution for big data storage.

[0003] DNA storage has significant advantages of high information storage density, low maintenance cost, strong stability and long life. As a new type of storage medium, DNA has the advantages of high storage density, long-term stability and low energy consumption. However, traditional storage solutions usually focus on high-density storage and ignore local stability and decoding robustness. DNA storage still faces many technical challenges. During DNA synthesis, amplification and sequencing, single nucleotides may have insertion, deletion and substitution errors (IDs). To solve these problems, many studies have tried to reduce errors by constraining DNA encoding, that is, balancing GC content (guanine and cytosine content between 40%-60%) or controlling homopolymer length (less than 4nt). Goldman et al. developed a 3-base rotation encoding scheme to control the homopolymer length to less than 2nt, but could not ensure GC balance. Grass et al. associated finite fields with base triplets and proposed an encoding algorithm that can avoid single-base repetitions with a length greater than 3. However, this algorithm did not solve the problem of regulating GC content under specific data structures. Mishra et al. proposed an effective method to compress digital data into DNA sequences using the minimum variance Huffman tree, combining compression coding with constraint coding, and the generated DNA sequences follow the GC content constraint and the homopolymer constraint of at most 1 nt. Blawat et al. used fixed mapping rules and optional mapping rules to interleave bytes to ensure that the homopolymer length is less than 4 nt. Erlich et al. proposed a fountain coding scheme, which packs binary numbers into droplets, and then performs Luby transformation and screening to obtain DNA sequences that meet both constraints. Wang et al. proposed a coding scheme that reduces redundancy by hiding addressing, thereby reducing the number of bases required for storage; the above coding scheme has good adaptability to indicators with low self-similarity and stable local GC content, reduces the error rate of DNA synthesis and sequencing, improves base utilization, and optimizes the reliability and efficiency of DNA storage. Ping et al. proposed a Yin-Yang coding and decoding system (YYC), which uses Yin and Yang coding rules to generate multiple coding schemes, and screens them to obtain DNA sequences that meet the constraints. However, the above schemes fail to take into account local GC balance and lack local stability. Zhao et al. proposed a De Bruijn Trim rotation graph encoding based on graph theory, which divides k-mers into De Bruijn Trim graphs and compresses the stored information according to the overlapping relationship to achieve high-density storage, which not only limits the length of homopolymers but also satisfies the local GC balance. However, its decoding relies only on error correction codes, and the fault tolerance rate is low.

[0004] In recent years, in order to improve the efficiency and accuracy of DNA storage data recovery, many studies have been devoted to error correction algorithms for substitution, insertion and deletion errors in DNA storage channels. For example, Yadi et al. developed an integrated pipeline to gradually correct IDs in multiple stages. The algorithm consists of a multiple sequence alignment (MSA) module, an iterative Burrows-Wheeler alignment (BWA) module and a post-processing module. The MSA module is used to obtain a consistent sequence, and the BWA module refines the sequence through some coding constraints to improve the reconstruction quality. Xue et al. constructed a special Levenshtein code that can correct single synchronization or substitution errors in codewords, but it is only suitable for small-scale data storage. Press et al. proposed a special convolutional code HEDGES, which uses a hash-based convolutional code as an internal code, which can directly correct the index, and the external RS block code reconstructs the fragments that the internal code cannot correct, but the recovery rate is significantly reduced when the error rate is greater than 10%. Song et al. constructed a DeBruijn graph (DBG) using multiple sequences and inferred the correct sequence through a graph search algorithm, but the recovery rate is significantly reduced when the error rate is greater than 10%. Zan et al. constructed a 6-base sequence codebook for English text and used a three-layer algorithm of Hamming distance, majority voting and spelling check to correct IDs. The data can be efficiently recovered at a low error rate, but the recovery rate is only 60% after the error rate reaches 10%. Moreover, this algorithm is only applicable to English text and has certain limitations. Welzel et al. proposed a coding scheme that supports variable-size coding sequences and can correct errors such as substitution, insertion and loss. Zhang et al. proposed a HMSA algorithm that combines Hamming distance error correction with multiple sequence alignment, which can correct mixed errors and has a certain universality for coding algorithms that rely on codebooks, but the error correction efficiency is low under large-scale data.

[0005] The error rate of second-generation synthesis and sequencing technology is about 1%-2%. Although the third-generation sequencing technology is faster and cheaper, its error rate is higher, about 10%-15%. Extreme GC content will also increase the error rate in DNA sequencing and synthesis. Currently, many storage solutions tend to be high-density, but ignore local GC stability and global robustness. As the amount of data stored increases, maintaining the integrity of data recovery requires more efficient error correction algorithms. Summary of the invention

[0006] In view of the technical problem that existing storage solutions ignore local GC stability and global robustness, the present invention proposes a storage method based on specific base sequences (DNA storage based on fixed base sequences, DNA-FBS), which encodes data into base sequences that meet constraints and have specific patterns, and uses the characteristics of mutant bases destroying specific patterns to perform synchronous error correction. When the non-continuous error rate is below 25%, more than 90% of the data can be recovered through only five sequencings. In addition, the present invention also proposes a multivariate Huffman rotation code for storing regular index sequences generated by DNA-FBS to improve storage density and maintain local stability. The present invention improves the error correction capability while ensuring that the local GC content is close to 50%, and optimizes the reliability and efficiency of DNA storage.

[0007] In order to achieve the above object, the technical solution of the present invention is implemented as follows: a highly robust DNA storage method based on a specific base sequence, the steps of which are as follows:

[0008] Step 1: Construct "0-1" mapping rules and specific base sequences;

[0009] Step 2: Use the "0-1" mapping rule and the specific base sequence to perform DNA encoding on the binary sequence to be encoded, and generate a rule index sequence and a data DNA sequence;

[0010] Step 3: compress the regular index sequence using multivariate Huffman coding, and encode the compressed sequence in combination with multivariate rotation coding to generate a DNA coding sequence;

[0011] Step 4: Sequence the data DNA sequence obtained in step 2 multiple times and compare it with the specific base sequence to detect and correct errors.

[0012] Preferably, the "0-1" mapping rule consists of a "0" rule and a "1" rule, wherein the "0" rule and the "1" rule respectively map each element in the binary set B = {00, 01, 10, 11} to a base in the DNA base set N = {A, C, G, T}, and satisfy two conditions: ① The "0" rule and the "1" rule are both one-to-one mappings, i.e., N 01 ≠N 02 ≠N 03 ≠N 04 , N 11 ≠N 12 ≠N 13 ≠N 14 ; Among them, N 01 、N 02 、N 03 、N 04 For the elements in the "0" rule, N 11 、N12 、N 13 、N 14 is an element in the "1" rule; ② The bases A / T and G / C are mutually exclusive;

[0013] According to the "0-1" mapping rule, bases A and T are regarded as one type of base, represented by P; bases G and C are regarded as another type of base, represented by O. The specific base sequence S mould P and O each account for half, and no more than three P or O can appear consecutively.

[0014] Preferably, the mathematical description of the "0-1" mapping rule is

[0015]

[0016] Where i = 0, 1, j = 1, 2, 3, 4;

[0017] There are 48 "0-1" mapping rules that meet the conditions:

[0018]

[0019]

[0020] The specific base sequence is a base sequence that appears periodically and repeatedly in a DNA storage sequence, and is composed of multiple identical sequences that satisfy specific constraints, wherein the specific constraints include a base combination pattern, a sequence repetition period, and a sequence structural stability; through any "0-1" mapping rule, a binary sequence is mapped into a base sequence composed of periodically repeated specific base sequences, thereby achieving efficient encoding and storage of data;

[0021] The mathematical description of a specific basis sequence of length l is:

[0022]

[0023] Among them, N 1 ,N 2 ,…,N l-1 ,N l Respectively represent l elements of a specific base sequence; half of them are A or T, and the other half are C or G;

[0024] The length l of the specific base sequence is an even number.

[0025] Preferably, the binary sequence to be encoded is C=(c 1 c 2 …c i-1 c i …c n-1 c n ), the i-th binary element ci ∈{0,1}, n is the length of the binary sequence C. Using the "0-1" mapping rule and a specific base sequence for DNA encoding is to map the binary sequence C to a specific base sequence that repeats and recurs. The implementation method is:

[0026] Step 1. Grouping: Group every two bits in the binary sequence C into groups and select a specific base sequence of a given length l;

[0027] Step 2. Determine the base: The first group c of the binary sequence C 1 and c 2 The dth group of binary sequence C is mapped to the dth position of the specific base sequence, and the corresponding base is determined according to the "0-1" mapping rule constraint, and the index number of the "0-1" mapping rule is recorded at the same time; by analogy, the dth group of binary sequence C is mapped to the dth position of the specific base sequence, and the corresponding base of the mapping is determined, and d <l;

[0028] Step 3. Cyclic mapping: The l+1th group of the binary sequence C is mapped to the first bit of the specific base sequence. When the mapping is successful, the index number of the "0-1" mapping rule is recorded. Similarly, every two bits in the binary sequence C are mapped and encoded to obtain the complete encoding result as the data DNA sequence, and the rule index sequence S is generated. rules .

[0029] Preferably, the multivariate Huffman coding is M-ary Huffman coding, and the regular index sequence is compressed into an M-ary sequence through M-ary Huffman coding. The multivariate rotation coding uses a rotation coding table to map the M-ary sequence into a base sequence that meets the constraints, and introduces RS code to ensure stable data recovery.

[0030] Preferably, the method for constructing the M-ary Huffman tree by M-ary Huffman coding is: rules Convert each t-bit binary number into decimal data and count the probability of each decimal data. The decimal data is called the source symbol. Assuming there are f source symbols, arrange them in order according to their probabilities. Give M source symbols with the smallest probabilities (a f-M-1 ,a f-M-2 ,…a f ) are assigned code words (0, 1, 2, ... M-1) respectively and the M source symbols are merged into a new symbol, and a new source with a length of (fM-1) is obtained, which is called the first source reduction; then the reduced new source is arranged again according to the probability from small to large and the source reduction process is repeated to obtain the M-ary Huffman tree with the shortest average code length;

[0031] The multivariate rotation coding is: for the obtained M-ary sequence, set the 0th code to be any one of the three bases, and look up the rotation coding table according to the 0th three bases to obtain the three bases mapped to the decimal number to be encoded in the first M-ary sequence, and look up the rotation coding table according to the first three bases to obtain the three bases mapped to the decimal number to be encoded in the second M-ary sequence; and so on, after rotating and mapping the M-ary sequences one by one, finally obtain the base sequence S 2 .

[0032] Preferably, the method for constructing the rotational coding table is as follows: under the condition of satisfying the homopolymer constraint and the GC content constraint, the base combination length in the rotational coding table is expanded to 3; there are 64 three-base combinations, and in order to satisfy the constraint that the maximum homopolymer length is less than 4, the combinations in which the first two bases are the same are removed, and the remaining 48 are used as rotational coding codewords, and the 48 three-base combinations are represented as rotational coding codewords H:(H 1 ,H 2 ,...H 47 ,H 48 ), and the rotation codeword is

[0033] <![CDATA[H 1 ]]> ATA <![CDATA[H 9 ]]> AGA <![CDATA[H 17 ]]> TCA <![CDATA[H 25 ]]> CAA <![CDATA[H 33 ]]> CGA <![CDATA[H 41 ]]> GTA <![CDATA[H 2 ]]> ATT <![CDATA[H 10 ]]> AGT <![CDATA[H 18 ]]> TCT <![CDATA[H 26 ]]> CAT <![CDATA[H 34 ]]> CGT <![CDATA[H 42 ]]> GTT <![CDATA[H 3 ]]> ATC <![CDATA[H 11 ]]> AGC <![CDATA[H 19 ]]> TCC <![CDATA[H 27 ]]> CAC <![CDATA[H 35 ]]> CGC <![CDATA[H 43 ]]> GTC <![CDATA[H 4 ]]> ATG <![CDATA[H 12 ]]> AGG <![CDATA[H 20 ]]> TCG <![CDATA[H 28 ]]> CAG <![CDATA[H 36 ]]> CGG <![CDATA[H 44 ]]> GTG <![CDATA[H 5 ]]> ACA <![CDATA[H 13 ]]> TAA <![CDATA[H 21 ]]> TGA <![CDATA[H 29 ]]> CTA <![CDATA[H 37 ]]> GAA <![CDATA[H 45 ]]> GCG <![CDATA[H 6 ]]> ACT <![CDATA[H 14 ]]> TAT <![CDATA[H 22 ]]> TGT <![CDATA[H 30 ]]> CTT <![CDATA[H 38 ]]> GAT <![CDATA[H 46 ]]> GCC <![CDATA[H 7 ]]> ACC <![CDATA[H 15 ]]> TAC <![CDATA[H 23 ]]> TGC <![CDATA[H 31 ]]> CTC <![CDATA[H 39 ]]> GAC <![CDATA[H 47 ]]> GCA <![CDATA[H 8 ]]> ACG <![CDATA[H 16 ]]> TAG <![CDATA[H 24 ]]> TGG <![CDATA[H 32 ]]> CTG <![CDATA[H 40 ]]> GAG <![CDATA[H 48 ]]> GCT

[0034] The positions of the 48 tri-bases are cyclically shifted to form a rotation coding matrix of size 48×48, thereby obtaining a rotation coding table. The 48-element rotation coding table is:

[0035]

[0036] Preferably, the error correction method is: first, the obtained data DNA sequence is sequenced multiple times to generate multiple read lengths, each read length is converted to OP to represent the base sequence of a specific base sequence and compared and aligned, so that the lengths of all read lengths are consistent; secondly, the weights of the bases at the corresponding positions of all read lengths are calculated, and the optimal base sequence to be decoded is obtained according to the weights; from the obtained rule index sequence, the corresponding decoding rule is selected in combination with the value of its corresponding position to restore the base sequence to be decoded, and finally the original binary sequence is obtained.

[0037] Preferably, in the process of data recovery, the base sequence S generated by the multi-element rotation encoding is 2 Divide into groups of three, and the default 0th base is H when encoding. 1 , when decoding the first three bases, the previous three bases are H 1 Look up the current rotation code table, and the value obtained is the code word corresponding to the first three bases. This process is iterated until the entire base sequence S is traversed. 2 ; According to the M-ary Huffman tree, the codeword sequence is decoded into a decimal sequence, and then the regular index sequence S is obtained rules ;

[0038] With the help of rule index sequence S rules The obtained base sequence S 1 The decoding method is: for the base sequence S 1 The bases in the sequence S are decoded bit by bit according to the rule index. rules The value at the corresponding position selects the corresponding decoding rule, specifically, in the decoding base sequence S 1 When the i-th base is found, the sequence S is indexed according to the rule rules The i-th bit value of is decoded by selecting the "0" rule or the "1" rule.

[0039] Preferably, the average code length is the number of bits required for each codeword in the sequence generated by M-ary Huffman coding, and M is selected as 48 according to the average code length;

[0040] In the encoding process, a specific base sequence with low error weight should be selected first; the length l of the specific base sequence is 8;

[0041] Before encoding the data to be stored, all mapping rules in the "0-1" mapping rule are used to encode the data to be stored, the occurrence frequency of the non-expected base sequence in the obtained base sequence is counted, and the mapping rule with the smallest frequency is selected.

[0042] Compared with the prior art, the beneficial effects of the present invention are as follows: a coding scheme based on a specific base sequence is proposed. First, 48 "0-1" mapping rules and specific base sequences are constructed to encode data into base sequences with a specific pattern, aiming to improve the error correction and error detection capabilities of the data. The storage of the rule index sequence is realized by combining the multivariate Huffman rotation coding, which not only improves the coding efficiency, but also maintains the local balance of the base sequence. For errors in the synthesis and sequencing process, the method of specific base sequence alignment is used to detect the error position, and the data is restored by combining multiple sequencing, which improves the efficiency of large-scale data recovery. The simulation results show that the local GC content of the encoded DNA sequence is between 45% and 55%, avoiding the generation of long homopolymers. At an error rate of 20%, only 5 sequencings are required to almost completely restore the data, and more errors can be tolerated by increasing the length of the specific base sequence and the number of sequencings. The technical solution of the present invention improves the robustness of storage under the condition of an average storage density of 1.33bit / nt, and provides a new idea for the reliable storage of large-scale DNA data. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0044] Figure 1 It is a flow chart of the present invention.

[0045] Figure 2 This is an example diagram of the encoding based on a specific base sequence of the present invention.

[0046] Figure 3 This is an example diagram of the 48-tuple Huffman rotation encoding of the present invention.

[0047] Figure 4 It is a schematic diagram of the error correction process based on a specific base sequence of the present invention.

[0048] Figure 5 Schematic diagram of the experimental analysis of the present invention, wherein (a) is the comparison of local GC content, (b) is the number of undesired base sequences in 48 regular mapping schemes, (c) is the error rate distribution of 12 3-mers, and (d) is the error weight corresponding to all specific base sequences that meet the coding constraints when the specific base sequence length is l=2, 4, and 8, and the x-axis 0-4 corresponds to l=4, the x-axis 5-6 corresponds to l=2, and the x-axis 7-62 corresponds to l=8.

[0049] Figure 6 1 is a comparison diagram of the error correction capabilities of different specific base sequences of the present invention.

[0050] Figure 7 This is an example diagram of the error correction capability when the length of the specific base sequence of the present invention is 8.

[0051] Figure 8 The specific base sequence length l and sequencing times N of the present invention under the continuous error rate of 20% and 40% r Comparison of error correction capabilities at different times, where (a) is 20% and (b) is 40%.

[0052] Fig. 9 This is a coding density comparison diagram when t takes different values ​​in the present invention. DETAILED DESCRIPTION

[0053] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0054] A highly robust DNA storage method based on specific base sequences proposes a coding scheme based on specific base sequences, which mainly includes the construction of "0-1" mapping rules, the design of specific base sequences, and the construction of multi-factor rotation coding matrices. First, the "0-1" mapping rule is constructed and the specific base sequence is designed; secondly, the data is encoded by DNA through the specific base sequence to generate a rule index sequence and a data DNA sequence. Finally, the Huffman technique is used to compress the rule index sequence, and combined with the multi-factor rotation coding matrix to generate a DNA sequence. This scheme detects and corrects errors by comparing the base sequence obtained by sequencing with the specific base sequence. The overall storage scheme is as follows: Figure 1 The specific implementation steps of the present invention are:

[0055] Step 1: Construct a "0-1" mapping rule to encode information into a DNA base sequence.

[0056] In order to meet biological constraints and improve sequence stability, a "0-1" mapping rule is constructed to encode information into a DNA base sequence. The "0-1" mapping rule consists of a "0" rule and a "1" rule, which maps each element in the binary set B = {00, 01, 10, 11} to a base in the DNA base set N = {A, C, G, T}. This process must meet two conditions: ① Both the "0" rule and the "1" rule are one-to-one mappings, that is, N 01 ≠N 02 ≠N 03 ≠N 04 , N 11 ≠N 12 ≠N 13 ≠N 14 ; Among them, N 01 、N 02 、N 03 、N 04 、N 11 、N 12 、N 13 、N 14 are all elements in the set N; ② The bases A / T and G / C are mutually exclusive, that is, A and T are one type of bases, G and C are another type of bases, and the two types of bases are mutually opposed. For example, if 00 is mapped to A or T in the "0" rule, then 00 must be mapped to G or C in the "1" rule. The mathematical description of the "0-1" mapping rule is shown in formula (1).

[0057]

[0058] Among them, i=0,1, j=1,2,3,4.

[0059] According to the above definition, there are 48 "0-1" mapping rules that meet these two conditions, as shown in Table 1. In the actual encoding process, one of them can be selected. At the same time, the probability of undesired base sequences appearing can be reduced by further screening these 48 mapping rules.

[0060] Table 1. 48 "0-1" mapping rules applied to specific base sequence encoding

[0061]

[0062] A specific base sequence is a base sequence that repeats periodically in a DNA storage sequence. It consists of multiple identical sequences that satisfy specific constraints. These specific constraints include the base combination pattern, the repetition period of the sequence, and the structural stability of the sequence. The design of the specific base sequence aims to improve the error correction and error detection capabilities of the data, because its periodic generation makes the detection and correction of errors more reliable and efficient. Through any "0-1" mapping rule, the binary sequence can be mapped to a base sequence composed of a specific base sequence, thereby achieving efficient encoding and storage of data. The mathematical description of a specific base sequence of length l is as follows:

[0063]

[0064] Among them, N 1 ,N 2 ,…,N l-1 ,N l They represent l elements of a specific base sequence, half of which are A or T and the other half are G or C.

[0065] To better balance the GC content, the length l of the specific base sequence is specified to be an even number. According to the "0-1" mapping rule, bases A and T are considered as one type of base, represented by P; bases G and C are considered as another type of base, represented by O. In order to keep the GC content of the generated base sequence balanced and avoid the generation of long homopolymers, the specific base sequence S is specified to be an even number. mould P and O each account for half, and no more than three P or O can appear consecutively. The number of specific base sequences of different lengths is shown in Table 2. For example, when the length l is 2, the specific base sequence S mould It can be "OP" or "PO". When the length l is 4, "OPOP" and "POPO" will appear, which are the same as when the length l is 2 and need to be discarded. The remaining four are "OPPO", "POOP", "OOPP", and "PPOO".

[0066] Table 2 The number of specific base sequences that meet the encoding constraints at different lengths

[0067] Specific base sequence length l 2 4 6 8 10 12 14 16 Number of available specific base sequences 2 4 18 56 192 592 1970 6282

[0068] Step 2: Use the "0-1" mapping rule and specific base sequence to perform DNA encoding on the binary sequence to be encoded, and generate a rule index sequence and a data DNA sequence.

[0069] For the binary sequence to be encoded C = (c 1 c 2 …c i-1 c i …c n-1 c n ), where the i-th binary element c i ∈{0,1}, n is the length of the binary sequence C, and the DNA encoding process is to map it into a specific base sequence that repeats the cycle. The specific description is as follows:

[0070] Step 1. Group mapping: Map every two bits in the binary sequence C as a group and select a specific base sequence of a given length l.

[0071] Step 2. Determine the base: The first group c of the binary sequence C 1 and c 2 Mapped to the first bit of a specific base sequence, determine the corresponding base according to the mapping rule constraints, and record the index number of the mapping rule; and by analogy, the dth group of the binary sequence C is mapped to the dth bit of a specific base sequence, and determine the corresponding base of the mapping.

[0072] Step 3. Cyclic mapping: Since the base sequence generated by the mapping is composed of repeated specific base sequences, the l+1th group of the binary sequence C is also mapped to the first bit of the specific base sequence. When the mapping is successful, the index number of this mapping rule is recorded. Similarly, every two bits in the binary sequence are mapped and encoded to obtain a complete encoding result, and a rule index sequence S is generated. rules .

[0073] To explain this encoding process in detail, the binary sequence "1001011001000101100111001101001001" and the specific base sequence "OOPOPPOP" are used as examples for mapping, and the encoding is formed into a base sequence of the "OOPOPPOPOOPOPPOPO..." pattern. The first mapping rule in Table 1 is selected. First, the binary sequence is grouped into groups of two digits each, and the first group "10" is mapped to the first bit O of the specific base sequence "OOPOPPOP", that is, the base G or C. According to the "0-1" mapping rule, only the "0" rule can map "10" to the base "G", and the rule index number "0" is recorded at the same time; the second group "01" is mapped to the second bit O of the specific base sequence "OOPOPPOP", and is mapped to "G" through the "1" rule, and the rule index number "1" is recorded. By analogy, the final coding sequence, i.e., the data DNA sequence, is "GGTGTAGTGGTCTTCAG". At the same time, the index numbers of the mapping rules are sequentially generated into the rule index sequence S rules is "01000010011110111", the specific process is as follows Figure 2 shown.

[0074] Step 3: Use multivariate Huffman rotation coding to compress the regular index sequence and combine it with the DNA coding sequence generated by multivariate rotation coding.

[0075] After the data to be stored is encoded based on a specific base sequence, a rule index sequence S of the "0-1" mapping rule will be generated. rules Given the high cost of DNA storage ($800 million / TB), an efficient, biologically-constrained, and robust encoding scheme is needed. Encoding rule index sequence S rules To this end, the rotation coding is combined with the multi-element Huffman coding, and a 48-element Huffman rotation coding scheme is proposed to encode the regular index sequence S rules to encode.

[0076] Rotational coding is one of the most important technologies highly respected in the field of data communication and storage. It can effectively control the length of homopolymers in DNA storage and meet homopolymer constraints. The ternary rotational coding proposed by Goldman determines the type of base to be mapped next based on the value of the previous base and the current ternary codeword, ensuring that the two adjacent bases in the generated base sequence are different and the maximum homopolymer length is 1. Due to the characteristics of rotational coding, the four bases in the generated base sequence are evenly distributed, so rotational coding can not only control the length of homopolymers but also effectively balance the GC content.

[0077] The basic base combination length in Zhao's rotational coding table is 2, and there are 16 types in total. In order to make the homopolymer length in the generated DNA sequence less than 4, the same base combination is removed, a total of 4 types, and the maximum rotational coding takes twelve elements. In order to further improve the coding efficiency, the base combination length in the rotational coding table is expanded to 3 under the conditions of satisfying the homopolymer constraint and GC content constraint. There are 64 three-base combinations in total. In order to meet the constraint that the maximum homopolymer length is less than 4, the combination with the same first two bases is removed, and the remaining 48 types can be used as rotational coding codewords. The 48 three-base combinations are represented as rotational coding codewords H: (H 1 ,H 2 ,...H 47 ,H 48 ), as shown in Table 3. By cyclically moving the positions of these 48 tri-bases, a rotation coding matrix of size 48×48 is formed, thereby obtaining a rotation coding table, as shown in Table 4. Since the number of "A" / "T" and "G" / "C" in these tri-base combinations is the same, all tri-base combinations in the generated DNA sequence are evenly distributed, satisfying the GC content and homopolymer constraints.

[0078] Table 3 Forty-eight-element rotation codewords

[0079] <![CDATA[H 1 ]]> ATA <![CDATA[H 9 ]]> AGA <![CDATA[H 17 ]]> TCA <![CDATA[H 25 ]]> CAA <![CDATA[H 33 ]]> CGA <![CDATA[H 41 ]]> GTA <![CDATA[H 2 ]]> ATT <![CDATA[H 10 ]]> AGT <![CDATA[H 18 ]]> TCT <![CDATA[H 26 ]]> CAT <![CDATA[H 34 ]]> CGT <![CDATA[H 42 ]]> GTT <![CDATA[H 3 ]]> ATC <![CDATA[H 11 ]]> AGC <![CDATA[H 19 ]]> TCC <![CDATA[H 27 ]]> CAC <![CDATA[H 35 ]]> CGC <![CDATA[H 43 ]]> GTC <![CDATA[H 4 ]]> ATG <![CDATA[H 12 ]]> AGG <![CDATA[H 20 ]]> TCG <![CDATA[H 28 ]]> CAG <![CDATA[H 36 ]]> CGG <![CDATA[H 44 ]]> GTG <![CDATA[H 5 ]]> ACA <![CDATA[H 13 ]]> TAA <![CDATA[H 21 ]]> TGA <![CDATA[H 29 ]]> CTA <![CDATA[H 37 ]]> GAA <![CDATA[H 45 ]]> GCG <![CDATA[H 6 ]]> ACT <![CDATA[H 14 ]]> TAT <![CDATA[H 22 ]]> TGT <![CDATA[H 30 ]]> CTT <![CDATA[H 38 ]]> GAT <![CDATA[H 46 ]]> GCC <![CDATA[H 7 ]]> ACC <![CDATA[H 15 ]]> TAC <![CDATA[H 23 ]]> TGC <![CDATA[H 31 ]]> CTC <![CDATA[H 39 ]]> GAC <![CDATA[H 47 ]]> GCA <![CDATA[H 8 ]]> ACG <![CDATA[H 16 ]]> TAG <![CDATA[H 24 ]]> TGG <![CDATA[H 32 ]]> CTG <![CDATA[H 40 ]]> GAG <![CDATA[H 48 ]]> GCT

[0080] Table 4 48-element rotation coding table

[0081]

[0082] In order to improve the coding efficiency, reduce the total number of code elements, and simplify the coding process, on the basis of satisfying the GC content and homopolymer constraints, the multivariate rotation coding is combined with the M-ary Huffman coding, and the M-ary Huffman rotation coding is proposed. Specifically, the regular index sequence is compressed into an M-ary sequence by M-ary Huffman coding, and the M-ary sequence is mapped into a base sequence that meets the constraints by using multivariate rotation coding, wherein the RS code is introduced, and the received codeword is substituted into the RS check matrix during decoding, and the decoder calculates the syndrome, and uses the syndrome and the known RS code to generate polynomials and construct a system of equations. By solving this system of equations, the decoder can determine the location of the error and the number of errors, thereby ensuring stable data recovery.

[0083] In M-ary Huffman rotation coding, the key is to construct an M-ary Huffman tree. Figure 3 As shown, the rule index sequence S rules Convert each t-bit binary number into decimal data and count the probability of each decimal symbol. The decimal symbols are collectively referred to as source symbols. Arrange the source symbols in ascending or descending order according to their probabilities, that is, p(a 1 )≥p(a2 )≥…≥p(a f ), where a 1 , a 2 …a f represents f source symbols. Given M source symbols with the smallest probability (a f-M-1 ,a f-M-2 ,…a f ) assigns codewords (0,1,2,…M-1) and merges the M source symbols into a new symbol, obtaining a new source with a length of (fM-1). This process is called the first source reduction. After that, the reduced sources are arranged again according to the probability from small to large and the source reduction process is repeated, and finally the M-ary Huffman tree with the shortest average code length can be obtained.

[0084] In the present invention, the maximum rotation code is 48 yuan, and the corresponding M in the Huffman code can be any integer value less than 49. The M value affects the Huffman coding efficiency, so a suitable M (M≤48) value should be selected in the coding process to improve the coding efficiency. The number of bits required for each codeword in the sequence generated by the M-ary Huffman code is called the average code length L avg , mathematically described as formula (3):

[0085]

[0086] Among them, p i represents the probability of the i-th source symbol appearing, l i Indicates the code length corresponding to the i-th source symbol.

[0087] Table 5 Average code length for encoding files of different sizes when M takes different values

[0088]

[0089] The average code length is an indicator to measure the coding efficiency. Generally, the shorter the average code length, the higher the coding efficiency, that is, fewer bits can be used to represent information, thereby achieving a better coding effect. In order to select a suitable M value, Huffman coding is performed on files of different sizes when M is 45, 46, 47 and 48 respectively, the packet length t is 8, and the corresponding average code length is counted. The results are shown in Table 5. After performing M-ary Huffman coding on files of sizes of 193kb, 510kb, and 769kb, it is found that as the M value increases, the average code length gradually decreases until the average code length is the smallest when M reaches 48, indicating that the coding efficiency is highest when M is 48. Therefore, the present invention takes M value as 48.

[0090] According to the constructed 48-element Huffman tree, the encoding result of each source symbol can be obtained. For the given information, that is, the regular index sequence S rulesAfter M-ary Huffman coding, a 48-ary symbol sequence is obtained. When performing rotation coding, let the 0th code be the three-base H 1 (For subsequent rotation encoding, the initial encoding base combination needs to be given, and the default 0th three bases are H 1 ).like Figure 3 As shown, for the 48-element sequence "30 16 0 10 18 0 8 8 1 37 39 13 0 5 24 31 3 ...", according to the 48-element rotation coding table mapping in Table 4, the three bases mapped by the current 48-element symbol depend on its own value and the previous three-base codeword. The 0th three bases are H 1 That is, when it is "ATA", the table shows that the first symbol to be encoded "30" should be mapped to H 31 That is, "CTC". When the current three bases are "CTC", the current symbol to be encoded "16" corresponds to H 47 That is, "GCA", so "16" is mapped to "GCA". After rotating and mapping the 48-tuple sequence one by one, the final base sequence S is obtained. 2 It is “CTCGCAGCAAGACACCACCGCGTCGTG…”

[0091] Step 4: Sequence the data DNA sequence obtained in step 2 multiple times and compare it with the specific base sequence to detect and correct errors.

[0092] In the process of DNA synthesis and sequencing, various types of errors (such as insertion, deletion, substitution) are inevitable, and most existing solutions rely on error correction codewords to improve data recovery capabilities. However, with the increase in the accuracy requirements for data recovery in the future, a more effective error correction scheme is needed. The coding scheme based on a specific base sequence proposed in the present invention converts data into a data DNA sequence that conforms to a certain rule. This sequence is composed of a combination of specific base sequences and changes periodically. In fact, the base type at each position is fixed and can only be A / T (P) or G / C (O). When an insertion / deletion error occurs, this fixed pattern will be broken, so that the position where the error occurs can be detected. By comparing it with a specific base sequence, the error can be further corrected. This error correction scheme that relies on coding characteristics is called error correction based on a specific base sequence.

[0093] For example, if the specific base sequence "OOPOPPOP" with a length of 8 is used for encoding, the correct data DNA sequence is "OOPOPPOPOOPOPPOP...". When an error occurs in this sequence, the length of the sequence and the bases at the corresponding position may change. Suppose the original sequence is "CGTCAAGT" and the specific base sequence is S mould =OOPOPPOP, the read length 2 after the introduction of the error is "CAGTCAAGT", which is converted to OP and represented as Sreads2 =OPOPOPPOP, compare and align the two to get "OPOPOPPOP", which means that an insertion error may occur in the second position. Delete the inserted "P" to get "OOPOPPOP", that is, "CGTCAAGT". At this time, the sequence has been aligned with the original sequence; read length 4 "CGTAAGT", OP is represented by S reads4 =OOPPPOP, and compare it with S mould Compare alignment, get OOP O PPOP, a deletion error may occur at the fourth position, thereby detecting the position where the insertion error occurred. Combining multiple sequencing, the consensus sequence obtained is "CGTCAAGT"; read length 5 is "CGTCAAGC", which is converted to OP and represented as S reads5 =OOPOPPOO, and compare it with the specific base sequence S mould Comparison revealed that a substitution error occurred at the eighth position, and combined with multiple sequencing attempts, the consensus sequence "CGTCAAGT" was obtained.

[0094] like Figure 4 As shown in Figure 1, the overall error correction process can be described as follows: first, the obtained data DNA sequence is sequenced multiple times to generate multiple read lengths, each read length is converted to OP representation and compared and aligned with a specific base sequence. At this time, the lengths of all read lengths are consistent. Secondly, the weights of various bases at the corresponding positions of all read lengths are calculated, and the bases with the largest weight at each position are output and combined to obtain the best base sequence to be decoded. Finally, the rule index sequence S is combined with the base sequence to obtain the optimal base sequence to be decoded. rules The original binary sequence is obtained by rotating the encoding table with forty-eight elements.

[0095] In the process of data recovery, we first need to obtain the regular index sequence S rules . The base sequence S encoded by the 48-tuple Huffman rotation 2 Divide into groups of three, and the default 0th base is H when encoding. 1 , so when decoding the first three bases, the previous three bases are H 1 Look up the current rotation code table, and the value obtained is the code word corresponding to the first three bases. This process is iterated until the entire base sequence S is traversed. 2 Finally, the codeword sequence is decoded into a decimal sequence according to the M-ary tree, and then the regular index sequence S is obtained. rules . Base sequence S 2 Base sequence S with a specific pattern decoded by RS 1 .

[0096] The second process of data recovery is to use the rule index sequence S rules The base sequence S obtained based on a specific base sequence encoding scheme 1Decoding, the two sequences have the same length, so for the base sequence S 1 The bases in the bitwise decoding need to refer to the regular index sequence S rules The value at the corresponding position selects the corresponding decoding rule, specifically, in decoding S 1 At the i-th base in the sequence, according to S rules The i-th bit value (0 or 1) selects 0 rule or 1 rule decoding.

[0097] The present invention will now simulate and analyze storage density, coding performance, error correction capability, etc., to verify the fault tolerance of the proposed error correction scheme and the performance of the generated DNA sequence.

[0098] The coding scheme proposed in the present invention shows excellent performance in terms of local GC content and maximum homopolymer length. When the length of the specific base sequence is 2, the maximum homopolymer length of the data sequence is 1, and the maximum homopolymer length of the regular index sequence is 3. The GC content of the present invention is maintained between 49% and 51%. In order to verify the local GC content of the coding scheme of the present invention, the same binary sequence is encoded using DNA Fountain and the present invention respectively, the first 1500nt base sequence is intercepted, and the GC content is counted every 15nt. The simulation analysis results are shown in Figure 5 As shown in (a), the local GC content of the coding scheme of the present invention is always maintained between 45% and 55%, while the local GC content of DNA Fountain varies greatly, with a maximum of 67% and a minimum of 20%. In longer DNA sequences, the local GC content has an important impact on the reliability of the overall coding. Therefore, while paying attention to the overall GC content, we should also pay attention to the changes in the local GC content to ensure the stability and reliability of the coding.

[0099] According to Organick's research, the probability of base errors depends not only on the base itself, but also on the types of adjacent nucleotides. Figure 5 As shown in (b), the 12 3-mers base combinations with the highest error rates are "GTC": 0.9%, "TGC": 0.85%, "CGC": 0.8%, "GTG": 0.8%, "GAC": 0.68%, "GCG": 0.67%, "AGA": 0.65%, "ATA": 0.67%, "CAC": 0.6%, "ACT": 0.55%, "TCT": 0.52%, "TAT": 0.44%. These 3-mers base combinations with high error rates are called undesirable base sequences. Their appearance will increase the difficulty of DNA synthesis and sequencing, and affect the stability and reliability of the overall DNA coding. Therefore, these undesirable base sequences should be avoided as much as possible during the coding process.

[0100] Table 6 OP forms and error weights corresponding to 12 high error rate 3-mers

[0101]

[0102] In the coding based on specific base sequences, the generation of base sequences depends first on the type of specific base sequences used, so the specific base sequence is the main factor affecting the probability of the occurrence of unexpected base sequences. Table 6 shows the results of converting these 12 3-mers base combinations with high error rates into the "OP" form, of which 4 3-mers can be converted to "OPO". When "OPO" appears in the specific base sequence, there are The probability of being a 3-mers base with a high error rate. Therefore, the five forms of "OPO", "POP", "OOO", "PPP" and "POO" should be avoided or reduced when selecting a specific base sequence. In order to quantify the impact of these five high error rate OP forms in a specific base sequence on the base sequence, the probability of each OP form 3-mer combination with a high error rate is set as the error weight w(s i ). Table 6 lists the error weights of “OPO”, “POP”, “OOO”, “PPP”, and “POO”, while the error weights of “PPO”, “OOP”, and “OPP” are 0. Suppose an OP form sequence of length L is S OP , and the error weight calculation formula is as follows:

[0103]

[0104] Among them, s i YesS OP The i-th 3-mers in w(s i ) is the error weight corresponding to the i-th 3-mers.

[0105] According to formula (4), the error weights of all specific base sequences that meet the coding constraints are calculated when the lengths are 2, 4, and 8. Since the calculation of the sequence error weight is related to its length L, the lengths of all specific base sequences are extended and unified to 10 for the convenience of comparison. Figure 5As shown in (d), when the length of the specific base sequence is l=2, the error weight of the sequence is the highest, reaching 3.5, so the probability of error in the sequence obtained by its encoding is the highest. When l=4, the error weight of all sequences is 0.25, so the encoding using a specific base sequence of length 4 has stronger robustness. When l=8, the error weights of the 56 specific base sequences are distributed differently, with the highest being 1.75 and the lowest being 0.75. The 8 specific base sequences with the lowest error weights are: "OOOPPOPP", "OOPPOPPO", "OPPOOOPP", "OPPOPPOO", "POOOPPOP", "POPPOOOP", "PPOOOPPO", and "PPOPPOOO". Therefore, in the encoding process, specific base sequences with low error weights should be given priority to improve the robustness of the encoding.

[0106] There are 48 "0-1" mapping rules proposed in the present invention, and the base sequences generated by each rule are different. In order to reduce the generation of unexpected base sequences, in addition to selecting specific base sequences with low error weights, it is also necessary to select appropriate "0-1" mapping rules. Before encoding the data to be stored, all mapping rules are used to encode the data, the frequency of occurrence of unexpected base sequences in the obtained base sequence is counted, and the mapping rule with the smallest frequency is selected to improve the stability of DNA coding and reduce the probability of errors. In order to verify this idea, a binary sequence is encoded using 48 mapping rules respectively, and the number of unexpected base sequences in each base sequence is counted. The results are as follows: Figure 5 As shown in (c), mapping rule 45 has the largest number of undesired base sequences, which is 407, while mapping rule 23 has the smallest number of undesired base sequences, which is 257. The simulation results show that by screening from 48 mapping rules, the frequency of undesired base sequences can be effectively reduced, thereby improving the reliability of coding.

[0107] The length l of the specific base sequence can be any even number. In order to verify the influence of specific base sequences of different lengths on the error correction capability, the present invention selects low error weight specific base sequences of lengths of 2, 4, 8, 12 and 16, encodes the same batch of binary sequences, and introduces 1%-49% mixed errors. The error correction scheme of the present invention is used for data recovery, and the number of sequencing is 5 times. The recovery rate calculation is based on the premise that the rule index sequence is completely correct, and all results are taken as the average of 50 experiments. The simulation results are shown in Figure 6As shown in the figure, when l is 2, the recovery rate is still above 80% when the error rate is 7%, but after the error rate reaches 10%, the recovery rate drops to 70%, and finally maintains at 50%. When l is 4, the error correction capability is significantly improved. The error correction capability is almost the same when the length is 8, 12, and 16. When the error rate reaches 25%, the recovery rate is still 90%. Therefore, as the length l increases, the error correction capability also increases. However, when the length l exceeds 8, the error correction capability remains at the same level. Considering the impact of the larger length l on the local GC content and homopolymer length of the code, it is believed that l=8 is the best choice, which ensures local stability without losing error correction capability.

[0108] The performance of the error correction scheme proposed in the present invention is not only related to the length of the specific base sequence, but also affected by the number of sequencing reads. In order to explore the relationship between the recovery rate and the number of reads, the present invention encodes a passage of text from "The Wandering Earth". The specific base sequence length is set to 8, and the specific base sequence is "OOPPOPPO", and a base sequence with a length of 10080nt is obtained. Various types of errors are randomly introduced into the base sequence, with an error rate of 1%-45%. Sequencing and error correction are simulated, and the number of sequencing times is 5-100. The simulation results are as follows Figure 7 As shown in the figure, within the error rate range of 1%-20%, the data can be fully recovered after at least five sequencings. As the error rate increases, the number of sequencings required to fully recover the data also increases. When the error rate reaches 40%, 100 sequencings are required to fully recover the data. The simulation results show that the error correction scheme of the present invention is extremely robust to non-continuous errors, especially when dealing with mixed errors.

[0109] When continuous errors occur in a sequence, especially insertions and deletions, the length of the sequence will change. The longer the length of the continuous errors, the greater the sequence change, and the difficulty of completely recovering the data will increase accordingly. The quaternary code error correction scheme proposed by Khuat TH et al. can correct up to 2 insertion or deletion errors. The present invention encodes a segment of characters in "The Wandering Earth". The specific base sequence length l is 8 and 16 respectively. Error correction simulation is performed when the error rate is 20% and 40%. The error form introduced is continuous. The simulation results are shown in Figure 2. Figure 8 As shown in (a) and (b), when the error rate is constant, the longer the continuous length of the error, the lower the data recovery rate. Increasing the length of the specific base sequence can improve the data recovery ability, but the effect is limited. In contrast, increasing the number of sequencing times is more effective in improving the data recovery ability. The simulation results show that even when the number of consecutive errors reaches 8 and the error rate is 40%, by increasing the length of the specific base sequence and the number of sequencing times, the data recovery rate can still be maintained at around 95%.

[0110] The efficiency of multivariate Huffman rotation coding is not only affected by the value of M, but also closely related to the value of the block length t. For a fixed M value, the choice of block length t directly affects the coding efficiency, so it is necessary to determine an appropriate t value to achieve efficient coding.

[0111] During the encoding process, the binary sequence is converted into a decimal sequence in groups of t bits. In the process of constructing the Huffman tree, there are 2 source symbols. t For a fixed M value, the present invention calculates the Huffman coding density for two jpg files with sizes of 193kb and 769kb and two txt files with sizes of 755kb and 957kb. The coding density of the four files varies with the value of the packet length t as shown in the following figure: Fig. 9 As shown in the figure, when the value of the packet length t is less than 18, the coding density reaches the highest value at t = 11. When the packet length t is greater than 18, the coding density increases with the increase of t, but at this time, too large a t value will cause a surge in the number of source symbols, which will cause the depth of the Huffman tree to increase sharply and the coding complexity to rise accordingly. Therefore, in the multivariate Huffman rotation coding, the pursuit of efficient coding is bound to be accompanied by high complexity. In the actual coding process, the appropriate t value should be selected according to the specific situation to balance the coding efficiency and coding complexity.

[0112] Table 7 lists the comparison between the technical solution of the present invention and the previous four classic solutions in terms of error correction algorithm, coding performance and error correction capability. The following is a specific analysis and discussion. Among them, the comparison method Zan comes from the literature [Zan, X., Yao, X., Xu, P., Chen, Z., Xie, L., Li, S., & Liu, W A hierarchical error correction strategy for text DNA storage. Interdisciplinary Sciences: Computational Life Sciences. 14(1), 141-150(2022).], the comparison method Press comes from the literature [Press, WH, Hawkins, JA, Jones Jr, SK, Schaub, JM, & Finkelstein, IJ HEDGES error-correcting code for DNA storage corrects indels and allows sequence constraints. Proceedings of the National Academy of Sciences. 117(31), 18489-18496(2020).], and the comparison method Song comes from the literature [Song, L., Geng, F., Gong, ZY. et al. Robust data storage in DNA by deBruijn graph-based de novo strand assembly.Nature Communications.13,5361(2022).], and the comparison method Xue comes from the literature [Xue, TB; Lau, FCM Construction of GC-Balanced DNA With Deletion / Insertion / Mutation Error Correction for DNA Storage System.IEEE Access.8,140972-140980(2020).].

[0113] Table 7 Comparison with previous methods

[0114]

[0115] The coding scheme of the present invention uses a specific base sequence with minimum error weight and multivariate rotation coding to make the DNA sequence fragment not only meet the GC content and homopolymer length constraints, but also avoid the generation of unexpected base sequences, thereby reducing the error rate in the DNA synthesis and sequencing process. In contrast, the other schemes in Table 7 do not have effective measures to avoid unexpected base sequences. The overall storage density of the present invention is about 1.33 bit / nt, which is higher than other schemes. The scheme of Song and Press introduced redundant RS codes [15, 16], while the present invention combines rotation coding with M-ary Huffman and introduces a compression algorithm to improve the coding density while increasing the coding time complexity. In the future, the regular index sequence S can be directly converted into a random number table. rules Convert to 48 yuan and utilize rotation coding mapping, thereby shorten the encoding and decoding time, although the coding density can be slightly reduced at this time, so need to select suitable mapping scheme according to specific circumstances.In the process of combining multivariate rotation coding with M-ary Huffman coding, the present invention only calculates the situation that the codeword length is 3 in the rotation coding table, when the codeword length is greater than 3, although the constraint of GC content and homopolymer length can also be maintained, the coding complexity increases, and the coding effect is not improved, therefore the situation that the codeword is greater than 3 is not considered.Although compression algorithms such as Lempel-Ziv-Welch, Gzip or run-length coding can be introduced before encoding, to further improve the coding density.But these compression algorithms are not specially designed for DNA storage, and compression may change the original information structure, causing the partial error fragment to affect the complete recovery of the overall data during data recovery.

[0116] The error correction scheme of the present invention can completely recover 40% of mixed errors when the number of sequencing times reaches 100 times, while the error correction capabilities of the schemes of Zan, Press, Song, Xue and others are 5%, 10%, 10% and 0.3% respectively. Xue established a single error correction system using improved Levenshtein code, but there are more errors under large-scale data storage. Song used the "multi-copy" feature to perform error-free reconstruction of chain sequences, but when errors increase, the chain reconstruction rate drops significantly. The HEDGES algorithm proposed by Press combines RS code to correct errors through the correlation between bits in the binary sequence, but the recovery rate drops significantly when the error rate is greater than 10%. Zan's three-layer error correction algorithm combines the six-base coding table with the Hamming distance, which can efficiently recover data at low error rates, but the recovery rate is only 60% after the error rate reaches 10%, and this algorithm is only applicable to English text, which has certain limitations.

[0117] In the scheme of the present invention, if the same type of base substitution occurs, for example, when the P-type base T mutates to A, the position where the substitution error occurs cannot be detected based on the specific base sequence. However, multiple sequencing compensates for the weak sensitivity of the scheme to substitution errors and enhances the overall error correction capability of the scheme.

[0118] The present invention proposes a coding scheme based on a specific base sequence, which fully considers the overall and local stability of the DNA sequence. The present invention uses 48 "0-1" mapping rules to encode any binary data into a specific form of DNA sequence, meeting the requirement that the local GC content is between 45% and 55%, and reduces the frequency of occurrence of undesirable base sequences by screening 48 mapping rules, thereby ensuring the reliability and diversity of storage. For the regular index sequence generated by the specific base sequence encoding, a rotation coding table with a codeword length of three and a size of 48×48 is constructed, combining the compression effect of rotation coding with Huffman coding, which not only improves the mapping potential, but also avoids the generation of long homopolymers, so that the maximum homopolymer length is not greater than 3.

[0119] With the development of DNA storage towards large-scale data storage, it has become a key issue to quickly and accurately recover DNA sequences with a large number of insertion, deletion, and substitution errors. The encoding scheme proposed in the present invention maps data into a base sequence of a specific pattern, which helps to detect the location where the error occurs and completely recover the data through multiple sequencing. The simulation results show that in a non-continuous error mode, when the error rate is below 25%, more than 90% of the data can be recovered through only five sequencings, and more errors can be tolerated as the number of sequencing times increases. Although continuous errors bring difficulties to data recovery, by increasing the length of specific base sequences or the number of sequencing times, complete decoding of data can also be guaranteed under high error rates and long continuous errors. The overall storage density of the scheme of the present invention is 1.33bit / nt, which sacrifices a certain storage density but obtains a high degree of robustness, providing a new idea for error correction in large-scale DNA data storage.

[0120] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A highly robust DNA storage method based on a specific base sequence, characterized in that: The steps are as follows: Step 1: Construct "0-1" mapping rules and specific base sequences; Step 2: Use the "0-1" mapping rule and the specific base sequence to perform DNA encoding on the binary sequence to be encoded, and generate a rule index sequence and a data DNA sequence; Step 3: compress the regular index sequence using multivariate Huffman coding, and encode the compressed sequence in combination with multivariate rotation coding to generate a DNA coding sequence; Step 4: Sequence the data DNA sequence obtained in step 2 multiple times and compare it with the specific base sequence to detect and correct errors.

2. The highly robust DNA storage method based on a specific base sequence according to claim 1, characterized in that: The "0-1" mapping rule consists of the "0" rule and the "1" rule. The "0" rule and the "1" rule respectively map each element in the binary set B = {00, 01, 10, 11} to a base in the DNA base set N = {A, C, G, T}, and meet two conditions: ① The "0" rule and the "1" rule are one-to-one mappings, that is, N 01 ≠N 02 ≠N 03 ≠N 04 , N 11 ≠N 12 ≠N 13 ≠N 14 ; Among them, N 01 、N 02 、N 03 、N 04 For the elements in the "0" rule, N 11 、N 12 、N 13 、N 14 is an element in the "1" rule; ② The bases A / T and G / C are mutually exclusive; According to the "0-1" mapping rule, bases A and T are regarded as one type of base, represented by P; bases G and C are regarded as another type of base, represented by O. The specific base sequence S mould P and O each account for half, and no more than three P or O can appear consecutively.

3. The highly robust DNA storage method based on a specific base sequence according to claim 2, characterized in that: The mathematical description of the "0-1" mapping rule is Where i = 0, 1, j = 1, 2, 3, 4; There are 48 "0-1" mapping rules that meet the conditions: The specific base sequence is a base sequence that appears periodically and repeatedly in a DNA storage sequence, and is composed of multiple identical sequences that satisfy specific constraints, wherein the specific constraints include a base combination pattern, a sequence repetition period, and a sequence structural stability; through any "0-1" mapping rule, a binary sequence is mapped into a base sequence composed of periodically repeated specific base sequences, thereby achieving efficient encoding and storage of data; The mathematical description of a specific basis sequence of length l is: Among them, N1, N2, ..., N l-1 ,N l Respectively represent l elements of a specific base sequence; half of them are A or T, and the other half are C or G; The length l of the specific base sequence is an even number.

4. The highly robust DNA storage method based on a specific base sequence according to any one of claims 1 to 3, characterized in that: The binary sequence to be encoded is C = (c1 c2 … c i-1 c i … c n-1 c n ), the i-th binary element c i ∈{0,1}, n is the length of the binary sequence C. Using the "0-1" mapping rule and a specific base sequence for DNA encoding is to map the binary sequence C to a specific base sequence that repeats and recurs. The implementation method is: Step 1. Grouping: Group every two bits in the binary sequence C into groups and select a specific base sequence of a given length l; Step 2. Determine the base: Map the first group c1 and c2 of the binary sequence C to the first bit of the specific base sequence, determine the corresponding base according to the "0-1" mapping rule constraint, and record the index number of the "0-1" mapping rule; and so on, the dth group of the binary sequence C is mapped to the dth bit of the specific base sequence, determine the corresponding base of the mapping, and d <l; Step 3. Cyclic mapping: The l+1th group of the binary sequence C is mapped to the first bit of the specific base sequence. When the mapping is successful, the index number of the "0-1" mapping rule is recorded. Similarly, every two bits in the binary sequence C are mapped and encoded to obtain the complete encoding result as the data DNA sequence, and the rule index sequence S is generated. rules .

5. The highly robust DNA storage method based on specific base sequences according to claim 4, characterized in that: The multivariate Huffman coding is an M-ary Huffman coding, which compresses the regular index sequence into an M-ary sequence through the M-ary Huffman coding, and the multivariate rotation coding uses a rotation coding table to map the M-ary sequence into a base sequence that meets the constraints.

6. The highly robust DNA storage method based on specific base sequences according to claim 5, characterized in that: The method of constructing the M-ary Huffman tree by M-ary Huffman coding is: rules Convert each t-bit binary number into decimal data and count the probability of each decimal data. The decimal data is called the source symbol. Assuming there are f source symbols, arrange them in order according to their probabilities. Give M source symbols with the smallest probabilities (a f-M-1 ,a f-M-2 ,…a f ) are assigned code words (0, 1, 2, ... M-1) respectively and the M source symbols are merged into a new symbol, and a new source with a length of (fM-1) is obtained, which is called the first source reduction; then the reduced new source is arranged again according to the probability from small to large and the source reduction process is repeated to obtain the M-ary Huffman tree with the shortest average code length; The multi-element rotation coding is: for the obtained M-ary sequence, set the 0th code to be any one of the three bases, and look up the rotation coding table according to the 0th three bases to obtain the three bases mapped by the decimal number to be encoded in the first M-ary sequence, and look up the rotation coding table according to the first three bases to obtain the three bases mapped by the decimal number to be encoded in the second M-ary sequence; and so on, after rotating and mapping the M-ary sequence one by one, finally obtain the base sequence S2.

7. The highly robust DNA storage method based on specific base sequences according to claim 6, characterized in that: The method for constructing the rotational coding table is as follows: under the condition of satisfying the homopolymer constraint and the GC content constraint, the base combination length in the rotational coding table is expanded to 3; there are 64 three-base combinations, and in order to satisfy the constraint that the maximum homopolymer length is less than 4, the combinations in which the first two bases are the same are removed, and the remaining 48 are used as rotational coding codewords, and the 48 three-base combinations are represented as rotational coding codewords H: (H1, H2, ... H 47 ,H 48 ), and the rotation codeword is The positions of the 48 tri-bases are cyclically shifted to form a rotation coding matrix of size 48×48, thereby obtaining a rotation coding table. The 48-element rotation coding table is:

8. The highly robust DNA storage method based on a specific base sequence according to any one of claims 5 to 7, characterized in that: The error correction method is as follows: first, the obtained data DNA sequence is sequenced multiple times to generate multiple read lengths, each read length is converted into an OP representation, compared with the base sequence of a specific base sequence and aligned, so that the lengths of all read lengths are consistent; secondly, the weights of the bases at the corresponding positions of all read lengths are calculated, and the best base sequence to be decoded is obtained according to the weights; the corresponding decoding rule is selected from the obtained rule index sequence in combination with the value of its corresponding position to restore the base sequence to be decoded, and finally the original binary sequence is obtained.

9. The highly robust DNA storage method based on specific base sequences according to claim 8, characterized in that: In the process of data recovery, the base sequence S2 generated by the multi-element rotation coding is divided into groups of three. When encoding, the default 0th three bases are H1. When decoding the first three bases, the current rotation coding table is searched according to the previous three bases as H1. The value obtained is the code word corresponding to the first three bases. This process is iterated until the entire base sequence S2 is traversed; the code word sequence is decoded into a decimal sequence according to the M-element Huffman tree, and then the regular index sequence S is obtained. rules ; With the help of rule index sequence S rules The method for decoding the base sequence S1 obtained by RS decoding is as follows: when decoding the bases in the base sequence S1 bit by bit, refer to the rule index sequence S rules The value at the corresponding position selects the corresponding decoding rule. Specifically, when decoding the i-th base in the base sequence S1, the rule index sequence S rules The i-th bit value of is decoded by selecting the "0" rule or the "1" rule.

10. The highly robust DNA storage method based on specific base sequences according to claim 9, characterized in that: The average code length is the number of bits required for each codeword in the sequence generated by the M-ary Huffman coding, and M is selected as 48 according to the average code length; In the encoding process, a specific base sequence with low error weight should be selected first; the length l of the specific base sequence is 8; Before encoding the data to be stored, all mapping rules in the "0-1" mapping rule are used to encode the data to be stored, the occurrence frequency of the non-expected base sequence in the obtained base sequence is counted, and the mapping rule with the smallest frequency is selected.

Images