annotate hbonds.py @ 0:ebeb72aa39a8 draft

planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 3ff06e3182c3a1546ea0a3b29e0d4383e12169e1

#!/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)