Dynamic programming algorithm-based mRNA sequence optimization method

Abstract

The invention discloses an mRNA sequence optimization method based on a dynamic programming algorithm, and relates to the field of mRNA sequence optimization. According to the method, related positive (promotion) negative (hindering) factors are respectively considered in the process of translating mRNA into a protein sequence after the mRNA enters a host cell, an optimal comprehensive effect is obtained by adopting a dynamic programming algorithm, the final protein translation yield of the mRNA is greatly improved, and the method has a good application prospect 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 an mRNA sequence optimization method 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 and high safety, and can be widely used in infectious disease vaccines, tumor immunity, protein replacement and other fields. However, due to the defect that mRNA is more easily degraded than DNA, in order to make it translate and produce more proteins, a more appropriate mRNA sequence design is necessary.

[0003] Among the common optimization technologies currently available 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 load, etc. These factors can be divided into two categories at the translation level: promotion or hindrance.

[0004] However, the current methods 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 the species. As can be seen from the calculation process of CAI, it only considers the performance of a single codon, and does not consider the connectivity of the context of the codon in the sequence. The GC content actually indicates the reduction of uracil (U), a core key base that can be cut. However, too high a GC content will cause the formation of low entropy repeat sequences in local areas. Low entropy repeat regions cause the sliding window effect of the polymerase in the in vitro transcription (IVT) step, and the quality inspection of the synthesized mRNA cannot be performed. Therefore, the optimization algorithm should not only achieve the reduction of uracil, but also prevent the local formation of low entropy repeat sequences. Recent reports have shown that the mRNA structure and the corresponding derivative calculated index (MFE) have very little blocking effect on the ribosome, and the main blocking effect on the ribosome is the abundance of tRNA corresponding to the codon. In the intracellular environment, the tertiary structure of mRNA is unstable and unpredictable. There may be degradation-sensitive structures recognized by protein kinase R (PKR), 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. A codon may have multiple synonymous codons that can be replaced. Taking 100 codons as an example, each codon has an average of three synonymous codons. If an algorithm that exhaustively enumerates all synonymous codons is used, the calculation amount is 3 100 , this number of calculations is unattainable.

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

[0006] The purpose of the present invention is to provide a mRNA optimization method.

[0007] In a first aspect of the present invention, a method for optimizing an mRNA sequence is provided, wherein the mRNA is used for expression in a target species or a target cell, 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 cell, and obtaining a codon set for replacement according to the codon occurrence frequency;

[0009] 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 cells;

[0010] S2) provides a sensitivity index for a dibase pair RNA endonuclease;

[0011] 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;

[0012] S4) providing an mRNA coding sequence to be optimized or an amino acid sequence encoded by the mRNA coding sequence, and 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 transfer 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, the maximum sequence score corresponding to the full-length mRNA coding sequence is obtained, thereby obtaining the optimal codon combination, that is, the optimized mRNA coding sequence,

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

[0015] In another preferred embodiment, in step S1, the target species CDS sequence is all protein-encoding genes in the target species; and the target cell CDS sequence is the top 1000 protein-encoding genes in the target cell type.

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

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

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

[0019] In another preferred embodiment, in step S2, the sensitivity index of the dibase pair RNA endonuclease represents the sensitivity of the phosphodiester bond between every two bases to the RNA endonuclease, and the higher the sensitivity, the lower the sensitivity index.

[0020] In another preferred example, the sensitivity index range of the two-base pair RNA endonuclease is [-10, 10], and the sensitivity index of each two-base combination is: 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.

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

[0022] In another preferred embodiment, in step S3, the calculation formula for a 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 ) scores for individual codons;

[0026] f(C j ) is codon C j The codon frequency;

[0027] D(C j ) is codon C j Internal rating of

[0028] w is the weight coefficient;

[0029] The calculation formula of 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] In the formula,

[0032] C j [1] C j [2] C j [3] represents the first, second, and third bases of codon Cj, respectively;

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

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

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

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

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

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

[0039] In another preferred embodiment, in step S4, the connection value of adjacent codons is assigned V(C j ,C k ) is calculated as follows:

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

[0041] In the formula,

[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 ith amino acid, the possible synonymous codons of the i-1th 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] In the formula,

[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 represents 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 codon frequency;

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

[0055] w is the weight coefficient.

[0056] In another preferred embodiment, 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 the score is traced back to obtain the optimal codon combination.

[0057] In another preferred embodiment, the maximum score T corresponding to the full-length mRNA coding sequence max Obtained by the following formula:

[0058]

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

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

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

[0062] 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:

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

[0064] In another preferred embodiment, the method further comprises the steps of:

[0065] S6) Eliminate 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.

[0066] In another preferred embodiment, the restriction endonuclease cleavage site is a sapI restriction endonuclease cleavage site.

[0067] In another preferred embodiment, 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 repeated connection of multiple di-bases or mono-bases.

[0068] In another preferred embodiment, the elimination is performed using an exhaustive method, and all possible synonymous codon substitutions are tested to see whether new negative factor sequences are introduced.

[0069] In another preferred embodiment, the mRNA coding sequence with the maximum sequence score is selected from the mRNA coding sequences from which the negative factor sequences have been eliminated, thereby obtaining an optimized mRNA coding sequence; the sequence score is calculated by the method described in step S4.

[0070] In another preferred embodiment, step S4 further comprises:

[0071] The mRNA coding sequence to be optimized or the amino acid sequence it encodes is preprocessed.

[0072] In another preferred embodiment, performing sequence preprocessing on the amino acid sequence comprises:

[0073] 1) Determine whether all the letters in the sequence are the abbreviations of amino acids;

[0074] 2) Determine whether the first letter of the sequence is methionine -M. If not, place M in the first place.

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

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

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

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

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

[0080] In another preferred embodiment, the preprocessing further includes determining the length of the CDS sequence, such as when the CDS sequence length is ≥ 2000 nt, selecting the mRNA of a shorter protein sequence encoded by the same gene as the input sequence.

[0081] In another preferred embodiment, the method further comprises: comparing the amino acid sequence translated from the CDS sequence before pretreatment with the amino acid sequence translated from the optimized mRNA coding sequence, thereby verifying the correctness of the optimized sequence.

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

[0083] The input module is used to input: target species and / or target cells, and mRNA sequences to be optimized or amino acid sequences encoded by them;

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

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

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

[0087] In a fourth aspect of the present invention, a computer-readable storage medium is provided, on which a computer program for implementing the method as 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, wherein 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-mentioned technical features of the present invention and the technical features specifically described below (such as embodiments) can be combined with each other to form a new or preferred technical solution. Due to space limitations, they will not be described one by one here. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] The following drawings are used to illustrate specific embodiments of the present invention and are not used to limit the scope of the present invention defined by the claims.

[0091] Figure 1 A schematic diagram of the system of the present invention is shown.

[0092] Figure 2 A flow chart of the optimization method of the present invention is shown.

[0093] Figure 3 A schematic diagram of the mRNA structure is shown. DETAILED DESCRIPTION

[0094] The inventors have developed a mRNA sequence optimization method based on a dynamic programming algorithm for the first time after extensive and in-depth research. The method of the present invention finds the optimal balance between mRNA anti-degradation and translation speed, the optimized mRNA has no repeated low information entropy sequence fragments, and there is no sapI restriction site in the sequence. Relative to 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 a full-length sequence. At the same time, the present invention uses a dynamic programming algorithm to find the optimal codon combination, which saves a lot of computing power and improves the computing speed relative to the traditional exhaustive method. On this basis, the present invention is completed.

[0095] Factors affecting sequence optimization

[0096] 1. Codon frequency

[0097] In organisms, DNA is transcribed to generate mRNA sequences. The translation region (coding sequence, CDS) of the mRNA sequence is composed of codons connected in sequence. Each codon is composed of three bases. Each codon corresponds to a unique amino acid. An amino acid may have multiple codons corresponding to it. Multiple codons corresponding to the same amino acid are called synonymous codons. The relative fitness (W) of a single codon that makes up the translation region (CDS sequence) of mRNA or the frequency of occurrence per thousand codons characterizes the usage frequency of each codon, which corresponds to the abundance of transfer RNA (tRNA) in the cell. The usage frequency corresponding to each codon is species-specific. During the synthesis of protein peptide chains, the codons on the mRNA bind to the tRNA molecules under the action of the ribosome, and the amino acids corresponding to the codons are sequentially connected to synthesize protein peptide chains. Different transfer RNAs have different abundances in cells, which affects the efficiency of amino acids binding to polypeptides.

[0098] 2. The third base of the codon

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

[0100] 3. RNA endonuclease

[0101] The structure of RNA is Figure 3As shown, the CDS region is located in the inner middle of the entire sequence, rather than at both ends. The degradation mechanism associated with it mainly occurs through the opening process of the phosphodiester bond by ribonuclease (RNA) endonucleases. According to investigations, the relevant RNA endonucleases include RNase A, RNase T1, RNase T2, and RNase L. RNAase 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 should be paid to RNase A, RNase T2, and RNaseL. The sequence characteristics of the cleavage of these three endonucleases are as follows: RNase A mainly cuts the phosphodiester bond after pyrimidine (C-cytosine, U-uracil); RNase T2 mainly cuts the phosphodiester bond between GU and AU; RNase L mainly cuts UU and UA. These three RNA endonucleases can also cut the connections of other base residues, but the efficiency is not high.

[0102] Dynamic Programming Algorithm

[0103] The present invention optimizes mRNA sequences based on a dynamic programming algorithm, which saves computing power and improves computing speed compared to the traditional exhaustive method. The dynamic programming algorithm of the present invention is implemented using a state transition equation, specifically including:

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

[0105] 2) State transfer: based on the selection of synonymous codons and their corresponding scores, the cumulative score is updated and the optimal path is recorded;

[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, all possible sequences are sorted according to the cumulative scores, and the one with the highest score is selected as the optimal solution.

[0108] Sequence Optimization Methods

[0109] The present invention provides an mRNA sequence optimization method based on a dynamic programming algorithm, which finds the optimal balance between anti-degradation and translation speed through the dynamic programming algorithm. The method of the present invention comprises the following steps:

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

[0111] S2: construct the sensitivity index of dibase pair ribonuclease (RNA). In the present invention, the degradation resistance of the sequence is defined by the sensitivity index of dibase pair ribonuclease (RNA). RNA endonucleases have different cutting characteristics, such as RNAaseA efficiently cuts at sensitive junctions such as pyrimidine (C, U) and RNAase T2 at GU, AU, etc., while RNAase L tends to efficiently cut UU and UA sequences. According to the upstream and downstream relationships of these enzymes in the intracellular reaction pathway, cutting efficiency, molecular similarity and other factors, the present invention defines the endonuclease sensitivity index for 16 dinucleotide combinations. In one embodiment, the custom assignment 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, and GG is 10.

[0112] S3: Sequence 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 the coding sequence (CDS) inside the messenger RNA (mRNA) is processed, replace thymine (T) with uracil (U), verify whether its length is a multiple of 3 and whether the base composition unit is ["A","U","C","G"]. The start codon is replaced with AUG, and the stop codon is replaced with UGA. If the translation sequence (CDS) of the mRNA to be optimized is too long (>=2knt), even after optimization design, the inevitable accumulation of negative factors in the sequence may cause it to degrade in the cell and lose its design value. Therefore, in the presence of shorter protein sequences, it is necessary to select mRNA with a shorter protein sequence encoded by the same gene from uniprot for design.

[0113] S4: Codon scoring. In order to translate and produce more proteins, the mRNA sequence needs two properties: 1) resistance to degradation, 2) high-speed translation by the ribosome; the property of resistance to degradation has been defined in step S2, and the main property of high-speed ribosome translation, the frequency of codon occurrence, has been calculated in S1. The entire sequence is composed of single codons connected in sequence, so calculating the sum of the two properties of steps S1+S2 within each codon is very necessary as a state transition 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. In summary, 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] represents the first, second, and third bases of codon Cj, respectively;

[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 by the sensitivity index of the above-mentioned two-base pair endonuclease.

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

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

[0122] Final codon scoring formula:

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

[0124] Symbol definitions:

[0125] S(C j ) scores for individual codons;

[0126] f(C j ) is codon C j The codon frequency;

[0127] D(C j ) is codon C j Internal rating of

[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 translation efficiency has a greater impact on the final yield, the weight coefficient is larger, and if the degradation rate has a greater impact, the weight coefficient is smaller. The weight coefficient w ranges from 0.1 to 1, preferably, w is selected from 0.1, 0.2, 0.3, 0.4, .05, 0.6, 0.7, 0.8, 0.9 or 1.

[0130] S5: Dynamic programming algorithm calculates the optimal sequence. The entire sequence consists of a single codon connected in sequence. Except for the start and stop codons, each codon connects 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 the 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 to ensure that the selection of each step 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, the following steps are included:

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

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

[0133] In the formula,

[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] Wherein, 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 Indicates the number of synonymous codons corresponding to the i-1th amino acid;

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

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

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

[0144] f(C ij ) is codon C ij The codon frequency;

[0145] D(C ij ) is codon C ij internal rating.

[0146] 3) Calculate the optimal score:

[0147]

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

[0149] 4) Optimal codon sequence backtracking

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

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

[0152] Through a well-defined scoring system and a dynamic programming model, this method can systematically select the optimal codon sequence based on the consideration of codon resistance to endonucleolytic degradation, codon usage frequency, and inter-codon linkage effects.

[0153] S6: Eliminate 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 cells used for in vitro transcription (IVT). In one embodiment, the restriction site includes a 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, the template DNA fragment will be cut and eventually degraded. When using synonymous codon replacement to eliminate the restriction site fragment, since there are few codons involved, all possible methods can be listed by exhaustive method, and the overall sequence of the synonymous replacement scheme finally selected is scored the highest according to the aforementioned scoring method.

[0158] 2) Low information entropy repetitive region sequences. If there are multiple repeated bases in the sequence, polymerase is required for mRNA synthesis and subsequent product sequencing verification. However, the fidelity of polymerase in repetitive base fragments is poor. Since the coverage range of the ribosome is 2 codons, if the repetitive region exceeds the length of two codons, it needs to be removed. "Repetitive region" refers to a fragment composed of multiple di-bases or one-base repeats. Similarly, the synonymous codon with the highest score is selected for replacement by 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 final optimized sequence, and verifying the correctness of the final optimized sequence.

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

[0161] 1) The present invention selects the factors that promote and hinder the process of mRNA translation into protein, and finally selects the frequency of codon occurrence, the third base of the codon, the connectivity of the adjacent upstream and downstream codons, and the elimination of restriction sites and low information entropy repeat regions as the main sequence optimization points. The optimized mRNA sequence U base content and the use of high-frequency codons reach the optimal balance, there is no repeated low information entropy sequence fragment, and there is no sapI restriction site 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 adopts a dynamic programming algorithm to find the optimal codon combination, which saves a lot of computing power and improves the computing speed compared with the traditional exhaustive method.

[0164] The present invention will be further described below in conjunction with specific examples. It should be understood that these examples are intended to illustrate the present invention only and are not intended to limit the scope of the present invention. The experimental methods in the following examples where specific conditions are not specified are usually performed 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 under conditions recommended by the manufacturer. Unless otherwise indicated, percentages and parts are weight percentages and weight parts.

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

[0166] The African swine fever protein I215L was transfected into HEK293T cells for translation and expression. HEK293T cells are human kidney-derived cells, so the host is human.

[0167] The I215L sequence is as follows, source:

[0168] https: / / www.ncbi.nlm.nih.gov / protein / AKO62822.1?report=fasta

[0169] >AKO62822.1 pI215L[African swine fever virus]

[0170] MVSRFLIAEYRHLIENPSENFKISVNENNITEWDVILRGPPDTLYEGGLFKAKVAFPPEYPYAPPKLTFTSEMWHPNIYPDGRLCISILHGDNAEEQGMTWSPAQKIDTI LLSVISLLNEPNPDSPANVDAAKSYRKYVYKEDLESYPMEVKKTVKKSLDECSPEDIEYFKNAASNVPPIPSDAYEDECEEMEDDTYILTYDDDEEEEDEEMDDE(SEQID NO:1)

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

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

[0173] Table 1 Codon 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 a sensitivity index for two-base pair ribonuclease (RNA) endonucleases. This index is mainly based on the known RNA endonucleases in eukaryotes, including RNAase A and RNase T2 families, and RNAase L cutting characteristics. RNAase A cuts efficiently at pyrimidine (C, U) and RNAase T2 cuts efficiently at sensitive junctions such as GU and AU, while RNAase L tends to cut UU and UA sequences efficiently. Other junction types are also cut by these endonucleases, but the efficiency is not as good as the above special sites.

[0176] There are 16 dinucleotide combinations. According to the upstream and downstream relationships of these enzymes in the intracellular reaction pathways, cleavage efficiency, molecular similarity and other factors, the custom value range here is [-10, 10], such as 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 dibase pair RNA endonucleases

[0178] Base pairs Assignment Base pairs Assignment Base pairs Assignment Base pairs 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: Basic sequence 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, M is placed first. If the coding sequence (CDS) inside messenger RNA (mRNA) is processed, replace thymine (T) in the sequence with uracil (U), verify whether its length is a multiple of 3 and whether the base unit is ["A","U","C","G"]. The start codon is replaced with AUG and the stop codon is replaced with UGA. In order to verify the correctness of the final optimized sequence, the input CDS sequence is converted into an amino acid sequence to be compared with the amino acid sequence corresponding to the final optimized sequence. If the translation sequence (CDS) of the mRNA to be optimized is too long (>=2knt), even after optimization design, the inevitable accumulation of negative factors in the sequence may cause it to degrade in the cell and lose its design value. Therefore, if there is a shorter protein sequence, it is necessary to select the mRNA with a shorter protein sequence encoded by the same gene from uniprot for design.

[0180] I215L is only 215 AA in length, and each amino acid letter corresponds correctly. The length does not exceed 2knt, which meets the requirements.

[0181] S4: Codon scoring. In order to translate and produce more proteins, the mRNA sequence needs 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 ribosome translation-the frequency of codon occurrence-is calculated in S1. The entire sequence is composed of single codons connected in sequence, so calculating the sum of the two properties of steps S1+S2 within each codon is very necessary as a state transition 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.

[0182] Table 3 Codon internal scoring

[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] Rare codons have been removed.

[0185] S5: Dynamic programming algorithm calculates the optimal sequence. The entire sequence is composed of a single codon connected in sequence. Except for the start and stop codons, each codon is connected to two codons upstream and downstream, that is, 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(Cj) of each codon, the connection assignment V(Ci,Ci+1) between codons is also an important factor affecting the overall score. The state transition equation is used to maximize the overall score.

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

[0187] AUGGUGAGCCGGUUCCUGAUCGCCGAGUACCGGCACCUGAUCGAGAACCCCAGCGAGAACUUCAAGAUCAGCGUGAACGAACAACAUCACCGAGUGGGACGUGAUCCUGCGGGG CCCCCCC GACACCCUCUACGAGGGCGGCCUCUUCAAGGCCAAGGUGGCCUU CCCCCCC GAGUACCCCUACG CCCCCCCC AAGCUGACCUUCACCAGCGAGAUGUGGCACCCCAACAUCUACCCCGACGGCCGGCUCUGCAUCAGCAUCCUGCACGGCGACAACGCCGAGGAGCAGGGCAUGACCUGGAGCCCCGCCCAGAAGAUCGACACCAUCCUGCUGAGCGUGAUCAGCCUGCUGAACGAGCCCAACCCCGACAGCCCCGCCAACGUGGACGCCGCCAAGAGCUACCGGAAGUACGUCUACAAG GAGGACCUGGAGAGCUACCCCAUGGAGGUGAAGAAGACCGUGAAGAAGCCUGGACGAGUGCAGCCCCGAGGACAUCGAGUACUUCAAGAACGCCGCCAGCAACGUGCCCCCCAUCCCCAGCGACGCCUACGAGGACGAGUGCGAGGAGAUGGAGGACGACACCUACAUCCUGACCUACGACGACGACGAGGAGGAGGAGGACGAGGAGAUGGACGACGAGUGA (SEQ ID NO:2)

[0188] S6: Eliminate negative factor fragment sequences. 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 underline in SEQ ID NO:2 shows the low information entropy repeat region sequence fragment, and the underline in SEQ ID NO:3 shows the corresponding sequence fragment after optimization.

[0189] AUGGUGAGCCGGUUCCUGAUCGCCGAGUACCGGCACCUGAUCGAGAACCCCAGCGAGAACUUCAAGAUCAGCGUGAACGAACAACAUCACCGAGUGGGACGUGAUCCUGCGGGG CCCACCC GACACCCUCUACGAGGGCGGCCUCUUCAAGGCCAAGGUGGCCUU CCCACCC GAGUACCCCUACG CCCCACCC AAGCUGACCUUCACCAGCGAGAUGUGGCACCCCAACAUCUACCCCGACGGCCGGCUCUGCAUCAGCAUCCUGCACGGCGACAACGCCGAGGAGCAGGGCAUGACCUGGAGCCCCGCCCAGAAGAUCGACACCAUCCUGCUGAGCGUGAUCAGCCUGCUGAACGAGCCCAACCCCGACAGCCCCGCCAACGUGGACGCCGCCAAGAGCUACCGGAAGUACGUCUACAAG GAGGACCUGGAGAGCUACCCCAUGGAGGUGAAGAAGACCGUGAAGAAGCCUGGACGAGUGCAGCCCCGAGGACAUCGAGUACUUCAAGAACGCCGCCAGCAACGUGCCCCCAAUCCCCAGCGACGCCUACGAGGACGAGUGCGAGGAGAUGGAGGACGACACCUACAUCCUGACCUACGACGACGACGAGGAGGAGGAGGACGAGGAGAUGGACGACGAGUGA (SEQ ID NO:3)

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

[0191] The African swine fever protein pB602L was transfected into HEK293T cells for translation and expression. HEK293T cells are human kidney-derived cells. The host is also human.

[0192] The pB602L sequence is as follows, source:

[0193] https: / / www.ncbi.nlm.nih.gov / protein / AAA65309.1?report=fasta

[0194] >AAA65309.1 pB602L[African swine fever virus]

[0195] (SEQ ID NO:4)

[0196] In this sequence optimization step, S1, S2, and S4 are consistent with the previous steps.

[0197] S3 run results: pB602L is 602 AA in length, and each amino acid letter corresponds correctly. The length of the corresponding nucleotide is 1806 nt, which does not exceed 2knt. The sequence meets the requirements.

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

[0199]

[0200] S6: Eliminate negative factor fragment sequences. 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. After optimization, these two negative factor fragment sequences do not exist.

[0201] All documents mentioned in the present invention are cited as references in this application, just as each document is cited as reference individually. In addition, it should be understood that after reading the above teachings of the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms also fall within the scope defined by the claims attached to this application.

Claims

1. A method for optimizing 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 according to 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 cells; 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 the mRNA coding sequence, and 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 transfer 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, that 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, characterized in that 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, characterized in that In step S3, the calculation formula for a single codon score is as follows: S(C j )=D(C j )+w·f(C j ) In the formula, S(C j ) scores for individual codons; f(C j ) is codon C j The codon frequency; 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]) In the formula, C j [1] C j [2] C j [3] represents the first, second, and third bases of codon Cj, respectively; p air_values(C j [1]C j [2]) represents the connection assignment between the first and second bases of codon Cj; pair_values(C j [2]C j [3]) represents the connection assignment between the second and third bases of codon Cj.

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

5.

5. The method according to claim 1, characterized in that 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: In the formula, 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 codon frequency; D(C ij ) is codon C ij Internal rating of w is the weight coefficient.

6. The method according to claim 5, characterized in that 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, characterized in that 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, characterized in that 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 whose repeat region exceeds the length of two codons, wherein the repeat region refers to a fragment formed by repeated connection of multiple dibases or monobases.

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 mRNA sequences to be optimized or amino acid sequences encoded by them; 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 of claim 1 is stored thereon.

Images