Finally, after two years of social distancing, we were able to continue the tradition of OPIGtreat – a 2-3 day escape to the countryside for a packed schedule of talks and fun.
This year, the lovely YHA Wilderhope Manor in Shropshire was chosen by Lewis, our trip organizer. With a hostel in the middle of nowhere, with no phone signal, this trip promised to be an exciting get-away from our plugged-in lives at the university.
Continue readingComputing is one of the only scientific fields which was once female-dominated. In the 30s and 40s, women made up the bulk of the workforce doing complex, tedious calculations in the fields including ballistics, astrophysics, aeronautics (think Hidden Figures) and code-breaking. Engineers themselves found that the female computers were far more reliable than themselves in doing such calculations [9]. As computing machines became available, there was no precedent set for the gender of a computer operator, and so the women previously doing the computing became the computer operators [10].
However, this was not to last. As computing became commercialised in the 50s, the skill required for computing work was starting to be recognised. As written in [1]:
“Software company System Development Corp. (SDC) contracted psychologists William Cannon and Dallis Perry to create an aptitude assessment for optimal programmers. Cannon and Perry interviewed 1,400 engineers — 1,200 of them men — and developed a “vocational interest scale,” a personality profile to predict the best potential programmers. Unsurprisingly given their male-dominated test group, Cannon and Perry’s assessment disproportionately identified men as the ideal candidates for engineering jobs. In particular, the test tended to eliminate extroverts and people who have empathy for others. Cannon and Perry’s paper concluded that typical programmers “don’t like people,” forming today’s now pervasive stereotype of a nerdy, anti-social coder.”
Continue readingMolecular dynamics (MD) simulations are a good way to explore the dynamical behaviour of a protein you might be interested in. One common problem is that they often have a relatively steep learning curve when using most MD engines.
What if you just want to run a simple, one-off simulation with no fancy enhanced sampling methods? OpenMM Setup is a useful tool for exactly this. It is built on the open-source OpenMM engine and provides an easy to install (via conda) GUI that can have you running a simulation in less than 5 minutes. Of course, running a simulation requires careful setting of parameters and being familiar with best practices and while this is beyond the scope of this post, there are many guides out there that can easily be found. Now on to the good stuff: using OpenMM Setup!
When you first run OpenMM Setup, you’ll be greeted by a browser window asking you to choose a structure to use. This can be a crystal structure or a model. Remember, sometimes these will have problems that need fixing like missing density or charged, non-physiological termini that would lead to artefacts, so visual inspection of the input is key! You can then choose the force field and water model you want to use, and tell OpenMM to do some cleaning up of the structure. Here I am running the simulation on hen egg-white lysozyme:
RDKit is very fussy when it comes to inputs in SDF format. Using the SDMolSupplier, we get a significant rate of failure even on curated datasets such as the PDBBind refined set. Pymol has no such scruples, and with that, I present a function which has proved invaluable to me over the course of my DPhil. For reasons I have never bothered to explore, using pymol to convert from sdf, into mol2 and back to sdf format again (adding in missing hydrogens along the way) will almost always make a molecule safe to import using RDKit:
from pathlib import Path
from pymol import cmd
def py_mollify(sdf, overwrite=False):
"""Use pymol to sanitise an SDF file for use in RDKit.
Arguments:
sdf: location of faulty sdf file
overwrite: whether or not to overwrite the original sdf. If False,
a new file will be written in the form <sdf_fname>_pymol.sdf
Returns:
Original sdf filename if overwrite == False, else the filename of the
sanitised output.
"""
sdf = Path(sdf).expanduser().resolve()
mol2_fname = str(sdf).replace('.sdf', '_pymol.mol2')
new_sdf_fname = sdf if overwrite else str(sdf).replace('.sdf', '_pymol.sdf')
cmd.load(str(sdf))
cmd.h_add('all')
cmd.save(mol2_fname)
cmd.reinitialize()
cmd.load(mol2_fname)
cmd.save(str(new_sdf_fname))
return new_sdf_fname
It is possible to get an installer for the crystallographer’s favourite molecular visualization tool for Windows machines, that is if you are willing to pay a fee. Fortunately, Christoph Gohlke has made available free, pre-compiled Windows versions of the latest PyMOL software, along with all of it’s requirements, it’s just not particularly straightforward to install. The PyMOLWiki offers a three-step guide on how to do this and I will break it down to make it somewhat clearer.
Download the Windows Installer (x-bit) for Python 3 from their website, x being your Windows architecture – 32 or 64.
Then, follow the instructions on how to install it. You can check if it has installed by running the following in PowerShell:
Continue readingI’m going to keep this one brief, because I am mid-confirmation-and-paper-writing madness. I have seen too many people – both beginners and seasoned veterans – wandering around their Linux filesystem blindfolded:
Whenever you want to see where you are, you have to execute pwd (present working directory), which will print your absolute location to stdout. If you have many terminals open at the same time, it is easy to lose track of where you are, and every other command becomes pwd; surely, I hear you cry, there has to be a better way!
Well, fear not! With a little tinkering with ~/.bashrc, we can display the working directory as part of the special PS1 environment variable, responsible for how your username and computer are displayed above. Putting the following at the top of ~/.bashrc
me=`id | awk -F\( '{print $2}' | awk -F\) '{print $1}'`
export PS1="`uname -n | /bin/sed 's/\..*//'`{$me}:\$PWD$ "… saving, and starting a new termanal window results in:
I haven’t used pwd in 3 years.
We have reached the era of design, not just ‘hunting’. Particularly exciting to me is the de novo design of proteins, which have a wide and ever increasing range of applications from therapeutics to consumer products, biomanufacturing to biomaterials. Protein design has been a) enabled by decades of research that contributed to our understanding of protein sequence, structure & function and b) accelerated by computational advances – capturing the information we have learned from proteins and representing it for computers and machine learning algorithms.
In this blog post, I will discuss three key methodological considerations for computational protein design:
Back-of-the-envelope calculations are one of our chief tools as scientists. When you spend most of your time wondering if your latest measurement is correct, having a tool to check if the numbers make sense is simply priceless. If you are lucky, a good estimate might just avoid a costly or laborious measurement — this is very common in disciplines like chemical engineering, which a friend described as “the art of estimating numbers and plugging them into some variation of Bernoulli’s continuity equation”. Unsurprisingly, these Fermi problems are now common interview questions at major consultancy and tech companies, and have even started to go viral.
Last week, I thought I would ask my biochemistry students to solve a back-of-the-envelope problem as part of their tutorial work. Disguised as an enzyme catalysis problem, I asked them to estimate the energy of a single hydrogen bond. Needless to say, they were puzzled. Some of them asked if I had forgotten to include some information in the problem sheet. For some reason, Fermi problems seem to be less common in chemistry and biology that they are in physics of engineering. Of course, estimating the energy of a hydrogen bond is in many ways much harder than guessing the number of ping pong balls that fit a Boeing 747. Nobody has seen a hydrogen bond in the flesh. And our minds struggle to grasp the vast numbers present at the molecular level. Nevertheless, guesstimates are incredibly useful
Continue readingAutoDock4 and AutoDock Vina are the most commonly used open-source software for protein-ligand docking. However, they both rely on a derivative of the “PDB” (Protein Data Base) file format: the “PDBQT” file (Protein Data Bank, Partial Charge (Q), & Atom Type (T)). In addition to the information contained in normal PDB files, PDBQT files have an additional column that lists the partial charge (Q) and the assigned AutoDock atom type (T) for each atom in the molecule. AutoDock atom types offer a more granular differentiation between atoms such as listing aliphatic carbons and aromatic carbons as separate AutoDock atom types.
The biggest drawback about the PDBQT format is that it does not encode for the bond order in molecules explicitly. Instead, the bond order is inferred based on the atom type, distance and angle to nearby atoms in the molecule. For normal sp3 carbons and molecules with mostly single bonds this system works fine, however, for more complex structures containing for example aromatic rings, conjugated systems and hypervalent atoms such as sulphur, the bond order is often not displayed correctly. This leads to issues downstream in the screening pipeline when molecules suddenly change their bond order or have to be discarded after docking because of impossible bond orders.
The solution to this problem is included in RDKit: The AssignBondOrdersFromTemplate function. All you need to do is load the original molecule used for docking as a template molecule and the docked pose PDBQT file into RDKIT as a PDB, without the bond order information. Then assign the original bond order from your template molecule. The following code snippet covers the necessary functions and should help you build a more accurate and reproducible protein-ligand docking pipeline:
#import RDKit AllChem from rdkit import Chem from rdkit.Chem import AllChem #load original molecule from smiles SMILES_STRING = "CCCCCCCCC" #the smiles string of your ligand template = Chem.MolFromSmiles(SMILES_STRING) #load the docked pose as a PDB file loc_of_docked_pose = "docked_pose_mol.pdb" #file location of the docked pose converted to PDB file docked_pose = AllChem.MolFromPDBFile(loc_of_docked_pose) #Assign the bond order to force correct valence newMol = AllChem.AssignBondOrdersFromTemplate(template, docked_pose) #Add Hydrogens if desired. "addCoords = True" makes sure the hydrogens are added in 3D. This does not take pH/pKa into account. newMol_H = Chem.AddHs(newMol, addCoords=True) #save your new correct molecule as a sdf file that encodes for bond orders correctly output_loc = "docked_pose_assigned_bond_order.sdf" #output file name Chem.MolToMolFile(newMol_H, output_loc)
We are all familiar with the idea of a correlation. In the broadest sense of the word, a correlation can refer to any kind of dependence between two variables. There are three widely used tests for correlation:
While these three measures give us plenty of options to work with, they do not work in all cases. Take for example the following variables, Y1 and Y2. These might be two variables that vary in a concerted manner.
Perhaps we suspect that a state change in Y1 leads to a state change in Y2 or vice versa and we want to measure the association between these variables. Using the three measures of correlation, we get the following results:
Continue reading