Mercurial > repos > chemteam > mdanalysis_extract_rmsd

comparison end-to-end.py @ 0:743bd6aa3c7a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 1b23e024af45cc0999d9142d07de6897d4189ec2"
author chemteam
date Mon, 24 Aug 2020 16:41:41 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:743bd6aa3c7a
1 #!/usr/bin/env python
2
3 import argparse
4 import itertools
5 import sys
6
7 import MDAnalysis as mda
8
9 import matplotlib
10 import matplotlib.pyplot as plt
11
12 import numpy as np
13 import numpy.linalg
14
15 matplotlib.use('Agg') # noqa
16
17
18 def parse_command_line(argv):
19 parser = argparse.ArgumentParser()
20 parser.add_argument('--itraj', help='input traj')
21 parser.add_argument('--istr', help='input str')
22 parser.add_argument('--itrajext', help='input traj ext')
23 parser.add_argument('--istrext', help='input str ext')
24 parser.add_argument('--isegid1', help='segid 1')
25 parser.add_argument('--ilabel', help='plot label')
26 parser.add_argument('--ititle1', help='plot title')
27 parser.add_argument('--output1', help='output1 - timeseries')
28 parser.add_argument('--o_plot', help='End to End plot')
29 return parser.parse_args()
30
31
32 args = parse_command_line(sys.argv)
33
34
35 u = mda.Universe(args.istr, args.itraj,
36 topology_format=args.istrext, format=args.itrajext)
37
38 ntermatoms = "(segid %s and name N)" % \
39 (args.isegid1)
40 ctermatoms = "(segid %s and name C)" % \
41 (args.isegid1)
42 # not sure how robust this selection really is
43 nterm = u.select_atoms(ntermatoms)[0] # first atom named N
44 cterm = u.select_atoms(ctermatoms)[-1] # takes the last atom named 'C'
45
46 enddist = []
47
48 for ts in u.trajectory: # iterate through all frames
49 r = cterm.position - nterm.position # e-to-e vector from atom positions
50 d = numpy.linalg.norm(r) # end-to-end distance
51 enddist.append((ts.frame, d))
52
53 enddist = np.array(enddist)
54
55
56 color = itertools.cycle(['r', 'b', 'gold'])
57
58 fig, axs = plt.subplots(1, 2, sharex=False, sharey=False, tight_layout=True)
59
60 params = {
61 'axes.labelsize': 8,
62 'legend.fontsize': 10,
63 'xtick.labelsize': 10,
64 'ytick.labelsize': 10,
65 'text.usetex': False,
66 'figure.figsize': [4.5, 4.5],
67 'figure.dpi': 300
68 }
69 plt.rcParams.update(params)
70
71 axs[0].plot(enddist[:, 0], enddist[:, 1], 'r-', lw=2, label=args.ilabel)
72 axs[0].set_xlabel("number of frames")
73 axs[0].set_ylabel(r"End to end distance ($\AA$)")
74 axs[0].legend()
75
76 n, bins, patches = axs[1].hist(enddist[:, 1], color=next(
77 color), label=args.ilabel, alpha=0.5, density=True, stacked=True)
78
79 axs[1].legend()
80 axs[1].set_ylabel('Density Normalised Frequency')
81 axs[1].set_xlabel(r'End to end distance ($\AA$)')
82 fig.suptitle(args.ititle1, fontsize=12, fontweight='bold')
83 fig.subplots_adjust(top=0.45)
84
85 print(
86 " \n".join(
87 [
88 'The End to End distance is measured between the following atoms:',
89 str(nterm),
90 str(cterm)]))
91
92 # svg is better but sticking with png for now
93 plt.savefig(args.o_plot, format='png')
94
95
96 np.savetxt(args.output1, enddist, delimiter='\t')