Normal Mode Algebra

This part shows how to use some handy features of Mode objects.

ANM Calculations

We will compare modes from two ANMs for the same protein, but everything applies to comparison of ANMs and PCAs (as long as they contain same number of atoms).

Let’s get started by getting ANM models for two related protein structures:

In [1]: from prody import *

In [2]: str1 = parsePDB('5uoj')

In [3]: str2 = parsePDB('1r39')

Find and align matching chains

In [4]: matches = matchChains(str1.protein, str2.protein)

In [5]: match = matches[0]

In [6]: ch1 = match[0]

In [7]: ch2 = match[1]

Minimize RMSD by superposing ch2 onto ch1:

In [8]: ch2, t = superpose(ch2, ch1)  # t is transformation, already applied to ch2

In [9]: calcRMSD(ch1, ch2)
Out[9]: 0.8271517534023315

Get ANM models for each chain

In [10]: anm1, ch1 = calcANM(ch1)

In [11]: anm2, ch2 = calcANM(ch2)

In [12]: anm1[0]
Out[12]: <Mode: 1 from ANM 5uoj>

Let’s rename these ANM instances, so that they print short:

In [13]: anm1.setTitle('5uoj_anm')

In [14]: anm2.setTitle('1r39_anm')

This is how they print now:

In [15]: anm1[0]
Out[15]: <Mode: 1 from ANM 5uoj_anm>

In [16]: anm2[0]
Out[16]: <Mode: 1 from ANM 1r39_anm>

Calculate overlap

We need NumPy in this part:

In [17]: from numpy import *

Multiplication of two Mode instances returns dot product of their eigenvectors. This dot product is the overlap or cosine correlation between modes.

Let’s calculate overlap for slowest modes:

In [18]: overlap = anm1[0] * anm2[0]

In [19]: overlap
Out[19]: -0.996773780512587

This shows that the overlap between these two modes is 0.98, which is not surprising since ANM modes come from structures of the same protein.

To compare multiple modes, convert a list of modes to a numpy.array():

In [20]: array(list(anm1[:3])) * array(list(anm2[:3]))
Out[20]: 
array([-0.996773780512587, -0.996017575839516, 0.9963391906808239],
      dtype=object)

This shows that slowest three modes are almost identical.

We could also generate a matrix of overlaps using numpy.outer():

In [21]: outer_product = outer(array(list(anm1[:3])), array(list(anm2[:3])))

In [22]: outer_product
Out[22]: 
array([[-0.996773780512587, -0.055565379414888975, 0.014704351225430195],
       [0.055699990428254346, -0.996017575839516, -0.02819580848734082],
       [0.015053866160526948, -0.028730152467655548, 0.9963391906808239]],
      dtype=object)

This could also be printed in a pretty table format using printOverlapTable():

In [23]: printOverlapTable(anm1[:3], anm2[:3])
Overlap Table
                      ANM 1r39_anm
                    #1     #2     #3
ANM 5uoj_anm #1   -1.00  -0.06  +0.01
ANM 5uoj_anm #2   +0.06  -1.00  -0.03
ANM 5uoj_anm #3   +0.02  -0.03  +1.00

Scaling

Mode instances can be scaled, but after this operation they will become Vector instances:

In [24]: anm1[0] * 10
Out[24]: <Vector: (Mode 1 from ANM 5uoj_anm)*10>

Linear combination

It is also possible to linearly combine normal modes:

In [25]: anm1[0] * 3 + anm1[1] + anm1[2] * 2
Out[25]: <Vector: (((Mode 1 from ANM 5uoj_anm)*3) + (Mode 2 from ANM 5uoj_anm)) + ((Mode 3 from ANM 5uoj_anm)*2)>

Or, we could use eigenvalues for linear combination:

In [26]: lincomb = anm1[0] * anm1[0].getEigval() + anm1[1] * anm1[1].getEigval()

It is the name of the Vector instance that keeps track of operations.

In [27]: lincomb.getTitle()
Out[27]: '((Mode 1 from ANM 5uoj_anm)*0.135300972474) + ((Mode 2 from ANM 5uoj_anm)*0.20786277222)'

Approximate a deformation vector

Let’s get the deformation vector between ch1 and ch2:

In [28]: defvec = calcDeformVector(ch1, ch2)

In [29]: abs(defvec)
Out[29]: 15.251924726841118

Let’s see how deformation projects onto ANM modes:

In [30]: array(list(anm1[:3])) * defvec
Out[30]: 
array([-5.205741393277132, 3.3492257559027094, -2.6041269714867803],
      dtype=object)

We can use these numbers to combine ANM modes:

In [31]: approximate_defvec = sum((array(list(anm1[:3])) * defvec) *
   ....:                          array(list(anm1[:3])))
   ....: 

In [32]: approximate_defvec
Out[32]: <Vector: ((-5.20574139328*(Mode 1 from ANM 5uoj_anm)) + (3.3492257559*(Mode 2 from ANM 5uoj_anm))) + (-2.60412697149*(Mode 3 from ANM 5uoj_anm))>

Let’s deform 1r39 chain along this approximate deformation vector and see how RMSD changes:

In [33]: ch2.setCoords(ch2.getCoords() - approximate_defvec.getArrayNx3())

In [34]: calcRMSD(ch1, ch2)
Out[34]: 0.7426555570303122

RMSD decreases from 0.89 A to 0.82 A.