Mercurial > repos > chemteam > mdanalysis_distance
changeset 3:489b25966bb9 draft
planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 3ff06e3182c3a1546ea0a3b29e0d4383e12169e1
| author | chemteam |
|---|---|
| date | Wed, 03 Apr 2019 15:46:56 -0400 |
| parents | 46892d756cec |
| children | 312f912de69d |
| files | dihedrals.py hbonds.py macros.xml pca_cosine.py ramachandran_plots.py test-data/Angle_Analysis_raw_data.tabular test-data/Dihedral_analysis_raw_data.tabular test-data/Distance_Analysis_raw_data.tabular test-data/Ramachandran_Plot_raw_data.tabular |
| diffstat | 9 files changed, 303 insertions(+), 36 deletions(-) [+] |
line wrap: on
line diff
--- a/dihedrals.py Sun Jan 13 03:22:41 2019 -0500 +++ b/dihedrals.py Wed Apr 03 15:46:56 2019 -0400 @@ -59,12 +59,13 @@ u = mda.Universe(args.ipdb, args.idcd, topology_format="PDB", format="DCD") data = np.array([(u.trajectory.frame, psi(u)) for ts in u.trajectory]) frame, psi = data.T +PSI = np.concatenate(psi, axis=0) -zip(frame, psi) +zip(frame, PSI) with open(args.output, 'w') as f: writer = csv.writer(f, delimiter='\t') - writer.writerows(zip(frame, psi)) + writer.writerows(zip(frame, PSI)) with open(args.output) as f: g = [xtmp.strip() for xtmp in f]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/hbonds.py Wed Apr 03 15:46:56 2019 -0400 @@ -0,0 +1,69 @@ +#!/usr/bin/env python + +import argparse +import csv +import sys + +import MDAnalysis.analysis.hbonds + +import pandas as pd + + +def parse_command_line(argv): + parser = argparse.ArgumentParser() + parser.add_argument('--idcd', help='input dcd') + parser.add_argument('--ipdb', help='input pdb') + parser.add_argument('--isegid1', help='segid 1') + parser.add_argument('--isegid2', help='segid 2') + parser.add_argument('--idistance', help='cutoff distance') + parser.add_argument('--iangle', help='ctoff angle') + parser.add_argument('--output', help='output') + parser.add_argument('--ofreq_output', help='frequency output') + parser.add_argument('--onumber_output', help='number of hbond output') + parser.add_argument('--otime_output', help='time steps output') + return parser.parse_args() + + +args = parse_command_line(sys.argv) + +selection1 = "segid %s" % args.isegid1 +selection2 = "segid %s" % args.isegid2 +distance = float(args.idistance) +angle = float(args.iangle) + +u = MDAnalysis.Universe( + args.ipdb, args.idcd, topology_format="PDB", format="DCD") + +h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis( + u, selection1, selection2, distance=distance, angle=angle) +h.run() +h.generate_table() + +df = pd.DataFrame.from_records(h.table) +df.to_csv(args.output, sep='\t') + +t1 = list(h.count_by_type()) +t2 = list(h.count_by_time()) +t3 = list(h.timesteps_by_type()) + +with open(args.ofreq_output, 'w') as f: + f.write("donor_index\tacceptor_index\t\ + donor_resname\tdonor_resid\tdonor_atom\t\ + hydrogen_atom\tacceptor_reansme\tacceptor_resid\t\ + acceptor_atom\tfrequency\n") + writer = csv.writer(f, delimiter='\t') + writer.writerows(t1) + + +with open(args.onumber_output, 'w') as f1: + f1.write("time_step\tno_of_h_bonds\n") + writer = csv.writer(f1, delimiter='\t') + writer.writerows(t2) + +with open(args.otime_output, 'w') as f2: + f2.write("donor_index\tacceptor_index\t\ + donor_resname\tdonor_resid\tdonor_atom\t\ + hydrogen_atom\tacceptor_reansme\tacceptor_resid\t\ + acceptor_atom\ttime_step\n") + writer = csv.writer(f2, delimiter='\t') + writer.writerows(t3)
--- a/macros.xml Sun Jan 13 03:22:41 2019 -0500 +++ b/macros.xml Wed Apr 03 15:46:56 2019 -0400 @@ -1,8 +1,9 @@ <macros> - <token name="@VERSION@">0.18</token> + <token name="@VERSION@">0.19</token> <xml name="requirements"> <requirements> - <requirement type="package" version="0.18.0">mdanalysis</requirement> + <requirement type="package" version="0.19.2">mdanalysis</requirement> + <yield/> </requirements> </xml> <xml name="analysis_inputs"> @@ -10,6 +11,18 @@ <param format="pdb" name="pdbin" type="data" label="pdb input"/> <yield/> </xml> + <xml name="sanitizer"> + <sanitizer invalid_char=""> + <valid initial="string.ascii_letters,string.digits"/> + </sanitizer> + <yield/> + </xml> + <xml name="sanitizer_resids"> + <sanitizer invalid_char=""> + <valid initial="string.digits"/> + </sanitizer> + <yield/> + </xml> <xml name="tests_inputs"> <param name="dcdin" value="test.dcd" /> <param name="pdbin" value="test.pdb" /> @@ -17,7 +30,9 @@ </xml> <xml name="citations"> <citations> + <citation type="doi">10.1093/bioinformatics/btz107</citation> <citation type="doi">10.1002/jcc.21787</citation> + <yield/> </citations> </xml> </macros>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pca_cosine.py Wed Apr 03 15:46:56 2019 -0400 @@ -0,0 +1,46 @@ +#!/usr/bin/env python + +import argparse +import csv +import sys + +import MDAnalysis as mda +import MDAnalysis.analysis.pca as pca + +import numpy as np + + +def parse_command_line(argv): + parser = argparse.ArgumentParser() + parser.add_argument('--idcd', help='input dcd') + parser.add_argument('--ipdb', help='input pdb') + parser.add_argument('--icomponents', help='number of principle components') + parser.add_argument('--iindex', help='index of the PC') + parser.add_argument('--output', help='output') + parser.add_argument('--cosout', help='cosine output') + return parser.parse_args() + + +args = parse_command_line(sys.argv) + +u = mda.Universe(args.ipdb, args.idcd, topology_format="PDB", format="DCD") + +components = int(args.icomponents) +pca_index = int(args.iindex) + +PSF_pca = pca.PCA(u, select='backbone') +PSF_pca.run() +n_pcs = np.where(PSF_pca.cumulated_variance > 0.95)[0][0] +atomgroup = u.select_atoms('backbone') + +pca_space = PSF_pca.transform(atomgroup, n_components=components) +cosine = mda.analysis.pca.cosine_content(pca_space, pca_index) + +PCA = list(pca_space) + +with open(args.output, 'w') as f: + writer = csv.writer(f, delimiter='\t') + writer.writerows(PCA) + +with open(args.cosout, 'w') as f1: + f1.write(str(cosine))
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ramachandran_plots.py Wed Apr 03 15:46:56 2019 -0400 @@ -0,0 +1,121 @@ +#!/usr/bin/env python + +import argparse +import csv +import sys +from collections import namedtuple + +import MDAnalysis as mda +from MDAnalysis.lib.distances import calc_dihedrals + +import matplotlib +matplotlib.use('Agg') # noqa +import matplotlib.pyplot as plt + +import numpy as np + +import seaborn as sns + + +def parse_command_line(argv): + parser = argparse.ArgumentParser() + parser.add_argument('--idcd', help='input dcd') + parser.add_argument('--ipdb', help='input pdb') + 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('--output', help='output') + parser.add_argument('--oramachandran_plot', help='dihedral plot') + return parser.parse_args() + + +args = parse_command_line(sys.argv) + +Dihedral = namedtuple( + 'Dihedral', ['atom1', 'atom2', 'atom3', 'atom4']) + +# 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) + + +def calc_torsion(dihedral): + """atom 1 -4 are valid atom selections. torsion in degrees is returned""" + A = u.select_atoms(dihedral.atom1).positions + B = u.select_atoms(dihedral.atom2).positions + C = u.select_atoms(dihedral.atom3).positions + D = u.select_atoms(dihedral.atom4).positions + + dihe = calc_dihedrals(A, B, C, D) + return np.rad2deg(dihe) + + +u = mda.Universe(args.ipdb, args.idcd, topology_format="PDB", format="DCD") + +phi_trajdata = np.array( + [(u.trajectory.frame, calc_torsion(dihe_phi)) for ts in u.trajectory]) +psi_trajdata = np.array( + [(u.trajectory.frame, calc_torsion(dihe_psi)) for ts in u.trajectory]) + +phi_frame, phi_series = phi_trajdata.T +psi_frame, psi_series = psi_trajdata.T + +phi_series = np.concatenate(phi_series, axis=0) +psi_series = np.concatenate(psi_series, axis=0) + +zip(phi_frame, phi_series, psi_series) + +with open(args.output, 'w') as f: + writer = csv.writer(f, delimiter='\t') + 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.ax_joint.set_xlim(-180, 180) + h.ax_joint.set_ylim(-180, 180) + plt.savefig(args.oramachandran_plot, format='png')
--- a/test-data/Angle_Analysis_raw_data.tabular Sun Jan 13 03:22:41 2019 -0500 +++ b/test-data/Angle_Analysis_raw_data.tabular Wed Apr 03 15:46:56 2019 -0400 @@ -1,15 +1,15 @@ -0.0 70.84918975830078 -1.0 70.97154998779297 -2.0 70.3709716796875 -3.0 70.12692260742188 -4.0 71.15946197509766 -5.0 71.91361999511719 -6.0 71.91268920898438 -7.0 71.97418212890625 -8.0 72.15215301513672 -9.0 72.30829620361328 -10.0 71.52912139892578 -11.0 81.56575775146484 -12.0 73.46330261230469 -13.0 67.04185485839844 -14.0 70.07105255126953 +0.0 70.84919180880273 +1.0 70.97155495136695 +2.0 70.37097279938007 +3.0 70.12692388497567 +4.0 71.15946406922332 +5.0 71.91362565497833 +6.0 71.91268728348935 +7.0 71.97418402125982 +8.0 72.1521609490865 +9.0 72.3083065219282 +10.0 71.5291313235259 +11.0 81.56576116363512 +12.0 73.46330915394758 +13.0 67.04185741201445 +14.0 70.0710549275429
--- a/test-data/Dihedral_analysis_raw_data.tabular Sun Jan 13 03:22:41 2019 -0500 +++ b/test-data/Dihedral_analysis_raw_data.tabular Wed Apr 03 15:46:56 2019 -0400 @@ -1,15 +1,15 @@ -0.0 -61.75347431106646 -1.0 -63.57940223575045 -2.0 -64.3468375822089 -3.0 -64.2932006449013 -4.0 -66.17815780069928 -5.0 -64.57418410935712 -6.0 -63.95470182210953 -7.0 -63.215898370961455 -8.0 -63.05227933072821 -9.0 -63.350881174296354 -10.0 -64.12889787645014 -11.0 -59.11099982465991 -12.0 -76.10678081593274 -13.0 -74.55530125465415 -14.0 -71.01850912317343 +0 -61.75347431106646 +1 -63.57940223575045 +2 -64.3468375822089 +3 -64.2932006449013 +4 -66.17815780069928 +5 -64.57418410935712 +6 -63.95470182210953 +7 -63.215898370961455 +8 -63.05227933072821 +9 -63.350881174296354 +10 -64.12889787645014 +11 -59.11099982465991 +12 -76.10678081593274 +13 -74.55530125465415 +14 -71.01850912317343
--- a/test-data/Distance_Analysis_raw_data.tabular Sun Jan 13 03:22:41 2019 -0500 +++ b/test-data/Distance_Analysis_raw_data.tabular Wed Apr 03 15:46:56 2019 -0400 @@ -1,9 +1,9 @@ 0 3.8900816 1 3.7835152 2 3.7693398 -3 3.7851017 +3 3.7851014 4 3.6331043 -5 3.6525044 +5 3.6525042 6 3.684632 7 3.7285593 8 3.7090275
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/Ramachandran_Plot_raw_data.tabular Wed Apr 03 15:46:56 2019 -0400 @@ -0,0 +1,15 @@ +0 -144.50369390751567 97.28315803510958 +1 -140.8182020539229 98.03359173250509 +2 -140.72994628351444 97.24540656140032 +3 -139.74848490172678 96.50804547826581 +4 -134.87886646306 95.11565955448816 +5 -133.72523437825183 98.27907111162767 +6 -132.95986510869454 98.00598700434969 +7 -132.83978048218802 97.90284983563569 +8 -132.446606581017 97.95362964393432 +9 -132.8171707044048 98.22639981540117 +10 -129.18708154562577 93.19031937895272 +11 -138.70826130463465 99.35902431554985 +12 -134.31372026825582 86.39732109628024 +13 -135.62845675103858 89.17557531169159 +14 -151.1966272020228 101.09732806451846