分子動力学でRNAのフォールディングを追跡しようとして失敗した記録
RNAの二次構造予測と言えば、RNAfoldとmfoldがバイオインフォマィテクス業界の鉄板だと思うけど、未だに、こいつらのアルゴリズムを理解してないわたしは、こいつら本当に正しいのか信用できないのだった。というわけで、分子動力学で、RNAの二次構造を見てみる。まぁ、分子動力学なら信用できるかというと、それは分からないが。練習問題として、適当な(pre-)miRNAを使う。MDとRNAfold/mfoldの結果が一致するなら、綺麗なステム・ループ構造が見られるはず。と言いたいけど、実際にフォールディングするまでに、どれくらいの時間かかるのか分からないのが問題。タンパク質の場合ミリ秒〜秒のオーダーでフォールディングすると言われるので、そこまでの計算は、そこらのマシンじゃほぼ無理
具体的なモノとしては、mir-206というのを選んだ。Pax3やPax7の3'-UTR領域をターゲットにして、筋細胞の増殖・分化に関わるらしいけど、今は素性は、特に意味はない
http://www.mirbase.org/cgi-bin/mirna_entry.pl?acc=MI0000490
最初は、pre-miRNA全長(90bp弱)でやろうと思ったら、うちの貧弱なマシンでは、存外計算時間がかかることがわかったので(1nsをやるのに推定一週間以上)、適当に両端をカットして、計算することにした。使用した配列は、AUCCCCAUAUGGAUUACUUUGCUAUGGAAUG。RNAFoldによる予測はこんなん
まず考えないといけないこととして、RNAを構成する原子の初期配置を決める必要がある。タンパクなら、結晶構造がある場合、適当にpdb取ってくればいいのだけど、RNAの場合は、塩基配列から初期位置を適当に決める必要がある。幸い、Ambertoolsのtleap(もしくは、xleap)というのを使えば、それは簡単にできる。Amberは無償では使えない(と思う)けど、Ambertoolsは誰でも、無償で使える。
AmbertoolsはAmberのHPで配布されている。現時点では、Ambertools13が最新らしい
http://ambermd.org/#AmberTools
Ambertoolsをビルドしたら、こんな感じ↓のテキストを書いて、genpdb.inという名前ででも保存してやる
source leaprc.ff99SB
mir206=sequence{RA5 RU RC RC RC RC RA RU RA RU RG RG RA RU RU RA RC RU RU RU RG RC RU RA RU RG RG RA RA RU RG3}
savepdb mir206 mir206.pdb
quitそれで、
tleap -f genpdb.in
とかやって、mir206.pdbができればOK。多分、直線状に塩基が並んでるだけのRNAができてる。pdbが出来れば、あとはGromacs tutorialの
Lysozyme in Water
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/
と同じようにやればいいはず(完)
まず、pdb2gmxでtopファイルとか作って、水を入れる。tleapで力場をamber99sbと指定しているので、合わせる
pdb2gmx -ignh -f mir206.pdb -o mir206_unbox.gro -p mir206.top -ff amber99sb -water spce #editconf -bt triclinic -d 0.85 -f mir206_unbox.gro -o mir206_newbox.gro editconf -box 14 10 8 -f mir206_unbox.gro -o mir206_newbox.gro genbox -cp mir206_newbox.gro -cs spc216.gro -o mir206_sol.gro -p mir206.top
適当にtriclinicでbox作ったら、RNAが一直線に並んでるせいで、やたら平べったいボックスができた(14x10x4くらいだった気がする)。何かバランスが悪い気がして、cubicにすると、無駄に水が増えるので、適当にボックスサイズを指定した。あと、系全体を中性にするためにイオンを適当に入れる必要がある。今は、Na+イオンを追加する。これは、一回.tpr作ってgenionコマンドでできる
grompp -p mir206.top -c mir206_sol.gro -f em.mdp -o mir206_sol.tpr echo 3 | genion -neutral -s mir206_sol.tpr -o mir206_all.pdb -pname NA -nname CL -p mir206.top
em.mdpは後で使うエネルギー最小化用のmdpファイルで、多分別に何でもいい。"echo 3 "は今の場合、水の一部をイオンに置き換えることを意味する
ここまで出来たら、あとは適当に.mdpファイルを用意して、以下の計算をやる
(1)エネルギー最小化(em.mdp)
(2)NVT平衡化/加熱(nvt.mdp)
(3)NPT平衡化(npt.mdp)
(4)プロダクションラン(md.mdp)
という4ステップ。mdpファイルは、Lysozyme in Waterに置いてあるmdpファイルを微修正して使った
#Energy minimization grompp -f em.mdp -c mir206_all.pdb -p mir206.top -o mir206_em.tpr mdrun -v -deffnm mir206_em #Heating grompp -f nvt.mdp -c mir206_em.gro -p mir206.top -o mir206_nvt.tpr mdrun -v -deffnm mir206_nvt #Equilibration grompp -f npt.mdp -c mir206_nvt.gro -t mir206_nvt.cpt -p mir206.top -o mir206_npt.tpr mdrun -v -deffnm mir206_npt #Production MD grompp -f md.mdp -c mir206_npt.gro -t mir206_npt.cpt -p mir206.top -o mir206_md.tpr mdrun -v -deffnm mir206_md
で、計算したので、解析しないといけない。とりあえず、見た目で変な計算結果になってないことを確認するのにVMDで眺める。
Download VMD:
http://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD
特に、何もしなくても、Gromacsの.groフォーマット読んで、.trrファイル読んでアニメーションしたりできる。水分子が邪魔なので、Graphics>Representationsメニューの真ん中あたりで、Selections選んで、Selected Atomsに"nucleic"とか入れてApplyボタン押すと、水とイオンが消える。まぁ、時間が足りないのか、全然フォールディングしてくれてない。
一応、水素結合の解析は、g_hbondコマンドでできるのだけど、今の場合は、見るからにだめっぽい。g_hbondの出力は、極めて不親切なので、hbmap.plとかいうのを使うと、分かりやすい出力に整形してくれる
echo 1 1 | g_hbond -noda -f mir206_md.xtc -s mir206_md.tpr -num hbnum.xvg -dist hbdist.xvg -hbn hbond.ndx -hbm hbmap.xpm perl hbmap.pl -s mir206_all.pdb -map hbmap.xpm -index hbond.ndx
というわけで、全然フォールディングしてくれないのだけど、計算時間が足らないのだとしても、1nsの計算でも半日くらいかかってるので、多少時間を増やしても、どうにもならない気がする。とりあえず、もうちょっとコンパクトな初期構造が欲しいと思って真空中で計算して見ることにした。系を中性にするために、イオンだけは入れておきたいけど、Gromacs付属のgenionは既存の原子をイオンに置き換えるツールなので、tleapでイオンを入れた。こんなん↓
source leaprc.ff99SB
mir206=sequence{RA5 RU RC RC RC RC RA RU RA RU RG RG RA RU RU RA RC RU RU RU RG RC RU RA RU RG RG RA RA RU RG3}
addions mir206 Na+ 0
savepdb mir206 mir206_amber.pdb
quittleapでイオンを追加した場合、Gromacsとイオンのresidue name/atom nameが違うので、適当に書き換える必要がある。あとは大体一緒。
tleap -f genpdb.in sed -e 's/Na+ Na+/Na NA/g' mir206_amber.pdb > mir206_gromacs.pdb pdb2gmx -ignh -f mir206_gromacs.pdb -o mir206_unbox.gro -p mir206_vac.top -ff amber99sb -water spce editconf -bt cubic -d 0.85 -f mir206_unbox.gro -o mir206_vac.pdb #Energy minimization grompp -p mir206_vac.top -c mir206_vac.pdb -f em.mdp -o mir206_vac_em.tpr mdrun -v -deffnm mir206_vac_em #Heating grompp -f nvt.mdp -c mir206_vac_em.gro -p mir206_vac.top -o mir206_vac_nvt.tpr mdrun -v -deffnm mir206_vac_nvt #Equilibration grompp -f npt.mdp -c mir206_vac_nvt.gro -t mir206_vac_nvt.cpt -p mir206_vac.top -o mir206_vac_npt.tpr mdrun -v -deffnm mir206_vac_npt
平衡化後の構造は、微妙にそれっぽい
この構造から始めて水を入れて再計算していくといいのかもしれない。一般的によい方法か分からないけど、とりあえず、tleapで生成されるよりコンパクトな初期構造ができれば、必要な水分子が減って計算時間が随分減らせる。tleapが適当にコンパクトな初期構造を生成してくれるようになれば、それが一番だと思う
#Solvation genbox -cp mir206_vac_md.gro -cs spc216.gro -o mir206_vac_sol.pdb -p mir206_vac.top #Energy minimization grompp -p mir206_vac.top -c mir206_vac_sol.pdb -f em.mdp -o mir206_em2.tpr mdrun -v -deffnm mir206_em2 #Heating grompp -f nvt.mdp -c mir206_em2.gro -p mir206_vac.top -o mir206_nvt2.tpr mdrun -v -deffnm mir206_nvt2 #Equilibration grompp -f npt.mdp -c mir206_nvt2.gro -t mir206_nvt2.cpt -p mir206_vac.top -o mir206_npt2.tpr mdrun -v -deffnm mir206_npt2 #Production MD grompp -f md.mdp -c mir206_npt2.gro -t mir206_npt2.cpt -p mir206_vac.top -o mir206_md2.tpr mdrun -v -deffnm mir206_md2
が、計算結果を見ても、全然それらしい水素結合が形成されてる気配はない。
とりあえず、ここまで殆ど何も調べず適当にやったのだけど
Why Can’t We Predict RNA Structure At Atomic Resolution?
http://daslab.stanford.edu/pdf/SripakdeevongBeauchampDas_Chapter_OnlinePDF.pdf
という文書を読むと"No existing algorithm can start with arbitrary RNA sequences and return the precise three-dimensional structures that ensure their biological function"と書いてあって、わたしのやろうとしたことは、現状無理っぽい。
けどまぁ、"In many respects, it now appears that RNA is easier to fold than other biopolymers. For example, unlike proteins, which typically require at least a dozen residues to form well-defined structures, the simplest RNAs with welldefined,recurrent structures are as small as eight residues (Jucker et al. 1996)."とも書いてあって、RNAならタンパクより簡単にフォールディング見れるんじゃねと思ったのは、的外れでもないらしい。ただ、上の計算で使ったRNAは30bpあるし、元のpre-miRNAは90bp近い。これらは、まだでかすぎるのかもしれない。というか、もっと簡単なpolyA/polyT混合した系とかから始めるべきだったかもしれない。。
1996年のJuckerという人たちの論文は
A network of heterogeneous hydrogen bonds in GNRA tetraloops.
http://www.ncbi.nlm.nih.gov/pubmed/9000624
らしい。2000年以降の論文を見ると、10bp前後の長さで計算してることが多いっぽいので、pre-miRNAの計算をやるのは、まだ難易度が高い
一応、使用したem.mdpとmd.mdp。nvt.mdpとnpt.mdpは、md.mdpから大体分かる
em.mdp
integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep = 0.002 ; Energy step size nsteps = 5000 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = grid ; Method to determine neighbor list (simple, grid) rlist = 2.0 ; Cut-off for making neighbor list (short range forces) rcoulomb = 2.0 ; Short-range electrostatic cut-off rvdw = 2.0 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions (yes/no)
こっちはmd.mdp
integrator = md ; leap-frog integrator nsteps = 500000 ; 2 * 500000 = 1000 ps, 1 ns dt = 0.002 ; 10 fs ; Output control nstxout = 1000 ; save coordinates every 0.2 ps nstvout = 1000 ; save velocities every 0.2 ps nstenergy = 1000 ; save energies every 0.2 ps nstlog = 1000 ; update log file every 0.2 ps ; Bond parameters continuation = yes ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs rlist = 1.0 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) rvdw = 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = RNA Water_and_ions tau_t = 0.1 0.1 ; time constant, in ps ref_t = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = no ; Velocity generation is off
