Mercurial > repos > chemteam > mdanalysis_distance
changeset 5:d540ea77b909 draft
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 12823af06b7926cc56aaf9b59cfea9f16d342b8c"
| author | chemteam |
|---|---|
| date | Thu, 06 Feb 2020 19:43:03 -0500 |
| parents | 312f912de69d |
| children | aa4090b50e7b |
| files | distance.py distance.xml distance_multiple.py distance_single.py test-data/list1.txt test-data/list2.txt |
| diffstat | 6 files changed, 208 insertions(+), 86 deletions(-) [+] |
line wrap: on
line diff
--- a/distance.py Mon Oct 07 12:52:40 2019 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -#!/usr/bin/env python - -import argparse -import sys - -import MDAnalysis as mda - -import matplotlib -matplotlib.use('Agg') # noqa -import matplotlib.pyplot as plt - -import numpy as np - - -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 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('--output', help='output') - parser.add_argument('--odistance_plot', help='odistance plot') - return parser.parse_args() - - -args = parse_command_line(sys.argv) - -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) - -u = mda.Universe(args.istr, args.itraj, - topology_format=args.istrext, format=args.itrajext) -x = u.select_atoms(atom1) -y = u.select_atoms(atom2) - -with open(args.output, 'w') as f: - for t in u.trajectory: - r = x.positions - y.positions - d = np.linalg.norm(r) - f.write(str(t.frame) + '\t ') - f.write(str(d) + '\n') - -with open(args.output) as f: - g = [xtmp.strip() for xtmp in f] - data = [tuple(map(float, xtmp.split())) for xtmp in g[0:]] - time = [xtmp[0] for xtmp in data] - distance = [xtmp[1] for xtmp in data] - plt.plot(time, distance) - plt.xlabel('Frame No.') - plt.ylabel(r'Distance ($\AA$)') - plt.savefig(args.odistance_plot, format='png')
--- a/distance.xml Mon Oct 07 12:52:40 2019 -0400 +++ b/distance.xml Thu Feb 06 19:43:03 2020 -0500 @@ -6,33 +6,57 @@ <expand macro="requirements" /> <command detect_errors="exit_code"> <![CDATA[ - python '$__tool_directory__/distance.py' - --itraj '$trajin' - --istr '$strin' - --itrajext '$trajin.ext' - --istrext '$strin.ext' - --isegid1 '$segid1' - --iresid1 '$resid1' - --iname1 '$name1' - --isegid2 '$segid2' - --iresid2 '$resid2' - --iname2 '$name2' - --output '$output' - --odistance_plot '$distance_plot' - 2>&1 + #if $inps.inps == 'one': + python '$__tool_directory__/distance_single.py' + --isegid1 '$inps.segid1' + --iresid1 '$inps.resid1' + --iname1 '$inps.name1' + --isegid2 '$inps.segid2' + --iresid2 '$inps.resid2' + --iname2 '$inps.name2' + --odistance_plot '$distance_plot' + 2>&1 + #elif $inps.inps == 'multiple': + python '$__tool_directory__/distance_multiple.py' + --list1 '$inps.list1' + --list2 '$inps.list2' + #end if + --itraj '$trajin' + --istr '$strin' + --itrajext '$trajin.ext' + --istrext '$strin.ext' + --output '$output' + $header + ]]></command> <inputs> <expand macro="analysis_inputs"/> - <param name="segid1" type="text" value="PRO" label="Segment ID of atom 1"/> - <param name="resid1" type="text" value="212" label="Residue ID of atom 1"/> - <param name="name1" type="text" value="OE2" label="Atom name of atom 1"/> - <param name="segid2" type="text" value="HET" label="Segment ID of atom 2"/> - <param name="resid2" type="text" value="3" label="Residue ID of atom 2"/> - <param name="name2" type="text" value="C1" label="Atom name of atom 2"/> + + <conditional name="inps"> + <param argument="inps" type="select" label="Number of pairwise distances to calculate?"> + <option value="one" selected="true">One</option> + <option value="multiple">Multiple</option> + </param> + <when value="one"> + <param name="segid1" type="text" value="PRO" label="Segment ID of atom 1"/> + <param name="resid1" type="text" value="212" label="Residue ID of atom 1"/> + <param name="name1" type="text" value="OE2" label="Atom name of atom 1"/> + <param name="segid2" type="text" value="HET" label="Segment ID of atom 2"/> + <param name="resid2" type="text" value="3" label="Residue ID of atom 2"/> + <param name="name2" type="text" value="C1" label="Atom name of atom 2"/> + </when> + <when value="multiple"> + <param name="list1" type="data" format="text" label="Selection groups (list 1)"/> + <param name="list2" type="data" format="text" label="Selection groups (list 2)"/> + </when> + </conditional> + <param name="header" type="boolean" truevalue="--header" falsevalue="" label="Include header in output file"/> </inputs> <outputs> - <data format="tabular" name="output" label="Distance Analysis raw data"/> - <data format="png" name="distance_plot" label="Distance Analysis Plot"/> + <data format="tabular" name="output" label="Distance Analysis raw data"/> + <data format="png" name="distance_plot" label="Distance Analysis Plot"> + <filter>inps == 'one'</filter> + </data> </outputs> <tests> <test> @@ -43,6 +67,7 @@ <param name="segid2" value="HET"/> <param name="resid2" value="3"/> <param name="name2" value="C1"/> + <param name="header" value="false"/> <output name="output" file="Distance_Analysis_raw_data.tabular" /> </test> <test> @@ -53,11 +78,27 @@ <param name="segid2" value="SYSTEM"/> <param name="resid2" value="3"/> <param name="name2" value="C1"/> + <param name="header" value="false"/> <output name="output"> <assert_contents> <has_n_columns n="2" /> - <has_line_matching expression="0\s+3.8.*" /> - <has_line_matching expression="14\s+3.2.*" /> + <has_line_matching expression="0\s+3.893.*" /> + <has_line_matching expression="14\s+3.262.*" /> + </assert_contents> + </output> + </test> + <test> + <expand macro="tests_inputs_gmx"/> + <param name="inps" value="multiple"/> + <param name="list1" value="list1.txt"/> + <param name="list2" value="list2.txt"/> + <param name="header" value="true"/> + <output name="output"> + <assert_contents> + <has_n_columns n="2" /> + <has_line_matching expression="Frame\s+resid_212_and_name_OE2-resid_3_and_name_C1" /> + <has_line_matching expression="0\s+3.893.*" /> + <has_line_matching expression="14\s+3.262.*" /> </assert_contents> </output> </test> @@ -67,7 +108,9 @@ **What it does** -This tool calculates and plots the distance between the two atoms. +This tool calculates and plots the distance between pairs of atoms. Two modes are +available: single pair mode, where distances are calculated between two specified +atoms, and multiple pair mode, where two lists of atoms need to be provided. _____ @@ -78,9 +121,27 @@ - Trajectory file (DCD). - PDB file. - - Segment IDs, Residue IDs and names of two atoms to calculate distances. + +In single pair mode, segment IDs, residue IDs and names of two atoms are selected. +Note that a MDAnalysis 'segment' is a larger organizational unit, for example one +protein or all the solvent molecules or simply the whole system. + +In multiple pair mode, two files need to be uploaded, each with one or more atom +groups defined using the MDAnalysis atom selection, each on a new line. For example: + +:: -Note that a MDAnalysis 'segment' is a larger organizational unit, for example one protein or all the solvent molecules or simply the whole system. + resid 163 + resid 56 + resid 12 and type N + +All possible distances between the two sets of atom groups will be calculated. For +example, if List 1 contains 5 atoms and List 2 contains 8 atoms, the output file will +contain 40 columns, each with the distance between one group in List 1 and one group +in List 2. + +Note that in multiple pair mode, if the group has multiple atoms, the center of mass +will be used for the calculation. _____
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/distance_multiple.py Thu Feb 06 19:43:03 2020 -0500 @@ -0,0 +1,56 @@ +import argparse +import sys + +import MDAnalysis as mda +from MDAnalysis.analysis import distances + +import numpy as np + + +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('--list1', help='list 2') + parser.add_argument('--list2', help='list 2') + parser.add_argument('--output', help='output') + parser.add_argument('--header', dest='header', action='store_true') + return parser.parse_args() + + +args = parse_command_line(sys.argv) + +u = mda.Universe(args.istr, args.itraj, + topology_format=args.istrext, format=args.itrajext) + +list1 = np.loadtxt(args.list1, dtype=str, delimiter="\t", ndmin=1) +list2 = np.loadtxt(args.list2, dtype=str, delimiter="\t", ndmin=1) + +sel1 = [u.select_atoms(selection) for selection in list1] +sel2 = [u.select_atoms(selection) for selection in list2] + +d = np.empty((u.trajectory.n_frames, list1.shape[0], list2.shape[0]),) + +for ts in u.trajectory: + c_o_m1 = np.array([selection.center_of_mass() for selection in sel1]) + c_o_m2 = np.array([selection.center_of_mass() for selection in sel2]) + distances.distance_array(c_o_m1, c_o_m2, result=d[ts.frame]) + +d = np.hstack(( + np.array(np.reshape(np.arange( + 0, d.shape[0]), (d.shape[0], 1)), dtype=int), # add column w frame + np.reshape(d, (d.shape[0], d.shape[1] * d.shape[2])) +)) + +if args.header: + header = 'Frame\t' + '\t'.join( + ['-'.join(pair) for pair in zip( + sum([[n, ] * len(list2) for n in list1], []), + list(list2) * len(list1),)]).replace(' ', '_') +else: + header = '' + +np.savetxt(args.output, d, header=header, comments='', + fmt=['%d'] + ['%f'] * (d.shape[1] - 1), delimiter='\t')
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/distance_single.py Thu Feb 06 19:43:03 2020 -0500 @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import argparse +import sys + +import MDAnalysis as mda + +import matplotlib +matplotlib.use('Agg') # noqa +import matplotlib.pyplot as plt + +import numpy as np + + +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 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('--output', help='output') + parser.add_argument('--odistance_plot', help='odistance plot') + parser.add_argument('--header', dest='header', action='store_true') + return parser.parse_args() + + +args = parse_command_line(sys.argv) + +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) + +u = mda.Universe(args.istr, args.itraj, + topology_format=args.istrext, format=args.itrajext) +x = u.select_atoms(atom1) +y = u.select_atoms(atom2) + +with open(args.output, 'w') as f: + if args.header: + f.write('Frame\tDistance') + for t in u.trajectory: + r = x.positions - y.positions + d = np.linalg.norm(r) + f.write(str(t.frame) + '\t ') + f.write(str(d) + '\n') + +with open(args.output) as f: + g = [xtmp.strip() for xtmp in f] + data = [tuple(map(float, xtmp.split())) for xtmp in g[0:]] + time = [xtmp[0] for xtmp in data] + distance = [xtmp[1] for xtmp in data] + plt.plot(time, distance) + plt.xlabel('Frame No.') + plt.ylabel(r'Distance ($\AA$)') + plt.savefig(args.odistance_plot, format='png')