Histograms are frequently used to visualize the distribution of a data set or to compare between multiple distributions. Python, via matplotlib.pyplot, contains convenient functions for plotting histograms; the default plots it generates, however, leave much to be desired in terms of visual appeal and clarity.
The two code blocks below generate histograms of two normally distributed sets using default matplotlib.pyplot.hist settings and then, in the second block, I add some lines to improve the data presentation. See the comments to determine what each individual line is doing.
Author Archives: Daniel Nissley
Writing “vector trajectories” with cpptraj
The program cpptraj, written by Daniel Roe (https://github.com/Amber-MD/cpptraj) and distributed Open Source with the AmberTools package (https://ambermd.org/AmberTools.php), is a powerful tool for analysis of molecular dynamics simulations. In addition to all of the expected analyses like Root Mean Square Deviation and native contacts, cpptraj also includes a suite of vector algebra functions.
While this vector algebra functionality is fairly well known and easy to find in the documentation, I think it is less well known that cpptraj can write trajectories of the computed vectors. These trajectories can then be loaded into Visual Molecular Dynamics (VMD) alongside the analysed trajectory and played as a movie. This functionality is a valuable tool for debugging your vector calculations to make sure they are doing precisely what you intend. It may also prove useful for generating visualizations of vectors alongside molecular structures for publications.
The cpptraj script below reads in an Amber parameter file and coordinate file and then calculates the angle between two planes.
parm 7bbg_fixed.prmtop
trajin 7bbg_fixed.rst7
vector v1 mask :65@NH1 :65@NH2
vector v2 mask :65@NH1 :65@NE
vector v3 mask :64@CA :66@CA
vector v4 mask :66@CA :68@CA
vectormath vec1 v1 vec2 v2 crossproduct name n1
vectormath vec1 v3 vec2 v4 crossproduct name n2
vectormath vec1 n1 vec2 n2 dotangle out 7bbg_ref_plane_angle.dat
The first plane is defined by two vectors in the plane of the guanidino group of a R65 residue (v1 and v2); the second plane is defined by two vectors between CA atoms of amino acids in the alpha helix containing R65 (v3 and v4). The first two vectormath calls determine the normal vectors to the planes and the final vectormath line computes the angle between the normal vectors. Taken together, these commands compute the angle between the arginine side chain and a plane passing through the CA atoms of the alpha helix. Let’s check that the vectors {v1, v2, v3, v4} are being computed correctly.
parm 7bbg_fixed.prmtop
trajin 7bbg_fixed.rst7
vector v1 mask :65@NH1 :65@NH2
vector v2 mask :65@NH1 :65@NE
vector v3 mask :64@CA :66@CA
vector v4 mask :66@CA :68@CA
run
writedata vectors.mol2 v1 v2 v3 v4 vectraj trajfmt mol2
The resulting vector trajectory vectors.mol2 can be loaded directly into VMD without a topology. Note that in this case we only analyzed a single frame, but you can run this same procedure on DCD files, too. This is what I get when I load the vectors into VMD alongside the structure:
I hope this vector trajectory functionality will be helpful to a few people who like to neurotically check their analyses like I do. You can download the example prmtop and rst7 files below. Note that you should rename them to remove the extra “.txt” file extension before attempting to use them for anything.
The information in this blog post is adapted from an Amber Archive post from Daniel Roe, dated 30-Oct-2018: http://archive.ambermd.org/201811/0058.html
Files for the example:
Coarse-grained models of antibody solutions
Various coarse-grained (CG) models have become increasingly common in studies of antibody-antibody interactions in solution. These models appear poised to enter development pipelines in the near future to help predict and understand how antibody-antibody interactions influence the suitability of a given monoclonal antibody (mAb) for mass production and delivery as an antibody therapy. This blog post is a non-exhaustive summary of some of the highlights I found during a recent literature search.
Continue readingDo you have cis peptide bonds in your simulation inputs?
People who run molecular simulations quickly become familiar with all of the things about a PDB file – missing residues, missing heavy atoms in residues, missing hydrogens, non-standard amino acids, multiple conformations, crystallization ligands, etc. – that might need to be fixed before setting up a simulation. This blog post is a reminder to check, after you have “fixed” your PDB, if you have accidentally introduced aberrant cis peptide bonds into your structure during rebuilding.
Continue readingMaking better plots with matplotlib.pyplot in Python3
The default plots made by Python’s matplotlib.pyplot module are almost always insufficient for publication. With a ~20 extra lines of code, however, you can generate high-quality plots suitable for inclusion in your next article.
Let’s start with code for a very default plot:
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(1)
d1 = np.random.normal(1.0, 0.1, 1000)
d2 = np.random.normal(3.0, 0.1, 1000)
xvals = np.arange(1, 1000+1, 1)
plt.plot(xvals, d1, label='data1')
plt.plot(xvals, d2, label='data2')
plt.legend(loc='best')
plt.xlabel('Time, ns')
plt.ylabel('RMSD, Angstroms')
plt.savefig('bad.png', dpi=300)
The result of this will be:
The fake data I generated for the plot look something like Root Mean Square Deviation (RMSD) versus time for a converged molecular dynamics simulation, so let’s pretend they are. There are a number of problems with this plot: it’s overall ugly, the color scheme is not very attractive and may not be color-blind friendly, the y-axis range of the data extends outside the range of the tick labels, etc.
We can easily convert this to a much better plot:
Continue readingMM(PB/GB)SA – a quick start guide
The MMPBSA.py program distributed Open Source in the AmberTools21 package is a powerful tool for end-point free energy calculations on molecular dynamics simulations. In its most simple application, MMPBSA.py is used to calculate the free energy difference between the bound and unbound states of a protein-ligand complex. In order to use it, however, you need to have an Amber-compliant trajectory file, which means you need to setup and run your simulation fairly carefully.
While the Amber Manual and the MMPBSA tutorial provide lots of helpful information, putting everything together into a full pipeline taking you from structure to a free energy is another story. The goal for this guide is to provide a schematic you can follow to get started. This guide assumes you are familiar with molecular dynamics simulations and the theory of MMPBSA.
The easiest way I have found to do this, using only Open Source software, is:
(1) Download your raw PDB file. If you are lucky and it contains a complete set of heavy atoms (excepting perhaps a terminal OXT here and there, which tleap will add for you in step 3) you are good to go.
(2) Use the H++ webserver to determine the protonation states of each residue and add hydrogens as needed. This webserver is particularly convenient because it will allow you to directly download a PQR file that you can use to generate your starting topology and coordinates. Note that you have various options to choose the pH and internal/external dielectric constants for the calculation.
(3) Use tleap to generate your topology (prmtop) and coordinate (mdcor) files for your simulations. Do not forget that you will need not only the prmtop for the solvated complex, but also a dry prmtop for each of the complex, receptor, and ligand. Load the PQR file from H++ and do not forget to set PBRadii *to the same value for all prmtops*. A typical tleap script for setting up your solvated complex would look something like:
Stochastic chemical kinetics – things randomly bumping into each other
In this blog post I describe the advantages of taking a stochastic view of chemical systems based on the work of D. T. Gillespie and subsequent publications. Gillespie presented his formalism for considering stochastic chemical kinetics, now referred to as the Gillespie Algorithm, in two papers published in 1976 and 1977 (Gillespie, D. T. J. Comp. Phys. 1976, Gillespie D. T. J. Phys. Chem. 1977) – if you want to see the full derivation for the Gillespie Algorithm along with many examples I recommend giving them both a read.
The essential question of chemical kinetics as stated by Gillespie is:
“If a fixed volume V contains a spatially uniform mixture of N chemical species which can inter-react through M specified chemical reaction channels, then given the numbers of molecules of each species present at some initial time, what will these molecular population levels be at any later time?”
Continue readingUnraveling the role of entanglement in protein misfolding
Proteins that fail to fold correctly may populate misfolded conformations with disparate structure and function. Misfolding is the focus of intense research interest due to its putative and confirmed role in various diseases, including neurodegenerative diseases such as Parkinson’s and Alzheimer’s Diseases as well as cystic fibrosis (PMID: 16689923).
Many open questions about protein misfolding remain to be answered. For example, how do misfolded proteins evade cellular quality control mechanisms like chaperones to remain soluble but non-functional for long timescales? How long do misfolded states persist on average? How widespread is misfolding? Experiments indicate that misfolding can even be caused by synonymous mutations that alter the speed of protein translation but not the sequence of the protein produced (PMID: 23417067), introducing the additional puzzle of how the protein maintains a “memory” of its translation kinetics after synthesis is complete.
A series of four recent preprints (Preprints 1, 2, 3, and 4, see below) suggests that these questions can be answered by the partitioning of proteins into long-lived self-entangled conformations that are structurally similar to the native state but with perturbed function. Simulation of the synthesis, termination, and post-translational dynamics of a large dataset of E. coli proteins suggests that misfolding and entanglement are widespread, with two thirds of proteins misfolding some of the time (Preprint 1). Many misfolded conformations may bypass proteostasis machinery to remain soluble but non-functional due to their structural similarity to the native state. Critically, entanglement is associated with particularly long-lived misfolded states based on simulated folding kinetics.
Coarse-grain and all-atom simulation results indicate that these misfolded conformations interact with chaperones like GroEL and HtpG to a similar extent as does the native state (Preprint 2). These results suggest an explanation for why some protein always fails to refold while remaining soluble, even in the presence of multiple folding chaperones – it remains trapped in entangled conformations that resemble the native state and thus fail to recruit chaperones.
Finally, simulations indicate that changes to the translation kinetics of oligoribonuclease introduced by synonymous mutations cause a large change in its probability of entanglement at the dimerization interface (Preprint 3). These entanglements localized at the interface alter its ability to dimerize even after synthesis is complete. These simulations provide a structural explanation for how translation kinetics can have a long-timescale influence on protein behavior.
Together, these preprints suggest that misfolding into entangled conformations is a widespread phenomenon that may provide a consistent explanation for many unanswered question in molecular biology. It should be noted that entanglement is not exclusive to other types of misfolding, such as domain swapping, that may contribute to misfolding in cells. Experimental validation of the existence of entangled conformations is a critical aspect of testing this hypothesis; for comparisons between simulation and experiment, see Preprint 4.
Preprint 1: https://www.biorxiv.org/content/10.1101/2021.08.18.456613v1
Preprint 2: https://www.biorxiv.org/content/10.1101/2021.08.18.456736v1
Preprint 3: https://www.biorxiv.org/content/10.1101/2021.10.26.465867v1
Preprint 4: https://www.biorxiv.org/content/10.1101/2021.08.18.456802v1
How fast can a protein fold?
A protein’s folding time is the time required for it to reach its unique folded state starting from its unfolded ensemble. Globular, cytosolic proteins can only attain their intended biological function once they have folded. This means that protein folding times, which typically exceed the timescales of enzymatic reactions that proteins carry out by several orders of magnitude, are critical to determining when proteins become functional. Many scientists have worked tirelessly over the years to measure protein folding times, determine their theoretical bounds, and understand how they fit into biology. Here, I focus on one of the more interesting questions to fall out of this field over the years: how fast can a protein fold? Note that this is a very different question than asking “how fast do proteins fold?”
Continue readingRibosome occupancy profiles are conserved between structurally and evolutionarily related yeast domains
Shameless plug for any OPIG blog readers to take a look at our recent publication in Bioinformatics. Consider giving it a read if the below summary grabs your attention.
Many proteins are now known to fold during their synthesis through the process known as co-translational folding. Translation is an inherently non-equilibrium process – one consequence of this fact is that the speed of translation can radically influence the ability of proteins to fold and function. In this paper we compare ribosome occupancy profiles between related domains in yeast to test the hypothesis that evolutionarily related proteins with similar native folds should tend to have similar translation speed profiles to preserve efficient co-translational folding. We find strong evidence in support of this hypothesis at the level of individual protein domains and across a set of 664 pairs of related domains for which we are able to compute high-quality ribosome occupancy profiles.
To find out more, view the Advance Article at Bioinformatics.