DNA storage method for composite code accompanying addressing coding

Abstract

The invention discloses a DNA (deoxyribonucleic acid) storage method of composite code accompanying addressing coding, which comprises the following steps: taking a composite code constructed by a plurality of pseudo-random sub-codes as an address sequence, carrying out accompanying coding and mapping on the address sequence and a data bit sequence of error correction coding, and combining to obtain a long DNA fragment or a short oligonucleotide molecule; during reading, the long fragment DNA is randomly broken or the oligonucleotide pool is read through a second-generation high-throughput sequencing technology and is demapped into two layers of bit sequences, then a plurality of pseudo-random sub-codes and damaged composite codes are utilized to carry out sliding correlation capture of a sub-code period, and a unique address is determined according to Chinese residual definite solution. And finally, consensus and decoding are carried out on the damaged data sequence. The method provided by the invention mainly solves the problems of high addressing difficulty, high assembly complexity and DNA breakage when long fragment DNA is randomly broken into short fragments in second-generation sequencing, supports rapid addressing identification of randomly intercepted short fragments in long fragment DNA shotgun method sequencing, and is also suitable for robust addressing of short fragments, tolerance fragment partial deletion, especially sequence head and tail deletion.

Description

Technical Field

[0001] The present invention relates to the technical field of DNA storage, and in particular to a DNA storage method with composite code accompanied by addressing coding. Background Art

[0002] The rapid development of information technology has accelerated the generation of data. According to the International Data Corporation (IDC), the total amount of global data will reach 175ZB in 2025. Entering the era of digital economy, data, as a new production factor, plays an increasingly important role in social production activities. According to the "Semiconductor Ten-Year Plan" released by the Semiconductor Industry Association (SIA) and the Semiconductor Research Corporation (SRC) in 2021, more than 60% of all stored data will become cold data that does not need to be accessed frequently. According to a report released by Horizons, the value of data has the characteristic of becoming important again as its value increases over time. With the development of artificial intelligence technology based on big data, the value of long-term data storage continues to highlight. Long-term data storage brings about the demand for long-term storage of massive data. However, the service life of existing data storage methods based on magnetism, light, and electricity generally does not exceed a few decades, the maintenance cost is high, and the density is difficult to increase, which makes it difficult to meet the explosive growth of data storage needs.

[0003] DNA data storage uses deoxyribonucleic acid (DNA) as a new data storage medium. It has the core advantages of high storage density, long usable time, and low maintenance cost. It is a very promising new mode for long-term stable storage of massive data, especially in the long-term storage of archived data. According to the "Semiconductor Synthetic Biology Roadmap" released by the National Institute of Standards and Technology (NIST) and SRC, DNA data storage has an improvement of about 7 orders of magnitude compared with the current high-capacity tapes with the highest storage density; according to the forecast of the Intelligence Advanced Research Projects Agency (IARPA), the use of DNA storage in EB-level data centers is expected to reduce power consumption from 200MW to 200KW. At the same time, DNA data storage can achieve offline storage, high security, and easy backup storage, and is expected to become a potential solution for long-term stable storage of massive archived data in the future. DNA data storage is a very cutting-edge emerging technology. In 2023, the UK National Research and Innovation Agency listed DNA data storage as one of the 50 most promising emerging technologies in the future, and it is one of the 8 technologies in the direction of "artificial intelligence, digital and computing technology". In 2023, the International Devices and Systems Roadmap released by the Institute of Electrical and Electronics Engineers (IEEE) will list DNA storage as one of the main storage media for massive data in the future.

[0004] In recent years, a series of breakthroughs have been made in proof-of-concept research on DNA data storage. The current mainstream DNA data storage modes mainly include short-fragment DNA storage mode and long-fragment DNA storage mode. The short-fragment DNA storage mode uses an oligonucleotide pool composed of a large number of short-chain DNA molecules (generally 100 to 300 bases) to store data, and usually uses second-generation high-throughput sequencing to read out data. The data readout of the long-fragment DNA storage mode can be achieved with the help of similar genome sequencing technologies, including second-generation high-throughput sequencing and third-generation nanopore sequencing.

[0005] Regardless of the mode, DNA data storage uses information to be dispersed into a large number of short DNA fragments for storage, and restoring the correct order of all fragments is the key to reliable data recovery. Therefore, how to achieve efficient addressing and fast readout of DNA stored data remains a key technical problem that needs to be solved in DNA data storage.

[0006] For the long-fragment DNA data storage mode, if the second-generation high-throughput sequencing is used for reading, similar to genome sequencing, the enriched DNA is usually randomly interrupted into a large number of short fragments by "shotgun" sequencing, and the original data sequence is reconstructed by a de novo assembly strategy. However, the interruption position and the length of the short DNA fragments generated by shotgun sequencing are random, which makes addressing difficult and requires a complex graph method based on overlapping relationships. In order to improve the accurate positioning of overlapping groups during de novo assembly and avoid gaps as much as possible, the sequencing coverage required for data recovery is usually high. At present, the software for assembling short fragments of second-generation sequencing is mainly based on the de Bruijin graph algorithm, and high sequencing coverage will lead to an increase in the number of k-mer nodes in the graph construction process, further increasing the computational complexity. The researchers used the assembly software Velvet and AbySS based on the de Bruijin graph to realize the de novo assembly of second-generation sequencing reads, and recovered the payload data with the help of the sequence alignment algorithm with the assistance of the guide sequence. The sequencing coverage required for error-free data recovery using these two assembly software exceeded 20×. The researchers used the minimum sequencing coverage required for assembly (usually at least 20×) to perform de novo assembly and recovery of long-fragment data DNA. The researchers used the graph-based short sequence assembly software SOAPdenovo2 for de novo assembly and combined it with overlapping group alignment and consensus to achieve data sequence reconstruction, but the processing complexity was very high.

[0007] For the short-fragment DNA data storage mode, that is, the regularly segmented short-fragment DNA oligonucleotide storage mode, a large number of oligonucleotide molecules have disordered characteristics, and the order of reading is random. Therefore, it is generally necessary to add an additional small base sequence to each oligonucleotide molecule as a unique index to identify the position of the data carried by the molecule in the entire data. At the same time, in the actual DNA storage and reading process, DNA molecules are at risk of degradation, resulting in quality degradation, and are very easy to hydrolyze in aqueous solution, resulting in frequent DNA chain breaks, which poses a threat to the reliability of DNA data storage. In response to the problem of DNA chain breaks, researchers proposed a sequence reconstruction and error correction algorithm based on the De Blein graph theory, and verified the high robustness of data recovery through accelerated aging experiments of sample breakage and degradation. However, the sequence reconstruction algorithm based on the De Blein graph generally relies on high sequencing coverage, has high computational complexity, and is difficult to adapt to the needs of rapid data readout. Summary of the invention

[0008] In order to overcome the shortcomings of the prior art, the present invention provides a DNA storage method with a composite code accompanied by addressing coding. The present invention realizes fast and robust recognition of large-scale, short-segment data DNA at any position with very low complexity, tolerates partial loss of the internal fragment, and is applicable to both short-segment and long-segment storage modes. See the following description for details:

[0009] A DNA storage method with composite code and addressing coding, the method comprising the following steps:

[0010] (1) At the data writing end, n pseudo-random sequences are selected as subcodes, where n is a positive integer, and a composite code is obtained as an address sequence through combinatorial logic construction. The user data is encoded using error correction coding to obtain an encoded data bit sequence. The composite code and the data bit sequence are then concomitantly encoded bit by bit to form a double-layer bit sequence. The data DNA sequence is then transcoded according to the mapping rule between bit pairs and bases. Finally, the data writing is completed through artificially synthesized DNA.

[0011] (2) At the data reading end, the second-generation high-throughput sequencing technology is used to obtain the sequencing reads of short DNA fragments. The sequencing reads are preprocessed and demapped according to the same mapping rules as the writing end to extract the damaged composite code and the damaged data sequence. Then, the known n subcodes are used to perform sliding correlation detection with the damaged composite code to obtain the correlation peak phase estimation of each subcode. Then, the Chinese remainder theorem is used to solve the problem and determine the unique position label of the read segment. Then, the in-phase autocorrelation value of the composite code is calculated and a threshold check is performed. The read segments that pass the check are kept as they are, while the read segments that fail the check are re-identified using the enhanced mode, and the insertion / deletion errors of the damaged data sequence are located and corrected. Finally, according to the identified position label, multiple copies of the same position sequence are used to perform a bit-by-bit majority voting decision to generate a consistent data sequence, which is sent to the decoder for decoding to restore the original user data.

[0012] The step (1) is specifically:

[0013] (1.1) Select n pseudo-random sequences as subcodes to form a subcode set {C1, C2, …, C n}, and then extend the period of each subcode to obtain a length of P = p1p2…p n The sequence is constructed by calculating the length of P = p1p2…p bit by bit according to the preset Boolean combination logic. n The composite code a is expressed in binary form; where C n represents subcode n, whose code length is p n , n ≥ 1 and is an integer;

[0014] (1.2) m linear block codes (N, K) are used to perform group error correction coding on user data to obtain a coded data bit sequence b, where m is the number of code words, m ≥ 1 and is an integer, N is the length of the coded code word, and K is the information bit length. Since the composite code period can be designed to be very long, a single composite code can carry multiple block code words, satisfying mN ≤ p1p2…p n ;

[0015] (1.3) The composite code a and the data bit sequence b are concomitantly encoded bit by bit to construct a double-layer bit sequence, and then converted into a DNA sequence according to the mapping rule between bit pairs and bases {00→A, 01→T, 10→G, 11→C}, and a long DNA fragment or a short oligonucleotide molecule is obtained by artificial synthesis technology;

[0016] The preset Boolean combination logic specifically includes: logical multiplication, modulo-two sum, and majority decision;

[0017] Among them, n pseudo-random sequences with short periods are selected as subcodes, which satisfy: the periods of the selected pseudo-random sequences are strictly mutually prime, and each of them has a sharp periodic autocorrelation characteristic;

[0018] Wherein, the step (1.3) is specifically as follows:

[0019] For long-fragment DNA data storage, the composite code is first truncated to obtain a shortened composite code of the same length as the data bit sequence, and then the shortened composite code is aligned with the data bit sequence to form a double-layer bit sequence, and directly mapped to obtain a long DNA fragment; for short-fragment oligonucleotide data storage, the data bit sequence is first divided into mN / l groups according to the length l, and the composite code is slidingly truncated according to the step length l to obtain mN / l shortened composite codes, and then the shortened composite code of length l is paired with the data bit sequence in sequence to construct mN / l groups of double-layer bit sequences, and after mapping, mN / l short-fragment oligonucleotides are finally generated; wherein l is a positive integer.

[0020] The step (2) is specifically:

[0021] (2.1) For the short-fragment oligonucleotide pool, the library is directly constructed and then read by the second-generation high-throughput sequencing technology to obtain the sequencing reads of the data DNA; for the long-fragment data DNA, high-throughput shotgun sequencing is used for sequencing readout. Specifically, the enriched DNA molecules are randomly interrupted and read by the second-generation high-throughput sequencing to obtain the sequencing reads;

[0022] (2.2) Preprocess the sequencing reads, demap them according to the same mapping rules as the writing end, and extract the damaged composite codes from the double-layer structure. With damaged data bit sequence Specifically, for the sequencing reads of short-fragment oligonucleotide molecules, the primer parts existing at both ends are trimmed off by alignment with the amplification primers, and the retained payload parts are demapped; for the reads obtained by long-fragment DNA shotgun sequencing, the sequencing reads are directly demapped;

[0023] (2.3) The damaged composite code The cross-correlation detector of each subcode is sent to the known subcode set {C1, C2, ..., C n} Perform sliding correlation capture, search for the position of the correlation peak of each subcode (i.e., the phase corresponding to the in-phase correlation), obtain the estimated phase of each subcode, and then determine the unique position label of the read segment through the Chinese residual theorem;

[0024] (2.4) Using the phase estimates of each subcode obtained above, regenerate an ideal composite code a′ with the same length and the same phase as the damaged composite code, calculate the autocorrelation value of the composite code, and retain the read segment if the autocorrelation value meets the preset conditions; for the read segment that fails the verification, use the method based on the edit distance alignment in the enhanced mode to re-identify it, determine the position label, and retain the read segment, and at the same time identify and correct the insertion and deletion error patterns according to the alignment path;

[0025] (2.5) Considering that there may be multiple copies of the data sequence at the same position, a majority vote is performed on the data sequence bit by bit to obtain a consistent data sequence. If the voting decision is not unique, the same erasure mark is used to fill the bit at that position to obtain a codeword, which is then sent to the decoder to correct the residual errors and restore the original user data.

[0026] The step (2.3) is specifically:

[0027] (3.1) Each subcode in the known subcode set is periodically extended, and each subcode is subjected to sliding correlation operation with the damaged composite code. The correlation window length is the same as the length of the damaged composite code, and the sliding step length is the subcode period. The cross-correlation function between subcode i and the damaged composite code can be expressed as:

[0028]

[0029] in, is the periodic extension sequence of subcode i, is the damaged composite code, τ is the subcode phase offset, L is the length of the correlation window, 1≤i≤n, and the periodic extension sequence With composite code It is represented by the symbol {+1, -1}, which satisfies the mapping relationship {+1→0, -1→1} with the binary sequence;

[0030] (3.2) According to the sliding correlation value between each subcode and the damaged composite code, the position of the correlation peak is searched to obtain the phase estimation of each subcode. The correlation peak of subcode i corresponds to the phase τ i It can be expressed as:

[0031]

[0032] Among them, argmax represents the position corresponding to the correlation peak;

[0033] (3.3) Using the estimated phases of each subcode obtained by the search, the unique position label of the corresponding read segment is obtained by the Chinese remainder theorem. The estimated phases of each subcode satisfy

[0034]

[0035] in, is the unique position label of the read segment obtained by solution, τ1, τ2, ..., τ n They correspond to the estimated phases of n subcodes respectively.

[0036] The step (2.4) is specifically:

[0037] (4.1) Based on the estimated phase of each subcode, regenerate a composite code with the same phase and the same length as the damaged one The same ideal composite code a′ is then used to calculate the damaged composite code The in-phase correlation value between the ideal composite code a′ and the obtained correlation value can be expressed as:

[0038]

[0039] Wherein, r and r′ are the damaged composite code and the ideal composite code respectively, both of which are represented by the symbol {+1, -1}, and satisfy the mapping relationship {+1→0, -1→1} with the binary form sequence;

[0040] (4.2) If the in-phase autocorrelation value of the composite code satisfies R r,r′ ≥ωL, the verification is passed and the current read segment is retained. ω represents the preset correlation threshold, 0≤ω≤1.

[0041] (4.3) For the reads that fail the above autocorrelation check, the damaged composite code is replaced by the edit distance comparison method in the enhanced mode. Compare with the known complete composite code a, search for the alignment position where the minimum edit distance is located, and use it as the position label of the read segment. At the same time, use the alignment path to locate the position where the insertion and deletion errors in the damaged composite code occur, and according to the characteristic that the positions of the insertion / deletion errors in the double-layer sequence structure are exactly the same, correct the insertion and deletion error patterns in the damaged data sequence; discard the position bits where the insertion errors occur, mark the position bits where the deletion errors occur as erased, and fill them with the determined markers; if the enhanced mode is not used, skip the current step.

[0042] Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects:

[0043] (1) The present invention makes full use of the unambiguous single-cycle characteristics of the composite code and the good cross-correlation characteristics between the subcode and the composite code, and uses the composite code generated by multiple short-cycle subcodes as the address sequence, and performs accompanying encoding with the original data sequence to form a DNA sequence; when reading data, a sliding correlation detection with a limited number of subcode cycles is performed, and then the actual position of the read segment is calculated using the Chinese remainder theorem, thereby quickly identifying the unique position label of the sequencing read segment, significantly improving the addressing efficiency of large-scale data storage read recovery.

[0044] (2) The present invention is applicable to the storage mode of long DNA fragments and short oligonucleotide pools, supports the robust recognition of DNA fragments at any position, and avoids the complex assembly operations required for sequence reconstruction of long fragments by shotgun sequencing; for the oligonucleotide pool with a fixed starting position, partial deletions are allowed within the fragment, thereby improving the readout efficiency of the DNA oligonucleotide pool. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 A system implementation block diagram of a DNA storage method with composite code accompanied by addressing coding provided by the present invention;

[0046] Figure 2 The present invention is a flowchart of implementing the composite code accompanying addressing encoding that constitutes the storage mode of the long DNA fragment or short oligonucleotide pool;

[0047] Figure 3 A flowchart of the read segment position label identification and data recovery for the second generation sequencing readout provided by the present invention;

[0048] Figure 4 A flow chart of composite code construction provided by the present invention;

[0049] Figure 5 A flowchart for realizing fast identification of read segments based on short-period subcode sliding correlation detection provided by the present invention;

[0050] Figure 6 This is a performance comparison diagram of the read segment recognition success rate and the available read segment ratio under different normalization thresholds in specific embodiment 1 of the present invention;

[0051] Figure 7 This is a performance diagram of the segment recognition success rate at different autocorrelation normalization thresholds after the enhancement mode is adopted in the specific embodiment 1 of the present invention;

[0052] Figure 8 This is a graph of error rate performance after merging multiple copies at different sequencing coverages in specific embodiment 1 of the present invention;

[0053] Fig. 9 This is a time complexity comparison diagram of specific embodiment 1 of the present invention;

[0054] Fig.10 This is a performance comparison diagram of the read segment recognition success rate and the available read segment ratio under different normalization thresholds in specific embodiment 2 of the present invention;

[0055] Fig.11 This is a graph of error rate performance after merging multiple copies at different sequencing coverages in specific embodiment 2 of the present invention;

[0056] Fig.12 This is the error-free recovery situation of 50 groups of repeated experiments at different sequencing coverages in specific embodiment 3 of the present invention. DETAILED DESCRIPTION

[0057] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention are described in further detail below.

[0058] In view of the problem of low sequence addressing efficiency in short-fragment and long-fragment DNA application scenarios, the present invention provides a new encoding method that is generally applicable to short-fragment DNA data storage and long-fragment DNA data storage, and a DNA storage method with composite code accompanying addressing encoding. This method constructs a DNA sequence by combining composite code segments formed by logical combination of multiple pseudo-random sequences with the data part, and realizes that DNA sequencing reads starting at any position (shotgun long-fragment sequencing mode) and starting at a fixed position (oligonucleotide pool sequencing readout mode) can be quickly and efficiently addressed during reading, and rapid data recovery can be achieved, and partial deletion of fixed-segment short-fragment DNA molecular sequences can be handled, thereby improving the readout efficiency.

[0059] The following describes the embodiments of the present invention in detail with reference to the accompanying drawings:

[0060] This embodiment introduces in detail a DNA storage method for composite code accompanied by addressing coding proposed by the present invention. For the overall system implementation block diagram, see Figure 1 , specifically including the following steps:

[0061] (1) Data writing end, see Figure 2 , select n pseudo-random sequences as subcodes, n is a positive integer, and obtain a composite code as an address sequence through combinatorial logic construction, encode the user data using error correction coding to obtain an encoded data bit sequence, and then perform bit-by-bit concomitant coding on the composite code and the data bit sequence to form a double-layer bit sequence, and according to the mapping rule between bit pairs and bases, transcode to obtain a data DNA sequence, and finally complete data writing through artificially synthesized DNA;

[0062] (2) Data readout terminal, see Figure 3, with the help of second-generation high-throughput sequencing technology, the sequencing reads of short-fragment DNA are obtained, the sequencing reads are preprocessed, and the mapping is performed according to the same mapping rules as the writing end to extract the damaged composite code and the damaged data sequence. Secondly, the known n subcodes are used to perform sliding correlation detection with the damaged composite code respectively to obtain the correlation peak phase estimation of each subcode. Then, the Chinese remainder theorem is used to solve the problem and determine the unique position label of the read segment. Then, the in-phase autocorrelation value of the composite code is calculated, and a threshold check is performed. The read segments that pass the check are kept as they are, while the read segments that fail the check are re-identified using the enhanced mode, and the insertion / deletion errors of the damaged data sequence are located and corrected. Finally, according to the identified position label, multiple copies of the same position sequence are used to perform bit-by-bit majority voting to generate a consistent data sequence, which is sent to the decoder for decoding to restore the original user data.

[0063] Wherein, step (1) is specifically as follows:

[0064] (1.1) Figure 4 As shown in the figure, n pseudo-random sequences are selected as subcodes to form a subcode set {C1, C2, ..., C n}, and then extend the period of each subcode to obtain a length of P = p1p2…p n The sequence is constructed by calculating the length of P = p1p2…p bit by bit according to the preset Boolean combination logic. n The composite code a is expressed in binary form; where C n represents subcode n, whose code length is p n , n ≥ 1 and is an integer;

[0065] (1.2) m linear block codes (N, K) are used to perform group error correction coding on user data to obtain a coded data bit sequence b, where m is the number of code words, m ≥ 1 and is an integer, N is the length of the coded code word, and K is the information bit length. Since the composite code period can be designed to be very long, a single composite code can carry multiple block code words, satisfying mN ≤ p1p2…p n ;

[0066] (1.3) The composite code a and the data bit sequence b are concomitantly encoded bit by bit to construct a double-layer bit sequence, and then converted into a DNA sequence according to the mapping rule between bit pairs and bases {00→A, 01→T, 10→G, 11→C}, and a long DNA fragment or a short oligonucleotide molecule is obtained by artificial synthesis technology;

[0067] The preset Boolean combination logic specifically includes: logical multiplication, modulo-two sum, and majority decision, etc.

[0068] Among them, n pseudo-random sequences with short periods are selected as subcodes, which satisfy: the periods of the selected pseudo-random sequences are strictly mutually prime, and each of them has a sharp periodic autocorrelation characteristic;

[0069] Among them, step (1.3) is specifically:

[0070] For long-fragment DNA data storage, the composite code is first truncated to obtain a shortened composite code of the same length as the data bit sequence, and then the shortened composite code is aligned with the data bit sequence to form a double-layer bit sequence, and directly mapped to obtain a long DNA fragment; for short-fragment oligonucleotide data storage, the data bit sequence is first divided into mN / l groups according to the length l, and the composite code is slidingly truncated according to the step length l to obtain mN / l shortened composite codes, and then the shortened composite code of length l is paired with the data bit sequence in sequence to construct mN / l groups of double-layer bit sequences, and after mapping, mN / l short-fragment oligonucleotides are finally generated; wherein l is a positive integer.

[0071] Wherein, step (2) is specifically as follows:

[0072] (2.1) For the short-fragment oligonucleotide pool, the library is directly constructed and then read by the second-generation high-throughput sequencing technology to obtain the sequencing reads of the data DNA; for the long-fragment data DNA, the high-throughput shotgun sequencing technology is used for sequencing readout. Specifically, the enriched DNA molecules are randomly interrupted and read by the second-generation high-throughput sequencing technology to obtain the sequencing reads;

[0073] (2.2) Preprocess the sequencing reads, demap them according to the same mapping rules as the writing end, and extract the damaged composite codes from the double-layer structure. With damaged data bit sequence Specifically, for the sequencing reads of short-fragment oligonucleotide molecules, the primer parts existing at both ends are trimmed off by alignment with the amplification primers, and the retained payload parts are demapped; for the reads obtained by long-fragment DNA shotgun sequencing, the sequencing reads are directly demapped;

[0074] (2.3) Figure 5 As shown, the damaged composite code The cross-correlation detector of each subcode is sent to the known subcode set {C1, C2, ..., C n} Perform sliding correlation capture, search for the position of the correlation peak of each subcode (i.e., the phase corresponding to the in-phase correlation), obtain the estimated phase of each subcode, and then determine the unique position label of the read segment through the Chinese residual theorem;

[0075] (2.4) Using the phase estimates of each subcode obtained above, regenerate an ideal composite code a′ with the same length and the same phase as the damaged composite code, calculate the autocorrelation value of the composite code, and retain the read segment if the autocorrelation value meets the preset conditions; for the read segment that fails the verification, use the method based on the edit distance alignment in the enhanced mode to re-identify it, determine the position label, and retain the read segment, and at the same time identify and correct the insertion and deletion error patterns according to the alignment path;

[0076] (2.5) Considering that there may be multiple copies of the data sequence at the same position, a majority vote is performed on the data sequence bit by bit to obtain a consistent data sequence. If the voting decision is not unique, the same erasure mark is used to fill the bit at that position to obtain a codeword, which is then sent to the decoder to correct the residual errors and restore the original user data.

[0077] Among them, step (2.3) is specifically:

[0078] (3.1) Each subcode in the known subcode set is periodically extended, and each subcode is subjected to sliding correlation operation with the damaged composite code. The correlation window length is the same as the length of the damaged composite code, and the sliding step length is the subcode period. The cross-correlation function between subcode i and the damaged composite code can be expressed as:

[0079]

[0080] in, is the periodic extension sequence of subcode i, is the damaged composite code, τ is the subcode phase offset, L is the length of the correlation window, 1≤i≤n, and the periodic extension sequence With composite code It is represented by the symbol {+1, -1}, which satisfies the mapping relationship {+1→0, -1→1} with the binary sequence;

[0081] (3.2) According to the sliding correlation value between each subcode and the damaged composite code, the position of the correlation peak is searched to obtain the phase estimation of each subcode. The correlation peak of subcode i corresponds to the phase τ i It can be expressed as:

[0082]

[0083] Among them, argmax represents the position corresponding to the correlation peak;

[0084] (3.3) Using the estimated phases of each subcode obtained by the search, the unique position label of the corresponding read segment is obtained by the Chinese remainder theorem. The estimated phases of each subcode satisfy:

[0085]

[0086] in, is the unique position label of the read segment obtained by solution, τ1, τ2, ..., τ n They correspond to the estimated phases of n subcodes respectively.

[0087] Among them, step (2.4) is specifically:

[0088] (4.1) Based on the estimated phase of each subcode, regenerate a composite code with the same phase and the same length as the damaged one The same ideal composite code a′ is then used to calculate the damaged composite code The in-phase correlation value between the ideal composite code a′ and the obtained correlation value can be expressed as:

[0089]

[0090] in, r′ is a damaged composite code and an ideal composite code, both of which are represented by symbols {+1, -1}, and satisfy the mapping relationship {+1→0, -1→1} with the binary sequence;

[0091] (4.2) If the in-phase autocorrelation value of the composite code satisfies If the verification is successful, the current read segment is retained. ω represents the preset correlation threshold, 0≤ω≤1.

[0092] (4.3) For the reads that fail the above autocorrelation check, the damaged composite code is replaced by the edit distance comparison method in the enhanced mode. Compare with the known complete composite code a, search for the alignment position where the minimum edit distance is located, and use it as the position label of the read segment. At the same time, use the alignment path to locate the position where the insertion and deletion errors in the damaged composite code occur, and according to the characteristic that the positions of the insertion / deletion errors in the double-layer sequence structure are exactly the same, correct the insertion and deletion error patterns in the damaged data sequence; specifically, discard the position bits where the insertion errors occur, mark the position bits where the deletion errors occur as erased, and fill them with the determined markers; if the enhanced mode is not used, skip the current step.

[0093] The present invention proposes a DNA storage method with composite code accompanied by addressing coding. The composite code uses the good correlation characteristics of the composite code and the distinguishable characteristics of any phase in a single cycle. The composite code constructed by multiple short-period pseudo-random subcodes is used as an address sequence, and is combined with a data sequence to form a double-layer sequence structure. After mapping, a DNA sequence is obtained. When reading data, combined with the Chinese remainder theorem, only a limited number of subcodes need to be used to perform sliding correlation detection of the subcode period to achieve rapid recognition of read segments, support robust recognition of DNA fragments at any position, and can tolerate partial deletions within the fragment. The method is generally applicable to the storage mode of long DNA fragments and short fragment oligonucleotide pools for second-generation sequencing readout.

[0094] Two specific embodiments are given below, and combined with the accompanying drawings, to illustrate the feasibility of the present invention and further understand the purpose, features and advantages of the invention. In the three specific embodiments given, the Boolean combination logic for constructing the composite code uses majority decision, and its expression can be expressed as in Represents the symbol sequence corresponding to the subcode n, which consists of {+1, -1} and satisfies the mapping relationship {-1→1, +1→0} with the bit sequence. The subcode number n satisfies n>1 and is an odd number.

[0095] Example 1

[0096] This embodiment uses a composite code generated by constructing 5 subcodes as the address sequence of the data DNA. By simulating the recognition performance of a sequencing read segment with a length of 130 bases, the low complexity and high robustness of the addressing method proposed in the present invention are verified. Five pseudo-random sequences with periods of 31, 35, 43, 47, and 59 are selected to form a subcode set, among which the subcode with a period of 31 is an m sequence, the subcode with a period of 35 is a double prime sequence, and the subcodes with periods of 43, 47, and 59 are Legendre sequences. Table 1 shows the bit sequence corresponding to each subcode.

[0097] Table 1 Five pseudo-random subcodes that constitute the composite code

[0098]

[0099] The specific implementation of the DNA storage method with composite code and addressing coding includes the following steps:

[0100] (1) The five selected subcodes are periodically extended to obtain period extension sequences of each subcode, each of which has a length of 129374315. Then, majority decision is used as the Boolean combination logic for constructing a composite code, and bitwise operations are performed to obtain a composite code with a period of P = 129374315, which is represented in binary form.

[0101] (2) Determine the starting phase of the composite code, then segment the data sequence and the composite code into groups of 130 bits in length, and obtain 995,187 shortened composite codes and data sequences, which are paired in pairs to form a double-layer bit sequence;

[0102] (3) According to the mapping rule {00→A, 01→T, 10→G, 11→C}, the bit pairs are converted into bases, and finally 995,187 DNA sequences of length L = 130 bases are obtained;

[0103] (4) Add random errors to each sequence to simulate and generate simulated sequencing reads, where the insertion (Ins.), deletion (Del.), and substitution (Sub.) error rates meet P Sub. =2P Ins. =2P Del. =0.0050, and Monte Carlo repeated tests were performed;

[0104] (5) directly demapping the sequencing reads to obtain the damaged composite code and data bit sequence;

[0105] (6) Using the five known subcodes, sliding correlation detection is performed with the damaged composite code, and a total of 31+35+43+47+59=215 subcode correlation operations are performed to obtain the estimated phase of each subcode correlation peak, and then the unique position label of the current read segment is obtained by the Chinese remainder theorem.

[0106] (7) Using the known subcode and the estimated phase of each identified subcode, an ideal composite code with a length of 130 bits is regenerated, and an in-phase autocorrelation operation is performed with the damaged composite code. The normalized autocorrelation value is compared with a preset threshold, and the threshold is set to 0.8. If the autocorrelation value is greater than or equal to the preset threshold, the position label recognition is valid, otherwise the read segment is discarded; in the enhanced mode, the read segments with autocorrelation values ​​less than the preset threshold are retained, and the damaged composite code is compared with the known complete composite code using a method based on edit distance comparison, and the alignment position corresponding to the minimum edit distance is searched as the position label, and the insertion and deletion error patterns are identified and corrected according to the alignment path;

[0107] (8) Based on the identified data location information, the damaged data sequence is consensus-generated and decoded to restore the original data.

[0108] Figure 6 The performance comparison of read segment recognition accuracy and available read segment ratio under different thresholds is given. The simulation results show that when the number of sub-codes is 5 and the verification threshold is set to 0.8, the read segment recognition accuracy can reach more than 98%, and the available read segment ratio is nearly 80%, which verifies that the method proposed in the present invention can achieve highly robust recognition in scenarios containing insertion / deletion errors.

[0109] Figure 7 The percentage of reads accurately identified in the enhanced mode is given, and simulation results show that when the threshold is set to 0.8 or above, the percentage of correctly identified reads in all reads increases to 97.69%, indicating that the enhanced mode of the method proposed in the present invention can effectively increase the number of valid reads that can be used while ensuring the accuracy of read identification.

[0110] Figure 8 The bit error rate performance of the method proposed in the present invention under different sequencing coverages is given. The simulation results show that after setting the threshold check, the bit error rate of the data after merging multiple copies can be greatly improved. When the enhanced mode is adopted, the error rate performance is further improved as the sequencing coverage increases, which is close to the method based on reference sequence alignment.

[0111] Fig. 9 A time complexity performance comparison is given, and 1,399,667 reads with a length of 130 bases are identified and the running time is counted. The results show that the processing speed of the method proposed in the present invention is about 230 times faster than that of the method based on reference sequence alignment, and the read segment recognition efficiency is significantly improved. After adopting the enhanced mode, the processing speed is also about 5 times faster than that of the method based on reference sequence alignment, which verifies that the implementation complexity of the method proposed in the present invention is low.

[0112] Example 2

[0113] This embodiment uses a composite code generated by constructing three subcodes as the address sequence of the data DNA. By simulating the recognition performance of a sequencing read segment with a length of 130 bases, it is further verified that the addressing method proposed in the present invention still has low complexity and high robustness when the number of subcodes is 3. Three pseudo-random sequences with periods of 499, 503, and 511 are selected to form a subcode set, among which the subcodes with periods of 499 and 503 are Legendre sequences, and the subcode with a period of 511 is an m sequence. Table 1 shows the bit sequence corresponding to each subcode.

[0114] Table 2 Three pseudo-random subcodes that constitute the composite code

[0115]

[0116]

[0117] Firstly, the three selected subcodes are periodically extended respectively, and the combinational logic of majority decision is adopted to construct a composite code with a period of P=128259467, which is expressed in binary form. Secondly, the data sequence and the composite code are segmented into groups of 130 bits in length, and 986611 shortened composite codes and data sequences are obtained respectively. They are paired in pairs in turn to form a double-layer bit sequence, and 986611 DNA fragments with a length of L=130 bases are obtained after mapping.

[0118] The sequencing readout process was further simulated by adding random errors to each sequence to generate simulated sequencing reads, where the insertion (Ins.), deletion (Del.) and substitution (Sub.) error rates met P Sub. =2P Ins. =2P Del. =0.0050.

[0119] The specific implementation of data reading refers to steps (5)-(8) of Example 1, except that this embodiment uses three known subcodes for sliding correlation detection.

[0120] Fig.10 The performance comparison of read segment recognition accuracy and available read segment ratio under different thresholds is given. The simulation results show that when the number of sub-codes is 3 and the verification threshold is set to 0.8, the read segment recognition accuracy can reach more than 98%, and the available read segment ratio is nearly 80%, which verifies that the method proposed in the present invention can still achieve highly robust recognition in scenarios containing insertion / deletion errors when the number of sub-codes is 3.

[0121] Fig.11 The bit error rate performance of the method proposed in the present invention under different sequencing coverages is given. The simulation results show that after setting the threshold check, when the sequencing coverage is about 5×, the error rate performance of the non-enhanced mode and the enhanced mode is basically consistent with the method based on reference sequence alignment, which verifies that the proposed method can obtain performance similar to that of the method based on reference sequence alignment when the number of subcodes is 3.

[0122] Example 3

[0123] This embodiment demonstrates the readout recovery of long-fragment DNA for second-generation sequencing. The five subcodes in Embodiment 1 are selected to form a subcode set, and the group error correction code adopts a binary LDPC (22680, 7560) code with a code rate R=1 / 3.

[0124] For the long-fragment DNA storage mode, the specific implementation of the composite code adjoint addressing encoding is:

[0125] The composite code construction process is the same as that in Example 1, and a composite code with a period of P=129374315 is obtained. The user data is encoded using a binary LDPC (22680, 7560) code to obtain a code word with a length of 22680 bits to form a data sequence. The starting phase of the composite code is determined and the composite code is shortened to obtain a shortened composite code with a length of 22680 bits. The shortened composite code is then formed with the data sequence into a double-layer bit sequence, and transcoded according to the mapping rule of bit pairs and bases {00→A, 01→T, 10→G, 11→C}, and finally a single long DNA fragment with a length of 22680 bases is obtained.

[0126] The second-generation genome sequencing data simulation software ART was used to simulate the random interruption and sequencing process of long DNA fragments, and the single-end 150 sequencing mode was set. The insertion, deletion and substitution error rates met the P Sub. =2P Ins. =2P Del. =0.0050, thereby obtaining a short-read sequencing read segment of 150 bases in length. The specific operation of read recovery is consistent with steps (5)-(8) in Example 1.

[0127] Fig.12 The error-free recovery conditions of 50 groups of repeated experiments at different sequencing coverages are given. The results show that the method proposed in the present invention can achieve error-free recovery at 2.6× sequencing coverage (a well-known term in the art). After adopting the enhanced mode, error-free recovery can be achieved at 2.2× sequencing coverage, and the performance is the same as that based on the reference sequence alignment method, which proves that the method proposed in the present invention can achieve completely error-free recovery of original data at lower sequencing coverage.

[0128] Those skilled in the art will appreciate that the accompanying drawing is only a schematic diagram of a preferred embodiment, and the serial numbers of the embodiments of the present invention are only for description and do not represent the advantages or disadvantages of the embodiments.

[0129] 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 DNA storage method with composite code and addressing encoding, characterized in that: The method comprises the following steps: (1) At the data writing end, n pseudo-random sequences are selected as subcodes, where n is a positive integer, and a composite code is obtained as an address sequence through combinatorial logic construction. The user data is encoded using error correction coding to obtain an encoded data bit sequence. The composite code and the data bit sequence are then concomitantly encoded bit by bit to form a double-layer bit sequence. The data DNA sequence is then transcoded according to the mapping rule between bit pairs and bases. Finally, the data writing is completed through artificially synthesized DNA. (2) At the data reading end, the sequencing reads of short DNA fragments are obtained by means of second-generation high-throughput sequencing. The sequencing reads are preprocessed and demapped according to the same mapping rules as the writing end to extract the damaged composite code and the damaged data sequence. Then, the known n subcodes are used to perform sliding correlation detection with the damaged composite code to obtain the correlation peak phase estimation of each subcode. Then, the Chinese remainder theorem is used to solve the problem and determine the unique position label of the read. Then, the in-phase autocorrelation value of the composite code is calculated and a threshold check is performed. The reads that pass the check are kept as they are, while the reads that fail the check are re-identified using the enhanced mode, and the insertion / deletion errors of the damaged data sequence are located and corrected. Finally, according to the identified position label, multiple copies of the same position sequence are used to perform a bit-by-bit majority voting decision to generate a consistent data sequence, which is sent to the decoder for decoding to restore the original user data.

2. The DNA storage method with composite code and addressing encoding according to claim 1, characterized in that: The step (1) is specifically: (1.1) Select n pseudo-random sequences as subcodes to form a subcode set {C1, C2, …, C n }, and then extend the period of each subcode to obtain a length of P = p1p2…p n The sequence is constructed by calculating the length of P = p1p2…p bit by bit according to the preset Boolean combination logic. n The composite code a is expressed in binary form; where C n represents subcode n, whose code length is p n , n ≥ 1 and is an integer; (1.2) Use m linear block codes (N, K) to perform group error correction coding on the user data to obtain the coded data bit sequence b; where m is the number of code words, m ≥ 1 and is an integer, N is the length of the coded code word, K is the information bit length, and satisfies mN ≤ p1p2…p n ; (1.3) The composite code a and the data bit sequence b are concomitantly encoded bit by bit to construct a double-layer bit sequence, and then converted into a DNA sequence according to the mapping rule between bit pairs and bases {00→A, 01→T, 10→G, 11→C}, and a long DNA fragment or a short oligonucleotide molecule is obtained by artificial synthesis; The preset Boolean combination logic specifically includes: logical multiplication, modulo-two sum, and majority decision; Among them, n pseudo-random sequences with short periods are selected as subcodes, which satisfy: the periods of the selected pseudo-random sequences are strictly mutually prime, and each of them has a sharp periodic autocorrelation characteristic; Wherein, the step (1.3) is specifically as follows: For long-fragment DNA data storage, the composite code is first truncated to obtain a shortened composite code of the same length as the data bit sequence, and then the shortened composite code is aligned with the data bit sequence to form a double-layer bit sequence, and directly mapped to obtain a long DNA fragment; for short-fragment oligonucleotide data storage, the data bit sequence is first divided into mN / l groups according to the length l, and the composite code is slidingly truncated according to the step length l to obtain mN / l shortened composite codes, and then the shortened composite code of length l is paired with the data bit sequence in sequence to construct mN / l groups of double-layer bit sequences, and after mapping, mN / l short-fragment oligonucleotides are finally generated; wherein l is a positive integer.

3. The DNA storage method with composite code and addressing encoding according to claim 1, characterized in that: The step (2) is specifically: (2.1) For short-fragment oligonucleotide pools, the library is directly constructed and then read by second-generation high-throughput sequencing to obtain sequencing reads of the data DNA; for long-fragment data DNA, high-throughput shotgun sequencing is used for sequencing reads; (2.2) Preprocess the sequencing reads, demap them according to the same mapping rules as the writing end, and extract the damaged composite codes from the double-layer structure. With damaged data bit sequence For sequencing reads of short-fragment oligonucleotide molecules, the primers at both ends are trimmed off by alignment with the amplification primers, and the retained payload is demapped; for reads obtained by long-fragment DNA shotgun sequencing, the sequencing reads are directly demapped; (2.3) The damaged composite code The cross-correlation detector of each subcode is sent to the known subcode set {C1, C2, ..., C n } Perform sliding correlation capture, search for the location of the correlation peak of each subcode, obtain the estimated phase of each subcode, and then determine the unique position label of the read segment through the Chinese residual theorem; (2.4) Using the phase estimates of each subcode obtained above, regenerate an ideal composite code a′ with the same length and the same phase as the damaged composite code, calculate the autocorrelation value of the composite code, and retain the read segment if the autocorrelation value meets the preset conditions; for the read segment that fails the verification, use the method based on the edit distance alignment in the enhanced mode to re-identify it, determine the position label, and retain the read segment, and at the same time identify and correct the insertion and deletion error patterns according to the alignment path; (2.5) Considering that there are multiple copies of the data sequence at the same position, a majority vote is performed on the data sequence bit by bit to obtain a consistent data sequence. If the voting decision is not unique, the same erasure mark is used to fill the bit at that position to obtain a codeword, which is then sent to the decoder to correct the residual errors and restore the original user data.

4. The DNA storage method with composite code and addressing encoding according to claim 3 is characterized in that: The step (2.3) is specifically: (3.1) Each subcode in the known subcode set is periodically extended, and each subcode is subjected to sliding correlation operation with the damaged composite code. The correlation window length is the same as the length of the damaged composite code, and the sliding step length is the subcode period. The cross-correlation function between subcode i and the damaged composite code is expressed as: in, is the periodic extension sequence of subcode i, is the damaged composite code, τ is the subcode phase offset, L is the length of the correlation window, 1≤i≤n, and the periodic extension sequence With composite code It is represented by the symbol {+1, -1}, which satisfies the mapping relationship {+1→0, -1→1} with the binary sequence; (3.2) According to the sliding correlation value between each subcode and the damaged composite code, the position of the correlation peak is searched to obtain the phase estimation of each subcode. The correlation peak of subcode i corresponds to the phase τ i It is expressed as: Among them, argmax represents the position corresponding to the correlation peak; (3.3) Using the estimated phases of each subcode obtained by the search, the unique position label of the corresponding read segment is obtained by the Chinese remainder theorem. The estimated phases of each subcode satisfy in, is the unique position label of the read segment obtained by solution, τ1, τ2, ..., τ n They correspond to the estimated phases of n subcodes respectively.

5. The DNA storage method with composite code and addressing encoding according to claim 3 is characterized in that: The step (2.4) is specifically: (4.1) Based on the estimated phase of each subcode, regenerate a composite code with the same phase and length as the damaged The same ideal composite code a′ is then used to calculate the damaged composite code The in-phase correlation value between the ideal composite code a′ and the obtained correlation value is expressed as: in, r′ is a damaged composite code and an ideal composite code, both of which are represented by symbols {+1, -1}, and satisfy the mapping relationship {+1→0, -1→1} with the binary sequence; (4.2) If the in-phase autocorrelation value of the composite code satisfies If the verification is successful, the current read segment is retained, where ω represents the preset correlation threshold, 0≤ω≤1; (4.3) For the reads that fail the above autocorrelation check, the damaged composite code is replaced by the edit distance comparison method in the enhanced mode. Compare with the known complete composite code a, search for the alignment position where the minimum edit distance is located, and use it as the position label of the read segment. At the same time, use the alignment path to locate the position where the insertion and deletion errors in the damaged composite code occur, and according to the characteristic that the positions of the insertion / deletion errors in the double-layer sequence structure are exactly the same, correct the insertion and deletion error patterns in the damaged data sequence; discard the position bits where the insertion errors occur, mark the position bits where the deletion errors occur as erased, and fill them with the determined markers; if the enhanced mode is not used, skip the current step.

Images