molPX Di-Ala example

 
Guillermo Pérez-Hernández  guille.perez@fu-berlin.de 

In this notebook we will be using a trajectory of Di-Ala-peptide to easily identify conformations in the Ramachandran plot.


In [4]:
from os.path import exists
import molpx
from matplotlib import pylab as plt
%matplotlib ipympl

import pyemma
import numpy as np

Start from files on disk


In [5]:
top = molpx._molpxdir(join='notebooks/data/ala2.pdb')
# What data do we  have?
if exists('/group/ag_cmb/scratch/gph82/Di-Ala-nbdata/ala2.dcd'):
    MD_trajfiles = ['/group/ag_cmb/scratch/gph82/Di-Ala-nbdata/ala2.dcd'] #long trajectory
elif exists('/home/guille/ala2.dcd'):
    MD_trajfiles = ['/home/guille/ala2.dcd'] # extra for Stralsund
else:
    MD_trajfiles = [molpx._molpxdir(join='notebooks/data/ala2.mini.xtc')] #short trajectory

Featurize to Ramachandran $(\phi,\psi)$-pairs with PyEMMA


In [6]:
feat = pyemma.coordinates.featurizer(top)
feat.add_backbone_torsions()
src  = pyemma.coordinates.source(MD_trajfiles, features=feat)
Y = src.get_output()


Visualize a FES and the molecular structures behind it

Execute the following cell and click either on the FES or on the slidebar


In [7]:
mpx_widget_box = molpx.visualize.FES(MD_trajfiles,
                                     top,
                                     Y,
                                     #proj_idxs=[1],
                                     nbins=50,       
                                     proj_labels=['$\phi$',
                                                  '$\psi$'],
                                     atom_selection="symbol != H",
                                     #n_overlays=5,
                                     #sticky=True, 
                                     #color_list='random'
                                    )
mpx_widget_box


Visualize trajectories, FES and molecular structures


In [8]:
from molpx import visualize, _linkutils
from imp import reload
reload(visualize)
reload(_linkutils)
mpl_wdg_box = molpx.visualize.traj(MD_trajfiles,  
                                   top, 
                                   Y,
                                   plot_FES = True,
                                   #dt = dt*1e-6, tunits='ms',
                                   max_frames=10000,
                                   proj_idxs=[0, 1],
                                   panel_height=2, 
                                   proj_labels=['$\phi$', '$\psi$']
                          )
mpl_wdg_box


Paths samples along the different projections (=axis)


In [10]:
paths_dict, idata = molpx.generate.projection_paths(MD_trajfiles, 
                                                    top, 
                                                    Y, 
                                                    n_points=50,
                                                    proj_idxs=[0,1],
                                                    n_projs=3,
                                                    proj_dim = 3, 
                                                    verbose=False, 
                                        )

In [11]:
# Choose the coordinate and the type of path
coord = 1
path_type = 'min_rmsd'
#path_type = 'min_disp'
igeom = paths_dict[coord][path_type]["geom"]
ipath = paths_dict[coord][path_type]["proj"]

# Choose the proj_idxs for the path and the FES 
# to be shown
proj_idxs = [0,1]

In [12]:
plt.ioff() # Turn of interactive plotting
plt.figure(figsize=(4,4))
h, (x,y) = np.histogramdd(np.vstack(Y)[:,proj_idxs], bins=50)
plt.contourf(x[:-1], y[:-1], -np.log(h.T), alpha=.50)
plt.ion()

linked_NGL_wdg, linked_ax_widget = molpx.visualize.sample(ipath[:,proj_idxs], 
                                                          igeom,  
                                                          plt.gca(), 
                                                          clear_lines=True,
                                                          n_smooth = 2, 
                                                          plot_path=True,  
                                                          #radius=True,
                                                         )
linked_NGL_wdg._set_size('4in', '4in')
from ipywidgets import HBox
HBox([linked_NGL_wdg, linked_ax_widget.canvas])


/home/guille/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log
  after removing the cwd from sys.path.

Let's do TICA and try to look a the correlations in a TICA analysis


In [13]:
feat = pyemma.coordinates.featurizer(top)
#feat.add_backbone_torsions(cossin=True)
feat.add_distances(feat.topology.select('symbol != H'))
src  = pyemma.coordinates.source(MD_trajfiles, features=feat)
tica = pyemma.coordinates.tica(src, lag=np.int(src.trajectory_lengths()/3000))
Y_tica = tica.get_output()


01-11-17 21:01:44 pyemma.coordinates.data.featurization.featurizer.MDFeaturizer[7] WARNING  The 1D arrays input for add_distances() have been sorted, and index duplicates have been eliminated.
Check the output of describe() to see the actual order of the features

In [15]:
mpx_wdg_box = molpx.visualize.FES(MD_trajfiles,
                                  top,                                                       
                                  Y_tica,                                                       
                                  n_overlays=5,
                                  atom_selection='backbone',
                                  #sticky=True,
                                  #color_list='rand'
          )
mpx_wdg_box



In [17]:
mpx_wdg_box = molpx.visualize.traj(MD_trajfiles,     
                                  top,                                                                                                                              
                                  Y_tica,                                        
                                  plot_FES = True,                                        
                                  #dt = dt*1e-6, tunits='ms',                                           
                                  max_frames=10000,
                                  #proj_idxs=[0,1],
                                  panel_height=2,                                              
                                  projection=tica
                                  )
mpx_wdg_box



In [20]:
paths_dict, idata = molpx.generate.projection_paths(MD_trajfiles, 
                                                    top, 
                                                    Y_tica, 
                                                    n_points=50,
                                                    proj_idxs=[0,1],
                                                    n_projs=2,
                                                    proj_dim = 2, 
                                                    verbose=False, 
                                        )

In [21]:
# Choose the coordinate and the type of path
coord = 0
path_type = 'min_rmsd'
#path_type = 'min_disp'
igeom = paths_dict[coord][path_type]["geom"]
ipath = paths_dict[coord][path_type]["proj"]

# Choose the proj_idxs for the path and the FES 
# to be shown
proj_idxs = [0,1]

In [22]:
plt.figure(figsize=(4,4))
h, (x,y) = np.histogramdd(np.vstack(Y_tica)[:,proj_idxs], bins=50)
plt.contourf(x[:-1], y[:-1], -np.log(h.T), alpha=.50)

linked_wdg, axes_widget = molpx.visualize.sample(ipath[:,proj_idxs], 
                                    igeom,
                                    plt.gca(), 
                                    clear_lines=True,
                                    n_smooth = 1, 
                                    plot_path=True,   
                                   )
# You can even choose to add the correlations a posteriori                                         
molpx.visualize.correlations(tica, widget=linked_wdg, proj_idxs=0)
linked_wdg.center_view()
linked_wdg


/home/guille/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log
  This is separate from the ipykernel package so we can avoid doing imports until
WARNING:root:DEPRECATED: Please use 'center' method

In [ ]: