The MDAnalysis.analysis.hole module in MDAnalysis can analyze trajectories of, say ion channels, and associate the per-frame pore radius profiles $R_\rho(\zeta)$ with an order parameter $\rho$ for each frame of the trajectory. It uses Oliver Smart's HOLE program.
In this example notebook we use a synthetic and short trajectory of the gramicidin A cation pore (gA). The trajectory was generated from the first elastic network mode of gA, generated by the elNémo server, using the gA structure 1GRM. (The trajectory is included as a test case trajectory MULTIPDB_HOLE
with the MDAnalysis unit tests and the X-ray structure is included as PDB_HOLE
.)
import MDAnalysis as mda
from MDAnalysis.analysis.hole import HOLEtraj
from MDAnalysis.analysis.rms import RMSD
from MDAnalysis.tests.datafiles import PDB_HOLE, MULTIPDB_HOLE
mda.start_logging()
/Volumes/Data/oliver/Biop/Projects/Methods/MDAnalysis/mdanalysis/testsuite/MDAnalysisTests/__init__.py:109: UserWarning: This call to matplotlib.use() has no effect because the backend has already been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot, or matplotlib.backends is imported for the first time. The backend was *originally* set to 'module://ipykernel.pylab.backend_inline' by the following code: File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/runpy.py", line 193, in _run_module_as_main "__main__", mod_spec) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/runpy.py", line 85, in _run_code exec(code, run_globals) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel_launcher.py", line 16, in <module> app.launch_new_instance() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/traitlets/config/application.py", line 658, in launch_instance app.start() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelapp.py", line 486, in start self.io_loop.start() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/ioloop.py", line 888, in start handler_func(fd_obj, events) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/stack_context.py", line 277, in null_wrapper return fn(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 450, in _handle_events self._handle_recv() File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 480, in _handle_recv self._run_callback(callback, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py", line 432, in _run_callback callback(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/tornado/stack_context.py", line 277, in null_wrapper return fn(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 283, in dispatcher return self.dispatch_shell(stream, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 233, in dispatch_shell handler(stream, idents, msg) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/kernelbase.py", line 399, in execute_request user_expressions, allow_stdin) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/ipkernel.py", line 208, in do_execute res = shell.run_cell(code, store_history=store_history, silent=silent) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/ipykernel/zmqshell.py", line 537, in run_cell return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2728, in run_cell interactivity=interactivity, compiler=compiler, result=result) File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2850, in run_ast_nodes if self.run_code(code, result): File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2910, in run_code exec(code_obj, self.user_global_ns, self.user_ns) File "<ipython-input-1-9dc852c9a767>", line 2, in <module> from MDAnalysis.analysis.hole import HOLEtraj File "/Volumes/Data/oliver/Biop/Projects/Methods/MDAnalysis/mdanalysis/package/MDAnalysis/analysis/hole.py", line 265, in <module> import matplotlib.pyplot as plt File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/matplotlib/pyplot.py", line 69, in <module> from matplotlib.backends import pylab_setup File "/Users/oliver/anaconda3/envs/mda_develop/lib/python3.6/site-packages/matplotlib/backends/__init__.py", line 14, in <module> line for line in traceback.format_stack() matplotlib.use('agg') MDAnalysis : INFO MDAnalysis 0.17.1-dev STARTED logging to 'MDAnalysis.log'
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')
We take the X-ray structure as the reference and we want to analyze the ENM-mode trajectory so we load them into two separate universes:
ref = mda.Universe(PDB_HOLE) # reference structure (needs MDAnalysis 0.17.1)
u = mda.Universe(MULTIPDB_HOLE) # trajectory
We calculate the all-atom RMSD from the reference structure as the order parameter $\rho$:
# calculate RMSD
R = RMSD(u, reference=ref, select="protein", weights="mass")
R.run()
MDAnalysis.analysis.rmsd: INFO RMS calculation for 260 atoms. MDAnalysis.analysis.base: INFO Starting preparation /Volumes/Data/oliver/Biop/Projects/Methods/MDAnalysis/mdanalysis/package/MDAnalysis/coordinates/base.py:821: UserWarning: Reader has no dt information, set to 1.0 ps warnings.warn("Reader has no dt information, set to 1.0 ps") MDAnalysis.analysis.base: INFO Finishing up
<MDAnalysis.analysis.rms.RMSD at 0x1148c5198>
The timeseries of the orderparameter is R.rmsd
(an N x 3 array, with each row (frame, time, RMSD)
):
frame, _, rho = R.rmsd.transpose()
For the trajectory from the ENM, time is meaningless so we just plot over frame number.
ax = plt.subplot(111)
ax.plot(frame, rho, '.-')
ax.set_xlabel("frame")
ax.set_ylabel(r"RMSD $\rho$ ($\AA$)");
Note that the ENM trajectory is symmetric to frame 5 (the original X-ray structure) and the mode deforms the structure so that the RMSD increases nearly linearly.
# HOLE analysis with order parameters
H = HOLEtraj(u, orderparameters=R.rmsd[:,2],
executable="~/hole2/exe/hole")
H.run()
MDAnalysis.analysis.hole: INFO HOLE analysis frame 0 (orderparameter 6.10501) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpgb079n8a.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpgb079n8a.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 1 (orderparameter 4.88398) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpcchs8ex1.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpcchs8ex1.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 2 (orderparameter 3.66304) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpj2x_ll39.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpj2x_ll39.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 3 (orderparameter 2.44202) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpau6yv34d.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpau6yv34d.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 4 (orderparameter 1.22101) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp65exmo_2.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp65exmo_2.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 5 (orderparameter 1.67286e-07) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpjan1v538.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpjan1v538.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 6 (orderparameter 1.221) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp7ywiy3t4.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp7ywiy3t4.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 7 (orderparameter 2.44202) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpb1_cfq5k.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpb1_cfq5k.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 8 (orderparameter 3.66303) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpqrz1nuyg.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpqrz1nuyg.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 9 (orderparameter 4.88398) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp3rp40h_v.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmp3rp40h_v.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames MDAnalysis.analysis.hole: INFO HOLE analysis frame 10 (orderparameter 6.10502) MDAnalysis.analysis.hole: INFO Setting up HOLE analysis for '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpa57uguid.pdb' MDAnalysis.analysis.hole: INFO Using radius file 'simple2.rad' MDAnalysis.analysis.hole: INFO HOLE will guess CPOINT MDAnalysis.analysis.hole: INFO HOLE will guess CVECT MDAnalysis.analysis.hole: INFO Starting HOLE on '/var/folders/ds/z7mbf7nx1gn31308g31t9tcc0000gp/T/tmpa57uguid.pdb' (trajectory: None) MDAnalysis.analysis.hole: INFO HOLE finished: output file 'hole.out' MDAnalysis.analysis.hole: INFO Collecting HOLE profiles for run with id 1 MDAnalysis.analysis.hole: INFO Run 1: Reading 1 HOLE profiles from 'hole.out' MDAnalysis.analysis.hole: INFO Collected HOLE radius profiles for 1 frames
As a quick check, plot the pore profiles $R_\rho(\zeta)$:
H.plot3D(rmax=2.5)
<matplotlib.axes._subplots.Axes3DSubplot at 0x116e53e48>
The 2D plot clearly shows how the pore profile opens up with increasing order parameter (i.e., increasing strength of the deformation along the mode):
ax = H.plot()
ax.set_ylim(1, 3)
ax.set_xlim(-15, 35)
plt.tight_layout()
Is the motion one that could possibly gate the pore, i.e., change the pore radius to a radius that makes more ions pass through? (Note: gA is already open so scientifically this is not a great question...)
We collect the minimum pore radius from each profile
$$ r(\rho) = \min_\zeta R_\rho(\zeta) $$and store it together with the corresponding order parameter:
r_rho = np.array([[rho, profile.radius.min()] for rho, profile in H])
The plot shows quantitatively that the pore opens up.
ax = plt.subplot(111)
ax.plot(r_rho[:, 0], r_rho[:, 1], lw=2)
ax.set_xlabel(r"order parameter RMSD $\rho$ ($\AA$)")
ax.set_ylabel(r"minimum HOLE pore radius $r$ ($\AA$)");
Comparison to the stacked pore profile plots it is clear that the constriction itself changes position. We can calculate the constriction point along the pore axis $$ \zeta_0(\rho) = \text{arg} \min_\zeta R_\rho(\zeta) $$
zeta0_rho = np.array([[rho, profile.rxncoord[np.argmin(profile.radius)]]
for rho, profile in H])
The plot clearly shows that for $\rho < 2.5$ A, the constriction is near the end but for larger distortions it is close to the center of the pore.
ax = plt.subplot(111)
ax.plot(zeta0_rho[:, 0], np.abs(zeta0_rho[:, 1]), lw=2)
ax.set_xlabel(r"order parameter RMSD $\rho$ ($\AA$)")
ax.set_ylabel(r"position of constriction $|\zeta_0|$ ($\AA$)");
(For the plot I symmetrized $\zeta_0$ using np.abs(zeta0)
because for small $\rho$, either of the two constriction sites near the ends is closed. However, the HOLE search procedure will typically not exactly reproduce the symmetry. For real work one might want to average over the two lowest values near the inside and the outside.)