Source code for MDAnalysis.analysis.contacts

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4## MDAnalysis --- https://www.mdanalysis.org# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors# (see the file AUTHORS for the full list of names)## Released under the GNU Public Licence, v2 or any higher version## Please cite your use of MDAnalysis in published work:## R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.# MDAnalysis: A Python package for the rapid analysis of molecular dynamics# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.# doi: 10.25080/majora-629e541a-00e## N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787#"""Native contacts analysis --- :mod:`MDAnalysis.analysis.contacts`================================================================This module contains classes to analyze native contacts *Q* over atrajectory. Native contacts of a conformation are contacts that existin a reference structure and in the conformation. Contacts in thereference structure are always defined as being closer then a distance`radius`. The fraction of native contacts for a conformation can becalculated in different ways. This module supports 3 different metricslisted below, as well as custom metrics.1. *Hard Cut*: To count as a contact the atoms *i* and *j* have to be at least as close as in the reference structure.2. *Soft Cut*: The atom pair *i* and *j* is assigned based on a soft potential that is 1 if the distance is 0, 1/2 if the distance is the same as in the reference and 0 for large distances. For the exact definition of the potential and parameters have a look at function :func:`soft_cut_q`.3. *Radius Cut*: To count as a contact the atoms *i* and *j* cannot be further apart than some distance `radius`.The "fraction of native contacts" *Q(t)* is a number between 0 and 1 andcalculated as the total number of native contacts for a given time framedivided by the total number of contacts in the reference structure.Examples for contact analysis-----------------------------One-dimensional contact analysis~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~As an example we analyze the opening ("unzipping") of salt bridgeswhen the AdK enzyme opens up; this is one of the example trajectoriesin MDAnalysis. :: import numpy as np import matplotlib.pyplot as plt import MDAnalysis as mda from MDAnalysis.analysis import contacts from MDAnalysis.tests.datafiles import PSF,DCD # example trajectory (transition of AdK from closed to open) u = mda.Universe(PSF,DCD) # crude definition of salt bridges as contacts between NH/NZ in ARG/LYS and # OE*/OD* in ASP/GLU. You might want to think a little bit harder about the # problem before using this for real work. sel_basic = "(resname ARG LYS) and (name NH* NZ)" sel_acidic = "(resname ASP GLU) and (name OE* OD*)" # reference groups (first frame of the trajectory, but you could also use a # separate PDB, eg crystal structure) acidic = u.select_atoms(sel_acidic) basic = u.select_atoms(sel_basic) # set up analysis of native contacts ("salt bridges"); salt bridges have a # distance <6 A ca1 = contacts.Contacts(u, selection=(sel_acidic, sel_basic), refgroup=(acidic, basic), radius=6.0) # iterate through trajectory and perform analysis of "native contacts" Q ca1.run() # print number of averave contacts average_contacts = np.mean(ca1.timeseries[:, 1]) print('average contacts = {}'.format(average_contacts)) # plot time series q(t) fig, ax = plt.subplots() ax.plot(ca1.timeseries[:, 0], ca1.timeseries[:, 1]) ax.set(xlabel='frame', ylabel='fraction of native contacts', title='Native Contacts, average = {:.2f}'.format(average_contacts)) fig.show()The first graph shows that when AdK opens, about 20% of the saltbridges that existed in the closed state disappear when the enzymeopens. They open in a step-wise fashion (made more clear by the movie`AdK_zipper_cartoon.avi`_)... _`AdK_zipper_cartoon.avi`: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2803350/bin/NIHMS150766-supplement-03.avi.. rubric:: NotesSuggested cutoff distances for different simulations* For all-atom simulations, cutoff = 4.5 Å* For coarse-grained simulations, cutoff = 6.0 ÅTwo-dimensional contact analysis (q1-q2)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Analyze a single DIMS transition of AdK between its closed and openconformation and plot the trajectory projected on q1-q2 [Franklin2007]_ :: import MDAnalysis as mda from MDAnalysis.analysis import contacts from MDAnalysisTests.datafiles import PSF, DCD u = mda.Universe(PSF, DCD) q1q2 = contacts.q1q2(u, 'name CA', radius=8) q1q2.run() f, ax = plt.subplots(1, 2, figsize=plt.figaspect(0.5)) ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 1], label='q1') ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 2], label='q2') ax[0].legend(loc='best') ax[1].plot(q1q2.timeseries[:, 1], q1q2.timeseries[:, 2], '.-') f.show()Compare the resulting pathway to the `MinActionPath result for AdK`_[Franklin2007]_... _MinActionPath result for AdK: http://lorentz.dynstr.pasteur.fr/joel/adenylate.phpWriting your own contact analysis~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~The :class:`Contacts` class has been designed to be extensible for your ownanalysis. As an example we will analyze when the acidic and basic groups of AdKare in contact which each other; this means that at least one of the contactsformed in the reference is closer than 2.5 Å.For this we define a new function to determine if any contact is closer than2.5 Å; this function must implement the API prescribed by :class:`Contacts`:: def is_any_closer(r, r0, dist=2.5): return np.any(r < dist)The first two parameters `r` and `r0` are provided by :class:`Contacts` when itcalls :func:`is_any_closer` while the others can be passed as keyword argsusing the `kwargs` parameter in :class:`Contacts`.Next we are creating an instance of the :class:`Contacts` class and use the:func:`is_any_closer` function as an argument to `method` and run the analysis:: # crude definition of salt bridges as contacts between NH/NZ in ARG/LYS and # OE*/OD* in ASP/GLU. You might want to think a little bit harder about the # problem before using this for real work. sel_basic = "(resname ARG LYS) and (name NH* NZ)" sel_acidic = "(resname ASP GLU) and (name OE* OD*)" # reference groups (first frame of the trajectory, but you could also use a # separate PDB, eg crystal structure) acidic = u.select_atoms(sel_acidic) basic = u.select_atoms(sel_basic) nc = contacts.Contacts(u, selection=(sel_acidic, sel_basic), method=is_any_closer, refgroup=(acidic, basic), kwargs={'dist': 2.5}) nc.run() bound = nc.timeseries[:, 1] frames = nc.timeseries[:, 0] f, ax = plt.subplots() ax.plot(frames, bound, '.') ax.set(xlabel='frame', ylabel='is Bound', ylim=(-0.1, 1.1)) f.show()Functions---------.. autofunction:: hard_cut_q.. autofunction:: soft_cut_q.. autofunction:: radius_cut_q.. autofunction:: contact_matrix.. autofunction:: q1q2Classes-------.. autoclass:: Contacts :members:"""from__future__importdivision,absolute_importfromsix.movesimportzipimportosimporterrnoimportwarningsimportbz2importnumpyasnpimportloggingimportMDAnalysisimportMDAnalysis.lib.distancesfromMDAnalysis.lib.utilimportopenany,deprecatefromMDAnalysis.analysis.distancesimportdistance_arrayfromMDAnalysis.core.groupsimportAtomGroupfrom.baseimportAnalysisBaselogger=logging.getLogger("MDAnalysis.analysis.contacts")

[docs]defhard_cut_q(r,cutoff):"""Calculate fraction of native contacts *Q* for a hard cut off. The cutoff can either be a float or a :class:`~numpy.ndarray` of the same shape as `r`. Parameters ---------- r : ndarray distance matrix cutoff : ndarray | float cut off value to count distances. Can either be a float of a ndarray of the same size as distances Returns ------- Q : float fraction of contacts """r=np.asarray(r)cutoff=np.asarray(cutoff)y=r<=cutoffreturny.sum()/r.size

[docs]defcontact_matrix(d,radius,out=None):"""calculate contacts from distance matrix Parameters ---------- d : array-like distance matrix radius : float distance below which a contact is formed. out: array (optional) If `out` is supplied as a pre-allocated array of the correct shape then it is filled instead of allocating a new one in order to increase performance. Returns ------- contacts : ndarray boolean array of formed contacts """ifoutisnotNone:out[:]=d<=radiuselse:out=d<=radiusreturnout

[docs]classContacts(AnalysisBase):"""Calculate contacts based observables. The standard methods used in this class calculate the fraction of native contacts *Q* from a trajectory. .. rubric:: Contact API By defining your own method it is possible to calculate other observables that only depend on the distances and a possible reference distance. The **Contact API** prescribes that this method must be a function with call signature ``func(r, r0, **kwargs)`` and must be provided in the keyword argument `method`. Attributes ---------- timeseries : list list containing *Q* for all refgroup pairs and analyzed frames """def__init__(self,u,selection,refgroup,method="hard_cut",radius=4.5,kwargs=None,**basekwargs):""" Parameters ---------- u : Universe trajectory selection : tuple(string, string) two contacting groups that change over time refgroup : tuple(AtomGroup, AtomGroup) two contacting atomgroups in their reference conformation. This can also be a list of tuples containing different atom groups radius : float, optional (4.5 Angstroms) radius within which contacts exist in refgroup method : string | callable (optional) Can either be one of ``['hard_cut' , 'soft_cut']`` or a callable with call signature ``func(r, r0, **kwargs)`` (the "Contacts API"). kwargs : dict, optional dictionary of additional kwargs passed to `method`. Check respective functions for reasonable values. verbose : bool (optional) Show detailed progress of the calculation if set to ``True``; the default is ``False``. """self.u=usuper(Contacts,self).__init__(self.u.trajectory,**basekwargs)ifmethod=='hard_cut':self.fraction_contacts=hard_cut_qelifmethod=='soft_cut':self.fraction_contacts=soft_cut_qelse:ifnotcallable(method):raiseValueError("method has to be callable")self.fraction_contacts=methodself.selection=selectionself.grA=u.select_atoms(selection[0])self.grB=u.select_atoms(selection[1])# contacts formed in referenceself.r0=[]self.initial_contacts=[]ifisinstance(refgroup[0],AtomGroup):refA,refB=refgroupself.r0.append(distance_array(refA.positions,refB.positions))self.initial_contacts.append(contact_matrix(self.r0[-1],radius))else:forrefA,refBinrefgroup:self.r0.append(distance_array(refA.positions,refB.positions))self.initial_contacts.append(contact_matrix(self.r0[-1],radius))self.fraction_kwargs=kwargsifkwargsisnotNoneelse{}self.timeseries=[]def_single_frame(self):# compute distance array for a framed=distance_array(self.grA.positions,self.grB.positions)y=np.empty(len(self.r0)+1)y[0]=self._ts.framefori,(initial_contacts,r0)inenumerate(zip(self.initial_contacts,self.r0)):# select only the contacts that were formed in the reference stater=d[initial_contacts]r0=r0[initial_contacts]y[i+1]=self.fraction_contacts(r,r0,**self.fraction_kwargs)iflen(y)==1:y=y[0]self.timeseries.append(y)def_conclude(self):self.timeseries=np.array(self.timeseries,dtype=float)@deprecate(release="0.19.0",remove="1.0.0")defsave(self,outfile):"""save contacts timeseries Parameters ---------- outfile : str file to save contacts """np.savetxt(outfile,self.timeseries,header="# q1 analysis\n",comments='')