Write your own native contacts analysis method¶
The contacts.Contacts class has been designed to be extensible for your own analysis. Here we demonstrate how to define a new method to use to determine the fraction of native contacts.
Last updated: June 29, 2020 with MDAnalysis 1.0.0
Minimum version of MDAnalysis: 1.0.0
Packages required:
MDAnalysisTests
See also
Fraction of native contacts over a trajectory (pre-defined metrics and a general introduction to native contacts analysis)
[1]:
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD
from MDAnalysis.analysis import contacts
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
Loading files¶
The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. ([BDPW09]) The trajectory DCD samples a transition from a closed to an open conformation.
[2]:
u = mda.Universe(PSF, DCD)
Background¶
Please see the Fraction of native contacts for an introduction to general native contacts analysis.
Defining salt bridges¶
We define salt bridges as contacts between NH/NZ in ARG/LYS and OE*/OD* in ASP/GLU. You may not want to use this definition for real work.
[3]:
sel_basic = "(resname ARG LYS) and (name NH* NZ)"
sel_acidic = "(resname ASP GLU) and (name OE* OD*)"
acidic = u.select_atoms(sel_acidic)
basic = u.select_atoms(sel_basic)
Define your own function¶
Any function you define must have r and r0 as its first and second arguments respectively, even if you don’t necessarily use them:
r: an array of distances between atoms at the current timer0: an array of distances between atoms in the reference
You can then define following arguments as keyword arguments.
In the function below, we calculate the fraction of native contacts that are less than radius, but greater than min_radius.
[4]:
def fraction_contacts_between(r, r0, radius=3.4, min_radius=2.5):
is_in_contact = (r < radius) & (r > min_radius) # array of bools
fraction = is_in_contact.sum()/r.size
return fraction
Then we pass fraction_contacts_between to the contacts.Contacts class. Keyword arguments for our custom function must be in the kwargs dictionary. Even though we define a radius keyword in my custom analysis function, it is not automatically passed from contacts.Contacts. We have to make sure that it is in kwargs.
[5]:
ca = contacts.Contacts(u,
select=(sel_acidic, sel_basic),
refgroup=(acidic, basic),
method=fraction_contacts_between,
radius=5.0,
kwargs={'radius': 5.0,
'min_radius': 2.4}).run()
One easy way to post-process results is to turn them into a dataframe.
[6]:
ca_df = pd.DataFrame(ca.timeseries,
columns=['Frame',
'Contacts from first frame'])
ca_df.head()
[6]:
| Frame | Contacts from first frame | |
|---|---|---|
| 0 | 0.0 | 1.000000 |
| 1 | 1.0 | 0.988764 |
| 2 | 2.0 | 0.943820 |
| 3 | 3.0 | 0.943820 |
| 4 | 4.0 | 0.943820 |
References¶
[1] Oliver Beckstein, Elizabeth J. Denning, Juan R. Perilla, and Thomas B. Woolf. Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of Open↔Closed Transitions. Journal of Molecular Biology, 394(1):160–176, November 2009. 00107. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164, doi:10.1016/j.jmb.2009.09.009.
[2] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. Proceedings of the 15th Python in Science Conference, pages 98–105, 2016. 00152. URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.
[3] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry, 32(10):2319–2327, July 2011. 00778. URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.