# 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: 1p38 and 1zz2

```In [4]: reference = parsePDB('1p38')

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, mobile)
```

`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 1p38 -> Chain A from 1zz2 from 1p38 (337 atoms)>

In [12]: mob_chain
Out[12]: <AtomMap: Chain A from 1zz2 -> Chain A from 1p38 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]: 96
```

This is the percentage of matched residues with respect to the longer chain. 1p38 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.93
```

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.86
```

## 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]: 34.196
```

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.8628014908695476
```

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.11, -0.52, -1.89, ...,  0.85, -0.18,  0.54])
```

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.42
```

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.42, -0.14,  0.49,  0.03, -0.17, -0.1 ])
```