A DNA storage coding optimization method with high robustness

Abstract

The present invention discloses a DNA storage coding optimization method with high robustness, which specifically comprises the following steps: constructing an optimal DNA coding sequence that meets constraints, firstly constructing a certain number of DNA sequences as an initial set, and evaluating and ranking the fitness of the initial set. Secondly, the obtained initial set is optimized by a gradient algorithm and a mutation strategy to obtain a sequence with higher fitness. Then, it is determined whether the sequence meets the constraints, and if so, it is added to the DNA coding set. Finally, the optimal DNA coding set is output. In addition to the three classic constraints, the present invention also proposes a new non-adjacent subsequence constraint, which can improve the quality of DNA coding. Therefore, the present invention can not only search for a relatively excellent number of DNA coding sets but also obtain a highly robust DNA coding set.

Description

Technical Field

[0001] The present invention relates to the field of coding design in DNA storage, and in particular to a DNA storage coding optimization method with high robustness. Background Art

[0002] In 1953, Watson and Crick published one of the most fundamental articles in the history of natural biology, revealing the structure of DNA molecules as carriers of genetic information. Since then, people have realized that the genetic information of organisms is stored in a linear sequence of four bases in DNA. In just ten years, many researchers have proposed the concept of storing specific information in DNA. However, because DNA synthesis and sequencing technology are still in their infancy, this concept has not been realized. In 1988, Davis et al. made the first attempt to build real DNA storage. He encoded the information of an image into a 28 base pair (bp) long DNA molecule and inserted it into Escherichia coli. After DNA sequencing and decoding, the original image was successfully restored. It was not until the first decade of the 21st century that the pioneering work of Church and Goldman opened the boom of DNA storage. At the same time, with the maturity of DNA synthesis and sequencing technology, it is a general trend that DNA storage technology will be widely used in the next few years. However, the current DNA encoding methods generally have the problem of low robustness. Summary of the invention

[0003] In view of the above problems in the prior art, the present application proposes a gradient optimization algorithm based on the Cauchy-Lévy mutation strategy and the constraint condition of no adjacent subsequences to construct a DNA coding sequence with higher robustness.

[0004] To achieve the above objectives, the technical solution of the present application is: a DNA storage coding optimization method with high robustness, comprising:

[0005] Step 1: Initialize the parameters, initial set and DNA encoding set required by the gradient algorithm based on the mutation strategy;

[0006] Step 2: Obtain the fitness value of the initial set, sort them according to the fitness value, and take the DNA sequence with the current best fitness as the optimal solution;

[0007] Step 3: Each DNA sequence in the initial set generates a random number r. When the random number r is greater than or equal to the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. When the random number r is less than the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. The DNA sequence is then updated through the local escape operator to obtain the fitness value of the corresponding DNA sequence.

[0008] Step 4: Using the updated DNA sequence as the input of the Cauchy mutation, the DNA sequence is perturbed by the Cauchy mutation operator to obtain the fitness value of the DNA sequence after mutation. If the fitness value after mutation is greater than the fitness value before mutation, the DNA sequence before mutation is replaced by the DNA sequence after mutation;

[0009] Step 5: Using the DNA sequence that has undergone the Cauchy mutation as the input of the Levy mutation, and using the Levy mutation operator to update the DNA sequence to obtain a fitness value after the mutation. If the fitness value after the mutation is greater than the fitness value before the mutation, the DNA sequence before the mutation is replaced by the DNA sequence after the mutation;

[0010] Step 6: Add the sequences in the candidate coding set that meet the constraints to the DNA coding set;

[0011] Step 7: Determine whether the maximum number of iterations has been reached. If so, proceed to step 8, otherwise return to step 2;

[0012] Step 8: Count all newly added DNA sequences and output the DNA coding set.

[0013] Furthermore, the gradient algorithm based on mutation strategy uses a gradient search rule to update the DNA sequence, as shown in formula (1);

[0014] x n m+1 =r a ×(r b ×X1 n m +(1-r b )×X2 n m )+(1-r a )×X3 n m (1)

[0015] Among them, r a 、r b Is a random number between 0 and 1; parameter X1 n m and X2 n m It is obtained using the following formula (2), parameter X3 n m It is obtained using the following formula (3);

[0016]

[0017] X3 n m =x nm -ρ1×(X2 n m -X1 n m ) (3)

[0018] Among them, x n m represents the current DNA sequence, randn is a random number that conforms to the normal distribution, ε is a number between 0 and 0.1, and x best and x worst In the iterative process, the DNA sequences with the best and worst fitness values ​​are obtained. rand represents a random number between 0 and 1, ρ1 and ρ2 are equilibrium parameters, and Δx is the difference value.

[0019] Furthermore, the balance parameters are obtained by the following formulas (4)-(6):

[0020] ρ=2×rand×α-α (4)

[0021]

[0022]

[0023] Among them, β min and β max They are 0.2 and 1.2 respectively, m is the current iteration number, and M is the total number of iterations.

[0024] Furthermore, the difference value is obtained by the following formulas (7)-(9):

[0025] Δx=rand(1:N)×|step| (7)

[0026]

[0027]

[0028] Among them, rand(1:N) is an N-dimensional random number, r1, r2, r3 and r4 (r1≠r2≠r3≠r4) are different random integers from 1 to N, x r1 m 、x r2 m 、x r3 m 、x r4 m are four randomly selected DNA sequences.

[0029] Furthermore, when the DNA sequence is updated again through the local escape operator, if the other random number is less than the threshold, formula (10) is executed, otherwise formula (11) is executed:

[0030] x n m+1 =x n m+1 +f1×(u1×x best -u2×x k m )+f2×ρ1×(u3×(X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2 (10)

[0031] x n m+1 =x best +f1×(u1×x best -u2×x k m )+f2×ρ1×(u3×(X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2 (11)

[0032] Among them, f1 is a uniform random number from -1 to 1, f2 is a random number that conforms to the normal distribution, and u1, u2, and u3 are three random numbers.

[0033] Furthermore, the updated DNA sequence is used as the input of the Cauchy mutation, and the DNA sequence is perturbed by the Cauchy mutation operator, specifically:

[0034] X new (i)=X(i)+X(i)×C(γ) i∈{1,...,N} (12)

[0035] Among them, X(i) is the updated DNA sequence, and C(γ) is the Cauchy mutation operator.

[0036] Furthermore, the Cauchy mutation operator is obtained according to the following formulas (13) and (14):

[0037]

[0038]

[0039] As a further step, the DNA sequence that has undergone Cauchy mutation is used as the input of Levy mutation, and the DNA sequence is updated using the Levy mutation operator.

[0040] X new (i)=X(i)+X(i)×L(β) i∈{1,...,N} (15)

[0041] Among them, X(i) is the DNA sequence that has undergone Cauchy mutation; L(β) is the random distribution number obtained from the Levy distribution, and the expression of the random distribution number is as follows:

[0042]

[0043]

[0044] in, is the value obtained by the sine function and the gamma function; μ and v are both random numbers that conform to the normal distribution, β is a constant, such as 1.5; Γ is the gamma function.

[0045] Furthermore, the constraint conditions include a full discontinuity constraint, a GC content constraint and a no adjacent subsequence constraint; wherein the full discontinuity constraint indicates that the same bases cannot appear continuously in a DNA sequence; the GC content constraint refers to the percentage of bases G and C in the entire DNA sequence in a DNA sequence; the no adjacent subsequence constraint indicates that the same subsequence cannot appear repeatedly in a DNA sequence; the Hamming distance constraint is used as the objective function, and the Hamming distance constraint is expressed as the number of different elements at the same position between two sequences.

[0046] The present invention adopts the above technical solution, and can achieve the following technical effects:

[0047] 1. The gradient algorithm based on mutation strategy mainly determines the update direction through Newton's method, and uses gradient search rules and local escape operators to explore the DNA sequence update domain. The gradient search rule uses a gradient-based method to enhance the algorithm's exploration ability and speed up convergence, so that better sequences can be obtained in the entire update space. The local escape operator is mainly to improve the efficiency of the algorithm in solving complex problems and escape from local optimality. These two main operators are conducive to the overall grasp of the data, allowing the DNA sequence to evolve in a favorable direction to obtain an excellent DNA coding set;

[0048] 2. Using the Cauchy mutation operator for random perturbation is conducive to increasing the diversity of DNA sequences, making the early iterative exploration range wider and more inclined to promising areas; Levy flight has a strong perturbation ability. It is a movement mode that alternates between high-frequency short-distance exploration and low-frequency long-distance exploration. It can not only expand the search range, but also enhance the ability of local search in specific areas. It can also avoid falling into the local optimum when searching for the optimal solution in a large range;

[0049] 3. Based on the problem that DNA sequences containing continuous repeated subsequences are more prone to errors during storage, a new non-adjacent subsequence constraint is proposed to limit the occurrence of such sequences, which can improve the robustness of the entire DNA coding set;

[0050] 4. The combination of gradient algorithm based on mutation strategy and constraint conditions can construct a DNA coding set with better quantity and higher robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 Flowchart of the optimization method for DNA storage encoding with high robustness. DETAILED DESCRIPTION

[0052] The following will be combined with the accompanying drawings in the present invention to clearly and completely describe the technical solutions in the implementation of the present invention. It can be understood that the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making creative work are within the protection scope of the present invention.

[0053] In addition to the three classic constraints (full discontinuity constraint, Hamming distance, GC content), the constraints involved in the present invention also include a newly proposed constraint (no adjacent subsequences). Among them, the full discontinuity constraint means that the same bases cannot appear continuously in a DNA sequence; the Hamming distance is expressed as the number of different elements at the same position of two sequences; the GC content refers to the percentage of bases in the entire DNA sequence occupied by bases G and C in a DNA sequence, and the GC content in this embodiment is 50%; the no adjacent subsequence constraint means that the same subsequence cannot appear repeatedly in a DNA sequence, and the present invention mainly focuses on the case where the number of subsequence bases is 2 or 3. The Hamming distance in the constraint is used as the objective function, and the other three are used as constraints.

[0054] Example 1

[0055] The embodiments of the present invention are implemented on the premise of the technical solution of the present invention, and detailed implementation methods and specific operation processes are given, but the protection scope of the present invention is not limited to the following embodiments. In the embodiments, the DNA encoding length n is 6, the Hamming distance constraint is d≥4, and the full discontinuity constraint, GC content constraint and no adjacent subsequence constraint are as described above.

[0056] Step 1: Initialize the parameters, initial set and DNA encoding set required by the gradient algorithm based on the mutation strategy;

[0057] Specifically, the DNA coding set can be 500 DNA coding sequences with a length of 6, and the 500 DNA coding sequences with a length of 6 are screened using the constraint condition to obtain the initial set. In this embodiment, a simulation experiment is performed using MATLAB, and the initial set is 29 obtained by using the GC content, full discontinuity, and no adjacent subsequence constraints.

[0058] Step 2: Obtain the fitness value of the initial set, sort them according to the fitness value, and take the DNA sequence with the current best fitness as the optimal solution;

[0059] Specifically, the fitness of the initial set is obtained using the fitness function and the fitness of the DNA sequences is sorted using the sort() function to select the DNA sequences with the best fitness and the worst fitness;

[0060] Step 3: Each DNA sequence in the initial set generates a random number r. When the random number r is greater than or equal to the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. When the random number r is less than the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. The DNA sequence is then updated through the local escape operator to obtain the fitness value of the corresponding DNA sequence.

[0061] Specifically, the 29 6-dimensional DNA sequences obtained are further optimized using the gradient algorithm. For each DNA sequence, a random number r is generated, and the size of the random number r and the threshold pr (pr can be 0.5) is determined:

[0062] 1) If r is greater than or equal to pr, only the gradient search rule is used to update the position, as shown in formula (1);

[0063] x n m+1 =r a ×(r b ×X1 n m +(1-r b )×X2 n m )+(1-ra )×X3 n m (1)

[0064] X1 n m and X2 n m It is obtained using the following formula (2), X3 n m It is obtained using the following formula (3).

[0065]

[0066] X3 n m =x n m -ρ1×(X2 n m -X1 n m ) (3)

[0067] where x n m represents the current DNA sequence, randn is a random number that conforms to the normal distribution, ε is a number between 0 and 0.1, and x best and x worst In the iterative process, the DNA sequences with the best and worst fitness values ​​are obtained. rand represents a random number between 0 and 1. ρ1 and ρ2 can be obtained by the following formulas (4)-(6), and Δx can be obtained by the following formulas (7)-(9).

[0068] ρ=2×rand×α-α (4)

[0069]

[0070]

[0071] where β min and β max They are 0.2 and 1.2 respectively, m is the current iteration number, and M is the total number of iterations.

[0072] Δx=rand(1:N)×|step| (7)

[0073]

[0074]

[0075] Where rand(1:N) is an N-dimensional random number, r1, r2, r3, and r4 (r1≠r2≠r3≠r4) are different random integers from 1 to N, and x r1 m 、x r2 m 、x r3 m 、x r4 m are four randomly selected DNA sequences.

[0076] 2) If r is less than pr, first use the gradient search rule to update the sequence, as shown above, and then use the local escape operator to update the sequence again. If the other random number is less than the threshold (the threshold can be 0.5), execute formula (10), otherwise execute formula (11);

[0077] x n m+1 =x n m+1 +f1×(u1×x best -u2×x k m )+f2×ρ1×(u3×(X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2 (10)

[0078] x n m+1 =x best +f1×(u1×x best -u2×x k m )+f2×ρ1×(u3×(X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2 (11)

[0079] Among them, f1 is a uniform random number from -1 to 1, f2 is a random number that conforms to the normal distribution, and u1, u2, and u3 are three random numbers.

[0080] Step 4: Using the updated DNA sequence as the input of the Cauchy mutation, the DNA sequence is perturbed by the Cauchy mutation operator to obtain the fitness value of the DNA sequence after mutation. If the fitness value after mutation is greater than the fitness value before mutation, the DNA sequence before mutation is replaced by the DNA sequence after mutation;

[0081] Specifically, the DNA sequence is disturbed by the Cauchy mutation operator, specifically:

[0082] X new (i)=X(i)+X(i)×C(γ) i∈{1,...,N} (12)

[0083] The Cauchy mutation operator C(γ) can be obtained according to the following formulas 13 and 14.

[0084]

[0085]

[0086] Step 5: Using the DNA sequence that has undergone the Cauchy mutation as the input of the Levy mutation, and using the Levy mutation operator to update the DNA sequence to obtain a fitness value after the mutation. If the fitness value after the mutation is greater than the fitness value before the mutation, the DNA sequence before the mutation is replaced by the DNA sequence after the mutation;

[0087] Specifically, the DNA sequence is updated using the Levy mutation operator, specifically:

[0088] X new (i)=X(i)+X(i)×L(β) i∈{1,...,N} (15)

[0089] Where L(β) is a random distribution number obtained from the Levy distribution. The expression of the Levy random number is as follows:

[0090]

[0091]

[0092] Step 6: Add the sequences in the initial set that meet the constraints to the DNA coding set;

[0093] Step 7: Determine whether the maximum number of iterations has been reached. If so, proceed to step 8, otherwise return to step 2;

[0094] Specifically, determine whether the maximum number of iterations, 1000 generations, is reached. If so, proceed to step 8, otherwise return to step 2;

[0095] Step 8: Count all newly added DNA sequences and output the DNA coding set.

[0096] The present invention proposes a highly robust DNA storage coding optimization method based on a gradient algorithm of a mutation strategy. First, the gradient algorithm is used to continuously perform iterative updates, and then the Cauchy and Levy mutation operators are used to continuously update the DNA sequence to the optimal DNA sequence, so as to obtain more sequences that meet the constraints of GC content, full discontinuity and no adjacent subsequences, and finally the obtained maximum DNA coding set is used as the output result. The present invention uses MATLAB to perform simulation experiments on the algorithm under the operating environment of Intel (R) CPU 3.6GHz, 10.0GB memory, and Windows 10. The experimental results show that the method results of this embodiment are better than the experimental results of other algorithms.

[0097] Table 1 Initial DNA sequence

[0098]

[0099] Table 2 DNA encoding set (n=6, d≥4)

[0100]

[0101] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. It should be pointed out that ordinary technicians in this technical field can make several improvements and modifications without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A DNA storage coding optimization method with high robustness, characterized in that: include: Step 1: Initialize the parameters, initial set and DNA encoding set required by the gradient algorithm based on the mutation strategy; Step 2: Obtain the fitness value of the initial set, sort them according to the fitness value, and take the DNA sequence with the current best fitness as the optimal solution; Step 3: Each DNA sequence in the initial set generates a random number r. When the random number r is greater than or equal to the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. When the random number r is less than the threshold pr, the gradient algorithm based on the mutation strategy uses the gradient search rule to update the DNA sequence. The DNA sequence is then updated through the local escape operator to obtain the fitness value of the corresponding DNA sequence. Step 4: Using the updated DNA sequence as the input of the Cauchy mutation, the DNA sequence is perturbed by the Cauchy mutation operator to obtain the fitness value of the DNA sequence after mutation. If the fitness value after mutation is greater than the fitness value before mutation, the DNA sequence before mutation is replaced by the DNA sequence after mutation; Step 5: Using the DNA sequence that has undergone the Cauchy mutation as the input of the Levy mutation, and using the Levy mutation operator to update the DNA sequence to obtain a fitness value after the mutation. If the fitness value after the mutation is greater than the fitness value before the mutation, the DNA sequence before the mutation is replaced by the DNA sequence after the mutation; Step 6: Add the sequences in the initial set that meet the constraints to the DNA coding set; Step 7: Determine whether the maximum number of iterations has been reached. If so, proceed to step 8, otherwise return to step 2; Step 8: Count all newly added DNA sequences and output the DNA coding set; The gradient algorithm based on mutation strategy uses gradient search rules to update the DNA sequence, as shown in formula (1); x n m+1 =r a ×(r b ×X1 n m +(1-r b )×X2 n m )+(1-r a )×X3 n m (1) Among them, r a 、r b Is a random number between 0 and 1; parameter X1 n m and X2 n m It is obtained using the following formula (2), parameter X3 n m It is obtained using the following formula (3); X3 n m =x n m -ρ1×(X2 n m -X1 n m ) (3) Among them, x n m represents the current DNA sequence, randn is a random number that conforms to the normal distribution, ε is a number between 0 and 0.1, and x best and x worst In the iterative process, the DNA sequences with the best and worst fitness values ​​are obtained, rand represents a random number between 0 and 1, ρ1 and ρ2 are equilibrium parameters, and Δx is the difference value; When the DNA sequence is updated again through the local escape operator, if the other random number is less than the threshold, formula (10) is executed, otherwise formula (11) is executed: x n m+1 =x n m+1 +f1 × (u1 × x best −u2×x k m ) + f2 × ρ1 × (u3 × (X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2(10) x n m+1 =x best +f1 × (u1 × x best −u2×x k m ) + f2 × ρ1 × (u3 × (X2 n m -X1 n m )+u2×(x r1 m -x r2 m )) / 2(11) Among them, f1 is a uniform random number from -1 to 1, f2 is a random number that conforms to the normal distribution, and u1, u2, and u3 are three random numbers; The constraint conditions include a full discontinuity constraint, a GC content constraint and a no adjacent subsequence constraint; wherein the full discontinuity constraint indicates that the same bases cannot appear continuously in a DNA sequence; the GC content constraint refers to the percentage of bases G and C in the entire DNA sequence in a DNA sequence; the no adjacent subsequence constraint indicates that the same subsequence cannot appear repeatedly in a DNA sequence; the Hamming distance constraint is used as the objective function, and the Hamming distance constraint is expressed as the number of different elements at the same position between two sequences.

2. The DNA storage coding optimization method with high robustness according to claim 1, characterized in that: The balance parameters are obtained by the following formulas (4)-(6): ρ=2×rand×α-α (4) Among them, β min and β max They are 0.2 and 1.2 respectively, m is the current iteration number, and M is the total number of iterations.

3. The DNA storage coding optimization method with high robustness according to claim 1, characterized in that: The difference value is obtained by the following formulas (7)-(9): Δx=rand(1:N)×|step| (7) Among them, rand(1:N) is an N-dimensional random number, r1, r2, r3 and r4 (r1≠r2≠r3≠r4) are different random integers from 1 to N, x r1 m 、x r2 m 、x r3 m 、x r4 m are four randomly selected DNA sequences.

4. The DNA storage coding optimization method with high robustness according to claim 1, characterized in that: The updated DNA sequence is used as the input of the Cauchy mutation, and the DNA sequence is perturbed by the Cauchy mutation operator, specifically: X new (i)=X(i)+X(i)×C(γ) i∈{1,...,N} (12) Among them, X(i) is the updated DNA sequence, and C(γ) is the Cauchy mutation operator.

5. The DNA storage coding optimization method with high robustness according to claim 4, characterized in that: The Cauchy mutation operator is obtained according to the following formulas (13) and (14): 。 6. The DNA storage coding optimization method with high robustness according to claim 1, characterized in that: The DNA sequence that has undergone Cauchy mutation is used as the input of Levy mutation, and the DNA sequence is updated using the Levy mutation operator. X new (i)=X(i)+X(i)×L(β) i∈{1,...,N} (15) Among them, X(i) is the DNA sequence that has undergone Cauchy mutation; L(β) is the random distribution number obtained from the Levy distribution, and the expression of the random distribution number is as follows: in, is the value obtained by the sine function and the gamma function; μ and v are both random numbers that conform to the normal distribution, β is a constant; Γ is the gamma function.

Images