Team:OUC-China/miniToe Family

Team OUC-China: Main

miniToe Family

1. Enzymes Mutation

1.1 The Four Keys in miniToe System

The wet lab members give us four important sites, Gln104, Tyr176, Phe155, His29, which play important roles in binding and cleavage in protein Csy4. Considering 20 kinds of amino acids, we have 80 mutants to explore and choose if we only have one site mutated.

Before designing the protein mutants, we first looking into the working process of miniToe structure to find the most important keys in our system.

Fig.1-1 The working process of miniToe system

All the reactions happened in our first system, miniToe, can be described chronologically by following five main steps[1]:

(1)The miniToe structure is produced and accumulated.
(2)The Csy4 is produced with IPTG induced.
(3)The Csy4 binds to the miniToe structure and form the Csy4-miniToe complex
(4)The Csy4 cleave the special site and divide the miniToe structure into two parts: the Csy4-crRNA complex and the mRNA of sfGFP.
(5)The sfGFP is produced.

From the description above, we can get four key problems in our system to make sure that our system can work successfully:

(1)Does the Csy4 dock correctly with the miniToe structure (hairpin)?
(2)How about the binding ability between the Csy4 and miniToe structure (hairpin)?
(3)How about the cleavage ability between the Csy4 and miniToe structure (hairpin)?
(4)Does cis-repressive RNA release from the RBS?

The most impressive way to explore four problems is to model our system at the atom level by molecular dynamics[2]. And there are lots of work in exploring the Csy4-RNA complex by molecular dynamics[3][4][5].

1.2 Molecular Dynamics

Molecular dynamics (MD)[6] is a computer simulation method for studying the physical movements of atoms and molecules. The atoms and molecules are allowing to interact for a fixed period of time, giving a view of the dynamic evolution of the system. In the most common version, the trajectories of atoms and molecules are determined by numerically solving Newton's equations of motion for a system of interacting particles, whose forces between the particles and their potential energies are often calculated using interatomic potentials or molecular mechanics force fields.

To a system which consists of molecule or atoms, the total energy of a system includes kinetic energy and potential energy，which can be described by the formula below:

E = E_{k i n} + U

where the E_kin denotes the kinetic energy and the U denotes the potential energy.

In a molecule system, the total potential energy can be calculated by adding the U_nb、 bond stretching potentials energy U_b、angle bending potentials energy U_θ、torsion angle potentials energy U_ϕ、out-of-plane potentials energy U_χ and some other cross effect U_cross together，which also can be described by the formula below:

U = U_{n b} + U_{b} + U_{θ} + U_{ϕ} + U_{χ} + U_{c r o s s}

This formula above also called the force field in molecular dynamics’ theory. There are many force field in the world that based on the statistical thermodynamics and empirical result. In the research of the protein and nucleic acid, the Amber force field is one of the best force field in the world. So we choose Amber as our force field and the formula of this field show below: \[\begin{array}{l} E = \sum\limits_{bond} {{K_b}{{({r_{ij}} - {r_0})}^2}} + \sum\limits_{angle} {{K_\theta }{{(\theta - {\theta _0})}^2}} + \sum\limits_{dihedral} {\frac{{{K_\phi }[1 + \cos (n\phi - {\phi _0})]}}{2}} \\ + \sum\limits_{impr} {\frac{{{K_\chi }[1 + \cos (n\chi - {\chi _0})]}}{2}} + \sum\limits_{nobond} {{\varepsilon _{ij}}\left[ {{{\left( {\frac{{R_{ij}^0}}{{{R_{ij}}}}} \right)}^{12}} - 2{{\left( {\frac{{R_{ij}^0}}{{{R_{ij}}}}} \right)}^6}} \right]} + \sum\limits_{nobond} {\frac{{{q_i}{q_j}}}{{{R_{ij}}}}} \end{array}\]

The items in the formula refers to bond stretching term、angle bending potentials、dihedral angle potentials、put of plane angle potentials、improper dihedral angle potentials、Van Der Waals interaction and Coulombic interaction terms in order.

Now considering the system which contains of consists of N molecule or atoms in the classic mechanics, the atom or molecule can be characterized as follow:

The position of atom i is r_i , the speed of atom i is v_i , the acceleration of atom i is a.

After the integral operation, we can get two formula of v_i and r_i :

{\overset{⇀}{v}}_{i} = {\overset{⇀}{v}}_{i}^{0} + {\overset{⇀}{a}}_{i} t

{\overset{⇀}{r}}_{i} = {\overset{⇀}{r}}_{i}^{0} + {\overset{⇀}{v}}_{i}^{0} + \frac{1}{2} {\overset{⇀}{a}}_{i} t^{2}

where the v_i ⁰ refers to initial speed and the r_i ⁰ refers to initial position.

According to the classic mechanics, the force applied to atoms is the negative gradient of potential energy:

{\overset{⇀}{F}}_{i} = - \nabla_{i} U = - (\frac{\partial}{\partial x} \overset{⇀}{i} + \frac{\partial}{\partial y} \overset{⇀}{j} + \frac{\partial}{\partial z} \overset{⇀}{k}) U

Due to Newton’s second law, the acceleration of atom i can also be described as:

{\overset{⇀}{a}}_{i} = \frac{{\overset{⇀}{F}}_{i}}{m_{i}}

Using the total formula above together, we can have the full process of the molecular dynamics described and the Fig.2-2 shows the flow chart of MD in program.

Fig.1-2 The flow chart of MD in program

1.3 Logic Line

When choosing the Csy4 mutants, we choose the four problems which discuss before as the four keys, and we choose the molecular dynamics as our main tools to look into our system in atom level. So what’s is the logic line between them?

We considering the mutants choosing-problem by following description:

What we know and proved by the experiment is that the wild-type Csy4 with the miniToe system is working well, which means that all the important key problems we discussion did not exist in the wild-type Csy4. The wild-type Csy4 can dock correctly with the miniToe structure and the Csy4 have a good ability to bind and cleave the miniToe structure, finally the crRNA release from the RBS. So we choose the wild-type Csy4 as a standard, and all the Csy4 mutant can check the four key problems by comparing to wild-type Csy4.

So the four key problems transform into two main problems: how to describe four key problems and compare between the wild-type and mutant in mathematics? That’s what we discuss in the following two parts.

1.4 How to describe four key problems in mathematical form

We are going to check the molecular dynamics of wild-type in the miniToe system to show you how to describe four key problems in mathematical form in this part.

For the first problem, we define a matrix called interaction matrix which can describe the interaction possibility between every amino acid of the protein and every nucleic acid of the hairpin, this interaction matrix can be calculated by the catRAPID graphic[7]. We submit the wild-type Csy4 and miniToe structure respectively to online service and it can return us the interaction matrix. The Fig.1-3 is the heat map of the interaction matrix for wild-type Csy4 in the hairpin region.

Fig.1-3 The heatmap of interaction matrix for wild-type Csy4.

And the rest three problems is solved by the molecular dynamics. The work of molecular dynamics is mostly based on the Jiří Šponer’s work[8], but still something different.

In order to explore the rest three key problems, we prepare two structure for the molecular dynamics. And the geometries of our miniToe system is based on the X-ray structure of Csy4/RNA complex with the cleaved RNA (PDB ID: 4AL5, resolution 2.0 A), it can be seen in the Fig.1-4.

Fig.1-4 the X-ray structure of Csy4/RNA complex

The first structure called precursor complex is prepared for the second and third problem: how about the binding ability and cleavage ability between the Csy4 and miniToe structure (hairpin)? The precursor complex consists of two part: wild-type Csy4 and miniToe structure that before cleaved. It describes the structure in the period the after the Csy4 binding to hairpin but didn’t cleave the hairpin in the special site. The Csy4 structure is coming from the X-ray structure we mentioned before while the miniToe structure is constructed totally by the rational model: we put the sequence into the mFold[9] to generate the secondary structure of hairpin then the tertiary structure is produced by the RNAComposer[10]. The molecular docking between Csy4 and miniToe structure is carried out by PatchDock[11]. The precursor complex of wild-type Csy4 can seem in the Fig.1-5

Fig.1-5. The precursor complex of wild-type Csy4

After getting the precursor complex, we begin to prepare for the simulation. Missing hydrogen atoms were added by PDBFixer[12], Forcefield we used is amber ff98SB[13]. The system is immersed in a rectangular TIP3P water box. After minimizing the energy of Protein/RNA system, we give some restriction to the RNA chain to make sure that the structure will not become an unreasonable structure when the temperature rises. All the reaction runs with the PBC condition under 300K and 1 atm in the NPT. The time step is 2 fs while the total simulation time is 50 ns. Gromacs[14] and OpenMM[15] are the most common software we used in Ubuntu 16.04 and Window10. The equipment we used is Intel i7 6700HQ with the NVIDIA GTX 960M 4G, it can simulate about 100-150 ns per day under the GPU acceleration. And the trajectories result is analyzed by Pymol[16] and MDAnalysis[17].

For the second problem, what we can get from the simulation data is protein binding free energy to describe the ability of binding. We use the data in 30-50 ns to calculate it, the 20 ns data in the beginning is aborted to make sure the structure is smooth when calculating. The result of binding free energy for wild-type Csy4 is

△ G_{b i n d i n g} = - 59154.9251 \overset{}{} k j / mol

.

For the third problem, what we can get from the simulation data is some significant distance of key interaction in the active site of Csy4 to describe the ability of cleavage. Jiří Šponer[8] points out some important key interactions of the active site including Ser148(OG)-G20(O2’)、Ser150(OG)-G20(O3’)、Ser151(OG)-G20(N2’). By exploring Jiří Šponer’s work, we finally choose the Ser151(OG)-G20(N2’) as our mathematical form in the third problem. The distance curve of Ser151(OG)-G20(N2’) for wild-type Csy4 can be seen in Fig.1-6. We get the similar result comparing to Jiří Šponer’s workp[8].

Fig.1-6. The distance of Ser151(OG)-G20(N2’) in wild-type Csy4

The second structure called product structure is prepared for the fourth problem: does crRNA release from the RBS. The product complex consists of two part: wild-type Csy4 and miniToe structure that after cleaved. It describes the structure in the period the after the Csy4 binding and cleaving the hairpin in the special site. The Csy4 structure is coming from the X-ray structure we mentioned before while the miniToe structure which is cleaved constructed totally by the rational model: we put two RNA sequence into the SimRNAweb[18] to finish the molecular docking of two chains RNA and generate the tertiary structure. And the molecular docking between Csy4 and miniToe structure is carried out by PatchDock[11]. The product complex of wild-type Csy4 can seem in the Fig.1-7

Fig.1-7 The product complex of wild-type Csy4

We also explore the product complex by molecular dynamics follow the protocol mentioned before, but this time we only set the restriction to RBS chain while the crRNA chain is free in moving.

For the fourth problem, what we can get from the simulation data is the RMSD describing the structure movement for the crRNA chain to be the mathematical form. We can see the RMSD in Fig.1-8. The RMSD is unstable which give an explanation to experiment that crRNA is release from RBS.

Fig.1-8 The RMSD of the product complex of wild-type Csy4

1.5 How to compare four mathematical forms between wild-type and mutants

We are going to check the molecular dynamics of the miniToe system with the mutant Csy4 to show you how to compare the four mathematical forms we choose before between the wild-type Csy4 and Csy4 mutants. In the following description, we give an example using the mutant Q104A to show you how to make comparing.

We first use the SwissModel[19] to generate the tertiary structure with the template Csy4 (PDB ID: 4AL5, resolution 2.0 A). And using the molecular dynamics, we can get the four mathematical form showing in the Fig.1-9.

Fig.1-9 The four key problems in mathematical forms for Csy4-Q104A

Now we divided the four curves into two kind of data: the matrix and the numerical value. The interaction matrix and the curve can be regard as matrix because the curve is discrete, and the binding free energy is just an numerical value.

For the matrix we can use Euclidean distance to describe the difference between two matric:

D (p, q_{W T}) = \sqrt{\sum_{i}^{m} \sum_{j}^{n} {(p_{i, j} - q^{W T}_{i, j})}^{2}}

For the free bind ing energy, we used the formula below to calculate the difference between the wild type and mutant[20]:

\ln (K_{d r e l}) = \ln (\frac{K_{d W T}}{K_{d M U T}}) = △ △ G_{b i n d i n g}

According to description above, we define four value used to compare four key problems between mutant and wild-type:

D_{1} (int t e r a c t i o n \overset{}{} m a t r i x)

,

\ln (K_{d r e l})

,

D_{3} (S e r 151 - G 20 \overset{}{} c u r v e)

,

D_{4} (R M S D)

.

For the mutant Q104A, the four is showing in the following chart

Csy4-Mutant	$D_{1}$	$\ln (K_{d r e l})$	$D_{3}$	$D_{4}$
Q104A	0.483	2483	9.48	30.82

1.6 Conclusion and Result

From the discussion in 1.4 and 1.5, we have the tools to evaluate the four key problems in mathematical form and find the method to compare between mutant and wild-type. The next step is to find out the mutant which fit our needs.

In order to save the computer resource, we only choose the top 10 mutants in D₁ to explore the molecular dynamics. And we mainly focus on the D₃ which describe the ability of cleavage while the D₄ and ln(K_rel) is an alternative value being considered. And for the D₃ we need limitation in the right border to make sure we will not choose the mutant which is totally inactive in cleavage. Luckily, we find that the Csy4-H29A is an inactive protein which has been proved[21]. So we choose 13.41, which is the D₃ value for Csy4-H29A as the right border.
According to all the things we have discussed, the five mutants were chosen in the following table:

Csy4-Mutant	$D_{1}$	$\ln (K_{d r e l})$	$D_{3}$	$D_{4}$
WT	0	0	0	0
Q104A	0.483	2483	9.48	30.82
Y176F	0.592	-382	11.61	40.62
F155A	0.233	-1627	13.41	35.71
H29A	0.173	833	15.29	316.22

2. Hairpin Mutation

2.1 The Large Mutant Library

Starting with the hairpin we want to design. We know that the hairpin used in miniToe structure is coming from the Repeat Area in CRISPR type I-F system, it can be recognized and cleaved by Csy4, and it also called DR in CRISPR system. The Fig.2-1 shows a stander CRISPR array. The yellow line between two DR is guide RNA（gRNA）which is caught from foreign DNA[22].

Fig.2-1 the structure of CRISPR Array

When design the hairpin mutations, it is impossible to use the MD method to find out which hairpin can fit our need because the mutation library is to large. In our miniTioe structure in Fig.2-2, except for the two important nucleic acids----G20 and C21, the mutant library is about. If we still use molecular dynamics to explore all the hairpin mutations, too much time and computer resource is wasted.

Fig.2-2 The structure of miniToe

2.2 Pre-processing Algorithm

Combining the bioinformatics and machine learning, we present an algorithm to pre-processing our big mutation library. Fig.2-3 is the flow chart of the pre-processing algorithm.

Fig.2-3 The flow chart of this pre-processing algorithm

The first step of pre-processing algorithm is to find all the repeat area in the genome as training input. We download the genome as much as possible from NCBI. With the help of PILER-CR[23], which is an algorithm used to find the Repeat Area in bioinformatics, we can get the Repeat Area and CROSPR array from genome quickly. And we only focus on the Repeat Area whose length are 28 and 29 bp because there are some research showing that the length of Repeat Area in CRISPR type I-F system are 28 and 29 bp[24]. In the first step, we download about 5000 genomes and find out 119 Repeat Areas which are 28 and 29 bp. The Fig.2-4 shows the 48 Repeat Areas whose length is 28 bp we find.

Fig.2-4 The Repeat Area which is 28bp

The second step of pre-processing algorithm is to scoring the hairpin we get in the first step. We create a score called DR-Score to evaluate the quality of Repeat area comparing to the wild-type hairpin. The calculation method is below:

{DR-Score}^{i}_{P A} {=E}_{H} • \sum_{n} \frac{D R_{size}}{D R_{size} + N_{mismatcn}}

This formula can be divided into two parts. The first part is E_H , regularized Levenshetein Distance value between evaluated hairpin and wild-type hairpin, which is an index to describe sequence similarity in bioinformatics. The second part is possibility value,

\sum_{n} \frac{D R_{size}}{D R_{size} + N_{mismatcn}}

which describe the quality of Repeat area in the genome it belongs on its own. DR_size refers to the length of Repeat Area, which is 28 or 29 bp in this case. The N_mismatch refers to the number of mutated nucleic acids in Repeat Area in the CRISPR array comparing to the common one. The n refers to the times that the common Repeat Areas apereas in the CRISPR array.

Now take the CRISPR loci in Fig2-5 as an example to calculated the possibility value.

Fig.2-5 one CRISPR array comes from CRISPR database

The CRISPR array shows in Fig.2-5 has 9 Repeat Areas which contains 3 kinds:

DR1: ‘ACTGTACCATGCCTTACTTTGGATTCAAGGCAAAAC’

DR2: ‘ACTGTACCATGCCTGATTTTGGATTCGAGGCAAAAC’

DR3: ‘ACTGTACCATGCCTTACTTTGGATTCAAGTAAATCG’

The first Repeat area is most common while the rest has some mutation. The 3 possibility values are listed following:

D R 1 - P o s s i b i l i t y V a l u e = \sum_{n} \frac{D R_{size}}{D R_{size} + N_{mismatcn}} = \sum_{9} \frac{36}{36 + 0} = 9

D R 2 - P o s s i b i l i t y V a l u e = \sum_{n} \frac{D R_{size}}{D R_{size} + N_{mismatcn}} = \sum_{9} \frac{36}{36 + 3} \approx 8.31

D R 3 - P o s s i b i l i t y V a l u e = \sum_{n} \frac{D R_{size}}{D R_{size} + N_{mismatcn}} = \sum_{9} \frac{36}{36 + 5} \approx 7.90

It is quite easy to understand why the DR-Score will divide into two part. We need the Levenshtein Distance to describe the sequence similarity because not just the length of Repeat Area in CRISPR type I-F system is 28 and 29 bp[26], we need it to distinguish the hairpin which is totally faulted. And we also need the possibility value to describe the quality of hairpin we get, and it is constructed by the following assumption: the more you occur, the better you are.

And If the same Repeat Areas occur in j different species，we use the weight means to calculate the DR-Score. Just assume that there are n_i repeats in the in the i th species, the DR-Score can be calculated by the formula below:

D R - S c o r e = \sum_{i = 1}^{j} (\frac{n_{i}}{\sum_{m}^{j} n_{m}} D R - S c o r e_{i})

The third step of pre-processing algorithm is training the SVM model. The SVM (support vector machine) is a machine learning algorithm which is used to classification and regression[25]. It can construct the complex relationship between input and output with the help of kernel function. And it has been successfully used in predicting the strength of promoter. In this step, we choose the sequence of Repeat Area as input and the DR-Score as output.

Before we training the SVM model, we should first change the sequence into mathematical representation using the following method.
The original sequence data Seq='seq1seq2...seqn-1seqn' which is coded by 'A', 'G', 'C'and 'T' can be transformed into the matrix x by the formula below:

x = {\begin{cases} [1, 0, 0, 0] \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} s e q_{i} = A \\ [0, 1, 0, 0] \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} s e q_{i} = G \\ [0, 0, 1, 0] \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} s e q_{i} = C \\ [0, 0, 0, 1] \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} \overset{}{} s e q_{i} = T \end{cases}

For example, the ‘AGCTA’ can be transformed into the matrix: [1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,1,0,0,0] .

After processing the sequence data, we divide the data into two part: training data and test data, then try to training the model.

The complex relationship between the Repeats sequence matrix x and its DR-Score y can be mapped by an SVM regression function y=f(x). In order to achieve our goals, the SVM model is constructed followed Vapnik et al[25]:

The performance of the SVM model is evaluated by two Squared correlation coefficient (R²)：

R^{2} = \frac{{(n \sum_{i = 1}^{n} f (x_{i}) y_{i} - \sum_{i = 1}^{n} f (x_{i}) \sum_{i = 1}^{n} y_{i})}^{2}}{(n \sum_{i = 1}^{n} f {(x_{i})}^{2} - {(\sum_{i = 1}^{n} f (x_{i}))}^{2}) (n \sum_{i = 1}^{n} y_{i} - {(\sum_{i = 1}^{n} y_{i})}^{2})}

Where the f(x_i) and y_i are the values of the prediction and of DR-score of the Repeat. Many kernel functions have been tested including polynomial function, the sigmoid function, and Gaussian radial basis function (RBF). The RBF is the best choice of our model finally. Fig.2-6 shows the training result by giving the R² of training data and test data which both have high R².

Fig.2-6 The training result

The last step of pre-processing algorithm is to evaluated the big mutation library using the SVM model we train in the third step before. And by the model we can give a score to the hairpin then

2.3 Re-Check and Result

After we pre-processing the big mutant library, we choose 10 hairpins whose score is high but also in gradient, which can make sure that we can achieve the goal that have different regulation level. Then we re-check the four key problems to
All in all, we finally choose five hairpin mutants in the following chart:

Hairpin-Mutant	$D R - S core$
miniToe1	76.6306
miniToe2	65.6278
miniToe3	66.7160
miniToe4	62.5537
miniToe5	52.9794

3.Comparing to the experiment

After designing the protein mutant and hairpin mutants, the wet lab members test the all the Csy4 mutants and hairpin mutants. The result can see in the Fig.3-1.

Fig.3-1 The experimental result of mutants

And we try to give a comparison between the special value we used before for evaluate the mutant and experimental result to check our model.

For the protein mutants, we give a comparison between D₃ and experimental result. Fig.3-2 is the result.

Fig.3-2 The comparison between model and experiment for protein mutant

As we can see in the Fig.3-2, we can find the inner relationship between D₃ and experiment result: the D₃ value describe the cleavage ability between the wild-type and mutant. The higher D₃ value means that it will have an big weaker cleavage ability than the wild-type Csy4.

For the hairpin mutants, we give a comparison between DR-Score and experimental result. Fig.3-3 is the result.

Fig.3-3 The comparison between model and experiment for hairpin mutant

As we can see in the Fig.3-3 we can also can find the inner relationship between DR-Score and experiment result except for the miniToe 5. It is reasonable because the machine learning is quite sensitive to the data amounts and the R² is not 1 in our training result of SVM model.
After all, our wet lab member test 30 combinations of our Csy4 and hairpin. Fig.3-4 is the heatmap result of it.

Fig.3-4 The heatmap result of 30 combination

4.Future Work

1) The QM/MM molecular dynamics is used to design protein mutant more and more frequently[26], we are going to use it to explore more and more details in Csy4/RNA complex.

2) We are gonging to find more and more Repeat Area to increase the data amount for our pre-processing algorithm.

5.Reference

[1].Du P, Miao C, Lou Q, et al. Engineering Translational Activators with CRISPR-Cas System[J]. Acs Synthetic Biology, 2016, 5(1):74.

[2].Molecular Dynamics: , Lecture Notes in Physics, Volume 258. ISBN 978-3-540-16789-1. Springer-Verlag, 1986

[3].Tang Y, Nilsson L. Molecular dynamics simulations of the complex between human U1A protein and hairpin II of U1 small nuclear RNA and of free RNA in solution[J]. Biophysical Journal, 1999, 77(3):1284-1305.

[4].Reyes C M, Kollman P A. Structure and thermodynamics of RNA-protein binding: using molecular dynamics and free energy analyses to calculate the free energies of binding and conformational change[J]. Journal of Molecular Biology, 2000, 297(5):1145-1158.

[5].Jr M A, Nilsson L. Molecular dynamics simulations of nucleic acid-protein complexes.[J]. Curr Opin Struct Biol, 2008, 18(2):194-199.

[6].Karplus M. Molecular Dynamics Simulations of Biomolecules[J]. Acc Chem Res, 2002, 9(9):321-323.

[7].Agostini F, Zanzoni A, Klus P, et al. catRAPID omics: a web server for large-scale prediction of protein–RNA interactions[J]. Bioinformatics, 2013, 29(22):2928-2930.

[8].Estarellas C, Otyepka M, Koča J, et al. Molecular dynamic simulations of protein/RNA complexes: CRISPR/Csy4 endoribonuclease.[J]. Biochimica Et Biophysica Acta, 2015, 1850(5):1072-1090.

[9].M. Zuker. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 31 (13), 3406-3415, 2003.

[10].Popenda M, Szachniuk M, Antczak M, et al. Automated 3D structure composition for large RNAs[J]. Nucleic Acids Research, 2012, 40(14):e112-e112.

[11].Schneidmanduhovny D, Inbar Y, Nussinov R, et al. PatchDock and SymmDock: servers for rigid and symmetric docking[J]. Nucleic Acids Research, 2005, 33(Web Server issue):363-7.

[12].https://github.com/pandegroup/pdbfixe

[13].Samudravijaya K. Comparison of multiple Amber force fields and development of improved protein backbone parameters[J]. Proteins-structure Function & Bioinformatics, 2010, 65(3):712-725.

[14].Pronk S, Schulz R, Larsson P, et al. GROMACS 4.5[J]. Bioinformatics, 2013, 29(7):845-854.

[15].Eastman P, Pande V S. OpenMM: A Hardware Independent Framework for Molecular Simulations[M]. IEEE Educational Activities Department, 2010.

[16].Shringi R P. PyMol Software for 3D Visualization of Aligned Molecules[J]. Biomaterials, 2005, 26(1):63-72.

[17].Michaudagrawal N, Denning E J, Woolf T B, et al. MDAnalysis: a toolkit for the analysis of molecular dynamics simulations.[J]. Journal of Computational Chemistry, 2011, 32(10):2319-2327.

[18].Boniecki M J, Lach G, Dawson W K, et al. SimRNA: a coarse-grained method for RNA folding simulations and 3D structure prediction[J]. Nucleic Acids Research, 2016, 44(7):e63-e63.

[19].Schwede T, Kopp J, Guex N, et al. SWISS-MODEL: An automated protein homology-modeling server.[J]. Nucleic Acids Research, 2003, 31(13):3381-3385.

[20].Estimation of Relative Protein–RNA Binding Strengths from Fluctuations in the Bound State

[21].Lee H Y, Haurwitz R E, Apffel A, et al. RNA-protein analysis using a conditional CRISPR nuclease[J]. Proceedings of the National Academy of Sciences of the United States of America, 2013, 110(14):5416-5421.

[22].Przybilski R, Richter C, Gristwood T, et al. Csy4 is responsible for CRISPR RNA processing in Pectobacterium atrosepticum.[J]. Rna Biology, 2011, 8(3):517-528.

[23].Edgar R C. PILER-CR: Fast and accurate identification of CRISPR repeats[J]. Bmc Bioinformatics, 2007, 8(1):1-6.

[24].Crawley A B, Henriksen J R, Barrangou R. CRISPRdisco: An Automated Pipeline for the Discovery and Analysis of CRISPR-Cas Systems[J]. 2018, 1(2).

[25].Schölkopf B, Tsuda K, Vert J P. Support Vector Machine Applications in Computational Biology[C]// MIT Press, 2004:71-92.

[26].Wang X, Li R, Cui W, et al. QM/MM free energy Simulations of an efficient Gluten Hydrolase (Kuma030) Implicate for a Reactant-State Based Protein-Design Strategy for General Acid/Base Catalysis[J]. Scientific Reports, 2018, 8.