deeptime-ml / deeptime

Python library for analysis of time series data including dimensionality reduction, clustering, and Markov model estimation
https://deeptime-ml.github.io/
GNU Lesser General Public License v3.0
746 stars 82 forks source link

Extracting pdb structures #298

Open GSINGH006 opened 1 month ago

GSINGH006 commented 1 month ago

Hi, is there a method in Deeptime similar to Pyemma's " pyemma.coordinates.save_traj" which can be used to extract pdb structures representing macrostates after coarse-graining using PCCA+

clonker commented 1 month ago

There is no direct MD dependency in deeptime so in that sense, no, but you can use deeptime to estimate your models and extract the indices, put these into pyemma :)

GSINGH006 commented 1 month ago

Hi @clonker I am using the following code for the same

topfile = "/home/lab/Downloads/penta_traj/pentapeptide-impl-solv.pdb" traj_list=["/home/lab/Downloads/penta_traj/pentapeptide-00-500ns-impl-solv.xtc","/home/lab/Downloads/penta_traj/pentapeptide-01-500ns-impl-solv.xtc"] import pyemma import deeptime as dt import numpy as np import matplotlib.pyplot as plt import mdshare plt.matplotlib.rcParams.update({'font.size': 16}) import numpy as np feat1= pyemma.coordinates.featurizer(topfile) tor=feat1.add_backbone_torsions(cossin=True) y= pyemma.coordinates.load(traj_list, feat1) from deeptime.decomposition import TICA tica_estimator = TICA(lagtime=1,dim=2) tica = tica_estimator.fit(y).fetch_model() tics=tica.transform(y) from tqdm.notebook import tqdm from deeptime.clustering import KMeans kmeans_estimator = KMeans(n_clusters=30, progress=tqdm) clustering = kmeans_estimator.fit(np.concatenate(tics)[::1]).fetch_model() fig, ax = plt.subplots() ax.plot(clustering.cluster_centers.T, 'k+') pyemma.plots.plot_free_energy(np.concatenate(tics).T, ax=ax) plt.savefig('energy_dist.png', dpi=300) dtrajs = [] for projected_trajectory in tics: dtrajs.append(clustering.transform(projected_trajectory)) **msm = MaximumLikelihoodMSM(lagtime=50).fitfetch(dtrajs) nstates = 2 cg = msm.pcca(nstates) for i, s in enumerate(cg.sets): print('π{} = {:f}'.format(i + 1, msm.stationary_distribution[s].sum())) pcca_dist1=cg.metastable_distributions indices = dt.markov.sample.compute_index_states(dtrajs) ind=dt.markov.sample.indices_by_distribution(indices,distributions=pcca_dist1,nsample=1) pyemma.coordinates.save_traj(traj_list, ind[0],'/home/lab/samples_pcca1.xtc' , top=topfile) Is this the correct way to extract the indices and use them for getting the representative structures? what does nsamples in this function exactly represent ? Thank you