PCA analysis on MD trajectory through cpptraj

-


PCA analysis through CPPTRAJ module of AmberTools:


>source $AMBERHOME/amber.sh

>$AMBERHOME/bin/cpptraj

>parm comp_solv.prmtop

>trajin comp_150ns.nc

>center

>autoimage

>rms first @CA mass

>average crdset avg_set

>createcrd traj_set

>run

>crdaction traj_set rms ref avg_set @CA

>crdaction traj_set matrix covar :116-229@CA name pro_covar out covmat.dat #writes covariance matrix based on the selection of residues

>runanalysis diagmatrix pro_covar out evecs.dat vecs 100 name myevecs nmwiz nmwizvecs 100 nmwizfile porcupine.nmd nmwizmask :116-229@CA #writes the information about 100 principal components (first 100 PCs) for that particular selection,  produces .nmd file that can be used to generate porcupine plots in normal mode wizard plugin of VMD.

>runanalysis modes eigenval name myevecs out evalfrac.dat gives output about the eigen values, fraction of contribution, cumulative contributions of  eigen vectors(or PCs) (in evalfrac file),

>runanalysis modes fluct out rmsfluct1-10.dat name myevecs beg 1 end 10 #to write the fluctuations of first 10 PCs

>runanalysis modes fluct out rmsfluct1-2.dat name myevecs beg 1 end 2 #to write the fluctuations of first two PCs

>runanalysis modes fluct out rmsfluct-mode1.dat name myevecs beg 1 #to write the fluctuations of first PC

>crdaction traj_set projection p1 modes myevecs beg 1 end 20 :116-229@CA out pca.dat #to write the projections of first 20 PCs

>crdaction traj_set projection p2 modes myevecs beg 1 end 2 :116-229@CA out pca12.dat #to write the projections of first 2 PCs

The PCA12.dat can be opened in xmgrace to see the projections of both PCs on x and y axis. The lines can be replaced with symbols such as circles to get the appearance of a scatter plot.

To generate the Free energy landscape (FEL) plots, u need to convert the PCA12.dat file i.e. the projection file to gibbs.txt. For this u need gromacs installed in your system and then follow the following steps.

gmx sham -f 2d_proj.xvg -ls gibbs.xpm -notime #for free energy, "sham" module uses PC1 and PC2 file

gmx xpm2ps -f gibbs.xpm  -o gibbs.eps -rainbow red

python xpm2txt.py -f gibbs.xpm -o gibbs.txt 

The link for xpm2txt.py and scripts to generate the images of different FEL plots are mentioned in our previous blogs: xpm2txt.py , Python script to plot 2D FEL, Python script to plot 3D FEL plots




For any queries regarding the script, kindly drop your message in the comments section below.

Please also provide your valuable feedback and stay tuned for more such blogs and scripts!!

Happy Computing!!

Comments

  1. AnonymousJanuary 11, 2024 at 5:33 AM

    I could not understand how we came to 2d_proj.xvg while using cpptraj we got PCA12.dat? xvg is only possible with gromacs?
    How can we now use the PCA12.dat to get FEL?

    ReplyDelete

Post a Comment

Popular posts from this blog

Calculation of RMSD values between two ligand poses using web servers/tools

-
#Web servers for calculation of RMSD values between a ligand poses.  Please note: One pose should be a reference pose and another should be a query pose. PLDbench:  The RMSD value between two molecules can be calculated using PLDbench. There are four modules in this webserver. Single: This module's functionality is limited to calculating RMSD for the 57 complexes utilized in the benchmarking study. Users must submit the docked pose of a ligand that corresponds to one of the 57 PDB-IDs offered in the drop down option in this module. It accepts three file formats: '.pdb', '.mol2', and '.sdf'. Users can choose the type of RMSD they want to compute. Standard heavy-atom RMSD, Hungarian (symmetry-corrected) heavy-atom RMSD, and Minimum-distance heavy-atom RMSD can all be calculated using this module. As a result, the Hungarian (symmetry-corrected) heavy-atom RMSD value is returned by default.  Link for the server:https://webs.iiitd.edu.in/raghava/pldbench/single.p...
Read more

AMBER TUTORIAL-2: How to simulate a protein- ligand system: basic steps for md simulation in amber (AMBER 16)

-
BLOG #3  Hi Guys,           This is a new tutorial on how to simulate a protein-ligand system: basic steps using AMBER. Before we begin, I hope you have read our previous post on AMBER TUTORIAL-1, since I will only be providing additional information with the ligand here, rather than repeating the previous instructions. #To parameterize the ligand file, we start with Antechamber and then parmchck.  These are a series of tools for creating files for organic molecules and some protein metal centers, which can then be read into LEaP. Junmei Wang created the Antechamber suite, which is intended to be used in conjunction with the general AMBER force field (GAFF) (gaff.dat). This is the package's most crucial program. It can convert a variety of files and assign atomic charges and types to them as well. Sqm (or, alternatively, mopac or divcon), atomtype, am1bcc, bondtype, espgen, respgen, and prepgen are among the programs that antechamber runs in respo...
Read more

Python script to plot 3D FEL plots

-
Image
  # To plot FEL in 3D use this python script and run as python2.7 fel3d.py. Here the file uses gibbs.txt. therefore before running the script make sure u have the gibbs.txt in the same folder of this script or else provide the location of the script.  I have attached some images from my research article using the same script for 3D FEL plot These plots are published in  Interplay among Structural Stability, Plasticity, and Energetics Determined by Conformational Attuning of Flexible Loops in PD-1 The 2D FEL scripts are also mentioned in the next blog:  Python script to plot 2D FEL Another point to consider is: when u r plotting FEL plots for multiple systems, make sure your color bar/gradient should be same in all the plots. So u can  change the values in vmax,  ax.set_zlim(0, 18),  ticks=range(0, 21, 2)) in the script below.  Like in my case the highest value was 18 in gibbs free energy (the z axis) in all the axis. therefore i mentioned 18 there...
Read more