This example is continued from Analysis. The aim of this part is to produce graphical comparison of experimental and theoretical data. We will reproduce the plots that was presented in our paper [AB09].

Load data

First, we load data saved in Calculations:

In [1]: from prody import *

In [2]: from pylab import *

In [3]: ion()
In [4]: pca = loadModel('p38_xray.pca.npz')

In [5]: anm = loadModel('1p38.anm.npz')

In [6]: ensemble = loadEnsemble('p38_X-ray.ens.npz')

In [7]: ref_chain = parsePDB('p38_ref_chain.pdb')

PCA - ANM overlap

In previous page, we compared PCA and ANM modes to get some numbers. In this case, we will use plotting functions to make similar comparisons:

In [8]: showOverlapTable(pca[:6], anm[:6]);

# Let's change the title of the figure
In [9]: title('PCA - ANM Overlap Table');

It is also possible to plot overlap of a single mode from one model with multiple modes from another:

In [10]: showOverlap(pca[0], anm);

In [11]: showCumulOverlap(pca[0], anm);
../../_images/ensemble_analysis_xray_overlap.png ../../_images/ensemble_analysis_xray_cumulative_overlap.png

Let’s also plot the cumulative overlap in the same figure:

Square fluctuations

In [12]: showSqFlucts(pca[:3]);

In [13]: showSqFlucts(anm[:3]);
../../_images/ensemble_analysis_xray_pca_sqflucts.png ../../_images/ensemble_analysis_xray_anm_sqflucts.png

Now let’s plot square fluctuations along PCA and ANM modes in the same plot:

In [14]: showScaledSqFlucts(pca[0], anm[2]);

In [15]: legend();
In [16]: showScaledSqFlucts(pca[1], anm[0]);

In [17]: legend();

In above example, ANM modes are scaled to have the same mean as PCA modes. Alternatively, we could plot normalized square fluctuations:

In [18]: showNormedSqFlucts(pca[0], anm[1]);

In [19]: legend();


Now we will project the ensemble onto PC 1 and 2 using showProjection():

In [20]: showProjection(ensemble, pca[:2]);

In [21]: axis([-0.8, 0.8, -0.8, 0.8]);

Now we will do a little more work, and get a colorful picture:

red unbound
blue inhibitor bound
yellow glucoside bound
purple peptide/protein bound
In [22]: color_list = ['blue', 'blue', 'blue', 'blue', 'blue', 'blue', 'blue',
   ....:               'blue', 'purple', 'purple', 'blue', 'blue', 'blue',
   ....:               'blue', 'blue', 'red', 'red', 'red', 'blue', 'blue',
   ....:               'blue', 'blue', 'blue','blue', 'blue', 'blue', 'blue',
   ....:               'blue', 'red', 'blue', 'blue','blue', 'blue', 'blue',
   ....:               'blue', 'blue', 'blue', 'blue', 'blue', 'blue', 'yellow',
   ....:               'yellow', 'yellow', 'yellow', 'blue', 'blue','blue',
   ....:               'blue', 'blue', 'blue', 'yellow', 'purple', 'purple',
   ....:               'blue', 'yellow', 'yellow', 'yellow', 'blue', 'yellow',
   ....:               'yellow', 'blue', 'blue', 'blue', 'blue', 'blue', 'blue',
   ....:               'blue', 'blue', 'blue', 'blue', 'blue', 'blue', 'blue',
   ....:               'blue', 'purple']

In [23]: color2label = {'red': 'Unbound', 'blue': 'Inhibitor bound',
   ....:                'yellow': 'Glucoside bound',
   ....:                'purple': 'Peptide/protein bound'}

In [24]: label_list = [color2label[color] for color in color_list]

In [25]: showProjection(ensemble, pca[:2], color=color_list,
   ....:                label=label_list);

In [26]: axis([-0.8, 0.8, -0.8, 0.8]);

In [27]: legend();

Now let’s project conformations onto 3d principal space and label conformations using text keyword argument and PDBEnsemble.getLabels() method:

In [28]: showProjection(ensemble, pca[:3], color=color_list, label=label_list,
   ....:                text=ensemble.getLabels(), fontsize=10);

The figure with all conformation labels is crowded, but in an interactive session you can zoom in and out to make text readable.


Finally, we will make a cross-projection plot using showCrossProjection(). We will pass scale='y' argument, which will scale the width of the projection along ANM mode:

In [29]: showCrossProjection(ensemble, pca[0], anm[2], scale="y",
   ....:                      color=color_list, label=label_list);

In [30]: plot([-0.8, 0.8], [-0.8, 0.8], 'k');

In [31]: axis([-0.8, 0.8, -0.8, 0.8]);

In [32]: legend(loc='upper left');
In [33]: showCrossProjection(ensemble, pca[1], anm[0], scale="y",
   ....:                     color=color_list, label=label_list);

In [34]: plot([-0.8, 0.8], [-0.8, 0.8], 'k');

In [35]: axis([-0.8, 0.8, -0.8, 0.8]);

It is also possible to find the correlation between these projections:

In [36]: pca_coords, anm_coords = calcCrossProjection(ensemble, pca[0], anm[2])

In [37]: print(np.corrcoef(pca_coords, anm_coords))
[[ 1.         -0.94621454]
 [-0.94621454  1.        ]]

This is going to print 0.95 for PC 1 and ANM mode 2 pair.

Finally, it is also possible to label conformations in cross projection plots too:

In [38]: showCrossProjection(ensemble, pca[1], anm[0], scale="y",
   ....:     color=color_list, label=label_list, text=ensemble.getLabels(),
   ....:     fontsize=10);

In [39]: plot([-0.8, 0.8], [-0.8, 0.8], 'k');

In [40]: axis([-0.8, 0.8, -0.8, 0.8]);