Simulation method and application of loop region and stem region of RNA (Ribonucleic Acid) pseudoknot

Abstract

The invention discloses a simulation method and application of a loop region and a step region of an RNA (Ribonucleic Acid) pseudoknot, and belongs to the field of bioinformatics. The thermodynamic and dynamic properties of the H-type false junction are studied through long-time high-temperature molecular dynamics simulation, and the thermodynamic and dynamic parameters of the base pair at the tail end of the stem region are obtained by utilizing the conversion balance between opening and closing of the base pair and analyzing a simulation track. And due to the existence of the loop region, the melting temperature of the stem region is reduced. However, the enthalpy change and the entropy change are increased due to the existence of the loop region. And under the same simulation temperature, the opening rate of the loop region is increased, and the closing rate is reduced. The method can effectively obtain the optimal solution of the problem, is high in prediction efficiency, and overcomes the defects of an existing method.

Description

Technical Field

[0001] The invention belongs to the field of bioinformatics and relates to a simulation method and application of a loop region and a stem region of an RNA pseudoknot. Background Art

[0002] RNA molecules perform a variety of biological functions in cells, from the transfer of genetic information, to serving as an intermediate in the transfer of DNA to proteins, and to enzyme catalysis. The three-dimensional structure of RNA molecules is largely determined by the interaction between base-paired nucleotide pairs. The secondary structure of an RNA molecule is a collection of base pairs in a three-dimensional structure. Therefore, the secondary structure can be used as a basis for predicting higher structural levels. The primary structural information of RNA can be obtained in large quantities through experiments, but obtaining the three-dimensional structure of RNA through experiments will make it costly and difficult, and not all three-dimensional structures can be measured experimentally.

[0003] With the continuous development of computers, people can use theoretical methods to predict the structure of RNA. The accuracy of the prediction depends on the assumptions of the folding structure and the potential energy model. In structure prediction, a common method is the maximum base pairing method, which is an algorithm developed based on the structural conservation of RNA and mainly relies on the sequence alignment effect. However, the main disadvantage of this method is the sequence alignment of the conserved region. Another method is the minimum free energy method, which mainly looks for a complementary structure in an energy-stable state. With continuous development and research, there are more and more methods for predicting RNA secondary structure.

[0004] Pseudoknot is a typical tertiary structure. Pseudoknot has different roles in biological functions, such as self-cleavage of introns, formation of catalytic cores of various ribozymes, telomerase, riboswitch and ribosomal frameshift. The biological function of pseudoknot is directly related to its folding stability and conformational changes. In RNA pseudoknot, stability and folding thermodynamics are mainly determined by the interaction between the loop region and the stem region. For example, pseudoknot involves several intermediates in the folding process. Studying how RNA pseudoknot structure folds and forms, misfolded intermediates may be formed during the folding process, so establishing a correct model is very important for predicting the structure of RNA pseudoknot.

[0005] In the prior art, the thermodynamic parameters of the classical RNA secondary structure were measured and estimated. Similarly, the thermodynamic parameters of some model structures of RNA pseudoknots were measured, indicating that the ion concentration in the solution is very important for the stability of the RNA pseudoknot. And the conformational transition involved in the folding process of the pseudoknot may include intermediate steps, mainly in the process of forming a hairpin. Therefore, it is very difficult to study the contribution of the loop region of the pseudoknot structure to the structural stability, especially the pseudoknot containing at least two loop structures. For the problem of pseudoknots, many algorithms have been developed to try to make accurate predictions. However, due to the compact structure of the RNA pseudoknot, non-bonded interactions between atoms will be formed in the pseudoknot, such as van der Waals interaction, electrostatic interaction, and volume repulsion. Therefore, it is inaccurate to predict the pseudoknot structure using thermodynamic parameters in the secondary structure. Summary of the invention

[0006] In view of the problems existing in the prior art, the present invention takes H-type RNA pseudoknot as the research basis, and provides a simulation method and application of the loop region and stem region of an RNA pseudoknot. The H-type RNA pseudoknot is composed of two helices forming a stem (stem) and two loops (loops) spanning the helical stem. The RNA pseudostructure is defined as the pairing between the structural loop with bases and other regions of the RNA. In addition, the secondary structure information of RNA is a key step in predicting its tertiary structure. The essence of predicting the secondary structure of RNA is to find out the pairing of bases in the stem region in the secondary structure and the thermodynamics of the loop region on the basis of the primary structure. The present invention uses all-atom dynamics simulation to study the influence of the loop region in the H-type pseudoknot on the thermodynamics and kinetic properties of the stem region.

[0007] In one aspect, the present invention provides a method for simulating the loop region and stem region of an RNA pseudoknot, which specifically comprises the following steps:

[0008] s1. Selection of simulation methods: All-atom force field and molecular dynamics simulation software are used to simulate RNA structure;

[0009] s2. Choice of RNA pseudoknot structure: H-type RNA pseudoknot;

[0010] s3. Construction of the simulation system: mutate the bases at the end of the H-type RNA pseudoknot stem region, balance the charge of the simulation system, set the ion concentration, and build the simulation system in a water box;

[0011] s4.Simulation parameter settings: handle electrostatic interactions and intermolecular interactions, constrain bond lengths, and calculate equations of motion.

[0012] Further, the crystal structure of the H-type RNA pseudoknot described in the present invention is PDB ID: 1hvu;

[0013] The PDB ID: 1hvu includes stem1 region, stem2 region, loop1 region and loop2 region;

[0014] The stem1 region contains 5 base stacks; the stem2 region contains 6 base stacks;

[0015] The loop1 region contains 2 nucleotides; the loop2 region contains 3 nucleotides.

[0016] Furthermore, the loop1 region in the PDB ID: 1hvu described in the present invention is located in the major groove formed by the stem2 region, and the loop2 region passes through the minor groove formed by the stem1 region.

[0017] Furthermore, in the method for simulating the loop region and stem region of an RNA pseudoknot described in the present invention, the construction of the simulation system includes: balancing the charge of the simulation system by adding Na + Make the simulated system charge neutral;

[0018] The set ion concentration includes the random addition of Na at the initial position of the crystal structure of PDB ID: 1hvu + and Cl - ;

[0019] The method of maintaining the water molecules in the simulation system completely rigid is to use the TIP3P water model in combination with the SETTLE algorithm;

[0020] The simulated water cube box is The simulated water cube box utilizes periodic boundary conditions to ensure that the simulated water cube box has three or four solvent layers in each direction.

[0021] The simulation system contains 26356 atoms.

[0022] Furthermore, in the method for simulating the loop region and stem region of an RNA pseudoknot of the invention, the simulation parameter settings include: the cutoff distance of the electrostatic interaction is

[0023] The cutoff distance for the intermolecular interaction is

[0024] The bond lengths of the solutes were updated every ten steps using a neighboring grid search method;

[0025] The time step of the equation of motion is 2fs and the coordinates are saved every 2ps.

[0026] The conformation of the simulation system at 290K is the initial structure.

[0027] On the other hand, the present invention provides an application of a method for simulating the loop region and stem region of an RNA pseudoknot in predicting the RNA pseudoknot structure.

[0028] Furthermore, the present invention predicts that the RNA pseudoknot structure includes a pair of nucleotides at the end of the open stem region and nucleotides in the loop region, and calculates the opening and closing probabilities of the base pairs in the stem region as well as the entropy change and enthalpy change between the open state and the closed state.

[0029] Furthermore, the present invention calculated the RMSD (root mean square deviation) and dihedral angle of two bases (U8 and A27) in the stem2 region. When the RMSD of base U8 is greater than and the dihedral angle is between -140° and 100°, or for base A27, the RMSD is greater than And the torsion angle is between about 50° and 100°. If one or both of these two conditions are met, the conformation is defined as the open state. In this state, one of the bases U8 and A27 is thrown into the solvent and the base pair and base stacking are destroyed.

[0030] Furthermore, the present invention calculates the thermodynamic parameters of the open state and the closed state and finds that the existence of the loop region reduces the melting temperature of the stem region, and the obtained enthalpy change becomes larger, and the entropy change also becomes larger.

[0031] Furthermore, the present invention calculated the kinetic parameters of the open state and the closed state and found that the presence of the loop region did not change the opening rate, but made the closing rate smaller.

[0032] Compared with the prior art, the technical solution provided by the present invention has at least the following beneficial effects or advantages:

[0033] The present invention studies the thermodynamic and kinetic properties of the H-type pseudoknot through long-term high-temperature molecular dynamics simulation, and obtains the thermodynamic and kinetic parameters of the base pairs at the end of the stem region by utilizing the conversion balance between the opening and closing of the base pairs, and by analyzing the simulation trajectory. The following results are obtained: the presence of the loop region reduces the melting temperature of the stem region. However, the presence of the loop region increases the enthalpy change and the entropy change. And at the same simulation temperature, the presence of the loop region increases the opening rate and decreases the closing rate. These conclusions show that the loop region has an impact on the thermodynamics and kinetics of the stem region. Therefore, the present invention provides reliable parameters for future structural predictions by studying the impact of the loop region on the stem region. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 The crystal structure of PDB ID: 1hvu.

[0035] Figure 2 is the relationship between base U8 and simulation time. a is the relationship between the RMSD of base 8 and simulation time; b is the relationship between the dihedral angle formed by base 8 and base 9 and simulation time; the simulation temperature is 360K, and the simulation time is 200-800ns.

[0036] Figure 3 is the change of RMSD and phosphate backbone torsion angle of base A27. Among them, a is the change of RMSD of base A27 with simulation time; b is the distribution of RMSD of base A27 during the whole simulation process; c is the distribution of phosphate backbone torsion angle of base A27 during the whole simulation process; simulation temperature is 360K; total simulation time is 10μs.

[0037] Figure 4 is the probability of the closed state in the simulation time at different simulation temperatures and ln(p o / p c ). Among them, a is the probability of the closed state during the simulation time at a simulation temperature of 360K to 400K; b is the ln(p o / p c ).

[0038] Figure 5 This is a bar chart of the average lifetime of the open and closed states at simulation temperatures of 360K to 400K.

[0039] Figure 6 This is a bar chart of the average lifetime distribution of the open and closed states at simulation temperatures of 360K to 400K.

[0040] Among them, a is the average lifetime distribution of the closed state at 360 K; b is the average lifetime distribution of the closed state at 370 K; c is the average lifetime distribution of the closed state at 380 K; d is the average lifetime distribution of the closed state at 390 K; e is the average lifetime distribution of the closed state at 400 K; f is the average lifetime distribution of the open state at 360 K; g is the average lifetime distribution of the open state at 370 K; h is the average lifetime distribution of the open state at 380 K; i is the average lifetime distribution of the open state at 390 K; j is the average lifetime distribution of the open state at 400 K. DETAILED DESCRIPTION

[0041] The technical scheme of the present invention is described below in conjunction with the embodiments, but the present invention is not limited to the following embodiments. The experimental methods and detection methods described in each embodiment are conventional methods unless otherwise specified; the reagents and materials described are all commercially available unless otherwise specified.

[0042] Example 1

[0043] This embodiment provides a method for studying the thermodynamic parameters and kinetic characteristics of the loop region on the stem region.

[0044] In this example, the Amber10 all-atom force field was selected during the simulation process. The simulation software used was Gromacs4.5.6, which was used to perform all molecular dynamics simulations. This example downloaded the crystal structure from the PDB library (PDB ID: 1hvu). Figure 1 As shown. The crystal structure of 1hvu contains two stem (S1 and S2) and two loop (L1 and L2) regions, which is the simplest H-type pseudoknot structure. The stem1 region has 5 base stacks, the stem2 region contains 6 base stacks, and loop1 and loop2 have 2 nucleotides and 3 nucleotides, respectively. Loop1 is in the major groove formed by the stem2 region, and loop2 passes through the minor groove formed by the stem1 region. In order to compare the thermodynamic and kinetic parameters with the 3'-AAUU-5' sequence when there is no loop, the U26A mutation in the stem2 region is used in this example. 26 Na+ molecules are added to the system. + The entire simulation system is charge neutral, and then Na is added at random initial positions. + and Cl - , and the concentration of ions in the field is 0.5 M. The TIP3P water model is used and the SETTLE algorithm is used to keep the water molecules completely rigid. The simulation water cube box is set to Periodic boundary conditions were used to ensure that there were three or four solvent layers in each direction. The final simulation system contained 26,356 atoms. The entire system was simulated at constant temperature and constant pressure using the random velocity and Parrinello-Rahman constant pressure algorithms, respectively. The Particle-Mesh Ewald method was used to deal with electrostatic interactions, where the cutoff distance for electrostatic interactions is The cutoff distance for the Lennard-Jones intermolecular interaction is also The bond length of the solute was constrained using the LINCS algorithm. A neighboring grid search method was used and updated every ten steps. The equations of motion were calculated using the Verlet algorithm. The time step was 2fs and the coordinates were saved every 2ps. At 290K, the system was energy minimized for 20ns and then the initial structure for the simulation was selected from the conformation. Due to the different lengths of the two loops, the melting temperatures of the bases in the stem region are also different. Therefore, in this example, stem2 and loop1 were first selected to study the influence of the loop region on the thermodynamic parameters and kinetic properties of the stem region. Except for the atoms on residues 8, residue 27, residue 6, and residue 7, the atoms of other residues on the nucleotide are fixed using additional forces (force constant: 1000kJ / mol·nm -2 ). In order to accelerate the conversion between the closed state and the open state, this embodiment uses high temperature simulation to simulate the entire system. The selected simulation temperatures are 360K, 370K, 380K, 390K and 400K, and each simulation temperature lasts for 10μs.

[0045] Example 2

[0046] This embodiment provides the movement trajectories of base U8 and base A27 at different simulation temperatures.

[0047] According to the motion trajectory of each simulated temperature, the base pair presents a transition between opening and closing. Therefore, in this embodiment, the RMSD (root mean square deviation) and dihedral angle (such as Figure 2 and Figure 3 ). For base U8 ( Figure 2 ), this embodiment calculates the dihedral angle C1'(i)-N1(i)-C1'(i+1)-N1(i+1) formed by the four atoms in the U8 base and the U9 base (such as Figure 2 As shown in b in FIG. 1 ). For base A27, this embodiment calculates the phosphate backbone torsion angle ζ of base A27, C3'(i)-O3'(i)-P(i+1)-O5'(i+1), where i represents the i-th nucleotide in the polynucleotide chain. This embodiment divides the simulated conformations into two categories: open state and closed state according to the calculated RMSD and the distribution of torsion angle ζ and dihedral angle.

[0048] The calculated RMSD of base U8 can be obtained (such as Figure 2 When the RMSD is about When RMSD is greater than When the base U8 flips into the solvent, and sometimes the base U8 falls into the major groove, and this state lasts for a long time.

[0049] The dihedral angle formed by base U8 and base U9 can be calculated through the simulation trajectory, C1'(i)-N1(i)-C1'(i+1)-N1(i+1) (such as Figure 2 When the central region of the dihedral angle is -160°, base U8 remains in the folded state and vibrates near the initial position. In this state, the RMSD corresponding to the conformation of base U8 is approximately When the dihedral angle is at 60° in the center, base U8 is thrown into the solvent and sometimes stuck in the groove for a long time. In this state, the RMSD of the corresponding conformation is greater than Base U8 is in a very active state.

[0050] According to the simulation trajectory, this embodiment also calculates the RMSD of base A27 and the phosphate backbone torsion angle C3'(i)-O3'(i)-P(i+1)-O5'(i+1). Figure 3 It can be seen that the RMSD is distributed in two regions, namely and When RMSD is about , base A27 remains in a folded state and vibrates around the start position. The other region is about In this state, base A27 deviates from its initial position and is thrown into the solvent, and base pairing and base stacking are destroyed. For the phosphate backbone torsion angle C3'(i)-O3'(i)-P(i+1)-O5'(i+1), base A27 is also distributed in two regions, approximately -75° and approximately 50°, respectively. When the torsion angle is approximately -75°, base A27 vibrates near its initial position. When the torsion angle is approximately 50°, base A27 is thrown into the solvent, and base pairing and base stacking are destroyed. When the RMSD is approximately When , the torsion angle corresponding to the conformation is about -75°. When RMSD is about When , the torsion angle corresponding to the conformation is about 50°. Therefore, the present invention divides the simulated conformation into two states: a closed state and an open state.

[0051] Closed state: When the RMSD of base U8 is about And the dihedral angle is about -80° to -140°. It is also required that for base A27, the RMSD is about And the torsion angle is between about -50° and -100°. In this state, base U8 and base A27 only vibrate at their initial positions, and the base pairing and base stacking between base U8 and base A27 are complete and not destroyed. The base pairing and base stacking of stem2 are complete, always maintaining 6bp.

[0052] Open state: when the RMSD of base U8 is greater than and the dihedral angle is between -140° and 100°, or for base A27, the RMSD is greater than And the torsion angle is between about 50° and 100°. If one or both of these two conditions are met, the conformation is defined as the open state. In this state, the base pair and base stacking are destroyed when one of the bases U8 and A27 is thrown into the solvent. And at least one base is thrown into the solvent. During the simulation, due to the influence of the loop region, the corresponding base in U8 can also be thrown into the solvent.

[0053] Example 3

[0054] This embodiment provides thermodynamic parameters and kinetic parameters of the open state and the closed state at different simulation temperatures.

[0055] At each simulation temperature, the probabilities of the open and closed states can be obtained as Where τ is the total simulation time τ = 10μs, Nc, No are the total number of occurrences of the corresponding conformations in the closed state and the open state respectively. c i , τ o i are the time periods that the conformation stays in the closed state and the open state for the i-th time, respectively. Figure 4 (a) shows the occupation probability of the closed state during the simulation at all simulation temperatures. The occupation probability of the closed state for all simulation temperatures up to 9600ns reaches a stable value, indicating that the transition balance between open and closed is achieved. To ensure that the simulation reaches the closed-open transition balance, each simulation temperature lasts until 10μs. Table 1 lists the probability of the open state and the closed state and the average lifetime at the simulation temperature.

[0056] Table 1 shows the average lifetime, number of times and probability of the closed and open states when the simulation time is 10μs and the simulation temperatures are 360K, 370K, 380K, 390K and 400K respectively.

[0057]

[0058] According to the statistics in Table 1, since the simulated conformations are divided into open and closed states, the free energy difference ΔG of the closure of U8 and A27 can be obtained by the probabilities of the closed and open states:

[0059]

[0060] The p in the formula o and p c are the probabilities of the open state and the closed state at the simulation temperature; k B is the Boltzmann constant.

[0061] like Figure 4 As shown in (b), It is linearly related to the inverse of temperature (1 / T), indicating the linear relationship between the free energy difference ΔG and the temperature T. ΔG = ΔH-TΔS, where ΔH and ΔS are the enthalpy and entropy changes of the transition between the closed state and the open state, respectively, and the thermodynamic parameters of the AU base pairs in the stem region of the pseudoknot are ΔH = -7.69 kcal / mol and ΔS = -21.46 eu.

[0062] The melting temperature at this time is 358K, while for the same terminal base pair without a pseudoknot, ΔH = -7.3 kcal / mol and ΔS = -18.5 eu, and the melting temperature at this time is 395K. Therefore, the existence of the loop region lowers the melting temperature of the stem region. And the enthalpy change obtained becomes larger, and the entropy change also becomes larger.

[0063] The average lifetime of the open and closed states can be obtained as follows: where τ ave is the average lifetime, N is the total number of base pairs formed by U8 and A27 in the stem region in the open or closed state, and τ i is the i-th lifetime of the conformation in the corresponding state. The average lifetime of the open and closed states at different simulation temperatures is as follows: Figure 5 As shown, Figure 6 The lifetime distribution of the open state and the closed state is shown. The average lifetime of the open state and the average lifetime of the closed state show different correlation behaviors with temperature. The simulated temperature has little effect on the average lifetime of the open state, and the open state is weakly correlated with the simulated temperature. The average lifetime of the closed state is greatly affected by temperature, and the closed state is strongly correlated with the simulated temperature.

[0064] For the open-closed two-state transition, the closing rate and opening rate of base 8 and base 27 in the stem region can be obtained from the average lifetime of the open state and the closed state as where τ o ave and τ c aveare the average lifetimes of the open state and the closed state, respectively. Relative to the average lifetimes of the open and closed states of a single base pair, it was found that at the same temperature, the average lifetime of the closed state became smaller, while the average lifetime of the open state became larger. This shows that the existence of the loop region did not change the opening rate, but made the closing rate smaller.

[0065] The embodiments described above are part of the embodiments of the present invention, but not all of them. The detailed description of the embodiments of the present invention is not intended to limit the scope of the invention claimed for protection, but only represents selected embodiments of the present invention. All other embodiments obtained without creative work and related deductions and substitutions made by ordinary technicians in the field under the conditions of the concept of the present invention belong to the scope of protection of the present invention.

Claims

1. A method for simulating the loop region and stem region of an RNA pseudoknot, characterized in that: The steps include: s1. Selection of simulation methods: All-atom force field and molecular dynamics simulation software are used to simulate RNA structure; s2. Choice of RNA pseudoknot structure: H-type RNA pseudoknot; s3. Construction of the simulation system: mutate the bases at the end of the H-type RNA pseudoknot stem region, balance the charge of the simulation system, set the ion concentration, and build the simulation system in a water box; s4.Simulation parameter settings: handle electrostatic interactions and intermolecular interactions, constrain bond lengths, and calculate equations of motion.

2. The simulation method according to claim 1, characterized in that: The crystal structure of the H-type RNA pseudoknot in s2 is PDB ID: 1hvu; The PDB ID: 1hvu includes stem1 region, stem2 region, loop1 region and loop2 region; The stem1 region contains 5 base stacks; the stem2 region contains 6 base stacks; The loop1 region contains 2 nucleotides; the loop2 region contains 3 nucleotides.

3. The simulation method according to claim 2, characterized in that: The loop1 region is located in the major groove formed by the stem2 region, and the loop2 region passes through the minor groove formed by the stem1 region.

4. The simulation method according to claim 1, characterized in that: The charge of the equilibrium simulation system in s3 is calculated by adding Na + Make the simulated system charge neutral; The set ion concentration includes the random addition of Na at the initial position of the crystal structure of PDB ID: 1hvu + and Cl - .

5. The simulation method according to claim 1, characterized in that: The water box in s3 uses periodic boundary conditions to ensure that there are three or four solvent layers in each direction of the water box.

6. The simulation method according to claim 1, characterized in that: The cutoff distance of the electrostatic interaction in s4 is The cutoff distance for the intermolecular interaction is 7. The simulation method according to claim 1, characterized in that: The bond length of the solute in s4 is updated every ten steps using a neighboring grid search method; The time step of the equation of motion is 2fs and the coordinates are saved every 2ps.

8. The simulation method according to claim 2, characterized in that: The conformation of the simulation system at 290K is the initial structure.

9. A simulation method for the loop region and stem region of an RNA pseudoknot and its application in predicting RNA pseudoknot structure.

10. The use according to claim 9, characterized in that: The predicted RNA pseudoknot structure includes a pair of nucleotides at the end of the open stem region and the nucleotides in the loop region, and the opening and closing probabilities of the base pairs in the stem region and the entropy change and enthalpy change between the open state and the closed state are counted.

Images