Connect with us

Bioinformatics Programming

Free_Energy_Landscape-MD: Python package to create Free Energy Landscape using PCA from GROMACS.

Dr. Muniba Faiza

Published

on

In molecular dynamics (MD) simulations, a free energy landscape (FEL) serves as a crucial tool for understanding the behavior of molecules and biomolecules over time. It is difficult to understand and plot a meaningful FEL and then extract the time frames at which the plot shows minima. In this article, we introduce a new Python package (Free_Energy_Landscape-MD) to generate an FEL based on principal component analysis (PCA) from MD simulation done by GROMACS [1].

Free_Energy_Landscape-MD Python Package

This package provides a 3D Free Energy Landscape (FEL) plot showing the minima points. It utilizes Principal Component Analysis (PCA) of the Molecular Dynamics (MD) trajectory to generate the FEL. Additionally, it offers specific time frames corresponding to the minima regions in the FEL. The package includes two main files:

  • FEL_Time_Frame_Extractor_CaseI.ipynb,
  • FEL_Time_Frame_Extractor_CaseII.ipynb,

along with other required input files such as the MD trajectory file and example output files. Before using the script to generate the FEL, PCA must be performed. The complete process, from performing PCA to extracting the time frames corresponding to the minima regions, is explained in the following sections.

Availability

The package is freely available on GitHub. There, you will find four directories and the scripts. These directories represent the different processes used to generate FEL and extract the time frames. They are organized to avoid confusion. However, you are not required to maintain your files in this specific structure. You can place all your input files in the same directory and run the scripts as you find suitable.

Usage

We will proceed in four fundamental steps:

1. Principal Component Analysis (PCA)

2. Calculation of Gibb’s Free Energy

3. Free Energy Landscape (FEL) generation and extraction of time frames

4. Conformational analysis

We will generate inputs for the FEL script in the first two steps. In the third step, we will generate the FEL and extract the time frames using the appropriate script (explained in step 3). Finally, in the last step, we will take snapshots of the time frames.

1. Principal Component Analysis (PCA)

For PCA, we will calculate the eigenvectors and eigenvalues. Eigenvectors are vectors associated with eigenvalues in linear algebra. In PCA, eigenvectors represent the directions along which the data varies the most. They are orthogonal (perpendicular) to each other and describe the principal axes of variation in the data. These vectors are useful for transforming the original data into a new coordinate system, where the dimensions are aligned with the directions of maximum variance. Using the gmx anaeig module, we will calculate the PC1 and PC2 and then merge them.

To do PCA of the MD trajectory obtained from GROMACS. Open a terminal and follow the steps given below:

1.1. Calculating the eigenvectors and eigenvalues using the gmx covar module.

$ gmx covar -s md_0_1.gro -f md_0_1_noPBC.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm

1.2. Calculating PC1 and PC2.

$ gmx anaeig -f md_0_1_noPBC.xtc -s md_0_1.gro -v eigenvectors.trr -last 1 -proj pc1.xvg

$ gmx anaeig -f md_0_1_noPBC.xtc -s md_0_1.gro -v eigenvectors.trr -first 2 -last 2 -proj pc2.xvg

1.3. Merging PC1 and PC2 in a single file.

$ paste pc1.xvg pc2.xvg | awk '{print $1, $2, $4}' > PC1PC2.xvg

2. Calculation of Gibb’s Free Energy

2.1. We will use the PC1PC2.xvg and gmx sham module to calculate Gibb’s Free Energy.

$ gmx sham -f PC1PC2.xvg -ls FES.xpm

This PC1PC2.xvg file contains time frames and the PC1 and PC2 values. We will use this file to extract time frames corresponding to the minima points.

2.2. Now we must convert the FES.xpm file to .dat format which we will use in free energy landscape (FEL) generation. For this, we will need a Python script, namely, xpm2txt.py. You can download it from here. Remember to use Python2.7 for the following step.

$ python2.7 xpm2txt.py -f FES.xpm -o fel.dat

This fel.dat file contains free energy values and the and coordinates. We will use this file as an input to generate the plot.

3. Free Energy Landscape (FEL) generation and extraction of time frames.

Case-I For long MD simulation without any noise and artifacts.

In this step, we will use the Python script FEL_Time_Frame_Extractor_CaseI.ipynb. You can download it and run it for your MD data. It will generate and save a 3D FEL plot (fel_with_minima) showing minima points along with a text file containing the minima points and their corresponding time frames (minima_time_frames.txt). By default, it is configured to display a maximum of 3 minima points, but you can adjust this setting in the script. After generating the plot, it will use the PC1PC2.xvg file to extract the time frames corresponding to the minima regions. The script searches for the time frames in the PC1PC2.xvg file by matching the minima x values or the minima y values.

NOTE: You might find only one minima point and may not find the time frames in the provided script on GitHub. This is because I ran a very short 1ns MD simulation using GROMACS [1] solely for demonstration purposes.

Case-II What if the script doesn’t provide any timeframe output matching the minima points?

In some cases, it does not show any time frame corresponding to the minima points. It can be due to many reasons including:

  • Firstly, check whether you have provided the correct input files.
  • Your MD simulation is too short. The PC1PC2.xvg file may not capture the full range of conformational states of the system, leading to missing minima points. This could happen if the simulation trajectory is biased or if certain conformations are under-sampled.
  • The resolution of the energy landscape plot and the data in the PC1PC2.xvg file might differ. Minima points that appear distinct in the energy landscape plot may not be resolved in the data from the PC1PC2.xvg file, especially if the sampling frequency or resolution is lower.
  • The system may have inherent complexity or dynamics that are not fully captured by the analysis. Certain states or transitions between states may be difficult to characterize accurately, leading to discrepancies between the energy landscape plot and the PC1PC2.xvg data.
  • Noise, artifacts, or inaccuracies in the data or analysis methods could also contribute to discrepancies. These issues may obscure or distort the identification of minima points in the PC1PC2.xvg data.

To address this issue, it is important to run a long MD for a sufficient amount of time and minimize artifacts during simulation. However, you can still find the minima points in the PC1PC2.xvg file but their precision may vary. They are the closest values to the minima points. For this, another similar script (FEL_Time_Frame_Extractor_CaseII.ipynb) is provided which will fetch you the corresponding time frames.

4. Conformational Analysis

In this last step, we will take snapshots of the protein/system at the extracted time frames from the MD trajectory.

$ gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o snapshot.pdb -dump 520

Visualize this PDB file in any molecular viewer such as Pymol [2]. It will show you the conformation. Repeat this step as needed to extract snapshots from different time frames in your trajectory. Adjust the parameters as necessary to extract the desired frames.

For more details on taking snapshots using GROMACS, read this article.


References

  1. Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., & Lindahl, E. (2015). GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX1, 19-25.
  2. The PyMOL Molecular Graphics System, Version 2.5, Schrödinger, LLC.

How to cite:

Faiza, M., 2024. Free_Energy_Landscape-MD: Python package to create Free Energy Landscape using PCA from GROMACS., 10(2): page 8-12. https://bioinformaticsreview.com/20240228/free_energy_landscape-md-python-package-to-create-free-energy-landscape-using-pca-from-gromacs/

Dr. Muniba is a Bioinformatician based in New Delhi, India. She has completed her PhD in Bioinformatics from South China University of Technology, Guangzhou, China. She has cutting edge knowledge of bioinformatics tools, algorithms, and drug designing. When she is not reading she is found enjoying with the family. Know more about Muniba

Bioinformatics Programming

How to copy and rename files simultaneously in same directory in Ubuntu (Linux)?

Dr. Muniba Faiza

Published

on

How to copy & rename files in same directory in Ubuntu (Linux)?

Copying and renaming files from one destination to another is easily accomplished by the cp command in Linux. However, it becomes more complex when we want to copy and rename all files in the same directory while assigning a serial number to each copy. In this article, we provide a shell script to copy and rename files present in a directory, adding a unique serial number to each copy.

(more…)

Continue Reading

Bioinformatics News

VS_Analysis: A Python package to perform post-virtual screening analysis

Dr. Muniba Faiza

Published

on

VS_Analysis: A Python package to perform post-virtual screening analysis

Virtual screening (VS) is a crucial aspect of bioinformatics. As you may already know, there are various tools available for this purpose, including both paid and freely accessible options such as Autodock Vina. Conducting virtual screening with Autodock Vina requires less effort than analyzing its results. However, the analysis process can be challenging due to the large number of output files generated. To address this, we offer a comprehensive Python package designed to automate the analysis of virtual screening results.

(more…)

Continue Reading

Bioinformatics Programming

vs_interaction_analysis.py: Python script to perform post-virtual screening analysis

Dr. Muniba Faiza

Published

on

vs_interaction_analysis.py: Python script to perform post-virtual screening analysis

Analyzing the results of virtual screening (VS) performed with Autodock Vina [1] can be challenging when done manually. In earlier instances, we supplied two scripts, namely vs_analysis.py [2,3] and vs_analysis_compounds.py [4]. This time, we have developed a new Python script to simplify the analysis of VS results.

(more…)

Continue Reading

Bioinformatics Programming

How to create a pie chart using Python?

Dr. Muniba Faiza

Published

on

How to create a pie chart using Python?

In this article. we are creating a pie chart of the docking score of five different compounds docked with the same protein. (more…)

Continue Reading

Bioinformatics Programming

How to make swarm boxplot?

Dr. Muniba Faiza

Published

on

How to make swarm boxplot?

With the new year, we are going to start with a very simple yet complicated topic (for beginners) in bioinformatics. In this tutorial, we provide a simple code to plot swarm boxplot using matplotlib and seaborn. (more…)

Continue Reading

Bioinformatics Programming

How to obtain ligand structures in PDB format from PDB ligand IDs?

Dr. Muniba Faiza

Published

on

How to obtain ligand structures in PDB format from PDB ligand IDs?

Previously, we provided a similar script to download ligand SMILES from PDB ligand IDs. In this article, we are downloading PDB ligand structures from their corresponding IDs. (more…)

Continue Reading

Bioinformatics Programming

How to obtain SMILES of ligands using PDB ligand IDs?

Dr. Muniba Faiza

Published

on

How to obtain SMILES of ligands using PDB ligand IDs?

Fetching SMILE strings for a given number of SDF files of chemical compounds is not such a trivial task. We can quickly obtain them using RDKit or OpenBabel. But what if you don’t have SDF files of ligands in the first place? All you have is Ligand IDs from PDB. If they are a few then you can think of downloading SDF files manually but still, it seems time-consuming, especially when you have multiple compounds to work with. Therefore, we provide a Python script that will read all Ligand IDs and fetch their SDF files, and will finally convert them into SMILE strings. (more…)

Continue Reading

Bioinformatics Programming

How to get secondary structure of multiple PDB files using DSSP in Python?

Dr. Muniba Faiza

Published

on

How to get secondary structure of multiple PDB files using DSSP in Python?

In this article, we will obtain the secondary structure of multiple PDB files present in a directory using DSSP [1]. You need to have DSSP installed on your system. (more…)

Continue Reading

Bioinformatics Programming

vs_analysis_compound.py: Python script to search for binding affinities based on compound names.

Dr. Muniba Faiza

Published

on

vs_analysis_compound.py: Python script to search for binding affinities based on compound names.

Previously, we have provided the vs_analysis.py script to analyze virtual screening (VS) results obtained from Autodock Vina. In this article, we have provided another script to search for binding affinity associated with a compound. (more…)

Continue Reading

Bioinformatics Programming

How to download files from an FTP server using Python?

Dr. Muniba Faiza

Published

on

How to download files from an FTP server using Python?

In this article, we provide a simple Python script to download files from an FTP server using Python. (more…)

Continue Reading

Bioinformatics Programming

How to convert the PDB file to PSF format?

Dr. Muniba Faiza

Published

on

How to convert the PDB file to PSF format?

VMD allows converting PDB to PSF format but sometimes it gives multiple errors. Therefore, in this article, we are going to convert PDB into PSF format using a different method. (more…)

Continue Reading

Bioinformatics Programming

smitostr.py: Python script to convert SMILES to structures.

Dr. Muniba Faiza

Published

on

smitostr.py: Python script to convert SMILES to structures.

As mentioned in some of our previous articles, RDKit provides a wide range of functions. In this article, we are using RDKit [1] to draw a molecular structure using SMILES. (more…)

Continue Reading

Bioinformatics Programming

How to preprocess data for clustering in MATLAB?

Dr. Muniba Faiza

Published

on

How to preprocess data for clustering in MATLAB?

Data preprocessing is a foremost and essential step in clustering based on machine learning methods. It removes noise and provides better results. In this article, we are going to discuss the steps involved in data preprocessing using MATLAB [1]. (more…)

Continue Reading

Bioinformatics Programming

How to calculate drug-likeness using RDKit?

Dr. Muniba Faiza

Published

on

How to calculate drug-likeness using RDKit?

RDKit [1] allows performing multiple functions on chemical compounds. One is the quantitative estimation of drug-likeness also known as QED properties. These properties include molecular weight (MW), octanol-water partition coefficient (ALOGP), number of hydrogen bond donors (HBD), number of hydrogen bond acceptors (HBA), polar surface area (PSA), number of rotatable bonds (ROTB), number of aromatic rings (AROM), structural alerts (ALERTS). (more…)

Continue Reading

Bioinformatics Programming

sdftosmi.py: Convert multiple ligands/compounds in SDF format to SMILES.

Dr. Muniba Faiza

Published

on

sdftosmi.py: Convert multiple ligands/compounds in SDF format to SMILES?

You can obtain SMILES of multiple compounds or ligands in an SDF file in one go. Here, we provide a simple Python script to do that. (more…)

Continue Reading

Bioinformatics Programming

tanimoto_similarities_one_vs_all.py – Python script to calculate Tanimoto Similarities of multiple compounds

Dr. Muniba Faiza

Published

on

tanimoto_similarities_one_vs_all.py – Python script to calculate Tanimoto Similarities of a compound with multiple compounds

We previously provided a Python script to calculate the Tanimoto similarities of multiple compounds against each other. In this article, we are providing another Python script to calculate the Tanimoto similarities of one compound with multiple compounds. (more…)

Continue Reading

Bioinformatics Programming

tanimoto_similarities.py: A Python script to calculate Tanimoto similarities of multiple compounds using RDKit.

Dr. Muniba Faiza

Published

on

tanimoto_similarities.py: A Python script to calculate Tanimoto similarities of multiple compounds using RDKit.

RDKit [1] is a very nice cheminformatics software. It allows us to perform a wide range of operations on chemical compounds/ ligands. We have provided a Python script to perform fingerprinting using Tanimoto similarity on multiple compounds using RDKit. (more…)

Continue Reading

Bioinformatics Programming

How to commit changes to GitHub repository using vs code?

Published

on

How to commit changes to GitHub repository using vs code?

In this article, we are providing a few commands that are used to commit changes to GitHub repositories using VS code terminal.

(more…)

Continue Reading

Bioinformatics Programming

Extracting first and last residue from helix file in DSSP format.

Dr. Muniba Faiza

Published

on

Extracting first and last residue from helix file in DSSP format.

Previously, we have provided a tutorial on using dssp_parser to extract all helices including long and short separately. Now, we have provided a new python script to find the first and last residue in each helix file. (more…)

Continue Reading

Bioinformatics Programming

How to extract x,y,z coordinates of atoms from PDB file?

Published

on

How to extract x,y,z coordinates of atoms from PDB file?

The x, y, and z coordinates of atoms are provided in the PDB file. One way to extract them is by using the Biopython package [1]. In this article, we will extract coordinates of C-alpha atoms for each residue from the PDB file using Biopython. (more…)

Continue Reading

LATEST ISSUE

ADVERT