Deformation Analysis

This example shows how to calculate the deformation vector describing the change between two structures of a protein. Two structures of the same protein in PDB format will be used. A Vector instance that contains the deformation vector describing the change in protein structure will be calculated. This object will be compared to ANM modes.

Parse structures

We start by importing everything from the ProDy package:

In [1]: from prody import *

In [2]: from matplotlib.pylab import *

In [3]: ion()

Let’s parse two p38 MAP Kinase structures: 5uoj and 1zz2

In [4]: reference = parsePDB('5uoj')

In [5]: mobile = parsePDB('1zz2')   # this is the one we want to superimpose

Match chains

ProDy offers the function matchChains() to find matching chains in two structures easily. We use it to find the chains for which we will calculate the deformation vector:

In [6]: matches = matchChains(reference.protein, mobile.protein)

matchChains() function returns a list. If there are no matching chains, list is empty, else the list contains a tuple for each pair of matching chains.

In [7]: len(matches)
Out[7]: 1

In [8]: match = matches[0]

There is only one match in this case. First item is a subset of atoms from the first structure (reference). Second item is a subset of atoms from the second structure (mobile).

In [9]: ref_chain = match[0]

In [10]: mob_chain = match[1]

Matched atoms are returned in AtomMap instances. We can get information on matched subset of atoms by entering the variable name:

In [11]: ref_chain
Out[11]: <AtomMap: Chain A from 5uoj -> Chain A from 1zz2 from 5uoj (337 atoms)>

In [12]: mob_chain
Out[12]: <AtomMap: Chain A from 1zz2 -> Chain A from 5uoj from 1zz2 (337 atoms)>

Both AtomMap instances refer to same number of atoms, and their name suggests how they were retrieved.

In addition, we can find out the sequence identity that the matched atoms (residues) share (third item in the tuple):

In [13]: match[2]
Out[13]: 99.40652818991099

The fourth item in the tuple shows the coverage of the matching:

In [14]: match[3]
Out[14]: 98

This is the percentage of matched residues with respect to the longer chain. 5uoj chain A contains 351 resiudes, 96% of it is 337 residues, which is the number of atoms in the returned atom maps.

RMSD and superpose

We calculate the RMSD using calcRMSD() function:

In [15]: calcRMSD(ref_chain, mob_chain).round(2)
Out[15]: 72.94

Let’s find the transformation that minimizes RMSD between these chains using calcTransformation() function:

In [16]: t = calcTransformation(mob_chain, ref_chain)

We apply this transformation to mobile structure (not to mob_chain, to preserve structures integrity).

In [17]: t.apply(mobile)
Out[17]: <AtomGroup: 1zz2 (2872 atoms)>

In [18]: calcRMSD(ref_chain, mob_chain).round(2)
Out[18]: 1.82

Deformation vector

Once matching chains are identified it is straightforward to calculate the deformation vector using calcDeformVector()

In [19]: defvec = calcDeformVector(ref_chain, mob_chain)

In [20]: abs(defvec).round(3)
Out[20]: 33.499

To show how RMSD and deformation vector are related, we can calculate RMSD from the magnitude of the deformation vector:

In [21]: (abs(defvec)**2 / len(ref_chain)) ** 0.5
Out[21]: 1.8248053285019694

Array of numbers for this deformation can be obtained as follows

In [22]: arr = defvec.getArray() # arr is a NumPy array

In [23]: arr.round(2)
Out[23]: array([-1.2 , -0.23, -2.05, ...,  0.9 , -0.17, -0.02])

Following yields the normalized deformation vector

In [24]: defvecnormed = defvec.getNormed()

In [25]: abs(defvecnormed)
Out[25]: 1.0

Compare with ANM modes

Let’s get ANM model for the reference chain using calcANM() (a shorthand function for ANM calculations):

In [26]: anm = calcANM(ref_chain)[0]

Calculate overlap between slowest ANM mode and the deformation vector

In [27]: (anm[0] * defvecnormed).round(2) # used normalized deformation vector
Out[27]: -0.41

We can do this for a set of ANM modes (slowest 6) as follows

In [28]: (array(list(anm[:6])) * defvecnormed).astype(float64).round(2)
Out[28]: array([-0.41,  0.21, -0.48, -0.06, -0.17, -0.08])