Markov state models
Current Models: Src https://github.com/choderalab/MSMs/blob/master/src/README.md
Markov state models
License: GNU General Public License v2.0
Markov state models
Current Models: Src https://github.com/choderalab/MSMs/blob/master/src/README.md
I guess we need to check this a bit more thoroughly in general re: ensembler models, but currently have a basic question:
If I plot the first frame of run0-clone0.h5 is that equivalent to plotting the starting structure?
Maybe there is a reason I am overlooking...
John relayed Frank's suggestion to use a subset of inter-residue contacts as input features for MSM-building: we would like to use just "the subset of contacts that change during a simulation." I've been looking into this a bit for modeling the data from 10468 (Abl).
tl;dr -- promising!
My biggest question is:
How should we choose the contact threshold?
To determine whether a given inter-residue contact "changes," we binarize the inter-residue distance as either > or <= a threshold distance, then select any residue pairs that aren't all 0s or all 1s throughout our simulations. PyEMMA has a default setting for this of 5 angstroms. However, the subset of contacts that will appear to be "changing" clearly depends on what threshold you use:
Note that the set of selected contacts at a threshold x+epsilon is not simply a superset of the contacts selected at a lower threshold x:
Here I've just set the threshold=1.0nm since that induces an acceptably small number of features (~4000) that seem to represent the extent of the protein, but I don't know how to do this in a principled way.
Comparing feature sets:
A feature set that does a good job at linearly representing the slow relaxation processes of a system will have a high GMRQ, i.e. a large sum of leading eigenvalues in the tICA solution.
I compared 3 feature sets in terms of their GMRQ for fixed-rank approximation (m=100) of the lag_time=100 propagator. I also computed the maximum cumulative sum of leading eigenvalues.
Feature set 1: Backbone torsions (all available: ‘phi’, ‘psi’, ‘omega’, ‘chi1’, ‘chi2’, ‘chi3’, ‘chi4’)
The sum of the first 100 eigenvalues = 55.7.
The sum of the first 542 eigenvalues = 76.3.
Feature set 2: Inter-residue distances between contacts that change (threshold=1.0nm)
(evaluated on a subset of 10 trajectories)
The sum of the first 100 eigenvalues = 61.9.
The sum of the first 1999 eigenvalues = 92.2.
Feature set 3: Kabsch-Sander hydrogen-bond energies
The sum of the first 100 eigenvalues = 40.4.
The sum of the first 468 eigenvalues = 48.5.
Visualizing the results
From the eigenvalue plots above we can see that the first 2 tICs represent only a tiny fraction of the kinetic behavior of the kinase, for any of these choices of input features-- projecting onto the first 2 tICs is a very lossy representation of the tICA model. Can we create a more trustworthy 2D picture?
Here I'm using the first 100 tICs as inputs into t-SNE, which produces a nonlinear 2D embedding that preserves local neighborhoods at the cost of introducing global distortions. This is useful if we want a picture of the cluster structure of a dataset, but means we can't read into the x and y axes as "components." Each dot represents a simulation frame, dots are colored by the order in which they were collected, and time-adjacent frames are connected by a thin grey line. Areas that are visited more than once should have a mix of colors in them. Areas that are color-homogenuous appear to have been visited during just one contiguous chunk of trajectory.
Note that the above pictures were computed using just 10 out of the ~500 trajectories available.
If we recompute the "Contact map --> top-100 tICs --> t-SNE" projection for all available simulation data, then we get the following picture:
This picture suggests that the majority of the available data samples transitions between large relevant states (the dense multi-colored regions in the map), and there may be a large number of negligibly occupied microstates that are visited ~once (the little "sprinkles" on the periphery). Caveat: It's possible that t-SNE is over-partitioning the dataset (i.e. hallucinating clusters)-- although I ran a quick sanity check that the k-nearest neighbor graph (for k=30) is mostly identical in the 2D embedding and in the 100D tICA model, there are some additional sanity checks I can run before we interpret these pictures too much.
Other comments
MDTraj offers several schemes for computing inter-residue distances. I've been using the scheme='ca'
(which returns the distance between their alpha carbons), since it's significantly cheaper than scheme='closest-heavy'
(which I presume is still cheaper than scheme='closest'
), but I don't know which of the available schemes is preferred here.
So I did simple experiment where I grabbed random selections of atom pairs, then calculated the tica eigenvectors. They're pretty robust even with a small subset of atom pairs. This suggests that we could probably do something very affordable here and get very reproducible results.
Here is the driver code:
import itertools
import numpy as np
import mdtraj as md
import msmbuilder as msmb
trj0 = md.load("./system.subset.pdb")
trajectories = [md.load("./Trajectories/trj%d.h5" % i) for i in range(5)]
atom_indices = np.arange(trj0.n_atoms)
pair_indices = np.array(list(itertools.combinations(atom_indices, 2)))
n_choose = 1000
n_trials = 10
for i in range(n_trials):
pair_indices_subset = pair_indices[np.random.choice(len(pair_indices), n_choose, replace=False)]
metric = msmb.metrics.AtomPairs(atom_pairs=pair_indices_subset)
tica = msmb.reduce.tICA(1, prep_metric=metric)
for trj in trajectories:
tica.train(trajectory=trj)
tica.solve()
print tica.vals[0:5]
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.