Mercurial > repos > chemteam > mdanalysis_ramachandran_protein

comparison ramachandran_plots.py @ 0:0f270722aca6 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:27:56 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0f270722aca6
1 #!/usr/bin/env python
2
3 import argparse
4 import csv
5 import sys
6 from collections import namedtuple
7
8 import MDAnalysis as mda
9 from MDAnalysis.lib.distances import calc_dihedrals
10
11 import matplotlib
12 import matplotlib.pyplot as plt
13 import matplotlib.ticker as ticker
14
15
16 import numpy as np
17
18 import seaborn as sns
19
20
21 import yaml
22
23 matplotlib.use('Agg') # noqa
24
25
26 def parse_command_line(argv):
27 parser = argparse.ArgumentParser()
28 parser.add_argument('--itraj', help='input traj')
29 parser.add_argument('--istr', help='input str')
30 parser.add_argument('--itrajext', help='input traj ext')
31 parser.add_argument('--istrext', help='input str ext')
32 parser.add_argument('--iyml', help='input in yml format')
33 parser.add_argument('--output', help='output')
34 parser.add_argument('--oramachandran_plot', help='dihedral plot')
35 return parser.parse_args()
36
37
38 args = parse_command_line(sys.argv)
39
40 with open(args.iyml) as file:
41 params = yaml.load(file, Loader=yaml.FullLoader)
42
43 Dihedral = namedtuple(
44 'Dihedral', ['atom1', 'atom2', 'atom3', 'atom4'])
45
46 for k, v in params.items():
47 for a in ['phi', 'psi']:
48 assert (a in v), "Key %s is missing in inputs: %s " % (a, k)
49 atoms = []
50 for b in ['atom1', 'atom2', 'atom3', 'atom4']:
51 assert (b in v[a]), "Key %s is missing in inputs: %s %s" % (
52 b, k, a)
53 for c in ['segid', 'resid', 'name']:
54 assert (c in v[a][b]), \
55 "Key %s is missing in inputs: %s %s %s " % (c, k, a, b)
56 atoms.append("(segid %s and resid %s and name %s)" %
57 (v[a][b]['segid'], v[a][b]['resid'], v[a][b]['name']))
58 print(atoms)
59 if a == 'phi':
60 dihe_phi = Dihedral(atoms[0], atoms[1], atoms[2], atoms[3])
61 if a == 'psi':
62 dihe_psi = Dihedral(atoms[0], atoms[1], atoms[2], atoms[3])
63
64 # order of dihedral atom is the crystallographic definition
65 # (see glycanstructure.org)
66
67 assert(dihe_phi), "phi dihedral doesn't exist"
68 assert(dihe_psi), "psi dihedral doesn't exist"
69
70
71 def calc_torsion(dihedral):
72 """atom 1 -4 are valid atom selections. torsion in degrees is returned"""
73 A = u.select_atoms(dihedral.atom1).positions
74 B = u.select_atoms(dihedral.atom2).positions
75 C = u.select_atoms(dihedral.atom3).positions
76 D = u.select_atoms(dihedral.atom4).positions
77
78 dihe = calc_dihedrals(A, B, C, D)
79 return np.rad2deg(dihe)
80
81
82 u = mda.Universe(args.istr, args.itraj,
83 topology_format=args.istrext, format=args.itrajext)
84
85 phi_trajdata = np.array(
86 [(u.trajectory.frame, calc_torsion(dihe_phi)) for ts in u.trajectory])
87 psi_trajdata = np.array(
88 [(u.trajectory.frame, calc_torsion(dihe_psi)) for ts in u.trajectory])
89
90 print(phi_trajdata, psi_trajdata)
91
92 phi_frame, phi_series = phi_trajdata.T
93 psi_frame, psi_series = psi_trajdata.T
94
95 phi_series = np.concatenate(phi_series, axis=0)
96 psi_series = np.concatenate(psi_series, axis=0)
97
98 zip(phi_frame, phi_series, psi_series)
99
100 with open(args.output, 'w') as f:
101 writer = csv.writer(f, delimiter='\t')
102 writer.writerows(zip(phi_frame, phi_series, psi_series))
103
104 with sns.axes_style("white"):
105 h = sns.jointplot(x=phi_series, y=psi_series,
106 kind="kde", space=0, legend=True)
107 h.set_axis_labels(r'$\phi$ (degrees)', r'$\psi$ (degrees)')
108 h.ax_joint.set_xlim(-180, 180)
109 h.ax_joint.set_ylim(-180, 180)
110 h.ax_joint.xaxis.set_major_locator(ticker.MultipleLocator(60))
111 h.ax_joint.yaxis.set_major_locator(ticker.MultipleLocator(60))
112 plt.savefig(args.oramachandran_plot, format='png', bbox_inches='tight')