DNA information storage method based on improved encoding table and HMSA

Abstract

The application provides a DNA information storage method based on an improved encoding table and HMSA, and steps are as follows: inspired by Meinolf Blawat, an improved mapping rule is used to construct an encoding table corresponding to ASCII code one by one; information to be stored is encoded into a base sequence according to the constructed encoding table; error correction is carried out by using HMSA; the length of reads and design sequence is compared, and three groups are divided, and error correction is carried out on each group; the minimum Hamming distance between each group of base fragments and the base fragments in the encoding table is calculated, the first group of base fragments with the minimum Hamming distance not being 0 is found, and preliminary correction is carried out; after the base sequence after preliminary correction is aligned and placed together, multi-sequence alignment method is used to carry out error correction again; and the corrected reads are divided every 5 base slices and decoded into stored information by using the encoding table. The application can not only effectively correct substitution errors in the DNA sequence, but also correct part of insertion and deletion errors in the DNA sequence.

Description

Technical Field

[0001] This invention relates to the technical field of DNA information storage, and more particularly to a DNA information storage method based on an improved coding table and HMSA. Background Technology

[0002] Information is the most important means of cultural transmission. Throughout human history, the advent of writing greatly promoted the progress of human civilization. However, we are now in an era of explosive growth in digital information. Based on current trends, it is estimated that by 2025, the global annual information generation will reach 175 ZB, equivalent to 491 EB of data generated daily. Faced with such a massive amount of information, traditional storage media such as optical discs, floppy disks, and USB flash drives cannot guarantee storage capacity and durability. Due to the exponential growth of digital information, silicon-based storage methods will not be able to meet global needs in the near future. DNA information storage stands out due to its potential advantages, including large storage capacity, environmental friendliness, long lifespan, high stability, low energy consumption, and low data maintenance costs.

[0003] DNA stands out among various storage media as a natural storage medium due to its numerous advantages. These include: 1. High stability and long lifespan: Under specific storage conditions, DNA has a half-life of approximately 521 years. If stored in silica, DNA can be preserved for over 2 million years. 2. High DNA storage density and large capacity: Theoretically, the capacity for storing digital information in DNA is 455 EB bytes per gram of single-stranded DNA (ssDNA). DNA information storage technology refers to the use of artificially synthesized adenine (A), guanine (G), thymine (T), and cytosine (C) bases for digital information storage. Due to the rapid development of DNA synthesis and amplification technologies (PCR) and sequencing technologies (Sanger sequencing, next-generation sequencing, third-generation sequencing, single-cell sequencing) in recent years, large-scale DNA information storage has become a groundbreaking storage technology in today's era of data explosion.

[0004] The concept of DNA information storage was first proposed by physicist Richard P. Feynman in 1959. In his book "There's plenty of room at the bottom (data storage): An Invitation to Enter a New Field of Physics," he described the future construction of objects (including organisms) similar to the microscopic world with the ability to store information. The process of storing and retrieving information from DNA mainly includes five parts: (1) Encoding part, which encodes the stored information into a DNA sequence. (2) Storage part, which synthesizes DNA information and stores it in a DNA pool. (3) Amplification part, which uses polymerase chain reaction (PCR) to amplify the target DNA. (4) Sequencing part, which uses sequencing technology to read the amplified DNA sequence information. (5) Decoding part, which decodes the DNA sequence information. The general flowchart is as follows: Figure 1 As shown.

[0005] The basic encoding rules used in the early stages (i.e., [00) A;01 C;10 G;11 During the synthesis and sequencing of DNA sequences, a large number of erroneous base sequences are generated. For example, Antkowiak et al. found that over 99.9% of the 15 million base sequences sequenced contained erroneous bases. Organick et al. found that over 88% of reads in third-generation sequencing methods contained insertion and deletion errors during synthesis and sequencing. Furthermore, large fragment losses can also occur during synthesis, amplification, and sequencing. Studies have shown that GC content and homopolymer length in DNA sequences are the main causes of these errors. Therefore, to reduce errors, most previous studies have introduced constraints when designing coding schemes, such as GC content between 40% and 60%; homopolymer length not exceeding 4; and minimizing the occurrence of secondary structures in DNA sequences. As early as 2012, Church et al. proposed a coding scheme that encodes 1 base per bit (A or C for 0, G or T for 1), encoding information in multiple ways to avoid extreme GC content, homopolymer quantity, and secondary structures. In 2013, Goldman et al. proposed a ternary rotation coding scheme, which cleverly removes homopolymers from DNA sequences by compressing information into ternary information and then using rotation coding. In 2015, Grass et al. combined the Galois domain of GF(47) with triplet bases, and proposed an encoding / decoding algorithm that avoids the problem of single-base repetitions longer than 3. Simultaneously, in terms of error correction, Grass et al. introduced error-correcting codes for the first time, adding horizontal and vertical Reed-Solomon codes to correct internal and inter-sequence errors in DNA sequences, respectively. In 2016, Blawat et al. proposed a forward error-correcting coding scheme that uses a dual-rule coding method to alternately encode bases at different positions. Finally, different selection rules are used to select the encoded base fragments to satisfy homopolymer constraints. In 2017, Erlich et al. proposed a DNA Fountain algorithm based on Luby Transform codes, which filters out DNA sequences that do not meet the constraints through a selection mechanism. Simultaneously, the error correction mechanism RS (Reed-Solomon) code was introduced. However, due to fundamental problems with the Luby Transform, fountain codes are at risk of decoding failure when processing specific binary sequences. In 2022, Zhi et al. proposed a Yin-Yang code based on the Yin-Yang philosophy of Chinese Taoism. This encoding method uses two rules to encode two binary digits into a nucleotide, constructing screening rules based on GC content, homopolymerity, and free energy to select DNA sequences that meet the constraints.

[0006] Reducing the occurrence of errors is only one aspect. Due to the current level of technology, errors are inevitable. Therefore, error correction mechanisms need to be added after the DNA sequence. Existing error correction mechanisms can be roughly divided into two categories in terms of redundancy. The first category is error correction methods that do not increase redundancy: most of them are multiple sequence alignment (MSA) or its extensions. These methods often require a large number of reads. Since there is a synchronization problem between sequences when insertion and deletion errors occur, this method also needs a high-performance clustering algorithm to solve this problem. For example, Xie and Philipp et al. used a clustering algorithm to align DNA sequences. Another example is that in 2021, Zan et al. improved on the idea of ​​MSA and set up a three-layer error correction method for text information. The second category of methods that increase redundancy can be divided into: (1) increasing physical redundancy, such as Goldman et al. using DNA sequences with 75% overlap for error correction, which increases redundancy by four times, and Bornholt et al. storing two DNA sequences after XORing, which increases redundancy by two times. (2) increasing logical redundancy, such as traditional error correction codes Hamming codes, RS codes, and LDPC codes. Traditional error-correcting codes can correct substitution errors in DNA sequences, but they are ineffective against insertion and deletion errors. Sequence insertions and deletions can lead to inconsistent DNA sequence lengths. When using RS codes for error correction, incorrectly long DNA sequences are typically removed. However, this approach significantly reduces the utilization rate of sequencing reads, resulting in substantial resource waste. Summary of the Invention

[0007] To standardize DNA encoding systems and ensure consistent character representation and transmission across different computers and communication devices, this invention, inspired by Brawat, proposes a DNA information storage method based on an improved encoding table and HMSA. This invention creates an encoding table that satisfies GC content and homopolymer constraints, with a one-to-one correspondence between the encoding table and ASCII codes, and an encoding density of 1.4 or 1.6 bits / nt. It addresses the issue of low read utilization in existing DNA information storage error correction methods. Furthermore, the HMSA error correction method is applicable to encoding schemes with all lookup table types.

[0008] To achieve the above objectives, the technical solution of the present invention is as follows: a DNA information storage method based on an improved coding table and HMSA, comprising the following steps:

[0009] Step 1: Improve the mapping rules of Meinolf Blawat and add filtering rules to construct an encoding table that corresponds one-to-one with ASCII codes;

[0010] Step 2: Encode the information to be stored into a base sequence according to the constructed encoding table;

[0011] Step 3: Error correction using HMSA: Compare the lengths of the Reads with the design sequence, divide the Reads into three groups, and correct each group separately: Slice the Reads into groups of 5 bases each, calculate the minimum Hamming distance between each group of base segments and the base segments in the coding table, find the first group of base segments with a non-zero minimum Hamming distance, and perform preliminary correction; After aligning the preliminarily corrected base sequences together, perform error correction again using multiple sequence alignment.

[0012] Step 4: Use an encoding table to divide the corrected Reads into 5-base slices and decode them into stored information.

[0013] Preferably, the rules for constructing the encoding table include rule A and rule B, wherein:

[0014] Rule A: {'00' 'C', '01' 'G', '10' 'T', '11' 'A'};

[0015] Rule B: {'00': ['AA', 'CC', 'GG', 'TT'], '01' ['AC', 'CG', 'GT', 'TA'], '10' ['AG', 'CT', 'GA', 'TC'], '11' ['AT', 'CA', 'GC', 'TG']}.

[0016] Preferably, when constructing the encoding table, {00000000, 00000001, 00000010, ..., 01111111} is used as input, and 128 8-bit binary numbers are encoded into their corresponding bases. The method is as follows: the first 6 bits of each 8-bit binary number are grouped in pairs, and the 7th and 8th bits are each grouped into two separate groups. Then, a set is used... To represent; use rule A to represent the sequence The binary information in the code is encoded into base information to obtain the bases. Use rule B to sequence The binary information in the code is encoded into base information to obtain the bases. The encoded base fragment is ,in, And use three filtering criteria to filter.

[0017] Preferably, the three screening conditions are:

[0018] Condition 1: Check if there are any three consecutive identical bases or the last two identical bases in the four base combinations after the encoding. If so, remove the base combination; if not, proceed to Condition 2.

[0019] Condition 2: Check if the GC content of the 5-position base combination is 40% or 60%. If so, put the result into the coding table of the corresponding order and proceed to Condition 3; otherwise, keep the result and proceed to Condition 3.

[0020] Condition 3: Check if the encoding table is empty. If it is empty, put the result of Condition 2 into the corresponding encoding table; if it is not empty, output the encoding table.

[0021] Preferably, the encoding table is as follows:

[0022] .

[0023] Preferably, the preliminary correction method is as follows: compare the length of Reads with the designed base length, and divide Reads into replacement error group, insertion error group and deletion error group; slice and group Reads in the replacement error group, insertion error group and deletion error group one by one, calculate the minimum Hamming distance between each group of base fragments and the base fragments in the encoding table, find the first base fragment group with a minimum Hamming distance of not 0, and correct the replacement error group, insertion error group and deletion error group respectively.

[0024] Preferably, when the length of the base sequence is equal to the designed base length of 150 nt and an error occurs during decoding, the base sequence is placed into the replacement error group;

[0025] The initial method for correcting and replacing the incorrect base sequences is as follows: the base sequence is sliced ​​into groups of 5 bases each, and each group of sliced ​​base fragments is then processed using... To represent, the base fragments in the coding table are represented by To represent; use lists Represents a base fragment Base fragments and coding tables A list of minimum Hamming distances between base segments in the nucleotide sequence;

[0026] Select the minimum Hamming distance Non-zero base fragments , base fragment Flip to calculate the minimum Hamming distance The corresponding base fragment j, until the list All values ​​in the value are 0, indicating a base fragment. All base fragments All in the encoding table middle;

[0027] Base fragments grouped by slice The merged DNA sequences yielded a preliminary corrected sequence.

[0028] Preferably, when the length of the Reads is greater than the designed DNA sequence length, i.e., 150 nt, and an error occurs during decoding, the Reads are placed in the insertion error group;

[0029] The initial method for correcting DNA sequences in groups with insertion errors is as follows: The DNA sequence is sliced ​​and grouped into segments of 5 bases each, with the last segment having fewer than 5 bases temporarily counted as a single group; each segment is then processed using a setter. Represents a set Base fragments and coding tables Minimum Hamming distance between medium base fragments using a list This indicates that, when calculating the minimum Hamming distance, if the last set of bases... Fewer than 5 people, making Hamming's distance Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of;

[0030] From list Find the first non-zero Hamming distance The corresponding base fragment Started correcting errors and simultaneously found a way to calculate Hamming distance. At that time, the encoding table The corresponding base fragment Delete base fragments Neutral base fragments Different first bases; base fragments All base fragments Merge, then re-slice the merged base sequence into 5-position groups, repeating the above process to begin cyclic error correction until the following conditions are met: list All values ​​are 0, the length of the corrected base sequence is equal to 150 nt, and the base sequence correction is initially completed.

[0031] Preferably, when the length of the Reads is less than the designed base length of 150 nt and an error occurs during decoding, the Reads are placed in the error deletion group;

[0032] The method for initially correcting the base sequence in the insertion error group is as follows: First, the bases in the base sequence are grouped into slices of 5 positions each, with the last group having fewer than 5 positions counted as a single group; each group of base segments after slicing is used as a set... This indicates that the base fragments in the encoding table are... This indicates that after grouping the base sequence slices that have deletion errors, the minimum Hamming distance calculated for the deletion error groups is obtained using a list. This indicates that, when calculating the minimum Hamming distance, if the last set of bases... If the number of digits is less than 5, then the Hamming distance is... Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of;

[0033] From list Find the first non-zero Hamming distance The corresponding base fragment Start error correction and find the corresponding encoding table. medium base fragment , j, in the base fragment Inserting a base into the base fragment makes the base fragment The 5-base fragment containing the inserted base and the base fragment Hamming distance Decrease or equal to 0; reduce the set All the base fragments Merge, regroup the merged base sequences into 5-position slices, and then extract them from the list again. Find the first non-zero Hamming distance The corresponding base fragment Start correcting errors until the list is satisfied. When all values ​​are 0 and the length of the corrected base sequence is equal to 150 nt, the initial base sequence correction is complete. If the above conditions are not met, the initial correction is complete. If all values ​​are 0 and the length of the corrected base sequence is equal to 150 nt, then proceed to the next base sequence, and so on, until all base sequences have been processed.

[0034] Preferably, all DNA sequences initially corrected from the replacement error group, insertion error group, and deletion error group are put together, and multiple sequence alignment is used to align all corrected reads, selecting the base with the largest proportion of bases at each position.

[0035] After correction, the corrected Reads are sliced ​​into 5-base segments using an encoding table and then decoded into stored information.

[0036] When decoding a base sequence, base fragments not in the encoding table are decoded as the fixed 'GTAGC'.

[0037] Compared with existing technologies, the beneficial effects of this invention are as follows: The encoding table designed in this invention not only satisfies the GC content constraint and the maximum homopolymer constraint, but also achieves an encoding density of 1.4 bits / nt. In text information encoding and decoding, the encoding density can be increased to 1.6 bits / nt, directly skipping the binary conversion process and directly encoding the data into A / T / C / G sequences, reducing the complexity of encoding and decoding. Regarding error correction, an improved first-type error correction method that does not increase redundancy is chosen, namely, an improvement to MSA. This invention combines the minimum Hamming distance and multiple sequence alignment methods to propose a novel error correction method that can correct three types of errors: substitution, insertion, and deletion, greatly improving the utilization rate of sequencing reads and avoiding resource waste. It can solve the sequence synchronization problem and avoids the waste of a large number of sequencing reads without using clustering algorithms. Simulation results show that the correction method of the present invention can not only effectively correct substitution errors in DNA sequences, but also correct partial insertion and deletion errors in DNA sequences. When processing reads with a maximum length deviation within ±4 compared to the correct length, regardless of whether the three types of errors occur consecutively, the present invention only requires 10 reads to achieve a 100% read recovery rate. Attached Figure Description

[0038] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0039] Figure 1 A flowchart for storing existing DNA information.

[0040] Figure 2 This is a flowchart illustrating the construction and encoding of the encoding table for this invention.

[0041] Figure 3 This is a schematic diagram of the correction methods for three error types of the present invention, wherein (a) is replacing the error group, (b) is inserting the error group, and (c) is deleting the error group.

[0042] Figure 4 The DNA sequence simulated for this invention.

[0043] Figure 5 This is a schematic diagram of the GC content and homopolymer size of the DNA sequence of the present invention, wherein (a) is the GC content and (b) is the homopolymer size.

[0044] Figure 6The diagram shows a comparison of the error correction performance of the present invention. In the diagram, (a) shows the simulation results of processing erroneous reads using the MSA method under different simulated sequencing runs, changing the error rate of substitution errors in the reads; (b) shows the simulation results of processing reads with substitution errors using the present invention; (c) shows the simulation results of processing reads with insertion errors using the HMSA method under different simulated sequencing runs, changing the number of reads with insertion errors; and (d) shows the simulation results of processing reads with deletion errors using the method of the present invention under different simulated sequencing runs, changing the number of reads with deletion errors.

[0045] Figure 7 This invention improves the error correction performance of Reads that correct discontinuous errors.

[0046] Figure 8 The error correction performance diagram for the present invention is used to correct consecutive replacement, insertion or deletion errors (Reads), where (a) represents 3 consecutive reads, (b) represents 4 consecutive reads, (c) represents 5 consecutive reads, and (d) represents 6 consecutive reads. Detailed Implementation

[0047] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0048] A DNA information storage method based on an improved coding table and HMSA, comprising the following steps:

[0049] Step 1: Modify the mapping rules of Meinolf Blawat and add filtering rules to construct an encoding table that corresponds one-to-one with ASCII codes.

[0050] Inspired by Meinolf Blawat and incorporating ASCII, an encoding table corresponding one-to-one with ASCII was constructed. This encoding table satisfies both GC content and homopolymer constraints, and, like ASCII, standardizes the encoding system, ensuring consistency in character representation and transmission across different computers and communication devices. The encoding table was constructed using mapping rules and encoding methods similar to Meinolf Blawat's. However, to satisfy the GC content constraint of the encoded base sequence, the basic rules were modified. Two mapping rules were given to satisfy the GC content constraint of the encoded base sequence; these are Rules A and Rules B of this invention. Specifically, each binary group in Rule B corresponds to four possible base combinations.

[0051] Rule A: {'00' 'C', '01' 'G', '10' 'T', '11' 'A'}.

[0052] Rule B: {'00': ['AA', 'CC', 'GG', 'TT'], '01' ['AC', 'CG', 'GT', 'TA'], '10' ['AG', 'CT', 'GA', 'TC'], '11' ['AT', 'CA', 'GC', 'TG']}.

[0053] To associate with the 128-bit ASCII code, the stored information is first converted into binary information. Then, the binary information is sliced ​​into 7-bit segments. Finally, a "0" is added to the beginning of each segment, making all the 7-bit segments into 8-bit segments, all within the set {00000000, 00000001, 00000010, ..., 01111111}. This establishes a connection with the ASCII code table. Therefore, when constructing the encoding table, {00000000, 00000001, 00000010, ..., 01111111} is used as input, and these 128 8-bit binary segments are encoded into their corresponding bases. Taking one segment (an 8-bit binary segment) as an example, firstly, the first 6 bits of the 8-bit binary segment are grouped in pairs, and the 7th and 8th bits are each grouped into a separate group. To represent. The encoded base fragment is represented using To represent. Base Is it using rule A to The binary information in the code is encoded into base information. (Bases) Is it using rule B to The binary information is encoded into base information. Since rule B has four options, the encoded 5-bit base... There are therefore four possible combinations. To meet the requirements for homopolymer and GC content, this invention establishes three screening criteria, the specific steps of which are as follows: Figure 2 As shown. The filtering criteria are as follows:

[0054] Condition 1: Check if there are any three consecutive identical bases or the last two identical bases in the four base combinations after encoding. If so, remove the base combination; if not, proceed to Condition 2.

[0055] Condition 2: Check if the GC content of the 5-base combination is 40% or 60%. If yes, put the result into the encoding table of the corresponding order and proceed to Condition 3; if not, keep the result and proceed to Condition 3.

[0056] Condition 3: Check if the encoding table is empty. If it is empty, put the result of Condition 2 into the corresponding encoding table; if it is not empty, output the encoding table.

[0057] After applying three filtering criteria, the results can form two encoding tables. When encoding the same file using these two encoding tables, their GC contents are approximately 45% and 55% respectively, both satisfying the GC content constraint. In subsequent simulation experiments, this invention selected the encoding table with a GC content of approximately 55%, as shown in Table 1.

[0058] Step 2: Encode the information to be stored into a base sequence directly according to the constructed encoding table.

[0059] Generally, when encoding all types of files, the first step is to convert the input file into 8-bit binary. The second step is to segment the binary information into 7-bit segments. If the last segment is less than 7 bits, 0s are added to the end until it reaches 7 bits. The third step is to add a "0" to the beginning of each segment, so that all the segmented binary segments are within {00000000, 00000001, 00000010, ..., 01111111}, that is, within the encoding table. The fourth step is to use the ASCII-DNA encoding table to encode the binary information processed in the third step into bases. The encoding process is as follows: Figure 2 As shown in Table 1, the addition of a "0" before the 7 binary bits is equivalent to adding a redundant bit, resulting in an actual calculated encoding density of 1.39 bits / nt. Since an encoding table directly associated with ASCII codes is constructed, especially for text information, the table can be used to directly encode the text information into a base sequence. For example, the text information 'abc' can be directly encoded as 'GTCCGGTACGGTACT' using the encoding table. This skips the intermediate conversion of binary information, significantly reducing time and space complexity, and also increasing the encoding density. It is equivalent to encoding 8 bits of valid binary information with 5 base bits, resulting in an encoding density of 1.6 bits / base.

[0060] Generally, decoding is the reverse of encoding. For all types of storage files, the unified decoding method is as follows: First, the sequencing reads are sliced ​​into groups of five bases each, and decoded using an encoding table. Second, the resulting binary sequence is sliced ​​into groups of eight bits each, and the leading zeros of each group are removed. Third, all fragments are merged to obtain binary information, and finally, the binary information is converted into a storage file. Specifically, for text files, the binary conversion process can be skipped, and the encoding table can be used directly to decode the sliced ​​(five-base-per-bit) base fragments into text information.

[0061] Table 1 Coding Table

[0062]

[0063] Step 3: Error correction using HMSA: Compare the length of the Reads with the design sequence, divide the Reads into three groups, and correct each group separately: Slice the Reads into groups of 5 bases each, calculate the minimum Hamming distance between each group of base segments and the base segments in the coding table, find the first group of base segments with a non-zero minimum Hamming distance, and perform preliminary correction; Align the base sequences that have completed preliminary correction and put them together, and then use MSA to perform error correction again.

[0064] Step 4: Use the encoding table to decode the base fragments after the corrected Reads are divided into 5-base slices into text information.

[0065] Since the encoding method proposed in this invention is codebook-based, the base segments grouped into 5-bit slices in a correct Read should all be present in the encoding table. In other words, assuming no error occurs in the Read, the minimum Hamming distance calculated between the base segments grouped in the Read and the base segments in the encoding table should all be 0. If the minimum Hamming distance is not 0, it can be determined that an error has occurred in the Read, and the location of the error can be found. Therefore, utilizing this principle, this invention combines minimum Hamming distance and multiple sequence alignment (MSA) to develop a novel error correction method (HMSA).

[0066] The HMSA error correction method of this invention can be roughly summarized into four steps:

[0067] Step 1: By comparing the length of the sequencing reads with the designed base length, the reads are divided into three groups (replacement error group, insertion error group, and deletion error group).

[0068] Step 2: Slice the Reads in the three groups one by one, calculate the minimum Hamming distance between each group of base fragments and the base fragments in the encoding table, and find the first group where the minimum Hamming distance is not 0 (the location where the error occurred).

[0069] Step 3: Initial correction, correcting errors in the replacement, insertion, and deletion groups respectively.

[0070] Step 4: Align all the Reads that have completed the initial correction and put them together, then use MSA to perform error correction again.

[0071] The specific method is as follows:

[0072] Case 1: Correcting base sequences in the substitution error group. Reads in the substitution error group are defined as follows: when the length of a read is equal to the designed base length (150 nt) and an error occurs during decoding, the read is placed in the substitution error group. Since reads in the substitution error group do not have sequence synchronization issues, multiple sequence alignment (MSA) can be directly used. However, sequence alignment may require a large number of sequencing iterations, resulting in higher sequencing costs. Therefore, before sequence alignment, the minimum Hamming distance is calculated to process the reads, i.e., the error correction method HMSA proposed in this invention. This not only corrects some errors in the reads in advance, thereby reducing the number of sequencing iterations required for MSA, but also improves the error correction capability of MSA.

[0073] Taking any base sequence Read in the replacement group as an example, firstly, Read is sliced ​​into groups of 5 bases each. Then, each group of sliced ​​bases is... To represent, the base fragments in the coding table are represented by To indicate the use of a list To represent base fragment groups The base fragments in A list of minimum Hamming distances between base segments in the code. Select the minimum Hamming distance. Non-zero base fragments This involves replacing the position of the incorrect base. The base fragment... Flip to calculate the minimum Hamming distance The corresponding base fragment j. Until the list All values ​​in the value set are 0, that is All base fragments All In the middle. At this point, the slices are grouped. After merging, the DNA sequence correction is initially completed, specifically as follows: Figure 3 As shown in (a).

[0074] Case 2: Correcting DNA sequences in the insertion error group. Reads in the insertion error group are defined as those whose length exceeds the designed DNA sequence length (150 nt) and an error occurs during decoding. It's important to note that readings in the insertion error group may experience sequence synchronization issues when using the MSA error correction method. Specifically, bases at the same position may not align, preventing MSA from correcting erroneous bases. The following is a specific method using this invention to solve the sequence synchronization problem in the insertion error group.

[0075] Taking any base sequence Read from the insertion group as an example, the bases in the Read sequence are first sliced ​​into groups of 5 bits each (the last group has fewer than 5 bits, so it is temporarily counted as one group). Each group of base fragments after slicing is then processed using a set. This indicates that when correcting incorrectly inserted groups... Base fragments and coding tables Minimum Hamming distance between medium base fragments using a list To represent. When calculating the minimum Hamming distance, if the last set of bases is... If the number of digits is less than 5, the Hamming distance cannot be calculated. Temporarily ordered Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of . The following is the correction process, from the list Find the first non-zero Hamming distance The corresponding base fragment Start error correction (i.e., this 5-base fragment) Not here (in Chinese). Simultaneously, find a way to calculate the Hamming distance. hour, The corresponding base fragment Delete base fragments One of the bases ( and (Different first bases). After deletion, the base fragment... All base fragments Merge. Re-slice the merged base sequence into 5-position groups and repeat the above process to begin cyclic error correction until two conditions are met: 1. List. All values ​​in the array are 0 (meaning that after grouping the Read slices, all base fragments are in the array). (Middle). 2. The length of the corrected base sequence Read is equal to 150 nt. The base sequence Read has been initially corrected. If both conditions cannot be met simultaneously, then process the next base sequence Read. It should be noted that the above correction process may not achieve perfect error correction, because... Some base segments in the encoding table are very similar, which may indicate that after an insertion error occurs, the erroneous base segment is still present in the encoding table. The corrected or distorted base fragment is not the correct base fragment, but another base fragment in the coding table. For example, the correct fragment is... If an insertion error occurs: ,because In the encoding table, error correction starts with 'ACCAC'. During correction, 'A' will be deleted, and the correction will be completed. In this case, error correction for this segment is not achieved, but the erroneous base segments are corrected and added to the encoding table, which will not hinder the decoding of subsequent base segments. The process of correcting insertion errors is as follows: Figure 3 As shown in (b).

[0076] Case 3. Reads in the deletion error group are defined as follows: when the length of a base sequence Read is less than the designed base length (150 nt) and an error occurs during decoding, the base sequence Read is placed in the deletion error group. When using MSA error correction, the base sequence Reads in the deletion error group also suffer from sequence synchronization problems. To solve the synchronization problem, the following is a specific method for correcting Reads in the deletion error group using the HMSA of this invention.

[0077] Taking the deletion of a base sequence Read from an incorrect group as an example, similarly, the bases in the base sequence Read are first grouped into slices of 5 bits each (the last group has fewer than 5 bits, so it's temporarily counted as one group). Each group of base segments after slicing is then processed using a set. To represent, the base fragments in the coding table are represented by To represent this, after grouping the base sequence Read slices that have deletion errors, the minimum Hamming distance for the deletion error groups is still calculated using a list. To represent. When calculating the minimum Hamming distance, if the last set of bases is... If the number of digits is less than 5, the Hamming distance cannot be calculated. Temporarily ordered Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of . The following is the correction process: from the list Find the first non-zero Hamming distance The corresponding base fragment Start error correction. Simultaneously locate the corresponding encoding table. medium base fragment ( j), in the base fragment Insert a base into the middle (to make the base fragment...) The 5-base fragment containing the inserted base and the base fragment More similar, i.e., Hamming distance (When decreasing or equal to 0). All the base fragments Merge, regroup the merged Read into 5-bit slices, repeat the above process, and start from the list again. Find the first non-zero Hamming distance The corresponding base fragment Start correcting errors until the list is satisfied. When all values ​​in the read are 0 and the corrected read length is equal to 150 nt, the initial correction of the base sequence read is complete. If both conditions are not met, the next read is processed. This continues until all reads have been processed. It's important to note that the above correction process cannot guarantee the accuracy of the error correction. HMSA also faces similar issues when handling deletion errors as it does when handling insertion errors. For example, the correct segment is: After a deletion error occurs: After correction: In this case, the HMSA of this invention does not correct the error in this segment, but it corrects the erroneous base segment to the encoding table, without affecting the decoding of subsequent base segments. The process for correcting and deleting errors is as follows: Figure 3 As shown in (c).

[0078] Finally, all DNA sequences output from Case 1, Case 2, and Case 3 are grouped together, and multiple sequence alignment (MSA) is used to align all corrected reads. For each position, the base with the highest percentage of bases at that position is selected (in multiple sequencing runs, the base with the highest percentage at that position is most likely the correct base). Since this invention uses a codebook-type encoding method, after HMSA correction of the reads (i.e., after final multiple sequence alignment), base fragments may still not be in the encoding table, leading to errors during decoding. To prevent this problem, this invention incorporates an if conditional statement into the decoding algorithm. After correction, when decoding the base sequence reads, base fragments not in the encoding table are decoded as the fixed 'GTAGC'.

[0079] This invention uses an English version of the novel "The Wandering Earth" (approximately 368KB) as text material for simulation experiments. Figure 4The diagram shows the DNA fragment designed in this invention, where 150 nt of encoded data constitutes the payload. This novel was encoded into 14318 DNA sequences, each 150 nt in length. While 25 nt index bits could encode indices from 0 to 99999, choosing 25 nt would be wasteful of resources and reduce encoding density due to the 14318 DNA sequences. Therefore, 20 nt index bits (encoding indices from 0 to 9999) were chosen, allowing the text material to be stored in two test tubes to meet storage requirements. The encoding / decoding and subsequent performance analysis of this invention were simulated using Python.

[0080] In studying the characteristics of encoded DNA sequences, the mapping encoding method of the coding table was compared with the basic encoding method (00-A, 01-C, 10-G, 11-T) and six other encoding methods, including Church, Goldman, Grass, Blawat, DNA Fountain, and Yin-Yang Code. Under the premise of encoding the same text material, "The Wandering Earth" (approximately 368KB), the number of DNA sequences satisfying GC content constraints and the number of the largest homopolymers were compared after encoding using different methods. The comparison results are as follows: Figure 5 As shown.

[0081] from Figure 5 As can be seen from (a), the DNA sequences encoded using the encoding method proposed in this invention almost all have a GC content in the range of 40% to 60%. This is closely related to the characteristics of the encoding table constructed in this invention, where the GC content of the base fragments in the encoding table is almost always 40% or 60%. Regarding homopolymer constraints, as... Figure 5 As shown in (b), the maximum homopolymer size of all DNA sequences was controlled to be 3 nt. The encoded DNA sequences satisfied both the GC content constraint and the homopolymer constraint, indicating that the encoding table constructed using this invention can achieve a storage-friendly encoding mode, which is not only beneficial for DNA synthesis and polymerase chain reaction (PCR), but also can reduce the occurrence of erroneous bases to a certain extent.

[0082] To test the error correction performance of HMSA facing three different error groups, the simulated sequencing reads in this invention all contained errors of the same type and were non-contiguous. For reads with replacement errors, such as... Figure 6As shown in (a) and (b), the average error correction performance of MSA and the HMSA method of this invention in processing reads is compared. Firstly, for both MSA and HMSA, the error correction capability continuously improves with increasing sequencing runs. The comparison shows that the improvement effect of MSA is very significant, especially when the sequencing runs are only 5. The HMSA of this invention can still achieve a 100% recovery rate when processing reads with an error rate of 25%, while MSA, with 5 sequencing runs and a 100% recovery rate, can only process reads with an error rate of 20%. During simulation experiments, it was found that when HMSA corrects reads with insertion and deletion errors, the correction effect is almost 0% when the error rate is 8%. To collect more experimental data, the variable was changed from error rate to the number of errors, which also reflects the actual situation of simulated sequencing results. The performance of the HMSA of this invention in correcting reads with insertion errors is as follows: Figure 6 As shown in (c), the figure shows that with 5 sequencing runs, 2 insertion errors can be perfectly corrected. However, when dealing with reads that contain deletion errors, as shown in (c), the results are not as good. Figure 6 As shown in (d), HMSA error correction performance is not ideal when the number of sequencing runs is small.

[0083] Although this invention did not conduct wet-scale experiments on simulated DNA sequencing results, extensive literature review revealed that, in reality, excluding DNA sequences with large deletions, approximately 90% of the sequencing reads were within ±3 of the correct read length. Therefore, in the subsequent simulation experiments, measuring the number of errors in a single simulated sequencing result by the number of reads is more realistic. In the following simulation experiments, the simulated sequencing reads consisted of three random error groups. For example, when the sequencing count was 5, the simulated sequencing result contained 2 reads for substitution errors, 1 read for insertion errors, and 2 reads for deletion errors, with the number of reads being random. Figure 7 The result of each simulation experiment is the mode of the recovery rate obtained from 20 simulation experiments. Figure 7 The display shows the error correction performance of the HMSA of this invention when three erroneous groups appear in the Reads. As the number of sequencing runs increases, the error correction performance of HMSA continuously improves. Furthermore, during continuous simulation experiments, it was found that when the proportion of Reads in the erroneous groups is large, the error correction performance of HMSA is even better. This is consistent with the performance analysis of HMSA in this invention, indicating that HMSA has a high tolerance for errors when correcting erroneous replacement groups.

[0084] Previous studies have shown that consecutive insertion and deletion errors frequently occur in sequencing reads, and the longer the consecutive errors, the more difficult it is to correct them. For example, HEDGES can handle consecutive deletions but cannot tolerate more than two consecutive insertions. Therefore, further simulations were conducted to investigate the impact of consecutive errors of different lengths on the decoding performance of HMSA. Figure 8 The error correction performance of the HMSA of this invention is shown when the number of consecutive errors is 3, 4, 5, and 6. It is easy to see that when faced with consecutive errors, the HMSA's ability to correct them continuously increases with the number of sequencing runs. At 15 sequencing runs, the HMSA can perfectly recover reads with an error rate of approximately 10%.

[0085] In most previous studies, incorrectly lengthed reads were typically discarded due to the inability to resolve sequence synchronization issues caused by insertion and deletion errors, potentially resulting in significant cost and time losses in DNA storage. However, the HMSA error correction method significantly alleviates this problem, and its error correction performance is highly reliable even with consecutive errors.

[0086] Table 2 compares the encoding scheme, error correction method, error correction type, storage file type, and encoding density of this invention with those of other studies. Compared to many other studies, this invention encodes all file types while maintaining encoding density, provided that GC content and homopolymer constraints are met. Furthermore, in terms of error correction, the HMSA proposed in this invention can correct insertion and deletion errors.

[0087] Table 2 Comparison of the present invention with previous research methods

[0088]

[0089] Compared to current media, DNA has a density 6-7 orders of magnitude higher and a lifespan hundreds of times longer. However, the error problems inherent in DNA storage are the most difficult to solve. The encoding table proposed in this invention, which corresponds one-to-one with ASCII codes, satisfies GC content constraints and homopolymer constraints, reducing the occurrence of erroneous bases. It can encode all types of files (text, images, audio, and video) and uses a standardized encoding system with an encoding density of 1.4 bits / nt. In particular, for text files, direct mapping is possible, eliminating the need to convert the text file to binary during encoding and decoding, achieving an encoding density of 1.6 bits / nt. The error correction method HMSA proposed in this invention has many advantages compared to other error correction methods. First, it does not require any redundancy; second, its decoding process is linear, facilitating efficient information retrieval; and finally, it is a combination of MSA and minimum Hamming distance, inheriting all the advantages of MSA. Most importantly, HMSA solves the sequence synchronization problem that MSA cannot address, enabling it to correct not only substitution errors but also insertion and deletion errors. This allows for the effective correction of reads of varying lengths generated during sequencing, significantly improving read utilization and reducing sequencing costs. In this invention, the error-correcting capability of HMSA for DNA storage was comprehensively studied. Simulation results show that HMSA can achieve complete recovery at low error rates and low sequencing runs. At moderate error rates, the error-correcting performance of HMSA remains highly reliable with increasing sequencing runs. Furthermore, HMSA is a universal method applicable to encoding methods of any codebook type, providing a new direction for error correction for all future codebook encoding methods.

[0090] 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 principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for storing DNA information based on an improved coding table and HMSA, characterized in that, The steps are as follows: Step 1: Improve the mapping rules of Meinolf Blawat and add filtering rules to construct an encoding table that corresponds one-to-one with ASCII codes; When constructing the encoding table, {00000000, 00000001, 00000010, ..., 01111111} is taken as input. The 128 8-bit binary numbers are encoded into their corresponding bases. This is achieved by grouping the first 6 bits of each 8-bit binary number in pairs, and then grouping the 7th and 8th bits into separate groups. This is done using set... To represent; use rule A to represent the sequence The binary information in the code is encoded into base information to obtain the bases. Use rule B to sequence The binary information in the code is encoded into base information to obtain the bases. The encoded base fragment is ,in, And use three filtering criteria to filter; Step 2: Encode the information to be stored into a base sequence according to the constructed encoding table; Step 3: Error correction using HMSA: Compare the lengths of the Reads with the design sequence, divide the Reads into three groups, and correct each group separately: Slice the Reads into groups of 5 bases each, calculate the minimum Hamming distance between each group of base segments and the base segments in the coding table, find the first group of base segments with a non-zero minimum Hamming distance, and perform preliminary correction; After aligning the preliminarily corrected base sequences together, perform error correction again using multiple sequence alignment. Step 4: Use an encoding table to divide the corrected Reads into 5-base slices and decode them into stored information.

2. The DNA information storage method based on the improved coding table and HMSA according to claim 1, characterized in that, The rules for constructing the encoding table include rule A and rule B, where... Rule A: {'00' 'C', '01' 'G', '10' 'T', '11' 'A'}; TypeB:{'00': ['AA', 'CC', 'GG', 'TT'], ['AC', 'CG', 'GT', 'TA'], ['AG', 'CT', 'GA', 'TC'], ['AT', 'CA', 'GC', 'TG']}.

3. The DNA information storage method based on the improved coding table and HMSA according to claim 1, characterized in that, The three filtering criteria are: Condition 1: Check if there are any three consecutive identical bases or the last two identical bases in the four base combinations after the encoding. If so, remove the base combination; if not, proceed to Condition 2. Condition 2: Check if the GC content of the 5-position base combination is 40% or 60%. If so, put the result into the coding table of the corresponding order and proceed to Condition 3; otherwise, keep the result and proceed to Condition 3. Condition 3: Check if the encoding table is empty. If it is empty, put the retained result of Condition 2 into the corresponding encoding table. If there is no space, output the encoding table.

4. The DNA information storage method based on the improved coding table and HMSA according to claim 3, characterized in that, The encoding table is as follows: 。 5. The DNA information storage method based on the improved coding table and HMSA according to any one of claims 1-4, characterized in that, The preliminary correction method is as follows: compare the length of Reads with the designed base length, and divide Reads into replacement error group, insertion error group and deletion error group; slice and group Reads in the replacement error group, insertion error group and deletion error group one by one, calculate the minimum Hamming distance between each group of base fragments and the base fragments in the encoding table, find the first base fragment group with a minimum Hamming distance of not 0, and correct the replacement error group, insertion error group and deletion error group respectively.

6. The DNA information storage method based on the improved coding table and HMSA according to claim 5, characterized in that, When the length of the base sequence is equal to the designed base length of 150 nt and an error occurs during decoding, the base sequence is placed in the replacement error group. The initial method for correcting and replacing incorrect base sequences is as follows: the base sequence is sliced ​​into groups of 5 bases each, and each group of sliced ​​base fragments is then processed using... To represent, the base fragments in the coding table are represented by To represent; use lists Represents a base fragment Base fragments and coding tables A list of minimum Hamming distances between base segments in the nucleotide sequence; Select the minimum Hamming distance Non-zero base fragments , base fragment Flip to calculate the minimum Hamming distance The corresponding base fragment j, until the list All values ​​in the value are 0, indicating a base fragment. All base fragments All in the encoding table middle; Base fragments grouped by slice The merged DNA sequences yielded a preliminary corrected sequence.

7. The DNA information storage method based on the improved coding table and HMSA according to claim 6, characterized in that, When the length of a Read is greater than the designed DNA sequence length (150 nt) and an error occurs during decoding, the Read is placed in the insertion error group. The initial method for correcting DNA sequences in groups with insertion errors is as follows: The DNA sequence is sliced ​​and grouped into segments of 5 bases each, with the last segment having fewer than 5 bases temporarily counted as a single group; each segment is then processed using a setter. Represents a set Base fragments and coding tables Minimum Hamming distance between medium base fragments using a list This indicates that, when calculating the minimum Hamming distance, if the last set of bases... Fewer than 5 people, making Hamming's distance Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of; From list Find the first non-zero Hamming distance The corresponding base fragment Started correcting errors and simultaneously found a way to calculate Hamming distance. At that time, the encoding table The corresponding base fragment Delete base fragments Neutral base fragments Different first bases; base fragments All base fragments Merge, then re-slice the merged base sequence into 5-position groups, repeating the above process to begin cyclic error correction until the following conditions are met: list All values ​​are 0, the length of the corrected base sequence is equal to 150 nt, and the base sequence correction is initially completed.

8. The DNA information storage method based on the improved coding table and HMSA according to claim 6 or 7, characterized in that, When the length of a Read is less than the designed base length of 150 nt and an error occurs during decoding, the Read is placed in the deletion error group. The method for initially correcting the base sequence in the insertion error group is as follows: First, the bases in the base sequence are grouped into slices of 5 positions each, with the last group having fewer than 5 positions counted as a single group; each group of base segments after slicing is used as a set... This indicates that the base fragments in the encoding table are... This indicates that after grouping the base sequence slices that have deletion errors, the minimum Hamming distance calculated for each deletion error group is obtained using a list. This indicates that, when calculating the minimum Hamming distance, if the last set of bases... If the number of digits is less than 5, then the Hamming distance is... Until the base fragments are continuously corrected When the number of bits is reduced to 5, the Hamming distance is calculated again. The truth value of; From list Find the first non-zero Hamming distance The corresponding base fragment Start error correction and find the corresponding encoding table. medium base fragment , j, in the base fragment Inserting a base into the base fragment makes the base fragment The 5-base fragment containing the inserted base and the base fragment Hamming distance Decrease or equal to 0; reduce the set All the base fragments Merge, regroup the merged base sequences into 5-position slices, and then extract them from the list again. Find the first non-zero Hamming distance The corresponding base fragment Start correcting errors until the list is satisfied. When all values ​​are 0 and the length of the corrected base sequence is equal to 150 nt, the initial base sequence correction is complete. If the above conditions are not met, the initial correction is complete. If all values ​​are 0 and the length of the corrected base sequence is equal to 150 nt, then proceed to the next base sequence, and so on, until all base sequences have been processed.

9. The DNA information storage method based on the improved coding table and HMSA according to claim 8, characterized in that, All DNA sequences initially corrected from the replacement error group, insertion error group, and deletion error group are put together, and multiple sequence alignment is used to align all corrected reads, selecting the base with the largest proportion of bases at each position. After correction, the corrected Reads are sliced ​​into 5-base segments using an encoding table and then decoded into stored information. When decoding a base sequence, base fragments not in the encoding table are decoded as the fixed 'GTAGC'.

Images