New Hybrid Modules

[jamesmkrieger]

This builds upon ClustENM to provide implementations for other hybrid methods.

It is tested with a new hybrid method built around adaptive ANM, which works successfully. The following test code produces a transition from 1akeA to 4akeA as confirmed by the plots below and visualisation of the output (not shown):

import matplotlib.pyplot as plt
plt.ion()

from prody import *

r = parsePDB('1akeA')
t = parsePDB('4akeA')

self = AdaptiveHybrid('test')
self.setAtoms(r, t)
self.run(n_gens=20)

plt.figure()
rmsd_from_start = self.getRMSDs()
plt.plot(rmsd_from_start)
plt.savefig('ake_adaptiveHybrid_RMSDs_A')
plt.close()

plt.figure()
rmsd_from_start = self.getRMSDsB()
plt.plot(rmsd_from_start)
plt.savefig('ake_adaptiveHybrid_RMSDs_B')
plt.close()

[jamesmkrieger]

We can now run ClustENM and then filter to only keep conformers that have the same or better CCC to a target map than a reference structure (this follows MDeNM-EMFit):

import matplotlib.pyplot as plt
plt.ion()

from prody import *

from TEMPy.protein.structure_blurrer import StructureBlurrer

o = parsePDB('1gruA')
c = parsePDB('1gr5A')

proc = 2

if proc == 1:
    ens = ClustENM('test')
    ens.setAtoms(o)
    ens.run(n_gens=3, maxclust=(10, 30, 50))

    saveEnsemble(ens, '1gruA_clustenm.ens.npz')
    writePDB('1gruA_clustenm.pdb', ens)
else:
    ens = loadEnsemble('1gruA_clustenm.ens.npz')

    atoms = ens.getAtoms()
    amap = alignChains(atoms, o.ca)[0]
    ens.setAtoms(amap)

    mean_conf = amap.copy()
    mean_conf.setCoords(ens.getCoords())

    plt.figure()
    rmsd_from_start = ens.getRMSDs()
    plt.plot(rmsd_from_start)
    plt.savefig('1gruA_clustenm_RMSDs')
    plt.close()

    cmap = alignChains(c, mean_conf)[0]
    c_alg, T = superpose(cmap, mean_conf, weights=cmap.getFlags('mapped'))

    rmsd_from_cmap = [calcRMSD(coordset, cmap, weights=cmap.getFlags('mapped'))
                     for coordset in ens.getCoordsets()]
    plt.figure()
    plt.plot(rmsd_from_cmap)
    plt.savefig('1gruA_clustenm_RMSDs_from_1gr5A')
    plt.close()

    blurrer = StructureBlurrer()
    c_Structure = c.toTEMPyStructure()
    c_blurred = blurrer.gaussian_blur(c_Structure, 5,
                                    filename=c.getTitle()+'.mrc')
    c_blurred.write_to_MRC_file(c_blurred.filename)

    defvec_to_cmap = calcDeformVector(mean_conf, cmap,
                                      weights=cmap.getFlags('mapped'))
    writeNMD('1gruA_to_1gr5A.defvec.nmd', defvec_to_cmap, mean_conf)

    omap = alignChains(o, mean_conf)[0]
    o_alg, T = superpose(omap, mean_conf, weights=cmap.getFlags('mapped'))

    rmsd_from_omap = [calcRMSD(coordset, omap, weights=omap.getFlags('mapped'))
                     for coordset in ens.getCoordsets()]
    plt.figure()
    plt.plot(rmsd_from_omap)
    plt.savefig('1gruA_clustenm_RMSDs_from_1gruA')
    plt.close()

    pca = PCA('1gruA ClustENM')
    pca.buildCovariance(ens)
    pca.calcModes()

    saveModel(pca, '1gruA_ClustENM.pca.npz')
    writeNMD('1gruA_ClustENM.pca.nmd', pca, mean_conf)

    plt.figure()
    showProjection(ens, pca[:2], label='ClustENM')
    showProjection(ens[0], pca[:2], c='c', markersize=10, label='1gruA (start)')
    showProjection(defvec_to_cmap, pca[:2], c='r', markersize=10, label='1gr5A')
    plt.legend()
    plt.savefig('1gruA_clustenm_pca_proj2d')
    xlim = plt.xlim()
    ylim = plt.ylim()
    plt.close()

    ens.setAtoms(atoms)
    ens.filterConformers(ref=o, target=c_blurred)

    ens.setAtoms(amap)
    plt.figure()
    showProjection(ens, pca[:2], label='ClustENM')
    showProjection(ens[0], pca[:2], c='c', markersize=10, label='1gruA (start)')
    showProjection(defvec_to_cmap, pca[:2], c='r', markersize=10, label='1gr5A')
    plt.legend()
    plt.xlim(xlim)
    plt.ylim(ylim)
    plt.savefig('1gruA_clustenm_filtered_pca_proj2d')
    plt.close()

    ens.setAtoms(atoms)

This can be seen on the PCA projections:

[jamesmkrieger]

It was previously not adding the conformers to build the ensemble correctly because negative indexing needs to be different depending on whether there are an even or odd number of them.

[jamesmkrieger]

It could be a good idea to sub/super class AdaptiveHybrid and CoMD so we don't need to keep changing both of them.