Robust multi-sequence reconstruction method based on maximum a posteriori probability in DNA storage
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- TIANJIN UNIV
- Filing Date
- 2023-08-07
- Publication Date
- 2026-06-02
Smart Images
Figure CN117116359B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of DNA storage technology, and more specifically to a robust multiple sequence reconstruction method based on maximum posterior probability in DNA storage. Background Technology
[0002] In the era of big data, traditional storage media can no longer meet the current demand for storing massive amounts of data. Media such as magnetic tape, hard drives, and flash memory face severe challenges due to their short storage time, high energy consumption, and susceptibility to data loss. At the same time, deoxyribonucleotide (DNA) molecules, as a medium for storing biological information, are expected to become the next generation of storage technology and replace traditional storage media due to their advantages such as high storage density, long storage time, and low energy consumption.
[0003] The DNA storage workflow typically involves five basic steps: encoding binary information into a DNA sequence (encoding), writing into DNA (synthesis), storing DNA (storage), reading the DNA sequence (sequencing), and recovering the binary information (decoding), such as... Figure 1 As shown.
[0004] Limited by existing DNA storage system technologies, errors are inevitably introduced during DNA synthesis, storage, and sequencing. There are three main types of errors: insertion, deletion, and substitution. Insertion and deletion errors are collectively referred to as synchronization errors. Currently, the probability of an error per base in mainstream next-generation sequencing is 1%-2.5%. Furthermore, in practical DNA storage, PCR amplification is typically performed to effectively reconstruct the original binary information. This process amplifies each original DNA reference sequence into multiple reads with errors, storing these erroneous reads from different reference sequences together in an unordered manner. Therefore, when recovering binary information, clustering methods are needed to process the sequencing files to group reads from the same original DNA reference sequence into a cluster. Multiple sequence reconstruction aims to infer the original DNA reference sequence or original binary information from multiple erroneous reads within a cluster. However, limited by existing DNA storage system technologies, the reliability of multiple sequence reconstruction methods is affected not only by IDS errors within each sequence but also by outliers within the cluster when reading information from the sequencing file. Outlier sequences have a large edit distance from the original sequence, which may be due to the following four reasons:
[0005] (1) Sequences from other clusters. Before reconstruction, sequencing files are typically clustered to ensure that reads from the same original DNA reference sequence are accurately assigned to their respective clusters. This perfect clustering assumption ensures that subsequent sequence reconstruction steps focus only on error correction of base IDS within each receiving cluster. However, this assumption is not always valid. In the clustering results, it is unavoidable to find erroneous clustered sequences from other DNA reference sequences.
[0006] (2) Strand breaks and rearrangements. Long-term storage of DNA exposes the bases to hydrolysis and degradation, leading to DNA sequence breaks. At the same time, during sequencing, non-specific amplification issues in the PCR process can also cause DNA sequence rearrangements, i.e., strand breaks and rearrangements.
[0007] (3) The inverse complementary strand of the original DNA reference strand. When a DNA sequence is stored in an organism, the two strands form a double helix structure. Once the double-stranded DNA fragment is sequenced, these two strands are separated, and usually only one strand is selected for decoding. This process results in the inverse complementary strand of the selected original DNA reference sequence being present in the sequencing file.
[0008] (4) Random sequences added to ensure data security. When considering the security of data storage, random sequences are intentionally added for data encryption.
[0009] Current multiple sequence reconstruction algorithms are limited to error correction of base IDS errors within each error cluster, ignoring outliers. In other words, during reconstruction, error reads within the same cluster are considered to originate from the same original DNA reference sequence and are therefore treated equally without differentiation in quality and reliability. To meet practical DNA storage needs, a multiple sequence reconstruction method robust to both in-cluster outliers and conventional base IDS errors is urgently required. Considering the impact of outliers in an IDS channel model based on base units further complicates the multiple sequence reconstruction problem. Summary of the Invention
[0010] To meet the practical needs of DNA storage, the purpose of this invention is to provide a robust multiple sequence reconstruction method based on maximum posterior probability in DNA storage.
[0011] To achieve the objectives of this invention, the technical solution provided by this invention is as follows:
[0012] A robust multiple sequence reconstruction method based on maximum a posteriori probability in DNA storage includes the following steps:
[0013] Step 1: Decode each sequence within the cluster using the improved BCJR decoder;
[0014] Step 2: Convert each sequence within the cluster into its corresponding time series representation;
[0015] Step 3: Determine the weight of each sequence;
[0016] Step 4: Derive the MAP algorithm formula for robust multi-sequence decoding, and incorporate sequence weights into the joint posterior probability calculation to suppress the influence of outlier sequences.
[0017] Compared with the prior art, the beneficial effects of the present invention are as follows:
[0018] This invention improves the BCJR decoding algorithm by proposing a base IDS channel model, thereby overcoming the limitation of existing statistical inference-based sequence reconstruction methods, which are only suitable for binary error-correcting codes. Therefore, this invention can handle quaternary error correction in the form of bases in DNA sequences.
[0019] This invention obtains sequence weights by mapping relative sequence quality to a parametric hyperbolic tangent function and incorporates these sequence weights into the joint posterior probability calculation. Considering the different contributions of each sequence within a cluster, this results in an improved weighted MAP separable multi-sequence decoding algorithm. In clusters containing a certain proportion of outlier sequences, this invention exhibits better reconstruction performance and stronger anti-interference capabilities compared to other comparison algorithms, thus demonstrating robustness. Attached Figure Description
[0020] Figure 1 A flowchart illustrating the DNA storage process;
[0021] Figure 2 This is a schematic diagram of the IDS channel model used in this application;
[0022] Figure 3 This is a diagram illustrating the process of encoding and converting binary information and obtaining multiple noise reads through a DNA storage channel in this application.
[0023] Figure 4 This is a diagram of the base IDS channel model in this application;
[0024] Figure 5 This is a schematic diagram of the crystal lattice in this application;
[0025] Figure 6 A schematic diagram illustrating an example of DNA time series representation;
[0026] Figure 7 A schematic diagram illustrating the robustness of the weighting method to different outlier proportions in practice;
[0027] Figure 8 This is a diagram comparing the error rate and success rate of editing. Figure 1 ;
[0028] Figure 9 This is a diagram comparing the error rate and success rate of editing. Figure 2 ;
[0029] Figure 10 This diagram illustrates the comparison of editing error rates at different DNA sequence lengths. Detailed Implementation
[0030] 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 of ordinary skill in the art without creative effort are within the scope of protection of the present invention.
[0031] The focus of this invention is to address the problem of multi-sequence reconstruction when outlier sequences exist within a cluster. Specifically, for sequences encoded and transformed using [2,1,2] convolutional codes, after multiple passes through an IDS channel and the addition of outlier sequences, noisy sequence clusters are generated, simulating the clustering of sequencing data in real DNA storage. The goal of this invention is to establish a robust reconstruction algorithm based on maximum a posteriori probability decoding for these noisy sequence clusters, capable of resisting the influence of outlier sequences and accurately inferring the encoded original binary information.
[0032] To address the above challenges, this invention proposes a robust multi-sequence decoding algorithm suitable for DNA storage, used for multi-sequence decoding of [2,1,2] convolutional codes. This algorithm not only resists insertion / deletion errors within the sequence but also resists the influence of outliers on reconstruction. Specifically,
[0033] First, an improved (single-sequence) decoding algorithm for [2,1,2] convolutional code encoding is proposed. By establishing the conversion relationship between binary bit errors and base errors in the channel model, an improved BCJR decoding algorithm is proposed for channels with base insertion, deletion and replacement errors. The relationship between DNA storage and convolutional code decoding is established, which alleviates the problem that the BCJR decoding method cannot be directly applied to quaternary decoding.
[0034] Secondly, robust multi-sequence reconstruction is performed using the improved single-sequence decoding algorithm in (1). By incorporating sequence weights into the joint posterior probability calculation and considering the different contributions of each sequence within the cluster, an improved weighted maximum a posteriori probability separable multi-sequence decoding algorithm is formed. This achieves automatic suppression of outlier sequences while preserving the contributions of high-quality sequences, thereby improving the robustness of multi-sequence decoding.
[0035] This invention employs the widely used IDS (Insertion-Deletion-Substitution) channel model to simulate errors occurring in DNA storage, such as... Figure 2 As shown. This channel has p i p d p s The three parameters represent the probabilities of insertion, deletion, and substitution errors, respectively. The DNA sequence is transmitted symbol by symbol into the channel, and the final received sequence is an error read with IDS (Input / Deletion / Substitution) errors. Assume that at time i, the base symbol x... i When queuing to enter the channel, three scenarios can occur: (e.g., p) i An insertion error occurs with a probability of 0, in which case a random base will be inserted into the received sequence, x. i Still in the channel; with p d The probability of a deletion error occurring is high; the received sequence did not add any bases. i+1 Queuing to enter the channel; with p t =1-p i -p d The probability of being transmitted is p; during transmission, the symbol either... s The probability is that it is replaced with other bases and added to the receiver sequence, or in 1-p... s The probability of the transmitted DNA sequence being correctly transmitted is high. However, due to insertion and deletion errors, the transmitted DNA sequence may differ in length from the received sequence.
[0036] The encoding scheme used in this invention is based on [2, 1, 2] convolutional codes. The binary information m = (m1, m2, ..., m...). K ), m i ∈{0,1}, encoded by [n,k,m] convolutional code, the generated codeword is e=(e1,e2,...,e 2L ), 2L=(K+m)n / k. Here, a [2,1,2] convolutional code is used. The codeword e corresponds to the DNA sequence x=(x1,x2,...,x... L The conversion between ) is determined according to ∑={A:00,T:01,C:10,G:11}, x i ∈{A, T, C, G}. x passes through the base IDS channel independently multiple times, generating M error reads carrying IDS errors. Figure 3 It demonstrates the process of encoding and converting binary information, and obtaining multiple noise reads (clusters) through a DNA storage channel.
[0037] set up (The corresponding binary sequence is) This is a cluster read during reconstruction, which includes not only multiple error copies from the original DNA reference sequence x, but also a small subset of outlier sequences generated at different stages of DNA storage. Next, the main algorithm flow of this invention is described. For the received noisy sequence cluster, a robust reconstruction algorithm based on maximum a posteriori probability decoding is established. This algorithm not only resists insertion / deletion errors within the sequence but also resists the influence of outlier sequences, accurately inferring the encoded original binary information.
[0038] The method provided by this invention includes four steps:
[0039] Step 1: Decode each sequence within the cluster using the improved BCJR decoder.
[0040] Step 2: Convert each sequence within the cluster into its corresponding time series representation.
[0041] Step 3: Determine the weight of each sequence.
[0042] Step 4: Derive the MAP algorithm formula for robust multi-sequence decoding, and incorporate sequence weights into the joint posterior probability calculation to suppress the influence of outlier sequences.
[0043] Preferably, step 1: Decoding each sequence within the cluster using an improved BCJR decoder adapted to the base IDS channel, including the following:
[0044] Step 1 improves the convolutional code BCJR (single sequence) decoder to adapt to the DNA storage channel feature model. The improved decoder is then used to decode each sequence within the cluster. That is, Step 1 is the single sequence y j The bit-by-bit maximum a posteriori probability decoder will output m at the i-th position of the original message. i discrete probability distribution p(m) i |y j ), (m i =0,1). The output of step 1 will be applied in step 4.
[0045] Consider a channel model suitable for DNA storage, where base IDS errors occur during storage. However, the current BCJR algorithm is designed for decoding binary convolutional codes and is not directly applicable to DNA channel models. This invention establishes a correlation between bit IDS errors in the binary received sequence and base IDS errors in the DNA received sequence. By imposing constraints on the decoding process, a maximum a posteriori probability decoding algorithm suitable for DNA storage channel characteristics is established.
[0046] Assume p = (p i ,p d ,p s )and These represent the IDS error probabilities of the DNA channel model and the binary channel model, respectively. This invention establishes... The connection with the base IDS channel model, such as Figure 4 As shown. Where e 2i-1 e 2i It corresponds to Figure 2 Chinese x i The binary representation of the error probabilities. The correspondence between the two IDS error probabilities is as follows:
[0047] (1) Insertion: The insertion of one base in DNA sequence x is equivalent to the insertion of two consecutive bits in binary sequence y. Therefore...
[0048] (2) Deletion: The deletion of one base in DNA sequence x is equivalent to the deletion of two consecutive bits in binary sequence y. Therefore,
[0049] (3) Substitution: A base substitution in DNA sequence x is equivalent to one or two consecutive bit substitutions in binary sequence y. Therefore,
[0050] For a received sequence y with an IDS error, the BCJR algorithm is used for decoding. This decoding algorithm is based on the Maximum A Posteriori (MAP) probability and compares the information bits m using the following formula. i The posterior probabilities of all possible values (i.e., a discrete probability distribution can be obtained for each bit of information):
[0051]
[0052] Where the joint probability p(y,m) i By marginalizing m i The position is obtained by using the states of all convolutional code shift registers. This position will be decoded to maximize the posterior probability m. i Values.
[0053] Introducing joint state variable σ i =(s i ,d i ), where s i The shift register state representing the convolutional code, the new hidden state variable "offset" d i For time i, e ni+1 Before entering the channel, e1,…,e ni The difference between the number of insertion errors and the number of deletion errors that have occurred. The joint probability can be expressed as:
[0054]
[0055] The summation takes into account the information bits m. i All possible state pairs can be taken. This formula can be further extended to the product of three parts:
[0056]
[0057] The above three parts are respectively denoted as: α i-1 (σ i-1 ),γ i (σ i-1 ,σ i ),β i (σ i The forward and backward recursion formulas are: The intermediate variable can be decomposed into: Where p(m) i ) is bit information m i The prior probability. Expression Efficient calculations can be performed using crystal lattice structures.
[0058] like Figure 5 As shown, in the crystal structure F n×μ In the diagram, the (n+1) rows correspond to a transmission sequence of length n. Each column corresponds to a received sequence of length μ. Horizontal movement in the lattice corresponds to a probability of The insertion event, vertical movement corresponds to a probability of The deletion event, moving along the diagonal, represents the transmission of one bit, if and If the corresponding element in the sequence matches, the transmission probability is . Otherwise, the transmission probability is In the binary model, q takes the value 2.
[0059] Based on the characteristic that errors in DNA channels occur at the base level, in the improved BCJR decoding, the offset vector d can only take even values at all times. The improved lattice eliminates impossible movement paths in DNA storage. Using F... r,l This represents the value in the r-th row and l-th column of the crystal lattice. For 2 ≤ r ≤ n, 2 ≤ l ≤ μ, the values of each element in the crystal lattice are calculated recursively:
[0060]
[0061] in:
[0062]
[0063] The initial values for lattice calculations are defined as follows:
[0064]
[0065] final, The above improvements establish the relationship between DNA storage and convolutional code decoding, alleviating the problem that the BCJR decoding method cannot be directly applied to quaternary decoding.
[0066] Preferably, step 2: converting each DNA sequence within the cluster into its corresponding time series representation, including the following:
[0067] Suppose x = (x1, x2, ..., x L Let s = (s1, s2, ..., sG) be a DNA sequence defined on {A, T, C, G}. It is transformed into a complex sequence (or vector sequence) using a one-to-one mapping defined in the complex plane. L Let ) be the corresponding complex sequence, where:
[0068]
[0069] satisfy:
[0070]
[0071] Here, i is the imaginary unit. To avoid symbol misuse, the phase sequence associated with the complex sequence is denoted as s, which is regarded as the time series representation of DNA. Figure 6 An example is given. The DNA sequence ACTAGTCAGG is converted into a complex sequence (0.5+0.87i,1.37+0.37i,1.87-0.5i,2.37+0.37i,3.24+0.87i,3.74,4.6-0.5i,5.1+0.37i,5.96+0.87i,6.83+1.37i), and its associated phase sequence s = (1.05,0.26,-0.26,0.15,0.26,0,-0.11,0.07,0.14,0.2) is the time series representation of this DNA sequence.
[0072] Preferably, step 3: using the DTW and DBA algorithms, assign weights to each sequence within the cluster, including the following:
[0073] Step 3 will use the important time series analysis algorithms DTW and DBA. Let's first briefly introduce these two classic algorithms.
[0074] Dynamic Time Warping (DTW) has long been considered the best approach for solving time series problems, often used as a similarity measure to quantify the dissimilarity between two series. DTW, based on the Levenshtein distance, allows for the global alignment of two time series with time offsets in a non-linear manner. Let A = (a1, ..., a... |A| ) and B = (b1,…,b |B| Given two time series, where |·| represents the series length. The cost of optimal alignment is recorded in a |A|×|B| matrix D, whose elements are obtained recursively:
[0075]
[0076] Among them, A i Represents the subsequence (a1,…,a) i ), where δ is the distance between sequence elements, such as Euclidean distance. The overall DTW distance is given by D(A |A| ),B |B| The following is given. Specifically, in optimal alignment, if a i With (b) j ,b j+1 If , ...) are aligned, then (b j ,b j+1 (,...) is called a i The correlation. The main steps of the DTW algorithm are summarized in Algorithm 1.
[0077]
[0078] DTW Barycenter Averaging (DBA) aims to find the average sequence of a given set of time series using the DTW algorithm. It is a global averaging method (i.e., insensitive to the order of the series) that iteratively optimizes the provisional average sequence to minimize the sum of squared DTW distances to the set of series. The process of one iteration of the DBA algorithm is as follows:
[0079] (1) Calculate the DTW distance between each sequence and the average sequence, and find the correlation between the coordinates of the average sequence and the coordinates of the sequence set;
[0080] (2) Update each coordinate of the average sequence to its associated centroid.
[0081] Assumption We want the average time series set, c = (c1, ..., c T () is the average sequence. In each iteration, the j-th coordinate of c is updated as follows:
[0082] cj =barycenter(assoc(c j )),
[0083] in Experimental results show that the DBA algorithm can achieve good results by randomly selecting a sequence from the sequence set as the initial averaging sequence. Furthermore, the DBA algorithm can be determined to have converged when the correlation of the averaging sequence stops updating or the sum of squared DTW distances from the averaging sequence to the sequence set stops decreasing, at which point iteration should be stopped.
[0084] Step 3 uses the time-series cluster of the DNA sequence obtained in Step 2, denoted as...
[0085] Calculate the DTW distance between any two time series within a cluster, denoted as . The DBA algorithm is used to estimate the average time series (cluster centroid) of the cluster, denoted as c. The DTW distance from each time series within the cluster to the cluster centroid c is calculated:
[0086]
[0087] Assume {d1,d2,…,d} M Sort in descending order. This represents the average distance from the sequence within a cluster to the cluster centroid. The relative quality of the sequences within the cluster is then sorted in ascending order as follows:
[0088]
[0089] in Clearly, sequences corresponding to larger values have better quality, and vice versa.
[0090] The following formula (2) is used to assign weights to each sequence, with the aim of suppressing the weights of outlier sequences to a smaller weight, while concentrating the weights of non-outlier sequences around 1. Given a cluster containing M DNA reads... Assume w j ∈(0,1) is the j-th sequence The weights are calculated as follows:
[0091]
[0092] Here, the parameter a in the hyperbolic tangent function tanh(ax) is calculated according to (3):
[0093]
[0094] Where ξ is a preset threshold close to 1, r is a preset maximum proportion of outlier sequences. In this invention, ξ = 0.98, corresponding to the k-th... * The weight of each read. This means that the weight of a reliable sequence (number ≈ (1-r)×M) will be greater than 0.98, while the weight of the next kth read will be less than 0.98. * Each sequence read (number ≈ r × M) will be assigned a lower weight.
[0095] The maximum outlier ratio *r* is a hyperparameter determined jointly by storage conditions, synthesis, and sequencing technologies. This weighting scheme is insensitive to this parameter; even if the actual outlier ratio in a cluster does not reach a predefined level, this method can still achieve a reasonable and discriminative weighting allocation. Figure 7 As shown. Assume a 60 nt DNA reference sequence enters the base IDS channel, generating a cluster containing M = 10 error reads, with a channel error rate of 2.5%. Replace the 0, 1, and 2 reads in each of the 10 reads with outlier sequences, generating three actual outlier proportions r. * Clusters with outlier ratios of 0, 0.1, and 0.2 are identified. The maximum outlier ratio is uniformly set to r = 0.2, allowing for the calculation of the corresponding parameter 'a' and the weights of each sequence within each cluster. Red markers represent the weights of outlier sequences within the cluster, while black markers represent the weights of error copies from the original reference sequence. It is observed that, under different actual outlier ratios, the proposed method consistently provides intra-cluster discriminative weight allocation; that is, the weights of outlier sequences are significantly suppressed to less than 0.5, while the weights of high-quality sequences are very close to 1.
[0096] Preferably, step 4: derives the MAP algorithm formula for robust multi-sequence decoding, incorporating sequence weights into the joint posterior probability calculation to suppress the influence of outlier sequences, including the following:
[0097] Step 4 uses each sequence y within the cluster output in Step 1. j The decoding result p(m) i |y j (This result can be represented as each original information bit m) i (discrete probability distribution on the data), and also uses the weights w of each sequence within the cluster output from step 3. j .
[0098] Based on the sequence weights in step 3, the posterior probability distribution of the single sequence obtained in step 1 is improved, resulting in formula (5). Combining the results of multi-sequence decoding, a MAP decoder that considers sequence weights is obtained, namely formula (6).
[0099] Theorem 1: (Boosted Posterior Probability Distribution): Suppose p(m) is the probability distribution of a discrete random variable m, where m takes q possible values, i.e. It is a real number. Define a binary function. for:
[0100]
[0101] So It has the following properties:
[0102]
[0103] (2) It is the probability distribution of the discrete random variable m;
[0104] (3) if So:
[0105]
[0106] (4)
[0107] (5)
[0108] Proof: Only property (2) will be proved. Consider a fixed w0∈[0,1]. First, by property (1), for all have Secondly, from the expression in (4), we have therefore, It is the probability distribution of the discrete random variable m.
[0109] Using (4), the following enhanced single-sequence posterior probability distribution is obtained:
[0110]
[0111] It is worth noting that when w j When =1, according to property (5), we have This means that the improved posterior probability distribution can be reduced to the original probability distribution. When w j When = 0, according to property (4), we will get m i Taking a discrete uniform distribution of all possible values means that there is no information available when decoding by comparing posterior probabilities.
[0112] This invention proposes a MAP decoder that considers sequence weights to solve the multi-sequence reconstruction problem:
[0113]
[0114] Among them, the binary information symbol m i Determined by hard decision, i.e., comparison m i The joint posterior probability of all possible values (if p(m) is satisfied)i =0|y1,…,y M )>p(m i =1|y1,…,y M If ), then the symbol will be decoded into m. i =0, and vice versa).
[0115] The underlying mechanism of equation (6) is analyzed below. Due to the prior probability distribution p(m) i The probability of a single sequence is fixed, and the decision depends only on the product of the posterior probabilities of multiple single sequences. Suppose y1 is an outlier sequence in the cluster. Since its weight w1 is very small, the function... The original probability distribution p(m) i |y1) is pulled toward a discrete uniform distribution (without changing the probability order among candidate values). This means that the influence of this sequence on the final decision will be weakened. The extreme case of w1=0 will lead to… This means that this sequence will not affect the final decoding result at all. Conversely, assume y2 is a high-quality sequence with weights close to 1. When applying the function... Afterwards, its original posterior probability distribution will only change slightly, which means that in equation (6), the decoding information obtained from this sequence is almost completely retained in the final decision. The extreme case when w2=1 will result in In other words, all decoding information from a single sequence is used in the reconstruction process. Essentially, the proposed MAP decoder distinguishes different reads based on their sequence quality, automatically suppressing the influence of sensitive outliers with small weights, while preserving the contribution of high-quality sequences.
[0116] Experiments and Results
[0117] The first five chapters of the English version of *The Little Prince* were used as experimental data, and each symbol was converted into a 7-bit binary representation using the standard ASCII code table to simulate data storage in real-world applications. This binary data was encoded using [2,1,2] convolutional codes and written into a short DNA sequence of 60 nt. This length is widely used in mainstream next-generation sequencing technologies. To simulate the DNA storage channel, the IDS error probability for each base was set to be equal, i.e., p i =p d =p s Furthermore, IDS errors were randomly introduced into the transmitted sequence. The error rate of traditional DNA storage systems is approximately 1%, therefore, the error rate for each base was set within the set {1%, 1.5%, 2%, 2.5%}. This was achieved using formula 1-(1-p... i (1-p) d (1-p) sThe probability of IDS errors can be calculated to be between 0.334% and 0.84%. Based on the formula relating base error rate and binary error rate described above, a numerical table of IDS error probabilities is provided, as shown in Table 1. These error probabilities will be used as prior information in the improved decoding algorithm; specifically, they will be substituted into the improved lattice calculation to accurately estimate and correct errors during the decoding process, thereby recovering the original sequence.
[0118] At each base IDS error rate, a dataset with 1000 distinct clusters was generated, each containing M error copies from a 60 nt long DNA reference sequence. To test the algorithm's robustness to outlier sequences, a certain percentage of the copy sequences within each cluster were replaced with outlier sequences, ranging from {0% to 10% to 20%}. Specific experimental parameters are shown in Table 2. In the experiment, the offset limit was set to... Where N = 2L is the binary codeword length. When the range of offset vector d is large, the summation term in the joint probability calculation will consider more possible cases, thus comprehensively estimating and correcting errors, but the computation time will increase. Conversely, when the range of offset vector values is small, some cases in the joint probability calculation are ignored, which may lead to a decrease in accuracy. Therefore, selecting an offset vector range that conforms to the channel characteristics is important, ensuring both accuracy and decoding efficiency. Based on the statistical regularity of the generated data, determining τ as 4 under different base error rates yields suitable offset vector values, thereby improving decoding efficiency while ensuring accuracy.
[0119] Table 1: Correspondence between base error rate and binary error rate
[0120]
[0121] Table 2: Parameters used in robust multiple sequence alignment experiments
[0122]
[0123]
[0124] The effectiveness of the proposed method is primarily evaluated through comparison with the MAP separable multiple sequence decoding algorithm, which treats each sequence in the cluster equally, contributing the posterior probability of each sequence indiscriminately to the joint posterior probability calculation. Simultaneously, it is compared with three reconstruction methods based on the BMA principle.
[0125] (1) MAP Separable Multi-Sequence Decoding Algorithm: Two methods are proposed for multi-sequence reconstruction. The first is a joint decoding algorithm, which uses the BCJR decoding algorithm to accurately calculate the posterior probability of each symbol in multi-dimensional space and joint grid. The other method is an approximation algorithm—the separable decoding algorithm—which estimates the posterior probability of each symbol in multi-sequences by combining the results of single-sequence decoding in an approximate manner. This method has lower computational complexity and is more practical. This invention uses the second method to decode [2,1,2] convolutional codes for performance comparison.
[0126] (2) DividerBMA: This is a variant of the BMA algorithm that corrects insertion and replacement errors in addition to deletion errors. The received sequence set is divided into three sub-clusters according to sequence length. A majority vote is performed on sequences without synchronization errors. Deletion and insertion errors are corrected in the shorter and longer sub-clusters, respectively.
[0127] (3) Hybrid: This is a hybrid algorithm. Specifically, if the number of sequences in a cluster is greater than 20, or the absolute value of the difference between the average sequence length in the cluster and the original sequence length is greater than 10% of the original sequence length, the Hybrid method will call other algorithms.
[0128] (4) Lookahead BMA: First, a majority vote is taken for each symbol. When a symbol in an error sequence does not match the current majority vote symbol, a lookahead window is used to check the next two (or more) symbols. Next, these symbols are compared with the majority vote symbol, and any deviations are classified as errors.
[0129] Decoding performance is evaluated using two criteria:
[0130] (1) Editing error rate:
[0131]
[0132] Edit distance (also known as Levenshtein distance) is a commonly used sequence alignment metric to describe the differences between two sequences of different lengths. It represents the minimum number of insertion, deletion, and substitution operations required to transform one sequence into another. This metric calculates the average edit distance of bases between the decoded sequence and the actual reference sequence.
[0133] (2) Success rate: The proportion of sequences that are completely decoded correctly out of the total number of transmitted sequences.
[0134]
[0135] This is a common standard for fully evaluating algorithm performance; a sequence only contributes to the success rate if it is perfectly reconstructed at every index position.
[0136] For a cluster size of M=5, the reconstruction performance of different multiple sequence reconstruction methods in terms of editing error rate and success rate was compared, and the results are shown in Tables 3 and 4, respectively. Similarly, for a cluster size of M=10, the experimental results of editing error rate and success rate are shown in Tables 5 and 6, respectively.
[0137] As shown in the table, compared with the main comparison object—the equally weighted MAP separable decoding algorithm—the proposed weighted MAP decoder improves decoding performance under both total IDS error rate and outlier proportion for a cluster size of M=10, achieving optimization in both metrics. When M=5, the decoding performance slightly lags behind the equally weighted method under the conditions of 1% / nt and 1.5% / nt total error rate and no outliers, with the maximum difference in edit error rate and success rate being 0.03% and 0.3%, respectively. Apart from these, the proposed method outperforms the equally weighted algorithm in both edit error rate and success rate under other conditions. The advantages of the proposed method become increasingly apparent with the increase of outlier proportion and total error rate per base, demonstrating its robustness to base IDS errors and outliers. Even under the most stringent condition of 2.5% error rate and 20% outlier proportion, the proposed method still achieves a low edit error rate. For M=5 and M=10, the proposed method achieved edit error rates of 1.69% and 0.41%, respectively, which surpasses the performance of the equal-weighted MAP decoder, which has edit error rates of 6.35% and 2.67%, respectively. Compared with BMA-based reconstruction methods, the proposed method outperforms DividerBMA and Hybrid in all cases. However, it consistently ranks second, behind Lookahead BMA, in terms of both total error rate and outlier sequence proportion. In particular, this variant of the BMA algorithm achieves perfect reconstruction in all cases when the number of sequences is sufficiently large, i.e., M=10. It is also noted that, in the presence of outliers, the reconstruction performance generally improves for all methods as the cluster size increases from M=5 to M=10.
[0138] Table 3: Comparison of editing error rates for different reconstruction methods when cluster size M=5
[0139]
[0140] Table 4: Comparison of success rates of different reconstruction methods when cluster size M=5
[0141]
[0142] Table 5: Comparison of editing error rates for different reconstruction methods when cluster size M = 10
[0143]
[0144]
[0145] Table 6: Comparison of success rates of different reconstruction methods when cluster size M=10
[0146]
[0147] Furthermore, under stringent conditions—a total base error rate of 2.5% and an outlier sequence ratio of 20%—the error rate and success rate of the proposed method and the equal-weighted MAP separable decoding algorithm are demonstrated when the cluster size increases to M=15. Figure 8 , Figure 9 As can be seen, the method of the present invention exhibits good reconstruction performance under different cluster sizes. Furthermore, when the cluster size is sufficiently large, a low editing error rate of only 0.1% can be achieved.
[0148] Due to limitations in current technology, excessively long DNA sequences cannot be synthesized during DNA storage. In mainstream next-generation sequencing technologies, DNA sequence lengths are typically within 200 nt, usually 60-180 nt. The impact of DNA sequences of different lengths on the reconstruction performance of this invention was investigated under the editing error rate standard. Experimental parameters are shown in Table 7, and experimental results are as follows: Figure 10 As shown, the editing error rate increases with increasing DNA sequence length. This is because the number of base errors increases with sequence length. For sequences shorter than 100 nt, decoding is almost perfectly accurate. However, for longer sequences, increasing the number of sequences can improve decoding accuracy, providing a solution for handling long DNA sequences in future practical applications.
[0149] Table 7: Parameters used in experiments with different sequence lengths
[0150]
[0151] Finally, it should be noted that the above embodiments are merely illustrative and explanatory of the present invention, and are not intended to limit the present invention to the scope of the described embodiments. Furthermore, those skilled in the art will understand that the present invention is not limited to the above embodiments, and many more variations and modifications can be made based on the teachings of the present invention, all of which fall within the scope of protection claimed by the present invention.
Claims
1. A robust multiple sequence reconstruction method based on maximum a posteriori probability in DNA storage, characterized in that, Includes the following steps: Step 1: Decode each sequence within the cluster using the improved BCJR decoder; the improved BCJR decoder in Step 1 is a single-sequence decoder. The bit-by-bit maximum a posteriori probability decoder, in the original message of the first bit... Each position will output information bits. discrete probability distribution , The output will be applied in step 4. By establishing the relationship between bit IDS errors in the binary received sequence and base IDS errors in the DNA received sequence, and by imposing constraints on the decoding process, a maximum a posteriori probability decoding algorithm suitable for the characteristics of DNA storage channels is established. The correspondence between the probabilities of the two types of IDS errors is then introduced into a joint state variable. ,in, For offset; based on the characteristic that errors in DNA channels occur in base units, in the improved BCJR decoder, the value of the offset vector can only be an even number at all times; the improved lattice removes impossible movement paths in DNA storage and recursively calculates each value of the lattice. Step 2: Convert each sequence within the cluster into its corresponding time series representation; Step 3: Determine the weight of each sequence; Step 4: Derive the MAP algorithm formula for robust multi-sequence decoding, incorporating sequence weights into the joint posterior probability calculation to suppress the influence of outlier sequences; Step 4 specifically includes using each sequence within the cluster output in Step 1. Decoding results The weights of each sequence within the cluster output in step 3. Based on the sequence weights in step 3, the posterior probability distribution of the single sequence obtained in step 1 is improved, resulting in the following formula: ； By combining the results of multi-sequence decoding, a MAP decoder that considers sequence weights is obtained, as follows: 。 2. The robust multiple sequence reconstruction method based on maximum a posteriori probability in DNA storage according to claim 1, characterized in that, Step 2 specifically includes the following: Assume It is defined in A DNA sequence; convert it into a complex sequence or vector sequence using a one-to-one mapping defined in the complex plane; let Let be the corresponding complex sequence, where: ； satisfy: ； in, It is the imaginary unit.
3. A robust multiple sequence reconstruction method based on maximum a posteriori probability in DNA storage according to claim 2, characterized in that, Step 3 specifically includes: using the DTW and DBA algorithms to assign weights to each sequence within the cluster.