Analyze Conformations

First, necessary imports:

In [1]: from prody import *

In [2]: from numpy import *

In [3]: from pylab import *

In [4]: ion()

In [5]: import os, glob

Parse conformations

Now, let’s read initial and refined conformations:

In [6]: initial = AtomGroup('p38 initial')

In [7]: refined = AtomGroup('p38 refined')
In [8]: for pdb in glob.glob('p38_ensemble/*pdb'):
   ...:     fn = os.path.splitext(os.path.split(pdb)[1])[0]
   ...:     opt = os.path.join('p38_optimize', fn + '.coor')
   ...:     parsePDB(pdb, ag=initial)
   ...:     parsePDB(opt, ag=refined)
   ...: 
In [9]: initial
Out[9]: <AtomGroup: p38 initial (5658 atoms; active #0 of 40 coordsets)>

In [10]: refined
Out[10]: <AtomGroup: p38 refined (5658 atoms; active #0 of 40 coordsets)>

Calculate RMSD change

We can plot RMSD change after refinement as follows:

In [11]: rmsd_ca = []

In [12]: rmsd_all = []

In [13]: initial_ca = initial.ca

In [14]: refined_ca = refined.ca

In [15]: for i in range(initial.numCoordsets()):
   ....:     initial.setACSIndex(i)
   ....:     refined.setACSIndex(i)
   ....:     initial_ca.setACSIndex(i)
   ....:     refined_ca.setACSIndex(i)
   ....:     rmsd_ca.append(calcRMSD(initial_ca, refined_ca))
   ....:     rmsd_all.append(calcRMSD(initial, refined))
   ....: 

In [16]: plot(rmsd_all, label='all');

In [17]: plot(rmsd_ca, label='ca');

In [18]: xlabel('Conformation index');

In [19]: ylabel('RMSD');

In [20]: legend();
../../_images/conformational_sampling_comparison.png

Select a diverse set

To select a diverse set of refined conformations, let’s calculate average RMSD for each conformation to all others:

In [21]: rmsd_mean = []

In [22]: for i in range(refined.numCoordsets()):
   ....:     refined.setACSIndex(i)
   ....:     alignCoordsets(refined)
   ....:     rmsd = calcRMSD(refined)
   ....:     rmsd_mean.append(rmsd.sum() / (len(rmsd) - 1))
   ....: 

In [23]: bar(arange(1, len(rmsd_mean) + 1), rmsd_mean);

In [24]: xlabel('Conformation index');

In [25]: ylabel('Mean RMSD');
../../_images/conformational_sampling_mean_rmsd.png

Let’s select conformations that are 1.2 Å away from other on average:

In [26]: selected = (array(rmsd_mean) >= 1.2).nonzero()[0]

In [27]: selected
Out[27]: 
array([ 0,  1,  2,  3,  4,  6,  7,  9, 12, 13, 15, 16, 17, 18, 20, 21, 23,
       24, 25, 26, 27, 28, 30, 31, 32, 33, 34, 35, 36, 37, 38])

In [28]: selection = refined[selected]

In [29]: selection
Out[29]: <Selection: 'index 0 to 4 6 ... to 28 30 to 38' from p38 refined (31 atoms; active #39 of 40 coordsets)>

Visualization

When you visualize the refined ensemble, you should see something similar to this:

../../_images/p38_sampling.png