# HG changeset patch
# User rnateam
# Date 1453801107 18000
# Node ID 0b9aab6aaebfcddd4b5d870af2594c84660b4fab
# Parent  303f6402a035645cf84611fbc6f9ea5097490de0
Uploaded 16cfcafe8b42055c5dd64e62c42b82b455027a40

diff -r 303f6402a035 -r 0b9aab6aaebf convert_bc_to_binary_RY.py
--- a/convert_bc_to_binary_RY.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/convert_bc_to_binary_RY.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,13 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from string import maketrans
+from sys import stdout
+from Bio import SeqIO
+from Bio.Seq import Seq
+from Bio.Alphabet import IUPAC
+
 tool_description = """
 Convert standard nucleotides to IUPAC nucleotide codes used for binary barcodes.
 
@@ -19,19 +27,6 @@
 Status: Testing
 """
 
-import argparse
-import logging
-from string import maketrans
-from sys import stdout
-from Bio import SeqIO
-from Bio.Seq import Seq
-from Bio.Alphabet import IUPAC
-
-# # avoid ugly python IOError when stdout output is piped into another program
-# # and then truncated (such as piping to head)
-# from signal import signal, SIGPIPE, SIG_DFL
-# signal(SIGPIPE, SIG_DFL)
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
diff -r 303f6402a035 -r 0b9aab6aaebf convert_bc_to_binary_RY.xml
--- a/convert_bc_to_binary_RY.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/convert_bc_to_binary_RY.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -5,7 +5,7 @@
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python convert_bc_to_binary_RY.py --version</version_command>
+  <version_command>python $__tool_directory__/convert_bc_to_binary_RY.py --version</version_command>
   <command interpreter="python"><![CDATA[
 convert_bc_to_binary_RY.py
 #if $positional_1 and $positional_1 is not None:
diff -r 303f6402a035 -r 0b9aab6aaebf coords2clnt.py
--- a/coords2clnt.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/coords2clnt.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,15 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from sys import stdout
+from pybedtools import BedTool
+from pybedtools.featurefuncs import five_prime
+# avoid ugly python IOError when stdout output is piped into another program
+# and then truncated (such as piping to head)
+from signal import signal, SIGPIPE, SIG_DFL
+signal(SIGPIPE, SIG_DFL)
+
 tool_description = """
 Given coordinates of the aligned reads, calculate positions of the crosslinked
 nucleotides. Crosslinked nts are assumed to be one nt upstream of the 5'-end of
@@ -25,16 +35,6 @@
 Status: Testing
 """
 
-import argparse
-import logging
-from sys import stdout
-from pybedtools import BedTool
-from pybedtools.featurefuncs import five_prime
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
diff -r 303f6402a035 -r 0b9aab6aaebf coords2clnt.xml
--- a/coords2clnt.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/coords2clnt.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -8,7 +8,7 @@
   <stdio>
     <exit_code level="fatal" range="1:"/>
   </stdio>
-  <version_command>python coords2clnt.py --version</version_command>
+  <version_command>python $__tool_directory__/coords2clnt.py --version</version_command>
   <command interpreter="python"><![CDATA[coords2clnt.py
 #if $positional_1 and $positional_1 is not None:
 $positional_1
diff -r 303f6402a035 -r 0b9aab6aaebf extract_aln_ends.py
--- a/extract_aln_ends.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/extract_aln_ends.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,17 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from sys import stdout
+from shutil import rmtree
+from tempfile import mkdtemp
+from pybedtools import BedTool
+import pysam
+# avoid ugly python IOError when stdout output is piped into another program
+# and then truncated (such as piping to head)
+from signal import signal, SIGPIPE, SIG_DFL
+signal(SIGPIPE, SIG_DFL)
+
 tool_description = """
 Extract alignment ends from sam file.
 
@@ -32,25 +44,12 @@
 Status: Development
 """
 
-import argparse
-import logging
-from sys import stdout
-from shutil import rmtree
-from tempfile import mkdtemp
-from pybedtools import BedTool
-import pysam
-
 
 class DefaultsRawDescriptionHelpFormatter(argparse.ArgumentDefaultsHelpFormatter,
                                           argparse.RawDescriptionHelpFormatter):
     # To join the behaviour of RawDescriptionHelpFormatter with that of ArgumentDefaultsHelpFormatter
     pass
 
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
diff -r 303f6402a035 -r 0b9aab6aaebf extract_aln_ends.xml
--- a/extract_aln_ends.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/extract_aln_ends.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -5,7 +5,7 @@
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python extract_aln_ends.py --version</version_command>
+  <version_command>python $__tool_directory__/extract_aln_ends.py --version</version_command>
   <command interpreter="python"><![CDATA[
 extract_aln_ends.py
 #if $positional_1 and $positional_1 is not None:
diff -r 303f6402a035 -r 0b9aab6aaebf extract_bcs.py
--- a/extract_bcs.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/extract_bcs.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,15 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+import re
+from sys import stdout
+from Bio.SeqIO.QualityIO import FastqGeneralIterator
+# avoid ugly python IOError when stdout output is piped into another program
+# and then truncated (such as piping to head)
+from signal import signal, SIGPIPE, SIG_DFL
+signal(SIGPIPE, SIG_DFL)
+
 tool_description = """
 Exract barcodes from a FASTQ file according to a user-specified pattern.
 
@@ -19,17 +29,6 @@
 Status: Testing
 """
 
-import argparse
-import logging
-import re
-from sys import stdout
-from Bio.SeqIO.QualityIO import FastqGeneralIterator
-
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
diff -r 303f6402a035 -r 0b9aab6aaebf extract_bcs.xml
--- a/extract_bcs.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/extract_bcs.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -5,7 +5,7 @@
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python extract_bcs.py --version</version_command>
+  <version_command>python $__tool_directory__/extract_bcs.py --version</version_command>
   <command interpreter="python"><![CDATA[
 extract_bcs.py
 #if $positional_1 and $positional_1 is not None:
diff -r 303f6402a035 -r 0b9aab6aaebf merge_pcr_duplicates.py
--- a/merge_pcr_duplicates.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/merge_pcr_duplicates.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,18 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from sys import stdout
+import pandas as pd
+from subprocess import check_call
+from shutil import rmtree
+from tempfile import mkdtemp
+from os.path import isfile
+# avoid ugly python IOError when stdout output is piped into another program
+# and then truncated (such as piping to head)
+from signal import signal, SIGPIPE, SIG_DFL
+signal(SIGPIPE, SIG_DFL)
+
 tool_description = """
 Merge PCR duplicates according to random barcode library.
 
@@ -28,24 +41,6 @@
 Status: Testing
 """
 
-import argparse
-import logging
-from sys import stdout
-from Bio import SeqIO
-import pandas as pd
-
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
-
-def fasta_tuple_generator(fasta_iterator):
-    "Yields id, sequence tuples given an iterator over Biopython SeqIO objects."
-    for record in input_seq_iterator:
-        yield (record.id, str(record.seq))
-
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
@@ -56,16 +51,11 @@
     help="Path to bed6 file containing alignments.")
 parser.add_argument(
     "bclib",
-    help="Path to fasta barcode library.")
+    help="Path to fastq barcode library.")
 # optional arguments
 parser.add_argument(
     "-o", "--outfile",
     help="Write results to this file.")
-parser.add_argument(
-    "--fasta-library",
-    dest="fasta_library",
-    action="store_true",
-    help="Read random barcode library as fasta format.")
 # misc arguments
 parser.add_argument(
     "-v", "--verbose",
@@ -78,7 +68,7 @@
 parser.add_argument(
     '--version',
     action='version',
-    version='0.1.0')
+    version='0.2.0')
 
 args = parser.parse_args()
 
@@ -96,50 +86,65 @@
     logging.info("  outfile: '{}'".format(args.outfile))
 logging.info("")
 
-# load barcode library into dictionary
-input_handle = open(args.bclib, "rU")
-if args.fasta_library:
-    input_seq_iterator = SeqIO.parse(input_handle, "fasta")
-else:
-    input_seq_iterator = SeqIO.parse(input_handle, "fastq")
-bcs = pd.DataFrame.from_records(
-    data=fasta_tuple_generator(input_seq_iterator),
-    columns=["read_id", "bc"])
+# see if alignments are empty and the tool can quit
+n_alns = sum(1 for line in open(args.alignments))
+if n_alns == 0:
+    logging.warning("WARNING: Working on empty set of alignments, writing empty output.")
+    eventalnout = (open(args.outfile, "w") if args.outfile is not None else stdout)
+    eventalnout.close()
+    exit(0)
+
+# check input filenames
+if not isfile(args.bclib):
+    raise Exception("ERROR: barcode library '{}' not found.")
+if not isfile(args.alignments):
+    raise Exception("ERROR: alignments '{}' not found.")
+
+try:
+    tmpdir = mkdtemp()
+    logging.debug("tmpdir: " + tmpdir)
 
-# load alignments
-alns = pd.read_csv(
-    args.alignments,
-    sep="\t",
-    names=["chrom", "start", "stop", "read_id", "score", "strand"])
+    # prepare barcode library
+    syscall1 = "cat " + args.bclib + " | awk 'BEGIN{OFS=\"\\t\"}NR%4==1{gsub(/^@/,\"\"); id=$1}NR%4==2{bc=$1}NR%4==3{print id,bc}' | sort -k1,1 > " + tmpdir + "/bclib.csv"
+    check_call(syscall1, shell=True)
+
+    # prepare alinments
+    syscall2 = "cat " + args.alignments + " | awk -F \"\\t\" 'BEGIN{OFS=\"\\t\"}{split($4, a, \" \"); $4 = a[1]; print}'| sort -k4,4 > " + tmpdir + "/alns.csv"
+    check_call(syscall2, shell=True)
 
-# keep id parts up to first whitespace
-alns["read_id"] = alns["read_id"].str.split(' ').str.get(0)
+    # join barcode library and alignments
+    syscall3 = "join -1 1 -2 4 " + tmpdir + "/bclib.csv " + tmpdir + "/alns.csv " + " | awk 'BEGIN{OFS=\"\\t\"}{print $3,$4,$5,$2,$6,$7}' > " + tmpdir + "/bcalib.csv"
+    check_call(syscall3, shell=True)
 
-# combine barcode library and alignments
-bcalib = pd.merge(
-    bcs, alns,
-    on="read_id",
-    how="inner",
-    sort=False)
+    # get alignments combined with barcodes
+    bcalib = pd.read_csv(
+        tmpdir + "/bcalib.csv",
+        sep="\t",
+        names=["chrom", "start", "stop", "bc", "score", "strand"])
+finally:
+    logging.debug("removed tmpdir: " + tmpdir)
+    rmtree(tmpdir)
+
+# fail if alignments given but combined library is empty
 if bcalib.empty:
     raise Exception("ERROR: no common entries for alignments and barcode library found. Please check your input files.")
-n_alns = len(alns.index)
+
+# warn if not all alignments could be assigned a barcode
 n_bcalib = len(bcalib.index)
 if n_bcalib < n_alns:
     logging.warning(
-        "{} of {} alignments could not be associated with a random barcode.".format(
-            n_alns - n_bcalib, n_alns))
+        "{} of {} alignments could not be associated with a random barcode.".format(n_alns - n_bcalib, n_alns))
 
 # remove entries with barcodes that has uncalled base N
 bcalib_cleaned = bcalib.drop(bcalib[bcalib.bc.str.contains("N")].index)
 n_bcalib_cleaned = len(bcalib_cleaned)
-if n_bcalib_cleaned < n_bcalib:
-    msg = "{} of {} alignments had random barcodes containing uncalled bases and were dropped.".format(
-        n_bcalib - n_bcalib_cleaned, n_bcalib)
-    if n_bcalib_cleaned < (0.8 * n_bcalib):
-        logging.warning(msg)
-    else:
-        logging.info(msg)
+# if n_bcalib_cleaned < n_bcalib:
+#     msg = "{} of {} alignments had random barcodes containing uncalled bases and were dropped.".format(
+#         n_bcalib - n_bcalib_cleaned, n_bcalib)
+#     if n_bcalib_cleaned < (0.8 * n_bcalib):
+#         logging.warning(msg)
+#     else:
+#         logging.info(msg)
 
 # count and merge pcr duplicates
 # grouping sorts by keys, so the ouput will be properly sorted
diff -r 303f6402a035 -r 0b9aab6aaebf merge_pcr_duplicates.xml
--- a/merge_pcr_duplicates.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/merge_pcr_duplicates.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -1,11 +1,15 @@
-<tool id="merge_pcr_duplicates.py" name="Merge PCR duplicates" version="0.1.0">
+<tool id="merge_pcr_duplicates.py" name="Merge PCR duplicates" version="0.2.0">
   <description>according to random barcode library.</description>
   <macros>
     <import>macros.xml</import>
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python merge_pcr_duplicates.py --version</version_command>
+  <requirements>
+      <requirement type="package" version="4.1.0">gnu_awk</requirement>
+      <requirement type="package" version="8.22">gnu_coreutils</requirement>
+  </requirements>
+  <version_command>python $__tool_directory__/merge_pcr_duplicates.py --version</version_command>
   <command interpreter="python"><![CDATA[merge_pcr_duplicates.py
 #if $positional_1 and $positional_1 is not None:
 $positional_1
@@ -18,7 +22,7 @@
 > $default]]></command>
   <inputs>
     <param area="false" label="bed6 file containing alignments." name="positional_1" type="data" format="bed"/>
-    <param area="false" label="fasta barcode library." name="positional_2" type="data" format="fastq"/>
+    <param area="false" label="fastaq barcode library." name="positional_2" type="data" format="fastq"/>
   </inputs>
   <outputs>
     <data format="bed" hidden="false" name="default"/>
diff -r 303f6402a035 -r 0b9aab6aaebf remove_tail.py
--- a/remove_tail.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/remove_tail.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,14 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from sys import stdout
+from Bio.SeqIO.QualityIO import FastqGeneralIterator
+# avoid ugly python IOError when stdout output is piped into another program
+# and then truncated (such as piping to head)
+from signal import signal, SIGPIPE, SIG_DFL
+signal(SIGPIPE, SIG_DFL)
+
 tool_description = """
 Remove a certain number of nucleotides from the 3'-tails of sequences in fastq
 format.
@@ -18,16 +27,6 @@
 Status: Testing
 """
 
-import argparse
-import logging
-from sys import stdout
-from Bio.SeqIO.QualityIO import FastqGeneralIterator
-
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
 # parse command line arguments
 parser = argparse.ArgumentParser(description=tool_description,
                                  epilog=epilog,
diff -r 303f6402a035 -r 0b9aab6aaebf remove_tail.xml
--- a/remove_tail.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/remove_tail.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -5,7 +5,7 @@
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python remove_tail.py --version</version_command>
+  <version_command>python $__tool_directory__/remove_tail.py --version</version_command>
   <command interpreter="python"><![CDATA[remove_tail.py
 #if $positional_1 and $positional_1 is not None:
 $positional_1
diff -r 303f6402a035 -r 0b9aab6aaebf rm_spurious_events.pl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rm_spurious_events.pl	Tue Jan 26 04:38:27 2016 -0500
@@ -0,0 +1,87 @@
+#!/usr/local/perl/bin/perl
+use feature ':5.10';
+use strict 'vars';
+use warnings;
+use Getopt::Long;
+use Pod::Usage;
+use List::Util qw/ min max /;
+use POSIX qw/ceil floor/;
+use File::Temp qw(tempdir);
+use File::Basename;
+
+=head1 NAME
+
+filterSquashedReads.pl --frac_max FLOAT
+
+=head1 SYNOPSIS
+
+this script reads sites after pcr duplicate removal.
+for all crosslinking sites sharing start and stop coordinates, the maximum number
+of squashed reads is determined.
+alignments having less than frac_max * max reads are discarded
+
+assumes bed entries to be sorted chr,strand,start,stop,score with score descending
+
+Options:
+
+    --frac_max  filter out alignments supported by less reads than this fraction of the maximum number of reads per position
+    -debug      enable debug output
+    -help       brief help message
+    -man        full documentation
+
+=head1 DESCRIPTION
+
+=cut
+
+###############################################################################
+# parse command line options
+###############################################################################
+my $frac_max = 0.1;
+my $help;
+my $man;
+my $debug;
+my $result = GetOptions (   "frac_max=f" => \$frac_max,
+                            "help"	=> \$help,
+							"man"	=> \$man,
+							"debug" => \$debug);
+pod2usage(-exitstatus => 1, -verbose => 1) if $help;
+pod2usage(-exitstatus => 0, -verbose => 2) if $man;
+($result) or pod2usage(2);
+
+###############################################################################
+# main
+###############################################################################
+my $current_chr = '';
+my $current_start = -1;
+my $current_stop = -1;
+my $current_strand = '';
+my $current_max = -1;
+my $current_threshold = -1;
+
+while (<>) {
+    my ($chr, $start, $stop, $id, $count, $strand) = split("\t");
+    # my ($count, undef, $start, $chr, $strand, $stop) = split("\t");
+
+    if ($current_start != $start or
+        $current_stop != $stop or
+        $current_chr ne $chr or
+        $current_strand ne $strand) {
+        # if this is the first occourence of these coordinates
+        # this must be the new maximum
+        # save new state
+        $current_start = $start;
+        $current_stop = $stop;
+        $current_chr = $chr;
+        $current_strand = $strand;
+        # print record and record maximum
+        $current_max = $count;
+        $current_threshold = $count*$frac_max;
+        print $_;
+        $debug and say "new threshold ${current_threshold} @ $chr $start $stop $strand $count";
+    } elsif ($count >= $current_threshold) {
+        # if it is not the first occourence, evaluate threshold and print if valid
+        print $_;
+    } else {
+        $debug and say "below threshold @ $chr $start $stop $strand $count";
+    }
+}
diff -r 303f6402a035 -r 0b9aab6aaebf rm_spurious_events.py
--- a/rm_spurious_events.py	Sat Dec 19 06:16:22 2015 -0500
+++ b/rm_spurious_events.py	Tue Jan 26 04:38:27 2016 -0500
@@ -1,5 +1,10 @@
 #!/usr/bin/env python
 
+import argparse
+import logging
+from subprocess import check_call
+import os
+
 tool_description = """
 Remove spurious events originating from errors in random sequence tags.
 
@@ -7,8 +12,6 @@
 of events the maximum number of PCR duplicates is determined. All events that
 are supported by less than 10 percent of this maximum count are removed.
 
-By default output is written to stdout.
-
 Input:
 * bed6 file containing crosslinking events with score field set to number of PCR
   duplicates
@@ -19,7 +22,7 @@
 
 Example usage:
 - remove spurious events from spurious.bed and write results to file cleaned.bed
-rm_spurious_events.py spurious.bed --out cleaned.bed
+rm_spurious_events.py spurious.bed --oufile cleaned.bed
 """
 
 epilog = """
@@ -30,89 +33,74 @@
 Status: Testing
 """
 
-import argparse
-import logging
-from sys import stdout
-import pandas as pd
-
 
 class DefaultsRawDescriptionHelpFormatter(argparse.ArgumentDefaultsHelpFormatter,
                                           argparse.RawDescriptionHelpFormatter):
     # To join the behaviour of RawDescriptionHelpFormatter with that of ArgumentDefaultsHelpFormatter
     pass
 
-# avoid ugly python IOError when stdout output is piped into another program
-# and then truncated (such as piping to head)
-from signal import signal, SIGPIPE, SIG_DFL
-signal(SIGPIPE, SIG_DFL)
-
-# parse command line arguments
-parser = argparse.ArgumentParser(description=tool_description,
-                                 epilog=epilog,
-                                 formatter_class=DefaultsRawDescriptionHelpFormatter)
-# positional arguments
-parser.add_argument(
-    "events",
-    help="Path to bed6 file containing alignments.")
-# optional arguments
-parser.add_argument(
-    "-o", "--outfile",
-    help="Write results to this file.")
-parser.add_argument(
-    "-t", "--threshold",
-    type=float,
-    default=0.1,
-    help="Threshold for spurious event removal."
-)
-# misc arguments
-parser.add_argument(
-    "-v", "--verbose",
-    help="Be verbose.",
-    action="store_true")
-parser.add_argument(
-    "-d", "--debug",
-    help="Print lots of debugging information",
-    action="store_true")
-parser.add_argument(
-    '--version',
-    action='version',
-    version='0.1.0')
 
-args = parser.parse_args()
-
-if args.debug:
-    logging.basicConfig(level=logging.DEBUG, format="%(asctime)s - %(filename)s - %(levelname)s - %(message)s")
-elif args.verbose:
-    logging.basicConfig(level=logging.INFO, format="%(filename)s - %(levelname)s - %(message)s")
-else:
-    logging.basicConfig(format="%(filename)s - %(levelname)s - %(message)s")
-logging.info("Parsed arguments:")
-logging.info("  alignments: '{}'".format(args.events))
-logging.info("  threshold: '{}'".format(args.threshold))
-if args.outfile:
-    logging.info("  outfile: enabled writing to file")
-    logging.info("  outfile: '{}'".format(args.outfile))
-logging.info("")
+def main():
+    # parse command line arguments
+    parser = argparse.ArgumentParser(description=tool_description,
+                                     epilog=epilog,
+                                     formatter_class=DefaultsRawDescriptionHelpFormatter)
+    # positional arguments
+    parser.add_argument(
+        "events",
+        help="Path to bed6 file containing alignments.")
+    # optional arguments
+    parser.add_argument(
+        "-o", "--outfile",
+        required=True,
+        help="Write results to this file.")
+    parser.add_argument(
+        "-t", "--threshold",
+        type=float,
+        default=0.1,
+        help="Threshold for spurious event removal."
+    )
+    # misc arguments
+    parser.add_argument(
+        "-v", "--verbose",
+        help="Be verbose.",
+        action="store_true")
+    parser.add_argument(
+        "-d", "--debug",
+        help="Print lots of debugging information",
+        action="store_true")
+    parser.add_argument(
+        '--version',
+        action='version',
+        version='0.1.0')
 
-# check threshold parameter value
-if args.threshold < 0 or args.threshold > 1:
-    raise ValueError("Threshold must be in [0,1].")
-
-# load alignments
-alns = pd.read_csv(
-    args.events,
-    sep="\t",
-    names=["chrom", "start", "stop", "read_id", "score", "strand"])
+    args = parser.parse_args()
 
-# remove all alignments that not enough PCR duplicates with respect to
-# the group maximum
-grouped = alns.groupby(['chrom', 'start', 'stop', 'strand'], group_keys=False)
-alns_cleaned = grouped.apply(lambda g: g[g["score"] >= args.threshold * g["score"].max()])
+    if args.debug:
+        logging.basicConfig(level=logging.DEBUG, format="%(asctime)s - %(filename)s - %(levelname)s - %(message)s")
+    elif args.verbose:
+        logging.basicConfig(level=logging.INFO, format="%(filename)s - %(levelname)s - %(message)s")
+    else:
+        logging.basicConfig(format="%(filename)s - %(levelname)s - %(message)s")
+    logging.info("Parsed arguments:")
+    logging.info("  alignments: '{}'".format(args.events))
+    logging.info("  threshold: '{}'".format(args.threshold))
+    if args.outfile:
+        logging.info("  outfile: enabled writing to file")
+        logging.info("  outfile: '{}'".format(args.outfile))
+    logging.info("")
 
-# write coordinates of crosslinking event alignments
-alns_cleaned_out = (open(args.outfile, "w") if args.outfile is not None else stdout)
-alns_cleaned.to_csv(
-    alns_cleaned_out,
-    columns=['chrom', 'start', 'stop', 'read_id', 'score', 'strand'],
-    sep="\t", index=False, header=False)
-alns_cleaned_out.close()
+    # check threshold parameter value
+    if args.threshold < 0 or args.threshold > 1:
+        raise ValueError("Threshold must be in [0,1].")
+
+    if not os.path.isfile(args.events):
+        raise Exception("ERROR: file '{}' not found.")
+
+    # prepare barcode library
+    syscall = "cat " + args.events + " | sort -k1,1V -k6,6 -k2,2n -k3,3 -k5,5nr | perl " + os.path.dirname(os.path.realpath(__file__)) + "/rm_spurious_events.pl --frac_max " + str(args.threshold) + "| sort -k1,1V -k2,2n -k3,3n -k6,6 -k4,4 -k5,5nr > " + args.outfile
+    check_call(syscall, shell=True)
+
+
+if __name__ == "__main__":
+    main()
diff -r 303f6402a035 -r 0b9aab6aaebf rm_spurious_events.xml
--- a/rm_spurious_events.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/rm_spurious_events.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -5,7 +5,11 @@
   </macros>
   <expand macro="requirements" />
   <expand macro="stdio" />
-  <version_command>python rm_spurious_events.py --version</version_command>
+  <requirements>
+    <requirement type="package" version="5.18">perl</requirement>
+    <requirement type="package" version="8.22">gnu_coreutils</requirement>
+  </requirements>
+  <version_command>python $__tool_directory__/rm_spurious_events.py --version</version_command>
   <command interpreter="python"><![CDATA[rm_spurious_events.py
 #if $positional_1 and $positional_1 is not None:
 $positional_1
@@ -14,7 +18,7 @@
 #if $threshold and $threshold is not None:
 --threshold $threshold
 #end if
-> $default]]></command>
+--outfile $default]]></command>
   <inputs>
     <param area="false" label="bed6 file containing alignments." name="positional_1" type="data" format="bed"/>
     <param help="(--threshold)" label="Threshold for spurious event removal." name="threshold" optional="true" type="float" value="0.1"/>
diff -r 303f6402a035 -r 0b9aab6aaebf tool_dependencies.xml
--- a/tool_dependencies.xml	Sat Dec 19 06:16:22 2015 -0500
+++ b/tool_dependencies.xml	Tue Jan 26 04:38:27 2016 -0500
@@ -7,18 +7,27 @@
         <repository name="package_biopython_1_65" owner="biopython"/>
     </package> -->
     <package name="pandas" version="0.16">
-        <repository changeset_revision="8766fb73e88f" name="package_python_2_7_pandas_0_16" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+        <repository changeset_revision="67abd81bd3ac" name="package_python_2_7_pandas_0_16" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
     <!-- <package name="pybedtools" version="0.7.4">
         <repository name="package_python_2_7_pybedtools_0_7_4" owner="iuc"/>
     </package> -->
     <package name="pybedtools" version="0.6.9">
-        <repository changeset_revision="7d07b5c08b8f" name="package_python_2_7_pybedtools_0_6_9" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+        <repository changeset_revision="c4641c3a869f" name="package_python_2_7_pybedtools_0_6_9" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
     <package name="pysam" version="0.8.3">
-        <repository changeset_revision="f88a5f4ac6c1" name="package_python_2_7_pysam_0_8_3" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+        <repository changeset_revision="7ac80143c68d" name="package_python_2_7_pysam_0_8_3" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
     <package name="flexbar" version="2.5">
         <repository changeset_revision="45d1b8373762" name="package_flexbar_2_5" owner="rnateam" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
+    <package name="gnu_awk" version="4.1.0">
+        <repository changeset_revision="0f0bdef2f686" name="package_gnu_awk_4_1_0" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+    </package>
+    <package name="gnu_coreutils" version="8.22">
+        <repository changeset_revision="8b60cf3e0c07" name="package_gnu_coreutils_8_22" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+    </package>
+    <package name="perl" version="5.18">
+        <repository changeset_revision="6f144bd786a8" name="package_perl_5_18" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+    </package>
 </tool_dependency>