For a long time crystallographers and subsequently the authors of AlphaFold2 had you believe that proteins are a static group of atoms written to a .pdb file. Turns out this was a HOAX. If you don’t want to miss out on the latest trend of working with dynamic structural ensembles of proteins this blog post is exactly right for you. MDAnalysis is a python package which as the name says was designed to analyse molecular dyanmics simulation and lets you work with trajectories of protein structures easily. <\/p>\n\n\n\n\n\n\n\n
Before we start install MDAnalysis through pip.<\/p>\n\n\n\n
pip install MDAnalysis<\/code><\/pre>\n\n\n\nThe main data structure MDAnalysis uses is the Universe<\/mark><\/code>. You can create a Universe<\/mark><\/code> by loading a trajectory file and a topology file. The trajectory file contains the trajectory of the protein structure as a list of coordinates. Typical formats are .xtc or .dcd. The topology file provides information on the identity of atoms and their connectivity. A pdb file can be used to this purpose.<\/p>\n\n\n\nimport MDAnalysis as mda\n\nu = mda.Universe('topology.pdb', 'trajectory.xtc')<\/code><\/pre>\n\n\n\nThe Universe<\/mark><\/code> now contains all the information of the protein trajectory and lets you perform operations on it.<\/p>\n\n\n\nThrough the Universe<\/mark><\/code> specific groups of residues and atoms and individual atoms can be accessed. Residues and atoms can be selected by indexing or a string selection.<\/p>\n\n\n\n# select all residues<\/strong>\nprint(u.residues)\n<ResidueGroup with 259 residues><\/mark>\n\n# select all C-alphas of alanine and glutamate residues<\/strong>\nprint(u.select_atoms('resname ALA or resname GLU and name CA'))\n<AtomGroup with 25 atoms><\/mark>\n\n# select first atom<\/strong>\nu.atoms[0]\n<Atom 1: N of type N of resname ASP, resid 1 and segid P000 and altLoc ><\/mark><\/code><\/pre>\n\n\n\nInformation can be accessed from ResidueGroup<\/mark><\/code>, AtomGroup<\/mark><\/code> and Atom<\/mark><\/code> objects.<\/p>\n\n\n\n# residue names<\/strong>\nu.residues[:3].resnames\narray(['ASP', 'ILE', 'GLN'], dtype=object)<\/mark>\n\n# atom coordinates in initial frame<\/strong>\nu.atoms[:3].positions\narray([[224.13002 , 30.02 , 729.97003 ],\n [224.34003 , 29.250002, 730.65 ],\n [223.12 , 30.260002, 730.04004 ],\n ...,\n [240.33002 , 61.730003, 749.9 ],\n [240.64001 , 62.640003, 749.09 ],\n [241.05002 , 60.720005, 750.13 ]], dtype=float32)<\/mark><\/code><\/pre>\n\n\n\nAccessing atom coordinates as above will return the positions in the first frame in the trajectory. If you want to access subsequent frames you need to work with the Universe.trajectory attribute. Information on frames is accessed by iterating through the trajectory or indexing it. At any stage the Universe<\/mark><\/code> will only contain information of a single timestep.<\/p>\n\n\n\n### radius of gyration over time<\/mark><\/strong>\nrgyr = []\ntime = []\nprotein = u.select_atoms('protein')\n\nfor ts in u.trajectory: # loop through individual frames in the trajectory<\/strong>\n time.append(u.trajectory.time)\n rgyr.append(protein.radius_of_gyration())\n\n\n### backbone RMSD against time<\/strong>\nfrom MDAnalysis.analysis import rms\n\nu.trajectory[0] # access first frame by indexing<\/strong>\nbb_0 = u.select_atoms('backbone').positions\n\ntime = []\nrmsd = []\nfor ts in u.trajectory:\n time.append(u.trajectory.time)\n bb_t = u.select_atoms('backbone').positions\n rmsd.append(rms.rmsd(bb_0, bb_t))<\/code><\/pre>\n\n\n\n