Repository 'mdanalysis_hbonds'
hg clone https://toolshed.g2.bx.psu.edu/repos/chemteam/mdanalysis_hbonds

Changeset 4:4c36f5ad2799 (2020-08-24)
Previous changeset 3:5efd0c95f97e (2020-05-20) Next changeset 5:dfda5e713926 (2020-10-28)
Commit message:
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 1b23e024af45cc0999d9142d07de6897d4189ec2"
modified:
hbonds.xml
macros.xml
ramachandran_plots.py
added:
end-to-end.py
extract_rmsd.py
ramachandran_auto_protein.py
ramachandran_auto_protein_html.j2
test-data/test.yml
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 end-to-end.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/end-to-end.py Mon Aug 24 16:53:30 2020 -0400
[
@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+
+import argparse
+import itertools
+import sys
+
+import MDAnalysis as mda
+
+import matplotlib
+import matplotlib.pyplot as plt
+
+import numpy as np
+import numpy.linalg
+
+matplotlib.use('Agg')  # noqa
+
+
+def parse_command_line(argv):
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--itraj', help='input traj')
+    parser.add_argument('--istr', help='input str')
+    parser.add_argument('--itrajext', help='input traj ext')
+    parser.add_argument('--istrext', help='input str ext')
+    parser.add_argument('--isegid1', help='segid 1')
+    parser.add_argument('--ilabel', help='plot label')
+    parser.add_argument('--ititle1', help='plot title')
+    parser.add_argument('--output1', help='output1 - timeseries')
+    parser.add_argument('--o_plot', help='End to End plot')
+    return parser.parse_args()
+
+
+args = parse_command_line(sys.argv)
+
+
+u = mda.Universe(args.istr, args.itraj,
+                 topology_format=args.istrext, format=args.itrajext)
+
+ntermatoms = "(segid %s and name N)" % \
+    (args.isegid1)
+ctermatoms = "(segid %s and name C)" % \
+    (args.isegid1)
+# not sure how robust this selection really is
+nterm = u.select_atoms(ntermatoms)[0]  # first atom named N
+cterm = u.select_atoms(ctermatoms)[-1]  # takes the last atom named 'C'
+
+enddist = []
+
+for ts in u.trajectory:  # iterate through all frames
+    r = cterm.position - nterm.position  # e-to-e vector from atom positions
+    d = numpy.linalg.norm(r)   # end-to-end distance
+    enddist.append((ts.frame, d))
+
+enddist = np.array(enddist)
+
+
+color = itertools.cycle(['r', 'b', 'gold'])
+
+fig, axs = plt.subplots(1, 2, sharex=False, sharey=False, tight_layout=True)
+
+params = {
+    'axes.labelsize': 8,
+    'legend.fontsize': 10,
+    'xtick.labelsize': 10,
+    'ytick.labelsize': 10,
+    'text.usetex': False,
+    'figure.figsize': [4.5, 4.5],
+    'figure.dpi': 300
+}
+plt.rcParams.update(params)
+
+axs[0].plot(enddist[:, 0], enddist[:, 1], 'r-', lw=2, label=args.ilabel)
+axs[0].set_xlabel("number of frames")
+axs[0].set_ylabel(r"End to end distance  ($\AA$)")
+axs[0].legend()
+
+n, bins, patches = axs[1].hist(enddist[:, 1], color=next(
+    color), label=args.ilabel, alpha=0.5, density=True, stacked=True)
+
+axs[1].legend()
+axs[1].set_ylabel('Density Normalised Frequency')
+axs[1].set_xlabel(r'End to end distance ($\AA$)')
+fig.suptitle(args.ititle1, fontsize=12, fontweight='bold')
+fig.subplots_adjust(top=0.45)
+
+print(
+    " \n".join(
+        [
+            'The End to End distance is measured between the following atoms:',
+            str(nterm),
+            str(cterm)]))
+
+# svg is better but sticking with png for now
+plt.savefig(args.o_plot, format='png')
+
+
+np.savetxt(args.output1, enddist, delimiter='\t')
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 extract_rmsd.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/extract_rmsd.py Mon Aug 24 16:53:30 2020 -0400
[
@@ -0,0 +1,129 @@
+import argparse
+import json
+
+import MDAnalysis as m
+from MDAnalysis.analysis import align, rms
+from MDAnalysis.analysis.base import AnalysisFromFunction
+from MDAnalysis.coordinates.memory import MemoryReader
+
+import numpy as np
+
+
+def calc_rmsd(str_files, traj_files, ref_str, str_format, traj_format,
+              ref_str_format, filepath_out, group, start, end, step,
+              fitting_atoms):
+    """
+    the function will cycle through range 0 to no_t and load all files found.
+
+    str_files: text file with filepaths for structures, one on each line
+    traj_files: text file with filepaths for trajectories, one on each line
+    ref_str: reference structure for fitting
+    filepath_in: directory where the files are located
+    filepath_out: pickle file where results (3D matrix) should be saved to
+
+    group: atoms for which RMSD should be calculated;
+        use the MDAnalysis selection language
+    fitting_atoms: atoms used for str alignment prior to RMSD calculation;
+        use the MDAnalysis selection language
+
+    start: first trajectory frame to calculate RMSD
+    end: last trajectory frame to calculate RMSD
+    step: how frequently frames are sampled between start and end; obviously,
+        the larger the step, the quicker the script finishes
+    """
+
+    # open list of files
+    with open(str_files) as f1, open(traj_files) as f2:
+        str_file_list = f1.read().strip().split('\n')
+        traj_file_list = f2.read().strip().split('\n')
+
+        if sum(1 for line in f1) != sum(1 for line in f2):
+            raise IOError('Number of structure and trajectory files unequal.')
+
+    no_t = len(traj_file_list)
+
+    data = np.zeros((no_t, no_t,
+                    int((end - start)/step + ((end - start) % step > 0))))
+
+    # load files
+    universes = {}
+
+    for traj in range(no_t):
+        mobile = m.Universe(str_file_list[traj], traj_file_list[traj],
+                            format=traj_format, topology_format=str_format)
+        ref = m.Universe(ref_str, topology_format=ref_str_format)
+
+        mobile.trajectory[-1]  # set mobile trajectory to last frame
+        ref.trajectory[0]  # set reference trajectory to first frame
+
+        # perform alignment
+        align.AlignTraj(mobile, ref, select=fitting_atoms,
+                        in_memory=True).run()
+
+        grp = mobile.select_atoms(group)
+        universes[traj] = m.core.universe.Merge(grp)  # create Universe w grp
+        coordinates = AnalysisFromFunction(lambda ag: ag.positions.copy(),
+                                           grp).run().results  # write to uv
+        universes[traj].load_new(coordinates, format=MemoryReader)
+
+    print("All trajs loaded by MDAnalysis")
+
+    # calculate differences
+    for traj1 in range(no_t):
+        print("Calculating differences for traj {}".format(traj1))
+        for traj2 in range(traj1):
+
+            u1 = universes[traj1]
+            u2 = universes[traj2]
+
+            l1 = u1.select_atoms(group)
+            l2 = u2.select_atoms(group)
+
+            rmsd = rms.RMSD(l1, l2)
+
+            rmsd.run()
+
+            data[traj1, traj2] = rmsd.rmsd[:, 2]
+            data[traj2, traj1] = rmsd.rmsd[:, 2]
+
+    with open(filepath_out, 'w') as f:
+        json.dump(data.tolist(), f, indent=4, sort_keys=True)
+
+    print("Done!")
+    return
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--trajs', required=True,
+                        help='File containing trajectory filepaths.')
+    parser.add_argument("--strs",
+                        help='File containing structure filepaths.')
+    parser.add_argument("--ref-str",
+                        help='File containing reference structure.')
+    parser.add_argument('--traj-format', required=True,
+                        help='Trajectory format.')
+    parser.add_argument("--str-format", help='Structure format.')
+    parser.add_argument("--ref-str-format",
+                        help='Reference structure format.')
+    parser.add_argument('-o', '--outfile',
+                        help="Path to the output JSON file")
+    parser.add_argument('--group', help="Atoms for which RMSD should be"
+                        "calculated in MDAnalysis selection language")
+    parser.add_argument('--fitting', help="Fitting atoms for alignment"
+                        "prior to RMSD calculation")
+    parser.add_argument('--start', type=int,
+                        help="First trajectory frame to calculate RMSD")
+    parser.add_argument('--end', type=int,
+                        help="Last trajectory frame to calculate RMSD")
+    parser.add_argument('--step', type=int,
+                        help="Frame sampling frequency for RMSD calculation")
+    args = parser.parse_args()
+
+    calc_rmsd(args.strs, args.trajs, args.ref_str, args.str_format,
+              args.traj_format, args.ref_str_format, args.outfile,
+              args.group, args.start, args.end, args.step, args.fitting)
+
+
+if __name__ == "__main__":
+    main()
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 hbonds.xml
--- a/hbonds.xml Wed May 20 13:05:00 2020 -0400
+++ b/hbonds.xml Mon Aug 24 16:53:30 2020 -0400
[
@@ -1,10 +1,11 @@
-<tool id="mdanalysis_hbonds" name="Hydrogen Bond Analysis" version="@VERSION@">
+<tool id="mdanalysis_hbonds" name="Hydrogen Bond Analysis" version="@TOOL_VERSION@+galaxy@GALAXY_VERSION@">
     <description>- analyze H-bonds between two segments</description>
     <macros>
         <import>macros.xml</import>
+        <token name="@GALAXY_VERSION@">0</token>
     </macros>   
     <expand macro="requirements">
-        <requirement type="package" version="0.25.1">pandas</requirement>
+        <requirement type="package" version="1.0.5">pandas</requirement>
     </expand>
     <command detect_errors="exit_code">
 <![CDATA[
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 macros.xml
--- a/macros.xml Wed May 20 13:05:00 2020 -0400
+++ b/macros.xml Mon Aug 24 16:53:30 2020 -0400
b
@@ -1,36 +1,36 @@
 <macros>
-    <token name="@VERSION@">0.20</token>
+    <token name="@TOOL_VERSION@">1.0.0</token>
     <xml name="requirements">
         <requirements>
-            <requirement type="package" version="0.20.1">mdanalysis</requirement>
+            <requirement type="package" version="@TOOL_VERSION@">mdanalysis</requirement>
             <yield/>
         </requirements>
     </xml>
     <xml name="analysis_inputs">
-        <param format="dcd,xtc" name="trajin" type="data" label="DCD/XTC trajectory input"/>
-        <param format="pdb,gro" name="strin" type="data" label="PDB/GRO input"/>
+        <param format="dcd,xtc" name="trajin" type="data" label="DCD/XTC trajectory input" />
+        <param format="pdb,gro" name="strin" type="data" label="PDB/GRO input" />
         <yield/>
     </xml>
     <xml name="sanitizer">
         <sanitizer invalid_char="">
-            <valid initial="string.ascii_letters,string.digits"/>
-         </sanitizer>
-         <yield/>
+            <valid initial="string.ascii_letters,string.digits" />
+        </sanitizer>
+        <yield/>
     </xml>
     <xml name="sanitizer_resids">
         <sanitizer invalid_char="">
-            <valid initial="string.digits"/>
+            <valid initial="string.digits" />
         </sanitizer>
         <yield/>
     </xml>
     <xml name="tests_inputs">
-        <param name="trajin" value="test.dcd" ftype="dcd"/>
-        <param name="strin" value="test.pdb" ftype="pdb"/>
+        <param name="trajin" value="test.dcd" ftype="dcd" />
+        <param name="strin" value="test.pdb" ftype="pdb" />
         <yield/>
     </xml>
     <xml name="tests_inputs_gmx">
-        <param name="trajin" value="test.xtc" ftype="xtc"/>
-        <param name="strin" value="test.gro" ftype="gro"/>
+        <param name="trajin" value="test.xtc" ftype="xtc" />
+        <param name="strin" value="test.gro" ftype="gro" />
         <yield/>
     </xml>
     <xml name="citations">
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 ramachandran_auto_protein.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ramachandran_auto_protein.py Mon Aug 24 16:53:30 2020 -0400
[
@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+
+import argparse
+import base64
+import importlib
+import sys
+
+import MDAnalysis as mda
+from MDAnalysis.analysis.dihedrals import Ramachandran
+
+import h5py
+
+from jinja2 import Environment, FileSystemLoader
+
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+import numpy as np
+import numpy.linalg
+
+import seaborn as sns
+
+
+matplotlib.use('Agg')  # noqa
+
+
+def parse_command_line(argv):
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--itraj', help='input traj')
+    parser.add_argument('--istr', help='input str')
+    parser.add_argument('--itrajext', help='input traj ext')
+    parser.add_argument('--istrext', help='input str ext')
+    parser.add_argument('--isegid1', help='segid 1')
+    parser.add_argument('--iresid1', help='resid start')
+    parser.add_argument('--iresid2', help='resid end')
+    parser.add_argument('--iresname', help='resname e.g. ALA')
+    parser.add_argument('--igroupby', help='groupby names or ids')
+    parser.add_argument('--itemplatepath', help='template path')
+    parser.add_argument('--o_plot1', help='MDA Ramachandran plot')
+    parser.add_argument('--o_plot2', help='Seaborn Ramachandran plot')
+    parser.add_argument('--o_data1', help='Timeseries in HDF5 format')
+    parser.add_argument('--o_html1', help='Html overview output of all plots')
+    return parser.parse_args()
+
+
+args = parse_command_line(sys.argv)
+
+currentpath = "."
+if args.itemplatepath is not None:
+    currentpath = args.itemplatepath
+
+
+u = mda.Universe(args.istr, args.itraj,
+                 topology_format=args.istrext, format=args.itrajext)
+selection = "(segid %s)" % \
+    (args.isegid1)
+
+if args.iresname is not None:
+    selection = "(segid %s and resname %s)" % \
+        (args.isegid1, args.iresname)
+
+if args.iresid1 is not None and args.iresid2 is not None:
+    assert(int(args.iresid1) > 0), "ResID numbering starts at 1 for this tool."
+    assert(int(args.iresid2) > 0), "ResID numbering starts at 1 for this tool."
+    assert(int(args.iresid2) > int(args.iresid1)
+           ), "ResID2 must be at least ResID1+1"
+    selection = "(segid %s and resid %s-%s)" % \
+        (args.isegid1, int(args.iresid1), int(args.iresid2))
+    if args.iresname is not None:
+        selection = "(segid %s and resid %s-%s and resname %s)" % \
+            (args.isegid1, int(args.iresid1), int(args.iresid2), args.iresname)
+
+r = u.select_atoms(selection)
+
+assert(r != u.select_atoms('name thiscannotpossiblyexist')
+       ), \
+               """The selection you specified returns an empty result.
+               Check segment names and residue ID's. Also check the
+                structure and trajectory file selected are the correct ones"""
+
+if args.igroupby is not None:
+    group_selections = {}  # dictionary of selections
+    if args.igroupby == 'name':
+        groupby = sorted(list(set(r.resnames)))
+        for e in groupby:
+            s = r & u.select_atoms("resname %s" % e)
+            this_sel = "%s and resname %s" % (selection, e)
+            group_selections[this_sel] = s
+    elif args.igroupby == 'id':
+        groupby = sorted(list(set(r.resids)))
+        for e in groupby:
+            s = r & u.select_atoms("resid %s" % e)
+            this_sel = "%s and resid %s" % (selection, e)
+            group_selections[this_sel] = s
+    else:
+        assert False, ("Invalid argument for igroupby. "
+                       "Only name and id are valid options.")
+
+
+def ramachandran_plot(atomgroup, selection, outputfile1, outputfile2,
+                      image_format='png'):
+    # plot standard mdanalysis and seaborn 2D with kde
+    R = Ramachandran(atomgroup).run()
+    fig, ax = plt.subplots(figsize=plt.figaspect(1))
+    R.plot(ax=ax, color='k', marker='.', ref=True)
+
+    a = R.angles.reshape(np.prod(R.angles.shape[:2]), 2)
+    # open hdf file
+    with h5py.File(args.o_data1, 'a') as f:
+        setname = "%s" % (selection)
+        f["/" + setname + "/ramachandran/phi"] = a[:, 0]
+        f["/" + setname + "/ramachandran/psi"] = a[:, 1]
+    plt.tight_layout()
+    # svg is better but sticking with png for now
+    plt.savefig(outputfile1, format=image_format)
+
+    sns.reset_defaults()
+    importlib.reload(plt)
+    importlib.reload(sns)
+    with sns.axes_style("white"):
+        h = sns.jointplot(x=a[:, 0], y=a[:, 1],
+                          kind="kde", space=0)
+        h.set_axis_labels(r'$\phi$ (deg)', r'$\psi$ (deg)')
+        h.ax_joint.set_xlim(-180, 180)
+        h.ax_joint.set_ylim(-180, 180)
+        h.ax_joint.xaxis.set_major_locator(ticker.MultipleLocator(60))
+        h.ax_joint.yaxis.set_major_locator(ticker.MultipleLocator(60))
+        plt.savefig(outputfile2, format=image_format, bbox_inches='tight')
+
+
+def get_base64_encoded_image(image_path):
+    """  encode image to string for use in html later"""
+    with open(image_path, "rb") as img_file:
+        return base64.b64encode(img_file.read()).decode('utf-8')
+
+
+plots = []
+if args.igroupby is not None:
+    for k, v in group_selections.items():
+        print(k, v)
+        try:
+            ramachandran_plot(v, str(k), "ramachandran1" +
+                              str(k), "ramachandran2" + str(k))
+            plots.append({'Name': "%s" % (k), 'plot1':
+                         get_base64_encoded_image("ramachandran1" + str(k)),
+                         'plot2': get_base64_encoded_image("ramachandran2"
+                          + str(k))})
+        except Exception as einstance:
+            print(type(einstance))
+            print(einstance.args)
+            print(einstance)
+
+ramachandran_plot(r, selection, args.o_plot1, args.o_plot2)
+plots.insert(0, {'Name': selection, 'plot1': get_base64_encoded_image(
+    args.o_plot1), 'plot2': get_base64_encoded_image(args.o_plot2)})
+
+template_environment = Environment(loader=FileSystemLoader(
+    currentpath), lstrip_blocks=True, trim_blocks=True)
+template = template_environment.get_template(
+    'ramachandran_auto_protein_html.j2')
+with open(args.o_html1, 'w+') as f:
+    f.write(template.render(title="Ramachandran Plots", plots=plots))
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 ramachandran_auto_protein_html.j2
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ramachandran_auto_protein_html.j2 Mon Aug 24 16:53:30 2020 -0400
[
@@ -0,0 +1,25 @@
+<html>
+
+<head>
+  <meta http-equiv="Content-Type" content="text/html; charset=UTF-8" />
+  <title>{{ title }}</title>
+</head>
+
+<body>
+  <table>
+    <tr>
+      <th>Selection</th>
+      <th>Ramachandran scatter plot</th>
+      <th>Ramachandran histogram </th>
+    </tr>
+    {% for plot in plots %}
+    <tr>
+      <td>{{ plot['Name'] }}</td>
+      <td style="vertical-align:center"><img src="data:image/png;base64,{{plot['plot1']}}" /> </td>
+      <td style="vertical-align:center"><img src="data:image/png;base64,{{plot['plot2']}}" /> </td>
+    </tr>
+    {% endfor %}
+  </table>
+</body>
+
+</html>
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 ramachandran_plots.py
--- a/ramachandran_plots.py Wed May 20 13:05:00 2020 -0400
+++ b/ramachandran_plots.py Mon Aug 24 16:53:30 2020 -0400
[
@@ -10,11 +10,16 @@
 
 import matplotlib
 import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
 
 import numpy as np
 
 import seaborn as sns
 
+
+import yaml
+
 matplotlib.use('Agg')  # noqa
 
 
@@ -24,30 +29,7 @@
     parser.add_argument('--istr', help='input str')
     parser.add_argument('--itrajext', help='input traj ext')
     parser.add_argument('--istrext', help='input str ext')
-    parser.add_argument('--isegid1', help='segid 1')
-    parser.add_argument('--iresid1', help='resid 1')
-    parser.add_argument('--iname1', help='name 1')
-    parser.add_argument('--isegid2', help='segid 2')
-    parser.add_argument('--iresid2', help='resid 2')
-    parser.add_argument('--iname2', help='name 2')
-    parser.add_argument('--isegid3', help='segid 3')
-    parser.add_argument('--iresid3', help='resid 3')
-    parser.add_argument('--iname3', help='name 3')
-    parser.add_argument('--isegid4', help='segid 4')
-    parser.add_argument('--iresid4', help='resid 4')
-    parser.add_argument('--iname4', help='name 4')
-    parser.add_argument('--isegid5', help='segid 1')
-    parser.add_argument('--iresid5', help='resid 1')
-    parser.add_argument('--iname5', help='name 1')
-    parser.add_argument('--isegid6', help='segid 2')
-    parser.add_argument('--iresid6', help='resid 2')
-    parser.add_argument('--iname6', help='name 2')
-    parser.add_argument('--isegid7', help='segid 3')
-    parser.add_argument('--iresid7', help='resid 3')
-    parser.add_argument('--iname7', help='name 3')
-    parser.add_argument('--isegid8', help='segid 4')
-    parser.add_argument('--iresid8', help='resid 4')
-    parser.add_argument('--iname8', help='name 4')
+    parser.add_argument('--iyml', help='input in yml format')
     parser.add_argument('--output', help='output')
     parser.add_argument('--oramachandran_plot', help='dihedral plot')
     return parser.parse_args()
@@ -55,35 +37,35 @@
 
 args = parse_command_line(sys.argv)
 
+with open(args.iyml) as file:
+    params = yaml.load(file, Loader=yaml.FullLoader)
+
 Dihedral = namedtuple(
     'Dihedral', ['atom1', 'atom2', 'atom3', 'atom4'])
 
+for k, v in params.items():
+    for a in ['phi', 'psi']:
+        assert (a in v), "Key %s is missing in inputs: %s " % (a, k)
+        atoms = []
+        for b in ['atom1', 'atom2', 'atom3', 'atom4']:
+            assert (b in v[a]), "Key %s is missing in inputs: %s %s" % (
+                b, k, a)
+            for c in ['segid', 'resid', 'name']:
+                assert (c in v[a][b]), \
+                        "Key %s is missing in inputs: %s %s %s " % (c, k, a, b)
+            atoms.append("(segid %s and resid %s and name %s)" %
+                         (v[a][b]['segid'], v[a][b]['resid'], v[a][b]['name']))
+        print(atoms)
+        if a == 'phi':
+            dihe_phi = Dihedral(atoms[0], atoms[1], atoms[2], atoms[3])
+        if a == 'psi':
+            dihe_psi = Dihedral(atoms[0], atoms[1], atoms[2], atoms[3])
+
 # order of dihedral atom is the crystallographic definition
 # (see glycanstructure.org)
 
-# phi
-atom1 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid1, args.iresid1, args.iname1)
-atom2 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid2, args.iresid2, args.iname2)
-atom3 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid3, args.iresid3, args.iname3)
-atom4 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid4, args.iresid4, args.iname4)
-
-dihe_phi = Dihedral(atom1, atom2, atom3, atom4)
-
-# psi
-atom1 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid5, args.iresid5, args.iname5)
-atom2 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid6, args.iresid6, args.iname6)
-atom3 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid7, args.iresid7, args.iname7)
-atom4 = "(segid %s and resid %s and name %s)" % \
-    (args.isegid8, args.iresid8, args.iname8)
-
-dihe_psi = Dihedral(atom1, atom2, atom3, atom4)
+assert(dihe_phi), "phi dihedral doesn't exist"
+assert(dihe_psi), "psi dihedral doesn't exist"
 
 
 def calc_torsion(dihedral):
@@ -120,8 +102,11 @@
     writer.writerows(zip(phi_frame, phi_series, psi_series))
 
 with sns.axes_style("white"):
-    h = sns.jointplot(x=phi_series, y=psi_series, kind="kde", legend=True)
-    h.set_axis_labels(r'$\Phi$ (degrees)', r'$\Psi$ (degrees)')
+    h = sns.jointplot(x=phi_series, y=psi_series,
+                      kind="kde", space=0, legend=True)
+    h.set_axis_labels(r'$\phi$ (degrees)', r'$\psi$ (degrees)')
     h.ax_joint.set_xlim(-180, 180)
     h.ax_joint.set_ylim(-180, 180)
-    plt.savefig(args.oramachandran_plot, format='png')
+    h.ax_joint.xaxis.set_major_locator(ticker.MultipleLocator(60))
+    h.ax_joint.yaxis.set_major_locator(ticker.MultipleLocator(60))
+    plt.savefig(args.oramachandran_plot, format='png', bbox_inches='tight')
b
diff -r 5efd0c95f97e -r 4c36f5ad2799 test-data/test.yml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/test.yml Mon Aug 24 16:53:30 2020 -0400
b
@@ -0,0 +1,36 @@
+ramachandran1:
+  phi:
+    atom1:
+     segid: HET
+     resid: 3
+     name: O5
+    atom2:
+     segid: HET
+     resid: 3
+     name: C1
+    atom3:
+     segid: HET
+     resid: 2
+     name: O4
+    atom4:
+     segid: HET
+     resid: 2
+     name: C4
+  psi:
+    atom1:
+     segid: HET
+     resid: 3
+     name: C1
+    atom2:
+     segid: HET
+     resid: 2
+     name: O4
+    atom3:
+     segid: HET
+     resid: 2
+     name: C4
+    atom4:
+     segid: HET
+     resid: 2
+     name: C5
+  comment: pick visually using VMD using labels. Go to labels, dihedral to see the information about resname resid and atomname and then lookup the segname for ach atom.