{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Write your own native contacts analysis method\n", "\n", "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.\n", "\n", "**Last updated:** December 2022 with MDAnalysis 2.4.0-dev0\n", "\n", "**Minimum version of MDAnalysis:** 1.0.0\n", "\n", "**Packages required:**\n", " \n", "* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n", "* MDAnalysisTests\n", "* [matplotlib](https://matplotlib.org)\n", "* [pandas](https://pandas.pydata.org)\n", "\n", "**See also**\n", "\n", "* [Fraction of native contacts over a trajectory](contacts_native_fraction.ipynb) (pre-defined metrics and a general introduction to native contacts analysis)\n", "* [Q1 vs Q2 contact analysis](contacts_q1q2.ipynb)\n", "* [Contact analysis: number of contacts within a cutoff](contacts_within_cutoff.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:00.563225Z", "iopub.status.busy": "2021-05-19T05:57:00.562481Z", "iopub.status.idle": "2021-05-19T05:57:01.656037Z", "shell.execute_reply": "2021-05-19T05:57:01.656427Z" } }, "outputs": [], "source": [ "import MDAnalysis as mda\n", "from MDAnalysis.tests.datafiles import PSF, DCD\n", "from MDAnalysis.analysis import contacts\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading files\n", "\n", "The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. (Beckstein *et al.*, 2009) The trajectory ``DCD`` samples a transition from a closed to an open conformation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:01.659952Z", "iopub.status.busy": "2021-05-19T05:57:01.659362Z", "iopub.status.idle": "2021-05-19T05:57:01.891750Z", "shell.execute_reply": "2021-05-19T05:57:01.892356Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/pbarletta/mambaforge/envs/mda-user-guide/lib/python3.9/site-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behaviour will be changed in 3.0 to be the same as other readers\n", " warnings.warn(\"DCDReader currently makes independent timesteps\"\n" ] } ], "source": [ "u = mda.Universe(PSF, DCD)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Background\n", "\n", "Please see the [Fraction of native contacts](contacts_native_fraction.ipynb#Background) for an introduction to general native contacts analysis." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining salt bridges" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:01.898428Z", "iopub.status.busy": "2021-05-19T05:57:01.897597Z", "iopub.status.idle": "2021-05-19T05:57:01.902430Z", "shell.execute_reply": "2021-05-19T05:57:01.902789Z" } }, "outputs": [], "source": [ "sel_basic = \"(resname ARG LYS) and (name NH* NZ)\"\n", "sel_acidic = \"(resname ASP GLU) and (name OE* OD*)\"\n", "acidic = u.select_atoms(sel_acidic)\n", "basic = u.select_atoms(sel_basic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define your own function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Any function you define *must* have `r` and `r0` as its first and second arguments respectively, even if you don't necessarily use them:\n", "\n", " - `r`: an array of distances between atoms at the current time\n", " - `r0`: an array of distances between atoms in the reference\n", "\n", "You can then define following arguments as keyword arguments.\n", "\n", "In the function below, we calculate the fraction of native contacts that are less than `radius`, but greater than `min_radius`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:01.907350Z", "iopub.status.busy": "2021-05-19T05:57:01.906760Z", "iopub.status.idle": "2021-05-19T05:57:01.908485Z", "shell.execute_reply": "2021-05-19T05:57:01.908901Z" } }, "outputs": [], "source": [ "def fraction_contacts_between(r, r0, radius=3.4, min_radius=2.5):\n", " is_in_contact = (r < radius) & (r > min_radius) # array of bools\n", " fraction = is_in_contact.sum()/r.size\n", " return fraction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:01.912994Z", "iopub.status.busy": "2021-05-19T05:57:01.912470Z", "iopub.status.idle": "2021-05-19T05:57:01.938317Z", "shell.execute_reply": "2021-05-19T05:57:01.938722Z" } }, "outputs": [], "source": [ "ca = contacts.Contacts(u, \n", " select=(sel_acidic, sel_basic),\n", " refgroup=(acidic, basic),\n", " method=fraction_contacts_between,\n", " radius=5.0,\n", " kwargs={'radius': 5.0,\n", " 'min_radius': 2.4}).run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One easy way to post-process results is to turn them into a dataframe." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-05-19T05:57:01.945491Z", "iopub.status.busy": "2021-05-19T05:57:01.944493Z", "iopub.status.idle": "2021-05-19T05:57:01.953817Z", "shell.execute_reply": "2021-05-19T05:57:01.954467Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", " | Frame | \n", "Contacts from first frame | \n", "
---|---|---|
0 | \n", "0.0 | \n", "1.000000 | \n", "
1 | \n", "1.0 | \n", "0.988764 | \n", "
2 | \n", "2.0 | \n", "0.943820 | \n", "
3 | \n", "3.0 | \n", "0.943820 | \n", "
4 | \n", "4.0 | \n", "0.943820 | \n", "