1 MiChroM - GROMACS Tutorial
This tutorial helps the user to run GROMACS Molecular Dynamics software using the MiChroM (Minimal Chromatin Model) parameters.
You can get the already files for MiChroM Gromacs simulation in MEGABASE page, or if you have your own chromatin sequence you can generate the input files in Gromacs Input File page, two input files are required. One of them should be the Chromatin sequence types, and the other should be the loops contact in Gromacs Input File page, two input files are required. One of them should be the Chromatin sequence types, and the other should be the loops contact.
For the chromatin sequence types, the file should have the format of 2 columns. The first column is the bead position in sequence and the second is the chromatin type. The types should follow the syntax A1, A2, B1, B2, B3, B4 and NA. In the Gromacs top file the chromatin types are named as ZA, OA, FB, SB, TB, LB and UN, respectively. The chromatin type group is presented as ChrA for types A1 and A2 , ChrB for types B1 , B2 , B3 and B4 and ChrU for the unknown type UN . A example of the chromatin type sequence file for chromosome 21 from GM12878 cell line can be found here.
For the loop pair interactions, the file also should have the format of 2 columns. The first column is the index of the bead i which interacts with the bead j described in second column. The loop interactions are described in \verb|[pairs]| section of Gromacs topfile. A example of the chromatin loops file for chromosome 21 from GM12878 cell line can be found here.
1.2 Files Description
It is expected nine different files if the code ran correctly, in which, five are table files *.xvg , two *.gro files and two *.top files. Examples of each output file also can be found in MEGABASE page, just pick a cell line.
For table files details please see the section (1.3 - Tables Details)
. Below are the descriptions of the *.gro and *.top files:
1.3 Tables Details
The MiChroM model uses some particular interactions which make necessary the usage tabulated potentials in Gromacs. For the details of the potentials, see the MiChroM page. To run the molecular dynamics simulation in Gromacs it will be necessary five different tables. Below are presented the name and the description of each one:
1.4 Running Gromacs
Here will be described the basic commands to allow the user to perform a chromatin molecular dynamics simulation using Gromacs. Click on the links to get the MDP files for the initial condition simulation MDP-init.mdp and for the state sampling simulation MDP-sim.mdp.
The first command creates an index file *.ndx which will be necessary forward. Here is the example for chromosome 21:
gmx make_ndx -f C21-init.gro -o C21-index
Chose the option 0 and then type q followed by Enter to quit.
1.4.1 Initial condition
The initial condition simulation starts the chromatin in an initial structure and collapses it inside the cell nucleu wall.
The grompp command is the Gromacs preprocessor which will check if everything is right to run the simulation and will create the *.tpr file.
gmx grompp -f MDP-init.mdp -n C21-index.ndx -p C21-init.top -c C21-init.gro
This command will generate a lot of warnings highlighted that the simulation has a "Duplicate atom index in angle_restraints". It is okay, the way of the implementation of the angle potential restraint needs that.
To run the simulation use the command:
-r C21-restraint.gro -o run-init.tpr -maxwarn -1
gmx mdrun -v -table table.xvg -tableb table_b0.xvg table_b1.xvg
-s run-init.tpr -c out.gro 1.4.2 Simulation
Once the file out.gro was generated correctly, the chromatin should be inside the cell nucleus wall. The simulation will assume the out.gro file as the initial structure to properly start the molecular dynamics simulation to sampling the chromatin states. Here it will be used the MDP-sim.mdp file, pay attention to the parameter of the simulation.
gmx grompp -f MDP-sim.mdp -n C21-index.ndx -p C21-sim.top
To run the long simulations use the command:
-c out.gro -o run-sim.tpr -maxwarn -1
gmx mdrun -v -table table.xvg -tablep tablep.xvg -tableb table_b0.xvg
table_b1.xvg table_b2.xvg -s run-sim.tpr -deffnm C211.5 How to analyze
Once the molecular dynamic simulation is finished, there is a large group of different analysis which could be done. Here it will be explained the steps to obtain the contact map analyzing the simulation trajectory .trr.
gmx trjconv -f C21.trr -s run-sim.tpr -o all-C21.gro
The output file is a *.gro file which contains the position of each in different simulation frames. To generate the contact map probability, the user can run the python code GAMinG-plot GAMinG plot described in section (1.6 - GAMinG-plot)
.
1.6 GAMinG-plot
The code GAMinG-plot.py reads the *.gro file a generates the contact map plot.
To run the code use the commands:
./GAMinG-plot.py or python2.7 GAMinG-plot.py. Use the flag -h or --help to get the description of the code options.
It is presented below an example of the usage of the code for chromosome 21.
./GAMinG-plot.py -f all-C21.gro -n C21 -plot yes
The code generates an output file *.dat: