Mercurial > repos > chemteam > mdanalysis_angle

changeset 3:ad135cf42274 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 3ff06e3182c3a1546ea0a3b29e0d4383e12169e1
author chemteam
date Wed, 03 Apr 2019 15:47:56 -0400
parents 4c93f7541218
children e7d0075052c9
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:24:24 2019 -0500
+++ b/dihedrals.py	Wed Apr 03 15:47: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:47: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:24:24 2019 -0500
+++ b/macros.xml	Wed Apr 03 15:47: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:47: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:47: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:24:24 2019 -0500
+++ b/test-data/Angle_Analysis_raw_data.tabular	Wed Apr 03 15:47: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:24:24 2019 -0500
+++ b/test-data/Dihedral_analysis_raw_data.tabular	Wed Apr 03 15:47: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:24:24 2019 -0500
+++ b/test-data/Distance_Analysis_raw_data.tabular	Wed Apr 03 15:47: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:47: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