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 executed: May 18, 2021 with MDAnalysis 1.1.1

Last updated: June 29, 2020 with MDAnalysis 1.0.0

Minimum version of MDAnalysis: 1.0.0

Packages required:

MDAnalysis ([MADWB11], [GLB+16])
MDAnalysisTests
matplotlib
pandas

See also

Fraction of native contacts over a trajectory (pre-defined metrics and a general introduction to native contacts analysis)
Q1 vs Q2 contact analysis
Contact analysis: number of contacts within a cutoff

[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 time
r0: 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

Plotting¶

We can plot directly from a dataframe (below), or you could use it with other plotting packags such as seaborn.

[7]:

ca_df.plot(x='Frame')

[7]:

<matplotlib.axes._subplots.AxesSubplot at 0x7f9304981690>

../../../_images/examples_analysis_distances_and_contacts_contacts_custom_16_1.png

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.