#!/usr/bin/python

import sys
from optparse import OptionParser
import os
import pysam

class CNVOutput(object):

    def __init__(self, chrom1, breakStart1, breakEnd1, chrom2, breakStart2, breakEnd2, ID, score, strand1, strand2, queryStart1, queryEnd1, queryStart2, queryEnd2, minNonOverlap, queryLength, qualScores, variant, variant_type, unaccounted_for_sequence, event_length,i ):
        self.chrom1 = chrom1
        self.breakStart1 = breakStart1 
        self.breakEnd1 = breakEnd1 
        self.chrom2 = chrom2
        self.breakStart2 = breakStart2
        self.breakEnd2 = breakEnd2
        self.ID = ID
        self.score = score
        self.strand1 = strand1
        self.strand2 = strand2
        self.queryStart1 =queryStart1
        self.queryEnd1 = queryEnd1
        self.queryStart2 = queryStart2 
        self.queryEnd2 = queryEnd2 
        self.minNonOverlap = minNonOverlap 
        self.queryLength = queryLength 
        self.qualScores = qualScores 
        self.variant = variant
        self.variant_type = variant_type 
        self.unaccounted_for_sequence = unaccounted_for_sequence 
        self.event_length = event_length
        self.event_id = self.chrom1+":"+self.breakEnd1 +"-" +self.breakStart2 + "_" + self.variant_type + "_" + str(i)

    def __str__(self):
        return '\t'.join(map(str, [self.chrom1, self.breakStart1, self.breakEnd1, self.chrom2, self.breakStart2, self.breakEnd2, self.ID, self.score, self.strand1, self.strand2, self.queryStart1, self.queryEnd1, self.queryStart2, self.queryEnd2  , self.minNonOverlap, self.queryLength, self.qualScores, self.variant_type, self.unaccounted_for_sequence, self.event_length,self.event_id]))


def _extract_breakpoints(region, output_prefix, slop=300, fasta_input_file="/Users/smilefreak/Programming/OpenSource/MyGitHub/long_read_pipeline/scripts/test/RM_scaffolds.fasta", second_event=None):
    """
        Extract the breakpoint sequence.
    """
    assembly = pysam.Fastafile(fasta_input_file)
    seq_fastq = assembly.fetch(region=region.ID)
    dirname = os.path.dirname(output_prefix)
    base_name_file = os.path.basename(fasta_input_file).split(".")[0]
    dirname = os.path.join(dirname, base_name_file) 
    reference = pysam.Fastafile("/data/james/assembly_alignment/pipeline_run/input_data/BY.fasta")
    ref_fasta = reference.fetch(region=region.chrom1)
    try:
        os.makedirs(os.path.join(dirname,"breakpoints"))
    except OSError:
        pass 
    with open(os.path.join(dirname,"breakpoints",region.event_id + ".brk1.fasta"), "w") as brk1_f:
        with open(os.path.join(dirname,"breakpoints",region.event_id + ".brk2.fasta"), "w") as brk2_f: 
            header =  region.event_id  
            seq = ""
            if region.strand1 == region.strand2: 
                start_read = region.queryEnd1  
                end_read = region.queryStart2
                if end_read >start_read:
                    #
                    # must be [   DUP    ]------[    DUP      ]
                    #
                    # We must have duplication situation, so just extract the sequence
                    # Could be no insertion.
                    if "duplication" in region.variant_type: 
                        h1 = "_" +region.chrom1 +":"+str(int(region.breakStart2)-slop + 1 ) +"-" + str(int(region.breakStart2)+slop) + "_strand_" + region.strand1 
                        h2 = "_" + region.chrom1+":"+str(int(region.breakStart1)-slop + 1) +"-" + str(int(region.breakStart1)+slop) + "_strand_" + region.strand1 
                    else: 
                        h1 = "_" +region.chrom1 +":"+str(int(region.breakStart1)-slop+ 1 ) +"-" + str(int(region.breakStart1)+slop) + "_strand_" + region.strand1 
                        h2 = "_" + region.chrom1+":"+str(int(region.breakStart2)-slop+ 1) +"-" + str(int(region.breakStart2)+slop) + "_strand_" + region.strand1            
                    if "deletion" in region.variant_type:
                        if int(region.unaccounted_for_sequence) < 1:
                            return
                    if int(end_read) - slop < 0:
                        end_idx = 0 
                        end_slop = slop - (slop - int(end_read)) 
                    else:
                        end_idx = int(end_read) - slop
                        end_slop = slop
                    if int(start_read) - slop < 0:
                        start_idx = 0 
                        start_slop = slop - (slop - int(start_read))
                    else:
                        start_idx = int(start_read) - slop
                        start_slop = slop
                    seq1 = seq_fastq[(start_idx):(int(start_read)+slop)]
                    seq2 = seq_fastq[(end_idx):(int(end_read)+slop)]
                    #print(seq)
                    if region.strand1 == "-":
                        seq1 = seq_fastq[(end_idx):(int(end_read)+slop)]
                        seq2 = seq_fastq[(start_idx):(int(start_read)+slop)]
                        seq1 = complement(seq1)
                        seq2 = complement(seq2)
                #print(seq)
                #print(seq)bww
                elif end_read <= start_read:
                    if "duplication" in region.variant_type: 
                        h2 = "_" +region.chrom1 +":"+str(int(region.breakStart2)-slop+ 1 ) +"-" + str(int(region.breakStart2)+slop) + "_strand_" + region.strand1 
                        h1 = "_" + region.chrom1+":"+str(int(region.breakStart1)-slop+ 1) +"-" + str(int(region.breakStart1)+slop) + "_strand_" + region.strand1 
                    else: 
                        h1= "_" +region.chrom1 +":"+str(int(region.breakStart1)-slop + 1 ) +"-" + str(int(region.breakStart1)+slop) + "_strand_" + region.strand1 
                        h2 = "_" + region.chrom1+":"+str(int(region.breakStart2)-slop + 1) +"-" + str(int(region.breakStart2)+slop) + "_strand_" + region.strand1 
                    if int(end_read) - slop < 0:
                        end_idx = 0 
                        end_slop = slop - (slop - int(end_read)) 
                    else:
                        end_idx = int(end_read) - slop
                        end_slop = slop
                    if int(start_read) - slop < 0:
                        start_idx = 0 
                        start_slop = slop - (slop - int(start_read))
                    else:
                        start_idx = int(start_read) - slop
                        start_slop = slop
                    seq1 = seq_fastq[(start_idx):(int(start_read)+slop)]
                    seq2 = seq_fastq[(end_idx):(int(end_read)+slop)]
                    #print(seq)
                    if region.strand1 == "-":
                        seq2= seq_fastq[(start_idx):(int(start_read)+slop)]
                        seq1 = seq_fastq[(end_idx):(int(end_read)+slop)]
                        seq1 = complement(seq1)
                        seq2 = complement(seq2)
                   #print(seq)
                # deal with inversions next
                brk1_f.write(">BRK1_" + header + h1 +"\n"+ seq1 + "\n") 
                brk2_f.write(">BRK2_" +  header + h2 +"\n"+ seq2 + "\n") 
                brk1_f.write(">"+header + h1 +"\n"+ seq1[:start_slop] + "\n") 
                brk2_f.write(">"+ header + h2 +"\n"+ seq2[:end_slop] + "\n") 
                brk1_f.write(">"+header + h1 +"\n"+ seq1[start_slop:] + "\n") 
                brk2_f.write(">"+header + h2 +"\n"+ seq2[end_slop:]+ "\n")
                # TODO: move messy code.
                if "deletion" in region.variant_type: 
                    with open(os.path.join(dirname,"event_fastas",region.event_id + ".deletion.fasta"), "w") as deletion_fasta:
                        deletion_fasta.write(">"+header + h1 + "\n" + ref_fasta[int(region.breakStart1):int(region.breakEnd2)]+"\n") 
            elif region.strand1 != region.strand2:
                start_read = region.queryStart2
                end_read = region.queryEnd2
                print(region)
                print(start_read)
                print(end_read)
                if end_read < start_read:
                    h2 = "_" +region.chrom1 +":"+str(int(region.breakStart1)-slop+ 1) +"-" + str(int(region.breakStart1)+slop) + "_strand1_" + region.strand1 +"_strand2_"+ region.strand2
                    h1 = "_" +region.chrom1 +":"+str(int(region.breakStart2)-slop+ 1) +"-" + str(int(region.breakStart2)+slop) + "_strand1_" + region.strand1 +"_strand2_"+ region.strand2
                    if int(end_read) - slop < 0:
                        end_idx = 0 
                        end_slop = slop - (slop - int(end_read)) 
                    else:
                        end_idx = int(end_read) - slop
                        end_slop = slop
                    if int(start_read) - slop < 0:
                        start_idx = 0 
                        start_slop = slop - (slop - int(start_read))
                    else:
                        start_idx = int(start_read) - slop
                        start_slop = slop
                    seq1 = seq_fastq[(start_idx):(int(start_read)+slop)]
                    seq2 = seq_fastq[(end_idx):(int(end_read)+slop)]
                    
                    if region.strand2 == "-":
                        seq1 = seq_fastq[(end_idx):(int(end_read)+slop)]
                        seq2 = seq_fastq[(start_idx):(int(start_read)+slop)]
                        seq1 = complement(seq1)
                        seq2 = complement(seq2)
                elif end_read >= start_read:
                    print("Always here")
                    h1 = "_" +region.chrom1 +":"+str(int(region.breakStart1)-slop+ 1) +"-" + str(int(region.breakStart1)+slop) + "_strand1_" + region.strand1 +"_strand2_"+ region.strand2
                    h2 = "_" +region.chrom1 +":"+str(int(region.breakStart2)-slop+ 1) +"-" + str(int(region.breakStart2)+slop) + "_strand1_" + region.strand1 +"_strand2_"+ region.strand2
                    print(h1)
                    print(h2)
                    print("HERE")
                    if int(end_read) - slop < 0:
                        end_idx = 0 
                        end_slop = slop - (slop - int(end_read)) 
                    else:
                        end_idx = int(end_read) - slop
                        end_slop = slop
                    if int(start_read) - slop < 0:
                        start_idx = 0 
                        start_slop = slop - (slop - int(start_read))
                    else:
                        start_idx = int(start_read) - slop
                        start_slop = slop
                    seq1 = seq_fastq[(start_idx):(int(start_read)+slop)]
                    seq2 = seq_fastq[(end_idx):(int(end_read)+slop)]
                    if region.strand2 == "-":
                        seq2= seq_fastq[(start_idx):(int(start_read)+slop)]
                        seq1 = seq_fastq[(end_idx):(int(end_read)+slop)]
                        seq1 = complement(seq1)
                        seq2 = complement(seq2)
                brk1_f.write(">BRK1_" + header + h1 +"\n"+ seq1 + "\n") 
                brk2_f.write(">BRK2_" +  header + h2 +"\n"+ seq2 + "\n") 
                brk1_f.write(">"+header + h1 +"\n"+ seq1[:start_slop] + "\n") 
                brk2_f.write(">"+ header + h2 +"\n"+ seq2[:end_slop] + "\n") 
                brk1_f.write(">"+header + h1 +"\n"+ seq1[start_slop:] + "\n") 
                brk2_f.write(">"+header + h2 +"\n"+ seq2[end_slop:]+ "\n")

def complement(sequence):
    sequence_updated = []    
    for char in sequence:
        if char == "G":
            sequence_updated.append("C")
        elif char == "C":
            sequence_updated.append("G")
        elif char == "A":
            sequence_updated.append("T")
        elif char =="T":
            sequence_updated.append("A")
        else:
            sequence_updated.append(char)
    return "".join(sequence_updated[::-1])

def processFasta(region, output_prefix, fasta_input_file="/Users/smilefreak/Programming/OpenSource/MyGitHub/long_read_pipeline/scripts/test/RM_scaffolds.fasta"):
    """
        Open Fasta input file and reconnect the sequence to the assembly.
    
    """
    assembly = pysam.Fastafile(fasta_input_file)
    seq_fastq = assembly.fetch(region=region.ID)
    dirname = os.path.dirname(output_prefix)
    base_name_file = os.path.basename(fasta_input_file).split(".")[0]
    dirname = os.path.join(dirname, base_name_file)
    reference = pysam.Fastafile("/data/james/assembly_alignment/pipeline_run/input_data/BY.fasta")
    ref_fasta = reference.fetch(region=region.chrom1)
    try:
        os.makedirs(os.path.join(dirname,"event_fastas"))
    except OSError:
        pass
    try:
        os.makedirs(os.path.join(dirname,"event_fastas", "novel"))
    except OSError:
        pass
    try:
        os.makedirs(os.path.join(dirname,"event_fastas", "reference"))
    except OSError:
        pass
    with open(os.path.join(dirname,"event_fastas","novel",region.event_id + ".fasta"), "w") as out_f1:
        with open(os.path.join(dirname,"event_fastas","reference",region.event_id + ".fasta"), "w") as out_f2:
            header = ">" + region.event_id  +"\n"
            seq = ""
            if region.strand1 != region.strand2:
                start_read = region.queryEnd1  
                end_read = region.queryStart2
            elif region.strand1 == region.strand2: 
                start_read = region.queryEnd1  
                end_read = region.queryStart2
                    # TODO: not ideal.
                if end_read > start_read:
                    #
                    # must be [   DUP    ]------[    DUP      ]
                    #
                    # We must have duplication situation, so just extract the sequence
                    # Could be no insertion.
                    seq = seq_fastq[start_read:end_read]
                    #print(seq)
                    if region.strand1 == "-":
                        seq = complement(seq)
                    #print(seq)
                else:
                    #
                    # must be [  DUP  ] 
                    #               [  DUP  ] 
                    # We must have some lost sequence relative to the reference, so we should just extract that sequence
                    seq = seq_fastq[end_read:start_read]
                    #print(seq)
                    if region.strand1 == "-":
                        seq = complement(seq)
                    #print(seq)
               # print(region.strand1)
            #if "duplication" in region.variant_type:
            out_f2.write(header + ref_fasta[int(region.breakStart1):int(region.breakEnd2)]+"\n") 
            #else:
            #    out_f2.write(">"+header + h1 + "\n" + ref_fasta[int(region.breakStart1):int(region.breakEnd2)]+"\n") 
                #print(region.strand1)

            if int(region.unaccounted_for_sequence) > 0:
                out_f1.write(header + seq + "\n")
            else: 
                out_f1.close()
                os.remove(os.path.join(dirname,"event_fastas","novel",region.event_id + ".fasta"))

def splitterToBreakpoint(contigFile,splitReadSlop,localRange, reference_genome, output_prefix, next_region=""):
#Define whether event is inv, del, or dup (classifying translocations as if they were local intrachromosomal events)
    if contigFile == "stdin":
        contigs = sys.stdin
    else:
        contigs = open(str(contigFile))
    inversion_out = open(os.path.join(output_prefix + ".inv"), "w")
    cnv_out = open(os.path.join(output_prefix +".cnv"),"w")
    header = "chrom1\tbreakstart1\tbreakend1\tchrom2\tbreakstart2\tbreakend2\tID\tscore\tstrand1\tstrand2\tquerystart1\tqueryend1\tquerystart2\tqueryend2\tminnonoverlap\tquerylength\tqualscores\tvariant_type\tunaccount_for_sequence\tevent_length\tevent_id\n"
    inversion_out.write(header)
    cnv_out.write(header)
    inversion_dict = {}
    for i, line in enumerate(contigs):
        split = line.split()    
        if (split[0] > split[3]): 
            chrom1 = split[3]
            start1 = int(split[4])
            end1 = int(split[5])
            chrom2 = split[0]
            start2 = int(split[1])
            end2 = int(split[2])
            ID = split[6]
            score = int(split[7])
            strand1 = split[9]
            strand2 = split[8]
            queryStart1 = int(split[12])
            queryEnd1 = int(split[13])
            queryStart2 = int(split[10])
            queryEnd2 = int(split[11])
            minNonOverlap = int(split[14])
            queryLength = int(split[15])
            if len(split) > 16: qualScores = split[16]
            else: qualScores = '*'

        else:
            chrom1 = split[0]
            start1 = int(split[1])
            end1 = int(split[2])
            chrom2 = split[3]
            start2 = int(split[4])
            end2 = int(split[5])
            ID = split[6]
            score = int(split[7])
            strand1 = split[8]
            strand2 = split[9]
            queryStart1 = int(split[10])
            queryEnd1 = int(split[11])
            queryStart2 = int(split[12])
            queryEnd2 = int(split[13])         
       	    minNonOverlap = int(split[14])
       	    queryLength = int(split[15])
            if len(split) > 16: qualScores = split[16]
            else: qualScores = '*'
        
        if (queryStart1 == queryStart2):
            raise StandardError("spitterToBreakpoint cannot handle bedpe with equal starting query offsets in " + ID);
        classification_read_map = start2 - end1 
        classification_query_distances = queryEnd1 - queryStart1 

        
        #
        #  start ================================ end   
        #                                       |           start ======================= end
        #
        #          
        if strand1 != strand2: #Classifying Inversions
            if chrom1  == chrom2 and abs(end2 - start1): 
                variant = "local_inversion"
                variant_type = "inversion"
            else:
                variant = "distant_inversion"
                variant_type = "inversion"
            # TODO Lengths are off for this one!!!
            classification_read_map = start2 - end1 
            event_length = abs(classification_read_map)
            clean_break = queryStart2 - queryEnd1
            unaccounted_for_sequence = clean_break
        elif(strand1 == "+"):
            classification_read_map = start2 - end1 
            event_length = abs(classification_read_map)
            if classification_read_map > 0:
                variant_class = ("Deletion")
                clean_break = queryStart2 - queryEnd1
                if clean_break < -50 or clean_break > 50:
                    variant_type = "complex_deletion" 
                else:
                    variant_type = "simple_deletion"
                unaccounted_for_sequence = clean_break
                
            elif classification_read_map < 0:
                clean_break = (queryStart2 - queryEnd1)
                #clean_break = classification_read_map - clean_break
                if clean_break < -50 or clean_break > 50:
                    variant_type = "complex_duplication" 
                else:
                    variant_type = "simple_duplication"
                unaccounted_for_sequence = clean_break
        elif(strand1 == "-"):
            classification_read_map = start1 - end2 
            event_length = abs(classification_read_map) 
            if classification_read_map > 0:
                variant_class = ("Deletion")
                clean_break = queryStart2 - queryEnd1 
                if clean_break < -50 or clean_break > 50:
                    variant_type = "complex_deletion" 
                else:
                    variant_type = "simple_deletion"
            elif classification_read_map < 0:
                variant_class = ("Duplication")
                clean_break = queryStart2 - queryEnd1 
                if clean_break < -50 or clean_break > 50:
                    variant_type = "complex_duplication" 
                else:
                    variant_type = "simple_duplication"
            unaccounted_for_sequence = clean_break
        # OLD STUff
        if strand1 != strand2: #Classifying Inversions
            if chrom1  == chrom2 and abs(end2 - start1) <= localRange:
                variant = "local_inversion"
            else:
                variant = "distant_inversion"
        elif ((strand1=="+" and strand2=="+" and queryStart1 < queryStart2) or (strand1=="-" and strand2=="-" and queryStart1 > queryStart2)): #Classifying Deletions
            if chrom1  == chrom2 and (end2 - start1) <= localRange:
                variant = "local_deletion"
            else:
                variant = "distant_deletion"
        elif ((strand1=="+" and strand2=="+" and queryStart1 > queryStart2) or (strand1=="-" and strand2=="-" and queryStart1 < queryStart2)): #Classifying Tandem Duplications
            if chrom1  == chrom2 and (end2 - start1) <= localRange:
                variant = "local_duplication"
            else:
                variant = "distant_duplication"
        else:
            print "ERROR: variant " + ID + " is lost in variant classification step. Contact Mitchell"
        # Below modifed by GGF to account for zero-based half-open coordinate system.
        # When a start takes its coordinate from an end, subtract one.
        # When an end takes its coordinate from a start, add one.       
        if variant == "local_inversion" or variant == "distant_inversion":
            breakStart1 = str(start2-1-splitReadSlop)
            breakEnd1 = str(start2+splitReadSlop)
            breakStart2= str(end2-1-splitReadSlop)
            breakEnd2 = str(end2+splitReadSlop)
        elif variant_type == "complex_deletion" or variant_type == "simple_deletion":
            if(strand1 == "+"):
                breakStart1 = str(end1-splitReadSlop)
                breakEnd1 = str(end1+1+splitReadSlop)
                breakStart2 = str(start2-1-splitReadSlop)
                breakEnd2 = str(start2+splitReadSlop)
            elif(strand1 == "-"):
                breakStart1 = str(end2-splitReadSlop)
                breakEnd1 = str(end2+1+splitReadSlop)
                breakStart2 = str(start1-1-splitReadSlop)
                breakEnd2 = str(start1+splitReadSlop)
        elif variant_type == "complex_duplication" or variant_type == "simple_duplication": 
            if unaccounted_for_sequence < 0:
                duplication_in_reference = unaccounted_for_sequence + event_length
                if duplication_in_reference < 25 and duplication_in_reference > -25:
                    variant_type = "misaligned_sv"
            if(strand1 == "+"):
                breakStart1 = str(start2-splitReadSlop)
                breakEnd1 = str(start2+1+splitReadSlop)
                breakStart2 = str(end1-splitReadSlop)
                breakEnd2 = str(end1+1+splitReadSlop)
            else:
                breakStart1 = str(start1-splitReadSlop)
                breakEnd1 = str(start1+1+splitReadSlop)
                breakStart2 = str(end2-splitReadSlop)
                breakEnd2 = str(end2+1+splitReadSlop)
        else:
            print "ERROR: " + ID + " is lost in breakpoint locus prediction step. Contact Mitchell"
        if chrom1 == chrom2:
            cnv_output = CNVOutput(chrom1, breakStart1, breakEnd1, chrom2, breakStart2, breakEnd2, ID, score, strand1, strand2, queryStart1, queryEnd1, queryStart2, queryEnd2, minNonOverlap, queryLength, qualScores, variant, variant_type, unaccounted_for_sequence, event_length, i)
            if variant_type != "inversion":
                processFasta(cnv_output, fasta_input_file=reference_genome, output_prefix=output_prefix)
                _extract_breakpoints(cnv_output, fasta_input_file=reference_genome, output_prefix=output_prefix)
            # For now ignore inversions
            if variant_type != "inversion":
                cnv_out.write(str(cnv_output) +"\n")
            elif variant_type == "inversion":
                try:
                    inversion_dict[ID+"_"+chrom1].append(cnv_output)
                except KeyError:
                    inversion_dict[ID+"_"+chrom1] = [cnv_output]
    for key, value in inversion_dict.items():
        if len(value) == 2:
            # Look for ABA arrangement.
            if(value[0].strand2 == value[1].strand1): 
                processFasta(value[0], fasta_input_file=reference_genome, output_prefix=output_prefix)
                _extract_breakpoints(value[0], fasta_input_file=reference_genome, output_prefix=output_prefix) 
                cnv_output = CNVOutput(chrom1, breakStart1, breakEnd1, chrom2, breakStart2, breakEnd2, ID, score, strand1, strand2, queryStart1, queryEnd1, queryStart2, queryEnd2, minNonOverlap, queryLength, qualScores, variant, variant_type, unaccounted_for_sequence, event_length, i)
            inversion_out.write(str(value[0]) + "\n")
    inversion_out.close()
    cnv_out.close()

def main():
    usage = "%prog -b <bedpe> -c <contigFile.bedpe> [options]\nVersion: 0.2\nAuthor: Mitchell L. Leibowitz\nEdited: Aug, 2011 by Colby Chiang (details below)\n \
\nEdits included:\n 1) When switching ordering of bedpe positioning (if start2>start1), also switch orientation \n 2) Makes sure all comparisons are integer/integer comparisons or \
string/string comparisons. \n 3) Allowed for bedpe files of larger lengths (17 fields; see -s help).\n\n \
Note: Add slop at your own risk.\n\nPrints modified bedpe so bedpe entry 1 is always less than bedpe entry 2\n\nModified by GGF on 2/13/2012 to account for zero-base half-open coordinates."
    parser = OptionParser(usage)
    parser.add_option("-i", dest="contigFile", metavar="FILE", help="BEDPE file containing the split contigs; requires 10 column bedpe plus 11-14 as query coordinates, field 15 as query length \
and field 16 minimum non overlap.")
    parser.add_option("-s", dest="splitReadSlop", metavar="INT", type="int", default=0, help="Bidirectional slop around the breakpoint [default = 1]")
    parser.add_option("-r","--reference", dest="reference_genome", default="/Users/smilefreak/Programming/OpenSource/MyGitHub/long_read_pipeline/scripts/test/RM_scaffolds.fasta")
    parser.add_option("-o","--output-prefix", dest="output_prefix")
    (opts, args) = parser.parse_args()

    if opts.contigFile is None:
        parser.print_help()
        print
    else:
        splitterToBreakpoint(opts.contigFile, opts.splitReadSlop, 10000, opts.reference_genome, opts.output_prefix)
   
if __name__ == "__main__":
        sys.exit(main())
