Mercurial > repos > chemteam > mdanalysis_angle

comparison ramachandran_plots.py @ 7:a0d210b9d287 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:32:30 -0400
parents 7c5fd4117a07
children
comparison
equal deleted inserted replaced
6:7c5fd4117a07 7:a0d210b9d287
8 import MDAnalysis as mda 8 import MDAnalysis as mda
9 from MDAnalysis.lib.distances import calc_dihedrals 9 from MDAnalysis.lib.distances import calc_dihedrals
10 10
11 import matplotlib 11 import matplotlib
12 import matplotlib.pyplot as plt 12 import matplotlib.pyplot as plt
13 import matplotlib.ticker as ticker
14
13 15
14 import numpy as np 16 import numpy as np
15 17
16 import seaborn as sns 18 import seaborn as sns
19
20
21 import yaml
17 22
18 matplotlib.use('Agg') # noqa 23 matplotlib.use('Agg') # noqa
19 24
20 25
21 def parse_command_line(argv): 26 def parse_command_line(argv):
22 parser = argparse.ArgumentParser() 27 parser = argparse.ArgumentParser()
23 parser.add_argument('--itraj', help='input traj') 28 parser.add_argument('--itraj', help='input traj')
24 parser.add_argument('--istr', help='input str') 29 parser.add_argument('--istr', help='input str')
25 parser.add_argument('--itrajext', help='input traj ext') 30 parser.add_argument('--itrajext', help='input traj ext')
26 parser.add_argument('--istrext', help='input str ext') 31 parser.add_argument('--istrext', help='input str ext')
27 parser.add_argument('--isegid1', help='segid 1') 32 parser.add_argument('--iyml', help='input in yml format')
28 parser.add_argument('--iresid1', help='resid 1')
29 parser.add_argument('--iname1', help='name 1')
30 parser.add_argument('--isegid2', help='segid 2')
31 parser.add_argument('--iresid2', help='resid 2')
32 parser.add_argument('--iname2', help='name 2')
33 parser.add_argument('--isegid3', help='segid 3')
34 parser.add_argument('--iresid3', help='resid 3')
35 parser.add_argument('--iname3', help='name 3')
36 parser.add_argument('--isegid4', help='segid 4')
37 parser.add_argument('--iresid4', help='resid 4')
38 parser.add_argument('--iname4', help='name 4')
39 parser.add_argument('--isegid5', help='segid 1')
40 parser.add_argument('--iresid5', help='resid 1')
41 parser.add_argument('--iname5', help='name 1')
42 parser.add_argument('--isegid6', help='segid 2')
43 parser.add_argument('--iresid6', help='resid 2')
44 parser.add_argument('--iname6', help='name 2')
45 parser.add_argument('--isegid7', help='segid 3')
46 parser.add_argument('--iresid7', help='resid 3')
47 parser.add_argument('--iname7', help='name 3')
48 parser.add_argument('--isegid8', help='segid 4')
49 parser.add_argument('--iresid8', help='resid 4')
50 parser.add_argument('--iname8', help='name 4')
51 parser.add_argument('--output', help='output') 33 parser.add_argument('--output', help='output')
52 parser.add_argument('--oramachandran_plot', help='dihedral plot') 34 parser.add_argument('--oramachandran_plot', help='dihedral plot')
53 return parser.parse_args() 35 return parser.parse_args()
54 36
55 37
56 args = parse_command_line(sys.argv) 38 args = parse_command_line(sys.argv)
57 39
40 with open(args.iyml) as file:
41 params = yaml.load(file, Loader=yaml.FullLoader)
42
58 Dihedral = namedtuple( 43 Dihedral = namedtuple(
59 'Dihedral', ['atom1', 'atom2', 'atom3', 'atom4']) 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])
60 63
61 # order of dihedral atom is the crystallographic definition 64 # order of dihedral atom is the crystallographic definition
62 # (see glycanstructure.org) 65 # (see glycanstructure.org)
63 66
64 # phi 67 assert(dihe_phi), "phi dihedral doesn't exist"
65 atom1 = "(segid %s and resid %s and name %s)" % \ 68 assert(dihe_psi), "psi dihedral doesn't exist"
66 (args.isegid1, args.iresid1, args.iname1)
67 atom2 = "(segid %s and resid %s and name %s)" % \
68 (args.isegid2, args.iresid2, args.iname2)
69 atom3 = "(segid %s and resid %s and name %s)" % \
70 (args.isegid3, args.iresid3, args.iname3)
71 atom4 = "(segid %s and resid %s and name %s)" % \
72 (args.isegid4, args.iresid4, args.iname4)
73
74 dihe_phi = Dihedral(atom1, atom2, atom3, atom4)
75
76 # psi
77 atom1 = "(segid %s and resid %s and name %s)" % \
78 (args.isegid5, args.iresid5, args.iname5)
79 atom2 = "(segid %s and resid %s and name %s)" % \
80 (args.isegid6, args.iresid6, args.iname6)
81 atom3 = "(segid %s and resid %s and name %s)" % \
82 (args.isegid7, args.iresid7, args.iname7)
83 atom4 = "(segid %s and resid %s and name %s)" % \
84 (args.isegid8, args.iresid8, args.iname8)
85
86 dihe_psi = Dihedral(atom1, atom2, atom3, atom4)
87 69
88 70
89 def calc_torsion(dihedral): 71 def calc_torsion(dihedral):
90 """atom 1 -4 are valid atom selections. torsion in degrees is returned""" 72 """atom 1 -4 are valid atom selections. torsion in degrees is returned"""
91 A = u.select_atoms(dihedral.atom1).positions 73 A = u.select_atoms(dihedral.atom1).positions
118 with open(args.output, 'w') as f: 100 with open(args.output, 'w') as f:
119 writer = csv.writer(f, delimiter='\t') 101 writer = csv.writer(f, delimiter='\t')
120 writer.writerows(zip(phi_frame, phi_series, psi_series)) 102 writer.writerows(zip(phi_frame, phi_series, psi_series))
121 103
122 with sns.axes_style("white"): 104 with sns.axes_style("white"):
123 h = sns.jointplot(x=phi_series, y=psi_series, kind="kde", legend=True) 105 h = sns.jointplot(x=phi_series, y=psi_series,
124 h.set_axis_labels(r'$\Phi$ (degrees)', r'$\Psi$ (degrees)') 106 kind="kde", space=0, legend=True)
107 h.set_axis_labels(r'$\phi$ (degrees)', r'$\psi$ (degrees)')
125 h.ax_joint.set_xlim(-180, 180) 108 h.ax_joint.set_xlim(-180, 180)
126 h.ax_joint.set_ylim(-180, 180) 109 h.ax_joint.set_ylim(-180, 180)
127 plt.savefig(args.oramachandran_plot, format='png') 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')