An mRNA sequence optimization method based on the dynamic programming algorithm

Abstract

The present invention discloses an mRNA sequence optimization method based on a dynamic programming algorithm, which relates to the field of mRNA sequence optimization. The present invention respectively considers the positive (promoting) and negative (hindering) factors related to the process of mRNA being translated into protein sequences after entering host cells, and uses the dynamic programming algorithm to obtain its optimal comprehensive effect, greatly improving the final protein translation yield of mRNA, and having good application prospects in the field of mRNA vaccines.

Description

Technical Field

[0001] The present invention relates to the fields of algorithms and molecular biology. Specifically, the present invention relates to a method for optimizing mRNA sequences based on a dynamic programming algorithm. Background Art

[0002] Compared with traditional attenuated inactivated vaccines, mRNA technology has the advantages of short R & D cycle, simple production process, strong immunogenicity, high safety, etc., and can be widely used in many fields such as infectious disease vaccines, tumor immunity, protein replacement, etc. However, due to the defect that mRNA is more easily degraded than DNA, in order to enable it to translate more proteins, a more appropriate mRNA sequence design is necessary.

[0003] In the existing common optimization technologies in the market, the optimization features considered for sequences include codon adaptation index (CAI), GC content, minimum free energy (MFE) + secondary structure, rare codons, low information entropy repeat region sequences, ribosome loading, etc. These factors can be divided into two categories at the translation level: promoting or hindering.

[0004] However, the methods currently adopted still need to be improved. For example, the codon adaptation index (CAI) is calculated by comparing the similarity between the codon usage frequency in a gene and the codon usage preference of highly expressed genes in that species. From the calculation process of CAI, it can be seen that it only considers the performance of a single codon and does not consider the connectivity of codons in the sequence context. The GC content actually indicates the reduction of the core key base - uracil (U) that can be cleaved. However, too high GC content will cause the formation of low information entropy repeat region sequences in local areas. The low information entropy repeat region causes the polymerase's sliding window effect in the in vitro transcription (IVT) step, and thus the synthesized mRNA cannot be quality inspected. Therefore, the optimization algorithm should not only reduce uracil but also prevent the formation of low information entropy repeat sequences locally. Recently, it has been reported that the mRNA structure and the corresponding derived calculation index (MFE) have a very small blocking effect on ribosomes, and the main factor that blocks ribosomes is the abundance of tRNAs corresponding to codons. In the intracellular environment, the tertiary structure of mRNA is unstable and changeable, and there may be degradation-sensitive structures recognized by protein kinase R (PKR), etc., but the final step of degradation is the cleavage and hydrolysis of RNA by endonucleases. In addition, the optimization technology also needs to consider the operation time. One codon may have multiple synonymous codons that can replace it. Taking 100 codons as an example, each codon has an average of three synonymous codons. If an algorithm of exhaustive search of all synonymous codons is adopted, the calculation amount is 3 100 , and this number of calculations is unattainable.

[0005] Therefore, there is a need in the art to develop an mRNA optimization method to maximize mRNA anti-degradation, improve translation and protein production efficiency, and increase computational speed. Summary of the Invention

[0006] The object of the present invention is to provide an mRNA optimization method.

[0007] In a first aspect of the present invention, there is provided an mRNA sequence optimization method, wherein the mRNA is for expression in a target species or target cells, and the method comprises the steps of:

[0008] S1) Providing the codon occurrence frequency in the coding sequence (CDS) of the target species or target cells, and obtaining a set of codons for replacement according to the codon occurrence frequency;

[0009] Wherein, the set of codons does not include rare codons, and the rare codons are codons with an occurrence frequency < 10 times per thousand codons in the CDS of the target species or target cells;

[0010] S2) Providing the sensitivity index of a dinucleotide RNA endonuclease;

[0011] S3) Calculating the single-codon score of each codon in the set of codons based on the codon occurrence frequency obtained in step S1 and the sensitivity index of the dinucleotide RNA endonuclease obtained in step S2;

[0012] S4) Providing the mRNA coding sequence to be optimized or its encoded amino acid sequence, for each amino acid site in the sequence, selecting all corresponding synonymous codons from the set of codons in step S1, and calculating the maximum sequence score corresponding to each synonymous codon at each amino acid site using a state transition equation based on the single-codon score obtained in step S3 and the connection assignment between adjacent codons;

[0013] S5) Based on the maximum sequence score corresponding to each synonymous codon obtained in step S4, obtaining the maximum sequence score corresponding to the full-length mRNA coding sequence, thereby obtaining the optimal codon combination, which is the optimized mRNA coding sequence,

[0014] Wherein, the optimal codon combination is the codon combination that maximizes the score of the full-length mRNA coding sequence.

[0015] In another preferred example, in step S1, the CDS sequence of the target species is all genes in the target species that can encode proteins; and the CDS sequence of the target cells is the top 1000 genes that can encode proteins with the highest expression levels in the target cell type.

[0016] In another preferred example, in step S1, the codon occurrence frequency is the number of times a certain codon appears per thousand codons in the CDS of the target species or target cells.

[0017] In another preferred example, in step S1, the codon set does not include the start codon AUG and the stop codon UGA.

[0018] In another preferred example, in step S2, the RNA endonuclease is selected from: RNAase A, RNase T2, RNAaseL, or a combination thereof.

[0019] In another preferred example, in step S2, the sensitivity index of the dinucleotide RNA endonuclease represents the sensitivity of the phosphodiester bond between every two bases to the RNA endonuclease. The higher the sensitivity, the lower the sensitivity index.

[0020] In another preferred example, the sensitivity index range of the dinucleotide RNA endonuclease is [-10, 10], and the sensitivity indices of each dinucleotide combination are: AA is 3, AU is -8, AC is 6, AG is 6, UA is -10, UU is -10, UC is 0, UG is 0, CA is 5, CU is 6, CC is 10, CG is 10, GA is 5, GU is -10, GC is 10, GG is 10.

[0021] In another preferred example, in step S3, the sensitivity index of the dinucleotide RNA endonuclease and the codon occurrence frequency are added through a weight parameter w, and the weight parameter w represents the contribution degree of the two components to the final translation yield.

[0022] In another preferred example, in step S3, the calculation formula for the single codon score is as follows:

[0023] S(C j ) = D(C j ) + w·f(C j )

[0024] In the formula,

[0025] S(C j ) is the single codon score;

[0026] f(C j ) is the codon occurrence frequency of codon C j ;

[0027] D(C j ) is the internal score of codon C j ;

[0028] w is the weight coefficient;

[0029] Among them, the calculation formula for the internal score is as follows:

[0030] D(C j ) = pair_values(C j [1]C j [2]) + pair_values(C j [2]C j [3])

[0031] Wherein,

[0032] C j [1], C j [2], C j [3] respectively represent the 1st, 2nd, and 3rd bases of the codon Cj;

[0033] pair_values(C j [1]C j [2]) represents the connection assignment between the 1st and 2nd bases of the codon Cj;

[0034] pair_values(C j [2]C j [3]) represents the connection assignment between the 2nd and 3rd bases of the codon Cj.

[0035] In another preferred example, the connection assignment is positively correlated with the sensitivity index of the dinucleotide RNA endonuclease.

[0036] In another preferred example, the range of the weight coefficient w is 0.1 - 1. Preferably, w is selected from 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, or 1.

[0037] In another preferred example, if C j [3] is G or C, then multiply the internal score D(C j ) by 1.5 times.

[0038] In another preferred example, in step S4, the initial score of the set sequence is 0.

[0039] In another preferred example, in step S4, the calculation formula for the connection assignment V(C j , C k ) between adjacent codons is as follows:

[0040] V(C j , C k ) = pair_values(C j [3]C k [1])

[0041] Wherein,

[0042] Cj 、C k denote the jth and kth codons of the mRNA sequence, respectively, where k = j + 1;

[0043] pair_values(C j [3]C k [1]) indicates codon C j The third base and codon C k The first base-to-base connection assignment.

[0044] In another preferred embodiment, in step S4, for the j-th synonymous codon of the i-th amino acid, the possible synonymous codons of the i-1-th amino acid are used as variables to calculate the maximum sequence score DP[i][j] that can be obtained.

[0045] In another preferred embodiment, in step S4, the maximum sequence score DP[i][j] corresponding to the j-th synonymous codon of the i-th amino acid is calculated as follows:

[0046]

[0047] Where,

[0048] DP[i-1][m] represents the maximum sequence score that the first i-1 amino acids can obtain when the i-1 amino acid uses its m-th synonymous codon;

[0049] k i-1 Indicates the number of synonymous codons corresponding to the i-1th amino acid;

[0050] C (i-1)m Indicates the mth synonymous codon corresponding to the i-1th amino acid;

[0051] C ij Indicates the jth synonymous codon corresponding to the i-th amino acid;

[0052] V(C (i-1)m ,C ij ) indicates codon C (i-1)m and C ij The adjacent codons between them are assigned values;

[0053] f(C ij ) is codon C ij The frequency of codon occurrence;

[0054] D(C ij ) is codon C ij Internal rating of

[0055] w is the weight coefficient.

[0056] In another preferred example, in step S5, the maximum score corresponding to the full-length mRNA coding sequence is selected from the sequence scores obtained in step S4, and each codon corresponding to this score is traced back to obtain the optimal codon combination.

[0057] In another preferred example, the maximum score Τ corresponding to the full-length mRNA coding sequence max is obtained by the following formula:

[0058]

[0059] where n is the number of amino acids in the full-length sequence.

[0060] In another preferred example, in step S5, the optimal codon combination is obtained by the following formula:

[0061] {C1, C2,..., C n} = Backtrack(j * );

[0062] where j* is the corresponding codon that maximizes the sequence score DP[n][j] among the full-length n amino acids, and it is obtained by the following formula:

[0063] j* = arg max j DP[n][j].

[0064] In another preferred example, the method further includes the step of:

[0065] S6) From the mRNA coding sequence obtained in step S5, eliminate the negative factor sequences through synonymous codon substitution, and the negative factor sequences include endonuclease cleavage sites and / or low information entropy repeat region sequences.

[0066] In another preferred example, the endonuclease cleavage site is the sapI cleavage site.

[0067] In another preferred example, the low information entropy repeat region sequence refers to a sequence fragment with a repeat region exceeding two codon lengths, and the repeat region is a fragment formed by multiple two-base or one-base repeats connected together.

[0068] In another preferred example, the elimination is performed by an exhaustive method, and it is checked whether new negative factor sequences are introduced for all possible synonymous codon substitutions.

[0069] In another preferred example, select the mRNA coding sequence with the maximum sequence score from the mRNA coding sequences from which the negative factor sequences have been eliminated to obtain the optimized mRNA coding sequence; the sequence score is calculated by the method described in step S4.

[0070] In another preferred example, step S4 further includes:

[0071] Preprocessing the mRNA coding sequence to be optimized or the amino acid sequence encoded thereby.

[0072] In another preferred example, sequence preprocessing of the amino acid sequence includes:

[0073] 1) Determining whether all letters in the sequence belong to the abbreviations of amino acids;

[0074] 2) Determining whether the first letter of the sequence is methionine - M, and if not, placing M at the first position.

[0075] In another preferred example, sequence preprocessing of the CDS sequence inside the mRNA includes one or more steps selected from the following group:

[0076] 1) Replacing thymine (T) in the sequence with uracil (U);

[0077] 2) Determining whether the sequence length is a multiple of 3;

[0078] 3) Determining whether the base composition units are A, U, C, and G;

[0079] 4) Replacing the start codon of the sequence with AUG and the stop codon with UGA.

[0080] In another preferred example, the preprocessing further includes determining the length of the CDS sequence. If the CDS sequence length ≥ 2000 nt, select the mRNA of a shorter protein sequence encoded by the same gene as the input sequence.

[0081] In another preferred example, the method further includes: comparing the amino acid sequence translated from the CDS sequence before preprocessing with the amino acid sequence translated from the optimized mRNA coding sequence to verify the correctness of the optimized sequence.

[0082] In the second aspect of the present invention, there is provided a device for mRNA sequence optimization, the device includes an input module, a processing module and an output module, wherein:

[0083] [[ID=3⑨The input module is used to input: the target species and / or target cells, and the mRNA sequence to be optimized or the amino acid sequence encoded thereby;

[0084] The processing module optimizes the mRNA sequence according to the method as described in the first aspect of the present invention;

[0085] The output module is used to output the optimized mRNA sequence.

[0086] In the third aspect of the present invention, a computer device is provided, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the method described in the first aspect of the present invention is implemented.

[0087] In the fourth aspect of the present invention, a computer-readable storage medium is provided, on which a computer program for implementing the method described in the first aspect of the present invention is stored.

[0088] In the fifth aspect of the present invention, an mRNA vaccine is provided, and the coding sequence of the mRNA is optimized by the method described in the first aspect of the present invention.

[0089] It should be understood that within the scope of the present invention, the above technical features of the present invention and the technical features specifically described below (such as in the embodiments) can be combined with each other to form new or preferred technical solutions. Due to space limitations, they will not be elaborated here one by one. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] The following drawings are used to illustrate the specific implementation of the present invention, rather than to limit the scope of the present invention defined by the claims.

[0091] Figure 1 Shows a schematic diagram of the system of the present invention.

[0092] Figure 2 Shows a flowchart of the optimization method of the present invention.

[0093] Figure 3 Shows a schematic diagram of the mRNA structure. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0094] Through extensive and in-depth research, the inventor of the present invention has developed for the first time an mRNA sequence optimization method based on the dynamic programming algorithm. The method of the present invention finds the optimal balance point between mRNA anti-degradation and translation speed. After optimization, the mRNA has no repetitive low-information-entropy sequence fragments, and there are no sapI restriction sites within the sequence. Compared with the existing codon adaptation index, the present invention comprehensively considers the connectivity within a single codon and the connectivity between adjacent codons, thereby providing a codon optimization method for the full-length sequence. At the same time, the present invention uses the dynamic programming algorithm to find the optimal codon combination, saving a large amount of computing power compared with the traditional exhaustive method and improving the operation speed. On this basis, the present invention is completed.

[0095] Influencing Factors of Sequence Optimization

[0096] 1. Codon Frequency

[0097] In an organism, DNA is transcribed into an mRNA sequence. The coding sequence (CDS) of the mRNA sequence is composed of codons connected in sequence. Each codon consists of three bases, and each codon corresponds to a unique amino acid. One amino acid may have multiple corresponding codons, and multiple codons corresponding to the same amino acid are called synonymous codons. The relative fitness (W) of a single codon or the frequency of occurrence per thousand codons in the mRNA coding sequence (CDS sequence) characterizes the usage frequency of each codon, corresponding to the abundance of transfer RNA (tRNA) in the cell. The usage frequency corresponding to each codon is species-specific. During the synthesis of a protein peptide chain, the codons on the mRNA bind to tRNA molecules under the action of ribosomes, and the amino acids corresponding to the codons are successively connected to synthesize the protein peptide chain. Different transfer RNAs have different abundances in the cell, thus affecting the efficiency of amino acids binding to the polypeptide.

[0098] 2. The third base of the codon

[0099] When the third base of the codon is G or C (GC3), the number of hydrogen bonds formed between the codon and the anticodon base on tRNA is 3, which is more stable than the 2 hydrogen bonds between A and T, helping to improve the specific binding of the codon and the anticodon, thus ensuring the accuracy of protein synthesis. In addition, both GC3 and the overall GC content of the sequence can predict a longer mRNA half-life.

[0100] 3. Ribonuclease (RNA)

[0101] The structural composition of RNA is as Figure 3As shown, the CDS region is located in the middle of the whole sequence, rather than at both ends. The related degradation mechanism mainly occurs through the opening of phosphodiester bonds by ribonucleases (RNAs). It is investigated that the related RNAses include RNase A, RNase T1, RNase T2 and RNase L. RNase T1 exists in fungi and bacteria, RNase A exists in vertebrates, RNase T2 exists in all organisms except archaea, and RNase L (ribonuclease L) is highly conserved in vertebrates, including mammals, birds, reptiles, amphibians and fish. Therefore, for the use of mRNA vaccines or related products, attention needs to be paid to RNase A, RNase T2 and RNase L. The sequence characteristics of the cleavage of these three endonucleases are as follows: RNase A mainly cleaves the phosphodiester bond after pyrimidines (C - cytosine, U - uracil); RNase T2 mainly cleaves the phosphodiester bond between GU and AU; RNase L mainly cleaves UU and UA. These three RNA endonucleases can also cleave the linkages of other base residues, but with low efficiency.

[0102] Dynamic programming algorithm

[0103] The present invention optimizes the mRNA sequence based on the dynamic programming algorithm, saving computing power and improving the operation speed compared with the traditional exhaustive method. The dynamic programming algorithm of the present invention is implemented using the state transition equation, specifically including:

[0104] 1) State definition: Each state represents a partial sequence constructed to a certain specific position (i.e., a certain specific amino acid or codon) and its cumulative score;

[0105] 2) State transition: According to the selection of synonymous codons and their corresponding scores, update the cumulative score and record the optimal path;

[0106] Among them, the initial state is the starting position of the sequence, and the cumulative score is zero.

[0107] Optimal sequence screening: After completing the dynamic programming process, sort all possible generated sequences according to the cumulative score, and select the one with the highest score as the optimal solution.

[0108] Sequence optimization method

[0109] The present invention provides a method for optimizing mRNA sequences based on the dynamic programming algorithm, and finds the optimal balance point between anti - degradation and translation speed through the dynamic programming algorithm. The method of the present invention includes the following steps:

[0110] S1: Calculate the occurrence frequency of codons. For mRNA vaccines known to be translated and expressed in certain specific cells, collect the RNAseq data of the target cells, calculate the expression profile, and obtain the top 1000 genes that can be translated into proteins with the highest expression levels. If there is no RNAseq data related to the target cells or the specific expressing cells are unknown, collect all the CDS sequences of the target species. For each codon, calculate the number of occurrences per thousand codons in all the CDS sequences. If the number of occurrences of a certain codon is less than 10, it is a rare codon, and this rare codon is not used as a synonymous substitution candidate. In particular, the start codon uses AUG, and the stop codon uses UGA.

[0111] S2: Construct the sensitivity index of dinucleotide ribonuclease (RNA) endonuclease. In the present invention, the anti-degradability of a sequence is defined by the sensitivity index of dinucleotide ribonuclease (RNA) endonuclease. RNA endonucleases have different cleavage characteristics. For example, RNase A efficiently cleaves at pyrimidines (C, U), and RNase T2 efficiently cleaves at sensitive junctions such as GU and AU, while RNase L tends to efficiently cleave UU and UA sequences. According to factors such as the upstream and downstream relationships, cleavage efficiency, and molecular similarity of these enzymes in the intracellular reaction pathway, the present invention defines the sensitivity index of 16 dinucleotide combinations. In one embodiment, the self-defined value range is [-10, 10], AA is 3, AU is -8, AC is 6, AG is 6, UA is -10, UU is -10, UC is 0, UG is 0, CA is 5, CU is 6, CC is 10, CG is 10, GA is 5, GU is -10, GC is 10, GG is 10.

[0112] S3: Sequence basic preprocessing. If the initial sequence to be processed is an amino acid sequence, verify whether all letters belong to the abbreviations of amino acids (["A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V"]). If the first letter is not methionine - M, add M at the first position. If it is the translation sequence (coding sequence, CDS) within messenger ribonucleic acid (mRNA), replace thymine (T) in the sequence with uracil (U), verify whether its length is a multiple of 3 and whether the base composition units are ["A","U","C","G"]. Replace the start codon with AUG and the stop codon with UGA. If the translation sequence (CDS) of the mRNA to be optimized is too long (>=2 knt), even after optimization design, the accumulation of inevitable negative factors within the sequence may lead to its degradation in the cell and loss of design value. Therefore, in the presence of a shorter protein sequence, it is necessary to select the mRNA of a shorter protein sequence encoded by the same gene from uniprot for design.

[0113] S4: Codon scoring. To translate and produce more proteins, the mRNA sequence requires two aspects of properties: 1) resistance to degradation, 2) high - speed translation by ribosomes; the property of resistance to degradation has been defined in step S2, and the main property of high - speed translation by ribosomes - the occurrence frequency of codons - has been calculated in S1. The whole sequence is composed of single codons connected in sequence. Therefore, it is very necessary to calculate the sum of the two properties of steps S1 + S2 within each codon for use as the state - transfer point in the subsequent dynamic programming algorithm. In addition, when the third position of the codon is G or C, the binding energy between the codon and the anticodon is stronger, which is conducive to improving the translation efficiency. Therefore, if the third position of the codon is G or C, the corresponding connection assignment is multiplied by 1.5 times. To sum up, the internal score D(C j ) is calculated as follows:

[0114] D(C j ) = pair_values(C j [1]C j [2]) + pair_values(C j [2]C j [3])

[0115] In the formula,

[0116] C j [1], C j [2], C j [3] respectively represent the 1st, 2nd, and 3rd bases of the codon Cj;

[0117] pair_values(C j [1]C j [2]) represents the connection assignment between the first and second bases of codon Cj;

[0118] pair_values(C j [2]C j [3]) represents the connection assignment between the second and third bases of codon Cj.

[0119] The connection assignment represents the connection strength between bases, which can be calculated through the sensitivity index of the above-mentioned dinucleotide endonuclease.

[0120] If C j [3] ∈ {G, C}, then <>

[0121] D(C j ) = 1.5 × D(C j )

[0122] Final scoring formula for codons:

[0123] S(C j ) = D(C j ) + w · f(C j )

[0124] Symbol definition:

[0125] S(C j ) is the score for a single codon;

[0126] f(C j ) is the occurrence frequency of codon C j ;

[0127] D(C j ) is the internal score of codon C j ;

[0128] w is the weight coefficient.

[0129] Among them, the weight coefficient w is used to balance the importance of anti-degradation ability and codon usage frequency. If the impact of translation efficiency on the final yield is greater, the weight coefficient is larger; if the degradation rate has a greater impact, the weight coefficient is smaller. The range of the weight coefficient w is 0.1 - 1. Preferably, w is selected from 0.1, 0.2, 0. = 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 or 1.

[0130] S5: Dynamic programming algorithm is used to calculate the optimal sequence. The entire sequence consists of a single codon connected in sequence. Except for the start and stop codons, each codon is connected to two codons upstream and downstream, so there are two connection assignments. Under the premise that the start codon (AUG) and the stop codon (UGA) are fixed, in order to avoid repeated calculations, only the upstream connection is calculated for each codon. Therefore, in addition to the internal score D (C j ) In addition, the connection between codons is assigned V(C i ,C i+1 ), is also an important factor affecting the overall score. If an exhaustive method is used to calculate each synonymous codon and then select the optimal path, the amount of calculation will be very large. Therefore, the present invention uses a dynamic programming algorithm to optimize the sequence. Dynamic programming takes both into account through the state transition equation, ensuring that each step of the selection not only optimizes the score of the current codon, but also takes into account the connection assignment with the previous codon, thereby maximizing the overall score. Specifically, it includes the following steps:

[0131] 1) Calculate the connection assignment of adjacent codons:

[0132] V(C j , C k )=pair_values(C j [3]C k [1])

[0133] Where,

[0134] C j 、C k denote the jth and kth codons of the mRNA sequence, respectively, where k = j + 1;

[0135] pair_values(C j [3]C k [1]) indicates codon C j The third base and codon C k The first base-to-base connection assignment.

[0136] 2) Dynamic programming state transfer:

[0137]

[0138] Where DP[i][j] represents the maximum sequence score that can be obtained when the i-th amino acid uses its j-th synonymous codon in the first i amino acid sequence;

[0139] DP[i-1][m] represents the maximum sequence score corresponding to the first i-1 amino acids when the i-1 amino acid uses its m-th synonymous codon;

[0140] K i-1 represents the number of synonymous codons corresponding to the (i - 1)-th amino acid;

[0141] C (i-1)m represents the m-th synonymous codon corresponding to the (i - 1)-th amino acid;

[0142] C ij represents the j-th synonymous codon corresponding to the i-th amino acid;

[0143] V(C (i-1)m ,C ij ) represents the adjacent codon connection assignment between codons C (i-1)m and C ij ;

[0144] f(C ij ) is the occurrence frequency of codon C ij ;

[0145] D(C ij ) is the internal score of codon C ij ;

[0146] 3) Calculate the optimal score:

[0147]

[0148] where n is the number of amino acids in the full-length sequence.

[0149] 4) Backtrack the optimal codon sequence

[0150] {C1, C2,..., C n} = Backtrack(j * )

[0151] where j* = arg max j DP[n][j]

[0152] By defining a clear scoring system and dynamic programming model, this method can systematically select the optimal codon sequence on the basis of considering the anti-endonucleolytic degradation ability of codons, codon usage frequency, and the connection effect between codons.

[0153] S6: Eliminate the negative factor fragment sequences. The algorithm of the present invention also includes eliminating small fragment sequences that affect in vitro transcription (IVT) production, including:

[0154] 1) Restriction sites. The algorithm of the present invention needs to eliminate restriction sites that may be affected by nucleases in the cells used for in vitro transcription (IVT). In one embodiment, the restriction site includes the sapI restriction site, and its sequence characteristics are as follows:

[0155] 5’...GCTC TTC(N)1▼...3’

[0156] 3’...CGAGAAG(N)4▲...5’

[0157] When this characteristic sequence fragment exists, it will cause the template DNA fragment to be cut and eventually degraded. When using synonymous codon replacement to eliminate the restriction site fragment, since the number of codons involved is small, all possible methods can be listed by the exhaustive method, and the overall sequence of the finally selected synonymous replacement scheme is scored highest according to the aforementioned scoring method.

[0158] 2) Low information entropy repeat region sequences. If multiple repeated bases appear in the sequence, both mRNA synthesis and subsequent product sequencing verification require polymerase. And the polymerase has poor fidelity in the repeated base fragment. Since the coverage range of the ribosome is 2 codons, if the repeat region exceeds two codon lengths, it needs to be removed. The "repeat region" refers to a fragment formed by the connection of multiple di-base or mono-base repeats. Similarly, the synonymous codons with the highest scores are selected for replacement by the exhaustive method.

[0159] The method of the present invention may further include: converting the input CDS sequence into an amino acid sequence to be compared with the amino acid sequence corresponding to the finally optimized sequence to verify the correctness of the finally optimized sequence.

[0160] The main advantages of the present invention include:

[0161] 1) The present invention selects the promoting factors and hindering factors that have obvious effects on the process of mRNA translation into protein. Finally, the codon occurrence frequency, the third base of the codon, the connectivity of adjacent upstream and downstream codons, and the elimination of restriction sites and low information entropy repeat region sequences are selected as the main sequence optimization points. The optimized mRNA sequence of the present invention reaches an optimal balance state between the U base content and the use of high-frequency codons, has no repeated low information entropy sequence fragments, and has no sapI restriction sites in the sequence.

[0162] 2) Compared with the existing codon adaptation index, the present invention comprehensively considers the connectivity within a single codon and the connectivity between adjacent codons, thereby providing a codon optimization method for the full-length sequence.

[0163] 3) The present invention uses the dynamic programming algorithm to find the optimal codon combination, saving a large amount of computing power compared with the traditional exhaustive method and improving the operation speed.

[0164] The present invention will be further described below in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and not to limit the scope of the present invention. The experimental methods without specific conditions noted in the following embodiments are generally carried out under conventional conditions, such as those described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or according to the conditions recommended by the manufacturer. Unless otherwise specified, percentages and parts are weight percentages and weight parts.

[0165] Example 1 Optimization of the mRNA sequence of classical swine fever protein I215L

[0166] The classical swine fever protein I215L was transfected into HEK293T for translational expression. HEK293T is a human renal cell, so the host is human.

[0167] The I215L sequence is as follows:

[0168] MVSRFLIAEYRHLIENPSENFKISVNENNITEWDVILRGPPDTLYEGGLFKAKVAFPPEYPYAPPKLTFTSEMWHPNIYPDGRLCISILHGDNAEEQGMTWSPAQKIDTILLSVISLLNEPNPDSPANVDAAKSYRKYVYKEDLESYPMEVKKTV

[0169] KKSLDECSPEDIEYFKNAASNVPPIPSDAYEDECEEMEDDTYILTYDDDEEEE

[0170] DEEMDDE(SEQ ID NO:1)

[0171] The design is divided into the following steps:

[0172] S1: Calculate the occurrence frequency of codons. Collect all the CDS sequences of humans. For each codon, calculate the number of occurrences per thousand codons in all the CDS sequences, as shown in Table 1 below. Rare codons in the table are not candidates for synonymous substitution during sequence optimization. The table does not include start and stop codons. The start codon uses AUG, and the stop codon uses UGA.

[0173] Table 1 Codon occurrence frequency

[0174] Codon Frequency Codon Frequency Codon Frequency Codon Frequency GCC 26 CGG 11 GAC 24 CAG 36 GCU 19 CGC 9 GAU 24 CAA 14 GCA 17 CGA 6 UGC 11 GAG 40 GCG 6 CGU 5 UGU 10 GAA 34 AGA 13 AAU 18 CAC 15 GGC 20 AGG 12 AAC 18 CAU 12 GGA 17 AUC 18 GGG 15 CUG 36 AGC 20 AUU 16 GGU 11 CUC 18 UCC 17 AUA 8 UUC 17 CUU 14 UCU 17 CUG 36 UUU 17 UUG 13 UCA 14 CUC 18 CCU 20 CUA 7 AGU 14 CUU 14 CCC 19 AAG 32 UCG 4 UUG 13 CCA 19 AAA 27 GUG 26 UUA 9 CCG 6 ACC 18 GUC 13 AAG 32 UGG 11 ACA 17 GUU 12 AAA 27 UAC 13 ACU 14 UAU 12 UGG 11 UAU 12 ACG 6 UAC 13 GUC 13 GUG 26 GUA 8 GUU 12

[0175] S2: Construct the sensitivity index of dinucleotide ribonucleic acid (RNA) endonuclease. This index is mainly custom-built based on the known RNA endonucleases in eukaryotes, including the RNAase A and RNase T2 families, as well as the cleavage characteristics of RNAase L. RNAase A efficiently cleaves at pyrimidines (C, U), and RNase T2 efficiently cleaves at sensitive junctions such as GU and AU. RNAase L tends to efficiently cleave UU and UA sequences. Other linkage types are also cleaved by these endonucleases, but with lower efficiency than the above special sites.

[0176] There are 16 dinucleotide combinations. According to factors such as the upstream and downstream relationships, cleavage efficiency, and molecular similarity of these enzymes in the intracellular reaction pathway, the custom-assigned value range here is [-10, 10]. For example, AA is 3, AU is -8, AC is 6, AG is 6, UA is -10, UU is -10, UC is 0, UG is 0, CA is 5, CU is 6, CC is 10, CG is 10, GA is 5, GU is -10, GC is 10, and GG is 10.

[0177] Table 2 Sensitivity index of dinucleotide RNA endonuclease

[0178] Base pair Assignment Base pair Assignment Base pair Assignment Base pair Assignment AA 3 CA 5 GA 5 TA -10 AC 6 CC 10 GC 10 TC 0 AG 6 CG 10 GG 10 TG 0 AT -8 CT 6 GT -10 TT -10

[0179] S3: Sequence-based preprocessing. If the initial sequence to be processed is an amino acid sequence, verify whether all letters belong to the abbreviations of amino acids (["A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V"]). If the first letter is not methionine - M, place M at the first position. If it is the translation sequence (coding sequence, CDS) within messenger ribonucleic acid (mRNA), replace thymine (T) in the sequence with uracil (U), verify whether its length is a multiple of 3 and whether the base composition units are ["A", "U", "C", "G"]. Replace the start codon with AUG and the stop codon with UGA. To verify the correctness of the finally optimized sequence, convert the input CDS sequence into an amino acid sequence to be compared with the amino acid sequence corresponding to the finally optimized sequence. If the translation sequence (CDS) of the mRNA to be optimized is too long (>= 2 knt), even after optimization design, the accumulation of inevitable negative factors within the sequence may lead to its degradation in the cell and loss of design value. Therefore, if there is a shorter protein sequence available, it is necessary to select the mRNA of the shorter protein sequence encoded by the same gene from uniprot for design.

[0180] The length of I215L is only 215 AAs, and each amino acid letter corresponds correctly. The length does not exceed 2 knt, meeting the requirements.

[0181] S4: Codon scoring. To translate and produce more proteins, the mRNA sequence requires two properties: 1) resistance to degradation, 2) high-speed translation by ribosomes; the property of resistance to degradation has been defined in step S2, and the main property of high-speed translation by ribosomes - the occurrence frequency of codons - was calculated in S1. The entire sequence is composed of individual codons connected in sequence. Therefore, it is necessary to calculate the sum of the two properties of steps S1+S2 within each codon for use as the state transition points in the subsequent dynamic programming algorithm. Additionally, when the third position of the codon is G or C, the binding energy between the codon and the anticodon is stronger, which is conducive to improving translation efficiency. Therefore, if the third position of the codon is G or C, the corresponding connection assignment is multiplied by 1.5 times. [[ID=??]]

[0182] Table 3 Internal Scoring of Codons

[0183] Codon Score Codon Score Codon Score Codon Score AAG 16.7 AAA 8.7 AAU -3.2 AAC 15.3 ACC 25.8 ACA 12.7 ACU 13.4 ACG 24.6 AGA 12.3 AGG 25.2 CGG 31.1 UCC 16.7 UCU 7.7 UCA 6.4 AGC 26 AUG -9.9 AUC -10.2 AUU -16.4 AGU -2.6 CAA 9.4 CAC 18 CAU -1.8 CAG 20.1 UUG -13.7 CCU 18 CCC 31.9 CCA 16.9 GAU -0.6 CUG 12.6 CUC 10.8 CUU -2.6 GCU 17.9 GAG 20.5 GAA 11.4 GAC 18.9 GGA 16.7 GCA 16.7 GGU 1.1 GCC 32.6 GUC -13.7 GGG 31.5 UGU -9 GGC 32 UAU -16.8 GUU -18.8 UAC -4.7 GUG -12.4 UUC -13.3 UGC 16.1 UGG 16.1 UUU -18.3

[0184] Among them, rare codons have been removed.

[0185] S5: Calculate the optimal sequence using the dynamic programming algorithm. The entire sequence is composed of individual codons connected in sequence. Except for the start and stop codons, each codon connects to the two codons upstream and downstream, that is, there are two connection assignments. On the premise that the start codon (AUG) and the stop codon (UGA) are already fixed, to avoid repeated calculations, only the upstream connection is calculated for each codon. Therefore, in addition to the internal score D(Cj) of each codon, the connection assignment V(Ci,Ci+1) between codons is also an important factor affecting the overall score. Use the state transition equation to maximize the overall score.

[0186] The result sequence obtained according to the above algorithm is as follows:

[0187] [[ID=2|1]]AUGGUGAGCCGGUUCCUGAUCGCCGAGUACCGGCACCUGAUCGAGAACCCCAGCGAGAACUUCAAGAUCAGCGUGAACGAGAACAACAUCACCGAGUGGGACGUGAUCCUGCGGGG CCCCCCC GACACCCUCUACGAGGGCGGCCUCUUCAAGGCCAAGGUGGCCUU CCCCCCC GAGUACCCCUACG CCCCCCCCAAGCUGACCUUCACCAGCGAGAUGUGGCACCCCAACAUCUACCCCGACGGCCGGCUCUGCAUCAGCAUCCUGCACGGCGACAACGCCGAGGAGCAGGGCAUGACCUGGAGCCCCGCCCAGAAGAUCGACACCAUCCUGCUGAGCGUGAUCAGCCUGCUGAACGAGCCCAACCCCGACAGCCCCGCCAACGUGGACGCCGCCAAGAGCUACCGGAAGUACGUCUACAAGGAGGACCUGGAGAGCUACCCCAUGGAGGUGAAGAAGACCGUGAAGAAGAGCCUGGACGAGUGCAGCCCCGAGGACAUCGAGUACUUCAAGAACGCCGCCAGCAACGUGCCCCCCAUCCCCAGCGACGCCUACGAGGACGAGUGCGAGGAGAUGGAGGACGACACCUACAUCCUGACCUACGACGACGACGAGGAGGAGGAGGACGAGGAGAUGGACGACGAGUGA(SEQ ID NO:2)

[0188] S6: Eliminate the negative factor fragment sequence. Mainly, small fragment sequences that affect in vitro transcription (IVT) production need to be eliminated. 1) sapI restriction site. 2) Low information entropy repeat region sequence. The low information entropy repeat region sequence fragment is underlined in SEQ ID NO:2, and the corresponding optimized sequence fragment is underlined in SEQ ID NO:3.

[0189] AUGGUGAGCCGGUUCCUGAUCGCCGAGUACCGGCACCUGAUCGAGAACCCCAGCGAGAACUUCAAGAUCAGCGUGAACGAGAACAACAUCACCGAGUGGGACGUGAUCCUGCGGGG CCCACCC GACACCCUCUACGAGGGCGGCCUCUUCAAGGCCAAGGUGGCCUU CCCACCC GAGUACCCCUACG CCCCACCCAAGCUGACCUUCACCAGCGAGAUGUGGCACCCCAACAUCUACCCCGACGGCCGGCUCUGCAUCAGCAUCCUGCACGGCGACAACGCCGAGGAGCAGGGCAUGACCUGGAGCCCCGCCCAGAAGAUCGACACCAUCCUGCUGAGCGUGAUCAGCCUGCUGAACGAGCCCAACCCCGACAGCCCCGCCAACGUGGACGCCGCCAAGAGCUACCGGAAGUACGUCUACAAGGAGGACCUGGAGAGCUACCCCAUGGAGGUGAAGAAGACCGUGAAGAAGAGCCUGGACGAGUGCAGCCCCGAGGACAUCGAGUACUUCAAGAACGCCGCCAGCAACGUGCCCCCAAUCCCCAGCGACGCCUACGAGGACGAGUGCGAGGAGAUGGAGGACGACACCUACAUCCUGACCUACGACGACGACGAGGAGGAGGAGGACGAGGAGAUGGACGACGAGUGA(SEQ ID NO:3)

[0190] Example 2 Optimization of the mRNA sequence of classical swine fever protein pB602L

[0191] The classical swine fever protein pB602L was transfected into HEK293T for translational expression. HEK293T is a human kidney-derived cell line. The host is also human.

[0192] The sequence of pB602L is as follows:

[0193] MAEFNIDELLKNVLEDPSTEISEETLKQLYQRTNPYKQFKNDSRVAFCSFTNLREQYIRRLIMTSFIGYVFKALQEWMPSYSKPTHTTKTLLSELITLVDTLKQETNDVPSESVVNTILSIADSCKTQTQKSKEAKTTIDSFLREHFVFDPNLHAQSAYTCASTCADTNVDTCASTCADTNVDTCASTCADTNVDTCASTCADTNVNTCASMCADTNVDTCASTCANTCASTEYTDLADPERIPLHIMQKTLNVPNELQADIDAITQTPQGYRAAAHILQNIELHQSIKHMLENPRAFKPILFNTKITRYLSQHIPPQDTFYKWNYYIEDNYEELRAATESIYPEKPDLEFAFIIYDVVDSSNQQKVDEFYYKYKDQIFSEVSSIQLGNWTLLGSFKANRERYNYFNQNNEIIKRILDRHEEDLKIGKEILRNTIYHKKAKNIQETGPDAPGLSIYNSTFHTDSGIKGLLSFKELKNLEKASGNIKKAREYDFIDDCEEKIKQLLSKENLTPDEESELIKTKKQLNNALEMLNVPDDTIRVDMWVNNNNKLEKEILYTKAEL(SEQ ID NO:4)

[0194] In the sequence optimization step, S1, S2, and S4 are the same as the previous steps.

[0195] The running result of S3: The length of pB602L is 602 amino acids (AA), and each amino acid letter is correct. The corresponding nucleotide length is 1806 nt, which does not exceed 2 knt. The sequence meets the requirements.

[0196] The sequence obtained through dynamic optimization in S5 is as follows:

[0197]

[0198] S6: Eliminate the negative factor fragment sequences. Mainly, the small fragment sequences that affect in vitro transcription (IVT) production need to be eliminated. 1) The sapI restriction site. 2) The low information entropy repeat region sequences. After optimization, these two negative factor fragment sequences do not exist.

[0199] All documents mentioned in this invention are cited herein for reference as if each individual document was cited for reference separately. In addition, it should be understood that after reading the above teachings of this invention, those skilled in the art can make various changes or modifications to this invention, and these equivalent forms also fall within the scope defined by the appended claims of this application.

Claims

1. A method for optimizing an mRNA sequence, wherein the mRNA is used for expression in a target species or target cell, characterized in that: The method comprises the steps of: S1) providing the codon occurrence frequency in the coding sequence (CDS) of the target species or target cell, and obtaining a codon set for replacement based on the codon occurrence frequency; Wherein, the codon set does not include rare codons, and the rare codons are codons that appear less than 10 times per thousand codons in the CDS of the target species or target cell; S2) provides a sensitivity index for a dibase pair RNA endonuclease; S3) calculating a single codon score for each codon in the codon set based on the codon occurrence frequency obtained in step S1 and the sensitivity index of the two-base pair RNA endonuclease obtained in step S2; S4) providing an mRNA coding sequence to be optimized or an amino acid sequence encoded by it, selecting all corresponding synonymous codons from the codon set described in step S1 for each amino acid site in the sequence, and calculating the maximum sequence score corresponding to each synonymous codon at each amino acid site using a state transition equation based on the single codon score obtained in step S3 and the connection assignment between adjacent codons; S5) Based on the maximum sequence score corresponding to each synonymous codon obtained in step S4, the maximum sequence score corresponding to the full-length mRNA coding sequence is obtained, thereby obtaining the optimal codon combination, which is the optimized mRNA coding sequence. The optimal codon combination is the codon combination that maximizes the score of the full-length mRNA coding sequence.

2. The method according to claim 1, wherein In step S2, the RNA endonuclease is selected from: RNAase A, RNase T2, RNAase L, or a combination thereof.

3. The method according to claim 1, wherein In step S3, the calculation formula for a single codon score is as follows: S(C j )=D(C j )+w·f(C j ) Where, S(C j ) scores for individual codons; f(C j ) is codon C j The frequency of codon occurrence; D(C j ) is codon C j Internal rating of w is the weight coefficient; The calculation formula of the internal score is as follows: D(C j )=pair_values(C j [1]C j [2])+pair_values(C j [2]C j [3]) Where, C j [1], C j [2], C j [3] respectively represent codon C j bases 1, 2, and 3; pair_values(C j [1]C j [2]) indicates codon C j The connection between the first and second bases is assigned; pair_values(C j [2]C j [3]) indicates codon C j The connection between the second and third bases is assigned.

4. The method according to claim 3, wherein Such as C j [3] is G or C, then the internal score is D(C j ) multiplied by 1.

5.

5. The method according to claim 1, wherein In step S4, the maximum sequence score DP[i][j] corresponding to the j-th synonymous codon of the i-th amino acid is calculated as follows: Where, DP[i-1][m] represents the maximum sequence score that the first i-1 amino acids can obtain when the i-1 amino acid uses its m-th synonymous codon; k i-1 Indicates the number of synonymous codons corresponding to the i-1th amino acid; C (i-1)m Indicates the mth synonymous codon corresponding to the i-1th amino acid; C ij Indicates the jth synonymous codon corresponding to the i-th amino acid; V(C (i-1)m ,C ij ) indicates codon C (i-1)m and C ij The adjacent codons between them are assigned values; f(C ij ) is codon C ij The frequency of codon occurrence; D(C ij ) is codon C ij Internal rating of w is the weight coefficient.

6. The method according to claim 5, wherein In step S5, the maximum score T corresponding to the full-length mRNA coding sequence max Obtained by the following formula: Where n is the number of amino acids in the full-length sequence.

7. The method according to claim 5, wherein In step S5, the optimal codon combination is obtained by the following formula: {C1,C2,...,C n }=Backtrack(j * ); Wherein, j* is the corresponding codon that maximizes the sequence score DP[n][j] in the full length of n amino acids, which is obtained by the following formula: j*=arg max j DP[n][j]。 8. The method according to claim 1, wherein The method further comprises the steps of: S6) Eliminating the negative factor sequence from the mRNA coding sequence obtained in step S5 by synonymous codon replacement, wherein the negative factor sequence includes an endonuclease cleavage site and / or a low information entropy repeat region sequence; The low information entropy repeat region sequence refers to a sequence fragment in which the repeat region exceeds the length of two codons, wherein the repeat region refers to a fragment formed by the repeated connection of multiple two bases or one base.

9. A device for mRNA sequence optimization, comprising an input module, a processing module, and an output module, wherein: The input module is used to input: target species and / or target cells, and the mRNA sequence to be optimized or the amino acid sequence encoded by it; The processing module performs mRNA sequence optimization according to the method of claim 1; The output module is used to output the optimized mRNA sequence.

10. A computer-readable storage medium, characterized in that A computer program for implementing the method according to claim 1 is stored thereon.

Images