An optimization method for sequence design

Abstract

The present invention discloses an optimization method for sequence design, which realizes fast generation of optimization results by introducing matrix tensors updated iteratively through backpropagation, has a fast convergence speed, and at the same time updates the neural network model through multiple rounds of iteration, uses the matrix tensors to guide the generation direction of the sequence, provides the gradients required for backpropagation of the matrix tensors for generating the sequence, and has higher optimization efficiency.

Description

Technical Field

[0001] The present invention relates to the technical field of gene sequence optimization, and more specifically, to an optimization method for sequence design. Background Art

[0002] Sequence optimization calculation is based on the interaction between sequence and function to design a suitable sequence and meet the optimization goal. Sequence optimization, such as the optimization of nucleotide and amino acid sequences, essentially expects to change several minimum constituent units in the sequence to enhance one or more specific capabilities of the sequence, such as enhancing the translation efficiency and stability of mRNA, enhancing the binding ability of proteins, etc.

[0003] In the prior art, there is a CDS (Coding Sequence) sequence corresponding to the optimal CAI (Codon Adaption Index) obtained based on codon adaptability calculation as the starting ancestor sequence of the genetic algorithm, and through multiple rounds of iteration to obtain a synonymous substitution CDS sequence with a lower minimum free energy; there is a method of designing eigenvalue such as base statistical characteristics, codon frequency, and segment free energy through feature engineering, then determining the protein expression label of the corresponding mRNA sample, establishing a relevant model to screen the corresponding features related to high / low expression, and performing selective directional mutation iteration on the 5'UTR sequence to obtain an optimized sequence (genetic algorithm); there is a method of splitting a long sequence into smaller units DFA, linearly optimizing the sequence from left to right based on the minimum free energy and CAI, and finally using beam search to reduce the influence of local optimum; there is a clustering algorithm (K-means / DBSCAN) that improves the genetic algorithm, and after designing a protein-DNA scoring function, more efficient iterative screening and optimization are realized; by replacing the 5' and 3' parts of the sequence and using known experience to design a cap structure, an RNA molecule that enhances a given sequence is synthesized.

[0004] The following problems still need to be solved in the above prior art:

[0005] (1) Most use genetic algorithms, with limited optimization efficiency. There is a large amount of randomness in each round of mutation, and the iterative convergence ability is poor;

[0006] (2) Some methods consider the optimization efficiency of long sequences and use the method of local optimization followed by global search to determine the result, but they cannot completely avoid falling into local optimal solutions;

[0007] (3) The reason why most methods use genetic algorithms for optimization is that there is a lack of an index calculation method with derivable results, and it is impossible to directly use the method of backpropagation by differentiation to update and optimize globally from the sequence;

[0008] (4) The disclosed method is based on empirical knowledge or existing data information, and there is no optimization method combined with experiments for unknown data results such as high-throughput experiment results. Summary of the Invention

[0009] The present invention provides an optimization method for sequence design, which solves the technical problem of limited optimization efficiency in the prior art when using genetic algorithms.

[0010] To solve the above technical problems, the technical solution of the present invention is as follows:

[0011] The present invention provides an optimization method for sequence design, including the following steps:

[0012] S1: Given the coding sequence of a preset gene as the original sequence, confirm the minimum unit of the original sequence and the number of types of the minimum unit, and encode the original sequence according to the minimum unit of the original sequence and the number of types of the minimum unit to obtain the encoding of the original sequence;

[0013] S2: According to a preset synonymous comparison table, perform random synonymous substitution on the minimum units at different positions in the original sequence to obtain a number of sequence variants, and encode the sequence variants to obtain the encoding of the sequence variants;

[0014] S3: Analyze all sequence variants to obtain target values respectively;

[0015] S4: Use the encoding of the sequence variants as input and the target values as output to train a neural network model;

[0016] S5: Establish a matrix tensor with a shape dimension of b*l*c according to the coding sequence of the current gene, where b is the number of generated sequence variants, l is the length of the original sequence, and c is the data dimension of the encoding of the sequence variants;

[0017] S6: Obtain the encoding of the input sequence according to the matrix tensor;

[0018] S7: Input the encoding of the input sequence into the neural network model to obtain the predicted target value;

[0019] S8: Calculate the loss function according to the predicted target value, and update the matrix tensor through backpropagation to make the loss function tend to 0. Count the number of iterations. If the number of iterations is less than the first preset number of times, return to step S6; if the number of iterations is not less than the first preset number of times, enter step S9;

[0020] S9: Sort the predicted target values of all the optimized sequences generated during the iteration process, select the top N optimized sequences, analyze them separately to obtain the target values, and then sort the top N optimized sequences according to the target values obtained from the re-analysis to obtain the optimized sequence corresponding to the optimal target value;

[0021] S10: If the number of loops is less than the second preset number, replace the original sequence with the optimized sequence corresponding to the optimal target value, and return to step S1; if the number of loops is not less than the second preset number, output the optimized sequence corresponding to the optimal target value in the current loop.

[0022] In the above technical means, by introducing the matrix tensor updated by backpropagation iteration, the optimization result can be quickly generated with a fast convergence speed. At the same time, by iteratively updating the neural network model for multiple rounds, the matrix tensor is used to guide the generation direction of the sequence, providing the gradient required for the backpropagation of the matrix tensor for the generated sequence, and the optimization efficiency is higher.

[0023] Further, in step S1, the encoding of the original sequence is in one-hot form or categorical number form.

[0024] Further, step S6 obtains the encoding of the input sequence according to the matrix tensor, including the following steps:

[0025] S6.1: Multiply the matrix tensor by the synonym replacement mask weight W(l*c), where the value W i is the selection in the synonym comparison table corresponding to the smallest unit of the i-th in the original sequence, with synonym replacement being 1 and non-synonym replacement being 0. After that, perform normalization conversion on the matrix tensor:

[0026]

[0027] In the formula, T i is the i-th dimension of the data dimension c in the third dimension of the matrix tensor T;

[0028] S6.2: Calculate the index value where the maximum value is located in sequence according to the data dimension c in the second dimension of the converted matrix tensor T, and encode the sequence variant corresponding to the maximum value to obtain the encoding of the input sequence.

[0029] Further, in step S6.2, a reparameterization operation is also included for the encoding of the input sequence.

[0030] Further, the loss function in step S8 includes:

[0031] Lm = -α / s

[0032] In the formula, Lm is the loss function, α is a constant, and s is the predicted target value.

[0033] Further, after sorting the first N optimized sequences according to the target value obtained from the re - analysis in step S9, the following steps are also included:

[0034] Select the optimized sequences with the top several target values before sorting, and count the selection situation Dis(l*c) of the minimum units at each position in each optimized sequence. The specific value Dis i is the selection probability of each minimum unit at the i - th position.

[0035] Further, in the loop of steps S1 to S9, except for the first loop, each subsequent loop samples and generates data for training the neural network model in the current loop based on the optimized sequence corresponding to the optimal target value obtained in the previous loop, all the generated optimized sequences, and the corresponding Dis(l*c).

[0036] Further, the process of sampling and generating data for training the neural network model in the current loop based on the optimized sequence corresponding to the optimal target value obtained in the previous loop, all the generated optimized sequences, and the corresponding Dis(l*c) includes:

[0037] Sampling and generating process:

[0038] Perform random synonymous substitution on the optimized sequence corresponding to the optimal target value obtained in the previous loop to generate sequence variants;

[0039] Select the top m optimized sequences with the lowest target values from all the optimized sequences generated in the previous loop;

[0040] According to the Dis(l*c) obtained in the previous loop, obtain the probability of random synonymous substitution for each position of the optimized sequence from Dis and perform distribution - based random substitution;

[0041] Combine the sequences generated in the sampling and generating process in a ratio of 2:4:4 for training the neural network model in the current loop.

[0042] Further, when the original sequence in step S1 is not limited to genes, in step S2, randomly substitute the minimum units at different positions in the original sequence, and in step S3, analyze the target values of all sequence variants generated through high - throughput experiments.

[0043] Further, when there are multiple target values, in step S3, all sequence variants are analyzed to obtain one of the target values respectively. In step S8, the loss function is calculated according to the predicted target value. It further includes inputting the encoding of the sequence variant into a pre-trained RPF prediction neural network model for prediction to obtain a second loss function. At the same time, other target values of the sequence variant are calculated to obtain corresponding loss functions. All the loss functions are added together according to different weights to obtain the total loss function, and the matrix tensor is updated by backpropagation using the total loss function.

[0044] Compared with the prior art, the beneficial effects of the technical solution of the present invention are as follows:

[0045] (1) Compared with the low efficiency of the genetic algorithm, the iterative optimization method based on backpropagation can quickly iteratively generate a sufficient number of optimization results with a fast convergence speed.

[0046] (2) Compared with the method of using local optimization to equivalent global optimization in some methods, in the process of sequence generation of the present invention, except for the necessary constraints, the generated space is not limited. The advantage of this is that there is a large enough generated space to provide optimization and iteration. Compared with local optimization, it is easier to find the theoretically global optimal solution.

[0047] (3) Based on the establishment of the neural network model, the present invention updates the neural network model through multiple rounds of iteration, and at the same time guides the generation direction of the sequence, provides the gradient required for backpropagation of the updatable matrix of the generated sequence, and can fit the performance of existing tools. If the existing tools are based on deep learning modeling and the input meets the conditions, it can be directly used to guide the optimization of the sequence with higher efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a schematic flowchart of an optimization method for sequence design provided by an embodiment of the present invention;

[0049] Figure 2 Provided by an embodiment of the present invention Figure 1 The flowchart of the method shown;

[0050] Figure 3 Based on two genes, Gluc and eGFP, provided by an embodiment of the present invention, using Figure 1 and Figure 2 The schematic diagram of the optimized MFE results of the method described and the existing method;

[0051] Figure 4 It is a schematic flowchart of the sampling generation process provided by an embodiment of the present invention;

[0052] Figure 5 It is a schematic flowchart of another optimization method for sequence design provided by an embodiment of the present invention;

[0053] Figure 6 provided for the embodiments of the present invention Figure 5 flow chart of the method shown

[0054] Figure 7 provided for the embodiments of the present invention using Figure 5 and Figure 6 schematic diagram of the result of the fitting trend line of the experimental measurement values by the method described

[0055] Figure 8 flow schematic diagram of the optimization method of the third sequence design provided for the embodiments of the present invention

[0056] Figure 9 provided for the embodiments of the present invention Figure 8 flow chart of the method shown Detailed implementation manners

[0057] The drawings are only for illustrative purposes and should not be construed as a limitation of this patent;

[0058] To better illustrate this embodiment, some components in the drawings are omitted, enlarged or reduced, and do not represent the size of the actual product;

[0059] For those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted.

[0060] The technical solutions of the present invention will be further described below with reference to the drawings and embodiments.

[0061] Embodiment

[0062] As Figure 1 and Figure 2 shown, Figure 1 is a flow schematic diagram of an optimization method for sequence design. The optimization method for sequence design uses the computable interesting characteristics in the bioinformatics sequence as the guiding direction, and generates a series of sequences optimized for one characteristic through an iterative optimization method, including the following steps:

[0063] S01: Given the coding sequence of a preset gene as the original sequence, confirm the minimum unit of the original sequence and the number of types of the minimum unit, and encode the original sequence according to the minimum unit of the original sequence and the number of types of the minimum unit to obtain the encoding of the original sequence;

[0064] In a specific embodiment, taking the optimization of the CDS (Coding sequence) of an mRNA sequence as an example, the goal is to obtain a sequence with the same functional expression and a sufficiently low minimum free energy (MFE). MFE is an index used to evaluate the stability of the mRNA coding region (CDS) and is an important parameter used in RNA secondary structure prediction to represent the lowest free energy of an RNA molecule, which can maintain the most stable structure even in different structural states. A lower MFE value indicates a more stable secondary structure of the RNA molecule, while a higher MFE value indicates a less stable secondary structure of the RNA molecule.

[0065] In step S01, the CDS sequence (ATGGGAATG…) of a specific gene is used as the original sequence. At this time, the minimum unit is the codon of the CDS. Each codon consists of three bases, and each base has 4 options (A / T / G / C). Then, there are 4*4*4 = 64 combinations of codons, and this CDS sequence can be encoded in one-hot form. Or in the form of category numbers [1, 0, 63, 4, …].

[0066] S02: According to a preset synonymous comparison table, randomly perform synonymous substitutions on the minimum units at different positions in the original sequence to obtain a number of sequence variants, and encode the sequence variants to obtain the encoding of the sequence variants.

[0067] In a specific embodiment, the rule for generating sequence variants is not to change the functional expression of the current sequence. Therefore, it is necessary to randomly replace the minimum unit codons with the same synonymous codons. According to existing knowledge, a synonymous codon comparison table can be determined to randomly perform synonymous substitutions on the codons at different positions in the original sequence.

[0068] S03: Analyze all sequence variants to obtain the target values respectively.

[0069] In a specific embodiment, the obtained sequence variants are analyzed using a tool software. The goal of this embodiment is to optimize the minimum free energy of the sequence. Therefore, the RNAFold tool is used to calculate the minimum free energy values of all generated sequence variants.

[0070] S04: Use the encoding of the sequence variants as the input and the target value as the output to train a neural network model.

[0071] In a specific embodiment, the input of the neural network model is set as the encoding of the sequence variants, and the output of the neural network model is the free energy value calculated using the RNAFold tool. The neural network model is trained to obtain the neural network model Mt.

[0072] S05: Establish a matrix tensor with a shape dimension of b*l*c according to the coding sequence of the current gene, where b is the number of generated sequence variants, l is the length of the original sequence, and c is the data dimension encoded by the sequence variant;

[0073] In a specific embodiment, when using one-hot encoding, c is 65, and when using other encodings, c is 1. This embodiment will be described by taking c = 65 as an example;

[0074] S06: Obtain the encoding of the input sequence according to the matrix tensor;

[0075] In a specific embodiment, in order to ensure that the matrix can generate synonymous substitution sequences of the CDS of the current gene, multiply the matrix tensor T by a synonymous substitution mask weight W(l*c), where l is the codon coding length of the original sequence and c is the number of codon coding types. The specific value of W i is the selection [0, 0, 1, 0, 1,..., 0] in the synonymous substitution list corresponding to the i-th codon in the original sequence. There are 65 selections, corresponding to the number of codon types. If it is a synonymous substitution, it is 1, and if it is not a synonymous substitution, it is 0. The selections at all positions i in the sequence are the synonymous substitution mask weight W. Then, perform a normalization conversion on the matrix tensor T of floating-point type:

[0076]

[0077] In the formula, T i is the i-th dimension of the data dimension in the third dimension c of the matrix tensor T;

[0078] Calculate the index value where the maximum value is located in sequence according to the data dimension c in the second dimension of the converted matrix tensor T, and encode the sequence variant corresponding to the maximum value to obtain the encoding of the input sequence:

[0079] T input = Onehot(argmax(T, axis = 1))

[0080] To avoid the problem of gradient loss caused by the non-differentiability of the Onehot and argmax operations, perform the following reparameterization operation on T input and add the information in the selection probability tensor T:

[0081] T input = T input - T.detach() + T

[0082] S07: Input the encoding of the input sequence into the neural network model to obtain the predicted target value s = Mt(T input );

[0083] S08: Calculate the loss function based on the predicted target value, and update the matrix tensor through backpropagation so that the loss function tends to 0. Count the number of iterations. If the number of iterations is less than the first preset number of times, return to step S06; if the number of iterations is not less than the first preset number of times, proceed to step S09;

[0084] In a specific embodiment, the loss function includes:

[0085] Lm = -α / s

[0086] In the formula, Lm is the loss function, α is a constant, which can be 50, 100, 1000, etc., and s is the predicted target value;

[0087] S09: Sort the predicted target values of all the optimization sequences generated during the iteration process, select the top N optimization sequences and analyze them again to obtain the target values respectively, and then sort the top N optimization sequences according to the target values obtained from the re-analysis to obtain the optimization sequence corresponding to the optimal target value;

[0088] In this embodiment, through backpropagation, by updating the tensor T, Lm tends to 0, and at the same time the MFE predicted value will also tend to a lower value. After iterating a certain number of rounds n, n*b synonymous substitution variant sequences of the current gene will be generated. Sort the variant sequences according to their predicted MFE values, select the top N optimal ones, and then use the rnafold tool to analyze them again to obtain the MFE value, and sort them again to obtain the sequence Seq_new with the lowest free energy;

[0089] S010: If the number of loops is less than the second preset number of times, replace the original sequence with the optimization sequence Seq_new corresponding to the optimal target value, and return to step S01; if the number of loops is not less than the second preset number of times, output the optimization sequence corresponding to the optimal target value in the current loop.

[0090] As Figure 3 shown, based on the two genes Gluc and eGFP, compare the results of optimizing MFE of the method in this embodiment with the existing method, Figure 3 In the left figure, the result of optimizing MFE of the Gluc gene is shown. The upper straight line is the MFE value before optimization, the lower straight line is the optimization result of the existing method, and the descending curve is the optimization result of the method in this embodiment, Figure 3 In the right figure, the result of optimizing MFE of the eGFP gene is shown. The upper straight line is the MFE value before optimization, the lower straight line is the optimization result of the existing method, and the descending curve is the optimization result of the method in this embodiment. It can be seen that after the number of iterations exceeds 5 times, the optimization result of the method in this embodiment is significantly better than that of the existing method.

[0091] This embodiment also provides a method for sampling from the iterative optimization results to accelerate the iterative process, such as Figure 4 shown, including:

[0092] (1) Taking the CDS optimization of mRNA as an example, according to the number of iterative rounds, if it is the first round of optimization, the sampling generation process of this round is: perform random synonymous substitution on each position of the original sequence according to the minimum unit codon, and the substitution probability dynamically changes from 0.1 to 0.5. After determining the substitution, the probability of substituting for a specific codon is the reciprocal of the number of optional substitutions of the current codon;

[0093] (2) If it is not the first round of optimization, the sampling generation process of this round is divided into three parts:

[0094] ① First, perform random synonymous substitution to generate variant sequences, the same as (1) above;

[0095] ② For the optimized variant sequences generated in the previous round of iteration, after sorting according to the MFE prediction values, select the lowest m ones;

[0096] ③ For the codon distribution Dis obtained in the previous round of iteration, according to this distribution, the probability of performing random synonymous substitution on each position of the sequence is obtained from Dis for distribution-based random substitution, and random generation adds the directionality of optimization;

[0097] (3) Combine the sequences of the above three parts in a ratio of 2:4:4, and train the target model Mt of the current round based on these variant sequences.

[0098] Such as Figure 5 and Figure 6 shown, Figure 5 is a flow schematic diagram of another sequence design optimization method. For sequence optimization that cannot be calculated by existing tool software and can only be verified through experiments, the method of the present invention can also be used to fit the results of high-throughput experiments in the form of a target model. Taking the optimization task of the expression level of a 110-base-long DNA sequence as an example, it includes the following steps:

[0099] S11: The original sequence does not limit the gene. Confirm the minimum unit of the original sequence and the number of types of the minimum unit, and encode the original sequence according to the minimum unit of the original sequence and the number of types of the minimum unit to obtain the encoding of the original sequence;

[0100] In a specific embodiment, the minimum unit of the DNA sequence is the base bp, and each base has 4 options (A / T / G / C). Taking the ATGC sequence as an example, this sequence can be encoded in one-hot form or in the form of type numbers [0, 1, 2, 3];

[0101] S12: Based on the sequence-based encoding, various variants are generated online, and 4 types of random substitutions are made for the minimum unit base thereof to obtain a large number of sequence variants;

[0102] S13: Perform a round of high-throughput experiments on the obtained sequences to analyze the expression levels of all sequence variants;

[0103] S14: Use the encoding of the sequence variant as the input and the expression level value as the output to train a neural network model;

[0104] In a specific embodiment, let the input of the neural network model be the encoding of the sequence variant, and the output of the neural network model be the expression level measured by high-throughput experiments, and train the neural network model to obtain the neural network model Mt;

[0105] S15: Establish a matrix tensor with a shape dimension of b*l*c according to the coding sequence of the current gene, where b is the number of generated sequence variants, l is the length of the original sequence, and c is the data dimension of the encoding of the sequence variant;

[0106] In a specific embodiment, when the base uses one-hot encoding, c is 4, and this embodiment is described by taking c = 65 as an example;

[0107] S16: Obtain the encoding of the input sequence according to the matrix tensor;

[0108] In a specific embodiment, multiply the tensor T by a selection template weight W(l*c), where l is the base encoding length of the original sequence and c is the number of base types of the codon. The specific value W i is the possible selection [1, 1, 1, 1] corresponding to the base of the i-th original sequence. There are 4 selections, corresponding to the number of base types. If the current position is selectable, it is 1, and if it is not selectable, it is 0. For example, if the current is [1, 1, 1, 1], it means that all four bases at the current position can be selected during optimization. If the current is [1, 0, 0, 0], it means that the base at the current position is fixed as A and cannot be changed during optimization. The selection of all positions i in the sequence is the selection template weight W. Then, perform a normalization conversion on the floating-point type matrix tensor T:

[0109]

[0110] In the formula, T i is the i-th dimension of the data dimension c in the third dimension of the matrix tensor T;

[0111] Calculate the index value where the maximum value is located in sequence according to the data dimension c in the second dimension of the converted matrix tensor T, and obtain the encoding of the input sequence by the sequence variant encoding corresponding to the maximum value:

[0112] T input = Onehot(argmax(T, axis = 1))

[0113] To avoid the problem of gradient loss caused by the non-differentiability of the Onehot and argmax operations, perform the following reparameterization operation on T input to incorporate the information in the selection probability tensor T:

[0114] T input = T input - T.detach() + T

[0115] S17: Input the encoding of the input sequence into the neural network model to obtain the predicted target value s = Mt(T input );

[0116] S18: Calculate the loss function based on the predicted target value s, and through backpropagation, update the matrix tensor to make the loss function tend to 0. Count the number of iterations. If the number of iterations is less than the first preset number, return to step S16; if the number of iterations is not less than the first preset number, proceed to step S19;

[0117] In a specific embodiment, the loss function includes:

[0118] Lm = -α / s

[0119] where Lm is the loss function, α is a constant, which can be 1, 5, 15, 20, etc., and s is the predicted target value;

[0120] S19: Sort the predicted target values of all the optimized sequences generated during the iteration process, select the top N optimized sequences and analyze them again to obtain the target values respectively, and then sort the top N optimized sequences according to the target values obtained from the re-analysis to obtain the optimized sequence corresponding to the optimal target value;

[0121] In a specific embodiment, through backpropagation, by updating the tensor T, Lm tends to 0, and at the same time, the expression quantity prediction value will also tend to a higher value. After iterating a certain number of rounds n, n * b sequences will be generated. Sort the sequences according to the expression quantity prediction value to find the sequence Seq_new with the highest expression prediction value;

[0122] S110: If the number of loops is less than the second preset number, replace the original sequence with the optimized sequence corresponding to the optimal target value, and return to step S11; if the number of loops is not less than the second preset number, output the optimized sequence corresponding to the optimal target value in the current loop;

[0123] In a preferred embodiment, while finding the sequence Seq_new with the highest expression prediction value, the optimal N sequences are selected, and then the base selection situation Dis(l*c) at each position of the sequence is counted. Here, l refers to the codon length of the sequence, and c is the number of base types. The specific value Dis i is the selection probability of each base at the i-th position. For example, [0.01, 0, 0.6, 0.399]. Except for the first round, in each subsequent round, the sampling method as Figure 4 shown is used. Based on the sequence Seq_new of the previous round, the n*b optimized sequences generated and the obtained codon distribution Dis are used for sampling to generate the data required for this round of training Mt, and sampling is performed from the iterative optimization results to accelerate the iterative process.

[0124] As Figure 7 shown, the continuously upward curve is the measured value of the high-throughput experiment, and the curve that first decreases and then increases is the predicted value of the neural network model. The trends of the two are close.

[0125] This embodiment is not limited to bioinformatics tools and known index calculation methods. It can also perform iterative fitting on certain indexes that need to be experimentally tested for unknown methods. Combining with high-throughput experiments, it can well fit the experimental measured values (the DNA expression level optimization experiment in the present invention) only after a few rounds of cycles.

[0126] As Figure 8 and Figure 9 shown, Figure 8 is a schematic flow chart of the optimization method for the third sequence design. While optimizing multiple objectives, taking the optimization of the CDS sequence of mRNA as an example, the optimization objectives are the RPF reflecting the mRNA expression level, the minimum free energy MFE and CSCG reflecting the mRNA stability. RPF (ribosome profiling) is an index used to evaluate the expression level of the coding region (CDS) of mRNA. It reflects the translation activity and protein expression level of this region by measuring the number of mRNA fragments bound to ribosomes; CSCG refers to the Codon Stability Coefficient, which is an index measuring the stability of different codons in mRNA. It is calculated by analyzing a large amount of genomic data and transcriptomic data and is used to describe the contribution degree of different codons to the stability of mRNA. The higher the value of CSCG, the higher the stability of the corresponding codon in mRNA and the smaller the impact on the degradation rate of mRNA. The optimization method for the sequence design of multiple objectives includes the following steps:

[0127] S21: Given the CDS sequence (ATGGGAATG…) of a specific gene as the original sequence, confirm the minimum unit of the original sequence and the number of types of the minimum unit, and encode the original sequence according to the minimum unit of the original sequence and the number of types of the minimum unit to obtain the encoding of the original sequence;

[0128] In a specific embodiment, the minimum unit of the CDS sequence is a codon. Each codon consists of three bases, and each base has 4 options (A / T / G / C). Then there are 4 * 4 * 4 = 64 combinations of codons, and this CDS sequence can be encoded in one - hot form Or in the form of type numbers [1, 0, 63, 4, …];

[0129] S22: According to a preset synonymous comparison table, perform random synonymous substitutions on the minimum units at different positions in the original sequence to obtain a number of sequence variants, and encode the sequence variants to obtain the encoding of the sequence variants;

[0130] In a specific embodiment, the rule for generating sequence variants is not to change the functional expression of the current sequence. Therefore, it is necessary to perform random synonymous substitutions of the same synonymous codons on its minimum unit codons. According to existing knowledge, a synonymous codon comparison table can be determined to perform random synonymous substitutions on the codons at different positions in the original sequence;

[0131] S23: Analyze all sequence variants to obtain target values respectively;

[0132] In a specific embodiment, the obtained sequence variants are analyzed using a tool software. The goal of this embodiment is to optimize the minimum free energy of the sequence. Therefore, the RNAFold tool is used to calculate the minimum free energy values of all generated sequence variants;

[0133] S24: Use the encoding of the sequence variants as the input and the target value as the output to train a neural network model;

[0134] In a specific embodiment, let the input of the neural network model be the encoding of the sequence variants, and the output of the neural network model be the free energy value calculated by the RNAFold tool. Train the neural network model to obtain the neural network model Mt;

[0135] S25: Establish a matrix tensor with a shape dimension of b * l * c according to the coding sequence of the current gene, where b is the number of generated sequence variants, l is the length of the original sequence, and c is the data dimension of the encoding of the sequence variants;

[0136] In a specific embodiment, when using one-hot encoding, c is 65, and when using other encodings, c is 1. This embodiment will be described by taking c = 65 as an example;

[0137] S26: Obtain the encoding of the input sequence according to the matrix tensor;

[0138] In a specific embodiment, in order to ensure that the matrix can generate the synonymous substitution sequence of the CDS of the current gene, we multiply the tensor T by a synonymous substitution mask weight W(l*c), where l is the codon encoding length of the original sequence, and c is the number of codon encoding types. The specific value W_i is the selection [0, 0, 1, 0, 1,..., 0] in the synonymous substitution list corresponding to the codon of the i-th original sequence. There are 65 selections, corresponding to the number of codon types. The selection for synonymous substitution is 1, and the selection for non-synonymous substitution is 0. The selection for all positions i in the sequence is the synonymous substitution mask weight W. Then, the floating-point type matrix tensor T is normalized and transformed:

[0139]

[0140] In the formula, T i is the i-th dimension of the data dimension in the third dimension c of the matrix tensor T;

[0141] Calculate the index value where the maximum value is located in sequence according to the data dimension c in the second dimension of the transformed matrix tensor T, and encode the sequence variant corresponding to the maximum value to obtain the encoding of the input sequence:

[0142] T input = Onehot(argmax(T, axis = 1))

[0143] To avoid the problem of gradient loss caused by the non-differentiability of the Onehot and argmax operations, the following reparameterization operation is performed on T input to add the information in the selection probability tensor T:

[0144] T input = T input - T.detach() + T

[0145] S27: Input the encoding of the input sequence into the neural network model to obtain the predicted target value Sm = Mt(T input ), input the encoding of the sequence variant into the pre-trained RPF prediction neural network model Mr for prediction to obtain the predicted value Sr = Mr(T input ), and simultaneously calculate other target values of the sequence variant; The RPF prediction neural network model Mr is an existing model;

[0146] S28: Calculate the loss function based on the predicted target values Sm, Sr, and other target values respectively, and update the matrix tensor through backpropagation to make the loss function tend to 0. Count the number of iterations. If the number of iterations is less than the first preset number of times, return to step S26; if the number of iterations is not less than the first preset number of times, proceed to step S29;

[0147] In a specific embodiment, the loss function includes:

[0148] Lm = -α / Sm

[0149] Lr = 1 / Sr

[0150] In the formula, Lm is the loss function, α is a constant, and can be 50, 100, 1000, etc. Sm is the predicted target value. Calculate the CSCG value of the current sequence based on the CSCG codon stability coefficient table, and at the same time design the corresponding loss function Lc. Specifically, establish the coefficient matrix W based on the codon stability coefficients statistically obtained before cscg , which is a 65*1 matrix representing the stability coefficient corresponding to each of the 65 codons, T inputi represents the i-th codon of the current sequence, and l represents the length of the codons in the sequence:

[0151]

[0152] Lc = 1 / Sc

[0153] Based on these three goals, design the overall loss:

[0154] L = ω m Lm + ω r Lr + ω c Lc

[0155] Among them, ω m , ω r , ω c are the weights of the three losses respectively, and different weights can be set customarily, such as 0.5, 0.25, 0.25, or can be adaptively updated based on the optimization iteration process;

[0156] S29: Sort the predicted target values of all the optimized sequences generated during the iteration process, select the top N optimized sequences to analyze and obtain the target values respectively, and then sort the top N optimized sequences according to the target values obtained from the re-analysis to obtain the optimized sequence corresponding to the optimal target value;

[0157] In a specific embodiment, through backpropagation, by updating the tensor T, Lm is made to tend to 0, and at the same time, the three target values will change in their respective better directions. After iterating a certain number of rounds n, n*b synonymous substitution variant sequences of the current gene will be generated. The variant sequences are sorted according to the predicted values of their MFE, the optimal N sequences are selected, and then the MFE values are obtained by analyzing them again using the rnafold tool. The N sequences are sorted again according to the MFE values, and the top n sequences are selected. The RPF values of the n sequences are sorted respectively to obtain the best n / 2 sequences. Then, the obtained sequences are sorted based on the CSCG values, and the sequence with the highest CSCG is selected as the optimized sequence Seq_new for this round;

[0158] S210: When the number of loops is less than the second preset number, the optimized sequence corresponding to the optimal target value replaces the original sequence, and step S21 is returned; when the number of loops is not less than the second preset number, the optimized sequence corresponding to the optimal target value in the current loop is output;

[0159] In a preferred embodiment, while finding the sequence Seq_new with the highest expression prediction value, based on the top 10 optimized sequences according to the CSCG value, the codon selection situation Dis(l*c) at each position of the sequence is statistically analyzed, where l refers to the codon length of the sequence and c is the number of codon types. The specific value Dis i is the probability of each codon being selected at the i-th position. For example, [0.01, 0, 0.6, 0, 3, …, 0.05]. Except for the first round, in each subsequent round, the sampling method shown in Figure 4 is used. Based on the sequence Seq_new of the previous round, the n*b optimized sequences generated and the obtained codon distribution Dis are used for sampling to generate the data required for training Mt in this round, and sampling is performed from the iterative optimization results to accelerate the iterative process.

[0160] This embodiment can combine multiple targets and optimize multiple targets of the sequence by setting weights or updating weights adaptively (in this embodiment, the multi-objective optimization of MFE, RPF, and CSCg enhances the stability and expression level of the CDS sequence corresponding to the gene respectively).

[0161] The terms describing the positional relationship in the drawings are only for illustrative purposes and should not be construed as a limitation of this patent;

[0162] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, rather than limitations on the implementation manners of the present invention. For those of ordinary skill in the art, other different forms of changes or alterations can be made based on the above description. It is not necessary and impossible to enumerate all implementation manners here. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the claims of the present invention.

Claims

1. A method for optimizing sequence design, characterized in that: The following steps are involved: S1: Given a coding sequence of a preset gene as an original sequence, confirming the minimum unit of the original sequence and the number of types of the minimum units, and encoding the original sequence according to the minimum unit of the original sequence and the number of types of the minimum units to obtain the encoding of the original sequence; S2: performing random synonymous replacement on the minimum units at different positions in the original sequence according to a preset synonymous comparison table to obtain a number of sequence variants, and encoding the sequence variants to obtain codes of the sequence variants; S3: Analyze all sequence variants and obtain target values ​​respectively; S4: using the encoding of the sequence variant as input and the target value as output, training to obtain a neural network model; S5: Establish a matrix tensor with a shape dimension of b*l*c according to the coding sequence of the current gene, where b is the number of generated sequence variants, l is the length of the original sequence, and c is the data dimension of the encoding of the sequence variant; S6: Obtaining the encoding of the input sequence according to the matrix tensor; S7: Inputting the encoding of the input sequence into the neural network model to obtain a predicted target value; S8: Calculate the loss function according to the predicted target value, and update the matrix tensor through back propagation so that the loss function tends to 0, count the number of iterations, and if the number of iterations is less than the first preset number, return to step S6; if the number of iterations is not less than the first preset number, proceed to step S9; S9: sorting the predicted target values ​​of all optimization sequences generated in the iteration process, selecting the first N optimization sequences and re-analyzing them to obtain target values ​​respectively, and then sorting the first N optimization sequences according to the target values ​​obtained by the re-analysis to obtain the optimization sequence corresponding to the optimal target value; S10: If the number of cycles is less than the second preset number, the optimization sequence corresponding to the optimal target value replaces the original sequence, and returns to step S1; if the number of cycles is not less than the second preset number, the optimization sequence corresponding to the optimal target value in the current cycle is output.

2. The method for optimizing sequence design according to claim 1, characterized in that: In step S1, the encoding of the original sequence is in onehot form or kind of digital form.

3. The method for optimizing sequence design according to claim 1, characterized in that: Step S6 obtains the encoding of the input sequence according to the matrix tensor, including the following steps: S6.1: Multiply the matrix tensor by the synonymous substitution mask weight W(l*c), value W i It is the selection in the synonym table corresponding to the i-th smallest unit in the original sequence. Synonymous replacement is 1, non-synonymous replacement is 0, and then the matrix tensor is normalized: Where, T i is the i-th dimension of the data dimension in the third dimension c of the matrix tensor T; S6.2: According to the data dimension c in the second dimension of the transformed matrix tensor T, the index value of the maximum value is calculated in sequence and the sequence variant corresponding to the maximum value is encoded to obtain the encoding of the input sequence.

4. The method for optimizing sequence design according to claim 3, characterized in that: Step S6.2 also includes re-parameterizing the encoding of the input sequence.

5. The method for optimizing sequence design according to claim 1, characterized in that: The loss function in step S8 includes: Lm=-α / s Where Lm is the loss function, α is a constant, and s is the predicted target value.

6. The method for optimizing sequence design according to claim 1, characterized in that: After the first N optimized sequences are sorted according to the target values ​​obtained by the re-analysis in step S9, the following steps are further included: Take the top optimization sequences of the target value, and count the minimum unit selection Dis(l*c) of each position in each optimization sequence. The specific value Dis i is the probability of each minimum unit at the i-th position being selected.

7. The method for optimizing sequence design according to claim 6, characterized in that: In the loop of steps S1 to S9, except for the first loop, each of the remaining loops is based on the optimization sequence corresponding to the optimal target value obtained in the previous loop, all generated optimization sequences and the corresponding Dis(l*c) to sample and generate data for training the neural network model in the current loop.

8. The method for optimizing sequence design according to claim 7, characterized in that: The method of sampling and generating data for training the neural network model in the current cycle based on the optimization sequence corresponding to the optimal target value obtained in the previous cycle, all generated optimization sequences and corresponding Dis(l*c) includes: Sampling generation process: Perform random synonymous substitution on the optimized sequence corresponding to the optimal target value obtained in the previous cycle to generate sequence variants; Select the first m optimization sequences with the lowest target value among all optimization sequences generated in the previous cycle; According to Dis(l*c) obtained in the previous cycle, the probability of random synonymous replacement for each bit of the optimized sequence is obtained from Dis to perform distribution-based random replacement; The sequences generated in the sampling generation process are combined in a ratio of 2:4:4 to train the neural network model in the current cycle.

9. The method for optimizing sequence design according to claim 8, characterized in that: When the original sequence does not limit the gene in step S1, the minimum units at different positions in the original sequence are randomly replaced in step S2, and the target values ​​of all sequence variants generated by high-throughput experimental analysis are obtained in step S3.

10. The method for optimizing sequence design according to claim 8, characterized in that: When there are multiple target values, all sequence variants are analyzed in step S3 to obtain one of the target values ​​respectively. In step S8, the loss function is calculated according to the predicted target value, and the encoding of the sequence variant is input into a pre-trained RPF prediction neural network model for prediction, and a second loss function is obtained. At the same time, other target values ​​of the sequence variant are calculated and corresponding loss functions are obtained. All loss functions are added together according to different weights to obtain a total loss function, and the total loss function is used for back propagation to update the matrix tensor.

Images