target_fragments.py

#!/bin/env python

import sys
import getopt
from add_fragment_to_bam import bed_reader


def main(hin=sys.stdin, hout=sys.stdout):
    """
    The main program loop
    """
    previous = None
    for current in bed_reader(hin):

        # first cycle and same chromosome
        if previous is None or current[0] != previous[0]:
            previous = (current[1], 0, 1, [""])

        # set the fragment coordinate
        start = previous[2]
        end = current[1]

        # skip fragments that overlap with previous
        # fragment (takes the maximal fragment)
        if start > end:
            continue

        # label the fragment
        label = "%s-%s" % (previous[3][0], current[3][0])

        # write the frament
        hout.write("%s\t%d\t%d\t%s\n" % (current[0], start, end, label))

        # prepare for next cycle
        previous = current


if __name__ == "__main__":

    def usage(message="", error=None):
        """
        Print the usage information
        """
        sys.stdout.write("""
Message %(message)s

Called as %(commandline)s

Usage:
    python %(tool)s -i [input bed file] -o [output bed file] -s [sequence] 

Options:    
    -i/--input      [file]  a bed file with the restriction sites
    -o/--output     [file]  a bed file with the restriction fragments            
    -h/--help       []       print this message

""" % {
    "message": message,
    "tool": sys.argv[0],
    "commandline": " ".join(sys.argv)
})
        #
        if error is not None:
            sys.exit(error)

    # main variables
    fin = sys.stdin
    fout = sys.stdout

    # option parsing
    shortopt = "i:o:h"
    longopt = ["input=", "output=", "help"]

    try:
        opts, args = getopt.getopt(sys.argv[1:], shortopt, longopt)
        for o, a in opts:
            if o in ("-i", "--input"):
                fin = open(a, "r")
            elif o in ("-o", "--output"):
                fout = open(a, "w")
            elif o in ('-h', '--help'):
                usage("Help was asked", error=0)
    except getopt.GetoptError as err:
        usage(str(err), error=2)

    # run the main program loop
    main(hin=fin, hout=fout)

    if not fin.closed and fin is not sys.stdin:
        fin.close()
    if not fout.closed and fout is not sys.stdout:
        fout.close()