#!/usr/bin/env python
#MIT License
#
#Copyright (c) 2019 Iñigo Prada Luengo
#
#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:
#
#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.


import argparse
import sys
from functools import partial
import os
import time
import pandas as pd
import pysam as ps
from extract_circle_SV_reads import readExtractor
from realigner import realignment
from bam2bam import bam2bam
from repeats import repeat
from utils import merge_final_output, filter_by_ratio, start_realign, start_simulate, mutate, insert_size_dist, write_clipped_read
from Coverage import coverage
import multiprocessing as mp
import pybedtools as bt
from simulations import sim_ecc_reads
import subprocess as sp
import glob
from tqdm import *
from circlemap.__version__ import __version__ as cm_version
import datetime

class circle_map:

    def __getpid__(self):

        pid = os.getpid()
        return (pid)

    def __init__(self):
        self.parser = argparse.ArgumentParser(
            description='Circle-Map',
            usage='''Circle-Map <subprogram> [options]

version=%s
contact= https://github.com/iprada/Circle-Map/issues

The Circle-Map suite

Commands:

    ReadExtractor   Extracts circular DNA read candidates
    Realign         Realign circular DNA read candidates
    bam2bam         Realign circular DNA read candidates and report them on a new BAM file
    Repeats         Identify circular DNA from repetitive regions
    Simulate        Simulate circular DNA

''' % cm_version)
        subparsers = self.parser.add_subparsers()

        self.readextractor = subparsers.add_parser(
            name="ReadExtractor",
            description='Extracts circular DNA read candidates',
            prog="Circle-Map ReadExtractor",
            usage='''Circle-Map ReadExtractor [options]'''

        )

        self.realigner = subparsers.add_parser(
            name="Realign",
            description='Realign circular DNA read candidates',
            prog="Circle-Map Realign",
            usage='''Circle-Map Realign [options]'''

        )

        self.repeats = subparsers.add_parser(
            name="Repeats",
            description='Identify circular DNA from repetitive regions',
            prog="Circle-Map Repeats",
            usage='''Circle-Map Repeats [options]'''

        )
        self.simulate = subparsers.add_parser(
            name="Simulate",
            description='Simulate eccDNA NGS datastes',
            prog="Circle-Map Reepeats",
            usage='''Circle-Map Simulate [options]'''

        )
        self.simulate = subparsers.add_parser(
            name="bam2bam",
            description='Realign the soft-clipped reads  and report',
            prog="Circle-Map bam2bam",
            usage='''Circle-Map bam2bam [options]'''

        )

        if len(sys.argv) <= 1:
            self.parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo argument given to Circle-Map"
                            "\nExiting\n")
            sys.exit(0)

        else:
            if sys.argv[1] == "ReadExtractor":


                self.subprogram = self.args_readextractor()
                self.args = self.subprogram.parse_args(sys.argv[2:])

                object = readExtractor(self.args.i, self.args.output, self.args.directory, self.args.quality,
                                        self.args.nodiscordant,
                                        self.args.nohardclipped, self.args.nosoftclipped, self.args.verbose,
                                        self.subprogram)
                object.extract_sv_circleReads()

            elif sys.argv[1] == "Realign":
                self.subprogram = self.args_realigner()
                self.args = self.subprogram.parse_args(sys.argv[2:])

                # get clusters
                splitted, sorted_bam, begin = start_realign(self.args.i, self.args.output, self.args.threads,
                                                            self.args.verbose, self.__getpid__(),
                                                            self.args.clustering_dist)




                sorted_bam.close()
                #get global insert size prior
                metrics = insert_size_dist(self.args.sample_size, self.args.insert_mapq, self.args.qbam)


                # pool based parallel of religment
                m = mp.Manager()

                lock = m.Lock()

                object = realignment(self.args.i, self.args.qbam, self.args.sbam, self.args.fasta,
                                        self.args.directory,
                                        self.args.mapq,
                                        self.args.insert_mapq, self.args.std, self.args.sample_size,
                                        self.args.gap_open,
                                        self.args.gap_ext, self.args.nhits, self.args.cut_off, self.args.min_sc,
                                        self.args.merge_fraction, self.args.interval_probability, self.args.output,
                                        self.args.threads, self.args.allele_frequency, lock, self.args.split,
                                        self.args.ratio, self.args.verbose, self.__getpid__(),
                                        self.args.edit_distance_fraction, self.args.remap_splits,
                                        self.args.only_discordants, self.args.split,
                                        self.args.split_quality, metrics,self.args.number_of_discordants)


                pool = mp.Pool(processes=self.args.threads)


                #progress bar
                with tqdm(total=len(splitted)) as pbar:
                    for i,exits in tqdm(enumerate(pool.imap_unordered(object.realign, splitted))):
                        pbar.update()
                        #kill if process returns 1,1
                        if exits == [1,1]:
                            pool.close()
                            pool.terminate()
                            pbar.close()
                            print(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S:"),
                                "An error happenend during execution. Exiting")
                            sys.exit()

                pbar.close()
                pool.close()
                pool.join()
                output = merge_final_output(self.args.sbam, self.args.output, begin, self.args.split,
                                            self.args.directory,
                                            self.args.merge_fraction, self.__getpid__())

                # compute coverage statistics
                if self.args.no_coverage == False:

                    coverage_object = coverage(self.args.sbam, output,
                                            self.args.bases, self.args.cmapq, self.args.extension,
                                            self.args.directory)

                    #Generator function for the coverage calculations
                    output = coverage_object.compute_coverage(coverage_object.get_wg_coverage())
                    filtered_output = filter_by_ratio(output, self.args.ratio)
                    filtered_output.to_csv(r'%s' % self.args.output, header=None, index=None, sep='\t', mode='w')

                else:
                    output.saveas("%s" % self.args.output)

            elif sys.argv[1] == "bam2bam":
                self.subprogram = self.args_bam2bam()
                self.args = self.subprogram.parse_args(sys.argv[2:])

                # get clusters
                splitted, sorted_bam, begin = start_realign(self.args.i, self.args.output, self.args.threads,
                                                            self.args.verbose, self.__getpid__(),
                                                            self.args.clustering_dist)


                #create output bam
                circle_sv_reads = ps.AlignmentFile(self.args.output, "wb", template=sorted_bam)

                sorted_bam.close()
                #get global insert size prior
                metrics = insert_size_dist(self.args.sample_size, self.args.insert_mapq, self.args.qbam)




                manager = mp.Manager()





                lock = manager.Lock()

                object = bam2bam(self.args.i,self.args.output,self.args.qbam, self.args.fasta,
                                self.args.directory,
                                self.args.mapq,
                                self.args.insert_mapq, self.args.std, self.args.sample_size,
                                self.args.gap_open,
                                self.args.gap_ext, self.args.nhits, self.args.cut_off, self.args.min_sc,
                                self.args.interval_probability,
                                self.args.threads, lock,
                                self.args.verbose, self.__getpid__(),
                                self.args.edit_distance_fraction, self.args.remap_splits,
                                self.args.only_discordants,
                                self.args.split_quality, metrics,manager)

                object.beta_version_warning()


                pool = mp.Pool(processes=self.args.threads)
                # create writer process

                writer_p =  mp.Process(target=object.listener_writer, args=(circle_sv_reads,))
                writer_p.daemon = True
                writer_p.start()
                #progress bar
                with tqdm(total=len(splitted)) as pbar:
                    for i,exits in tqdm(enumerate(pool.imap_unordered(object.realign,splitted))):
                        pbar.update()
                        # kill if process returns 1,1
                        if exits == [1, 1]:
                            pool.close()
                            pool.terminate()
                            pbar.close()
                            print(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S:"),
                                "An error happenend during execution. Exiting")
                            sys.exit()

                pbar.close()
                pool.close()
                pool.join()
                object.kill()
                writer_p.terminate()
                writer_p.join()

                circle_sv_reads.close()
                print("Done")


            elif sys.argv[1] == "Repeats":

                self.subprogram = self.args_repeats()
                self.args = self.subprogram.parse_args(sys.argv[2:])


                object = repeat(self.args.i, self.args.directory, self.args.mismatch, self.args.fraction,
                                self.args.read_number)

                bed = object.find_circles()

                coverage_object = coverage(self.args.i, bed,
                                            self.args.bases, self.args.cmapq, self.args.extension,
                                            self.args.directory)

                output = coverage_object.compute_coverage(coverage_object.get_wg_coverage())

                filtered_output = filter_by_ratio(output, self.args.ratio)
                filtered_output.to_csv(r'%s' % self.args.output, header=None, index=None, sep='\t', mode='w')



            elif sys.argv[1] == "Simulate":

                self.subprogram = self.args_simulate()
                self.args = self.subprogram.parse_args(sys.argv[2:])

                sim_pid = start_simulate(self.__getpid__())

                lock = mp.Lock()

                paired_end_fastq_1 = open("%s_1.fastq" % self.args.base_name, "w")
                paired_end_fastq_2 = open("%s_2.fastq" % self.args.base_name, "w")
                paired_end_fastq_1.close()
                paired_end_fastq_2.close()

                # mutate reference genome
                if self.args.variants == True:
                    mutate(self.args.g, sim_pid, self.args.Indels, self.args.substitution, self.args.java_memory)


                manager = mp.Manager()
                # Shared memory object
                circle_list = manager.list()
                skipped_circles = mp.Value('i', 0)
                correct_circles = mp.Value('i', 0)
                jobs = []
                # init the processes

                for i in range(self.args.processes):
                    p = mp.Process(target=sim_ecc_reads,
                                args=(self.args.g, self.args.read_length, self.args.directory,
                                        int(round(self.args.read_number / self.args.processes)),
                                        self.args.skip_region, self.args.base_name,
                                        self.args.mean_insert_size, self.args.error,
                                        self.args.mean_coverage, lock, i, circle_list,
                                        "%s_1.fastq" % self.args.base_name, "%s_2.fastq" % self.args.base_name,
                                        skipped_circles,
                                        correct_circles, self.args.insRate, self.args.insRate2, self.args.delRate,
                                        self.args.delRate2, sim_pid,))
                    jobs.append(p)
                    p.start()
                # kill the process
                for p in jobs:
                    p.join()
                print("Skipped %s circles, that overlapped the provided regions to exclude" % skipped_circles.value)
                print("Simulated %s circles across %s parallel processes" % (
                correct_circles.value, self.args.processes))
                print("Writting to disk bed file containing the simulated circle coordinates")

                bt.BedTool(list(circle_list)).saveas(self.args.output)

            else:
                self.parser.print_help()
                time.sleep(0.01)
                sys.stderr.write("\nWrong argument given to Circle-Map"
                                "\nExiting\n")
                sys.exit(0)

    def args_readextractor(self):

        parser = self.readextractor

        parser._action_groups.pop()
        required = parser.add_argument_group('required arguments')
        optional = parser.add_argument_group('optional arguments')
        # prefixing the argument with -- means it's optional
        # input and output

        required.add_argument('-i', metavar='', help="Input: query name sorted bam file")

        if "-i" in sys.argv:
            optional.add_argument('-o', '--output', metavar='',
                                    help="Ouput: Reads indicating circular DNA structural variants",
                                    default="circle_%s" % sys.argv[sys.argv.index("-i") + 1])

            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            # mapping quality cutoff

            optional.add_argument('-q', '--quality', type=int, metavar='',
                                    help="bwa-mem mapping quality cutoff. Default value 10",
                                    default=10)

            # read extraction options
            # extract discordant reads
            optional.add_argument('-nd', '--nodiscordant', help="Turn off discordant (R2F1 oriented) read extraction",
                                    action='store_true')

            # soft-clipped argument
            optional.add_argument('-nsc', '--nosoftclipped', help="Turn off soft-clipped read extraction",
                                    action='store_true')
            # extract hard-clippped reads
            optional.add_argument('-nhc', '--nohardclipped', help="Turn off hard-clipped read extraction",
                                    action='store_true')

            # verbose level

            optional.add_argument('-v', '--verbose', type=int, metavar='',
                                    help='Verbose level, 1=error,2=warning, 3=message',
                                    choices=[1, 2, 3], default=3)

        else:
            optional.add_argument('-o', '--output', metavar='',
                                    help="Ouput: Reads indicating circular DNA structural variants")

            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            # mapping quality cutoff
            optional.add_argument('-q', '--quality', type=int, metavar='',
                                    help="bwa-mem mapping quality cutoff. Default value 10",
                                    default=10)

            # read extraction options
            # extract discordant reads
            optional.add_argument('-nd', '--nodiscordant', help="Turn off discordant (R2F1 oriented) read extraction",
                                    action='store_true')

            # soft-clipped argument
            optional.add_argument('-nsc', '--nosoftclipped', help="Turn off soft-clipped read extraction",
                                    action='store_true')
            # extract hard-clippped reads
            optional.add_argument('-nhc', '--nohardclipped', help="Turn off hard-clipped read extraction",
                                    action='store_true')

            # verbose level

            optional.add_argument('-v', '--verbose', type=int, metavar='',
                                    help='Verbose level, 1=error,2=warning, 3=message. Default=3',
                                    choices=[1, 2, 3], default=3)

            parser.print_help()

            time.sleep(0.01)
            sys.stderr.write(
                "\nNo input or output input given to readExtractor, be sure that you are providing the flags'-i' and '-o'"
                "\nExiting\n")
            sys.exit(0)

        # parse the commands

        if len(sys.argv[2:]) == 0:
            parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo arguments given to read extractor. Exiting\n")
            sys.exit(0)

        return (parser)

    def args_realigner(self):
        parser = self.realigner

        # declare the different groups for the parser
        parser._action_groups.pop()
        io_options = parser.add_argument_group('Input/Output options')
        alignment_options = parser.add_argument_group('Alignment options')
        out_decision = parser.add_argument_group('eccDNA output options')
        i_size_estimate = parser.add_argument_group('Insert size estimation options')
        interval = parser.add_argument_group('Interval processing options')
        coverage_metrics = parser.add_argument_group('Coverage metrics options')
        running = parser.add_argument_group('Running options')

        io_options.add_argument('-i', metavar='',
                                help="Input: bam file containing the reads extracted by ReadExtractor")
        io_options.add_argument('-qbam', metavar='', help="Input: query name sorted bam file")
        io_options.add_argument('-sbam', metavar='', help="Input: coordinate sorted bam file")
        io_options.add_argument('-fasta', metavar='', help="Input: Reference genome fasta file")

        if "-i" and "-qbam" and "-fasta" in sys.argv:
            # output

            io_options.add_argument('-o', '--output', metavar='', help="Output filename",
                                    default="circle_%s.bed" % sys.argv[sys.argv.index("-i") + 1])

            # alignment
            alignment_options.add_argument('-n', '--nhits', type=int, metavar='',
                                            help="Number of realignment attempts. Default: 10",
                                            default=10)

            alignment_options.add_argument('-p', '--cut_off', type=float, metavar='',
                                            help="Probability cut-off for considering a soft-clipped as realigned: Default: 0.99",
                                            default=0.99)

            alignment_options.add_argument('-m', '--min_sc', type=float, metavar='',
                                            help="Minimum soft-clipped length to attempt the realignment. Default: 8",
                                            default=8)

            alignment_options.add_argument('-g', '--gap_open', type=int, metavar='',
                                            help="Gap open penalty in the position specific scoring matrix. Default: 5",
                                            default=5)

            alignment_options.add_argument('-e', '--gap_ext', type=int, metavar='',
                                            help="Gap extension penalty in the position specific scoring matrix. Default: 1",
                                            default=1)

            alignment_options.add_argument('-q', '--mapq', type=int, metavar='',
                                            help="Minimum mapping quality allowed in the supplementary alignments. Default: 20",
                                            default=20)

            alignment_options.add_argument('-d', '--edit_distance-fraction', type=float, metavar='',
                                            help="Maximum edit distance fraction allowed in the first realignment. Default (0.05)",
                                            default=0.05)

            alignment_options.add_argument('-Q', '--split_quality', type=float, metavar='',
                                            help="Minium split score to output an interval. Default (0.0)",
                                            default=0.0)

            alignment_options.add_argument('-R', '--remap_splits', help="Remap probabilistacally the split reads",
                                            action='store_true')

            # insert size

            i_size_estimate.add_argument('-iq', '--insert_mapq', type=int, metavar='',
                                            help="Mapq cutoff for stimating the insert size distribution. Default 60",
                                            default=60)

            i_size_estimate.add_argument('-sd', '--std', type=int, metavar='',
                                            help="Standard deviations of the insert size to extend the intervals. Default 5",
                                            default=4)

            i_size_estimate.add_argument('-s', '--sample_size', type=int, metavar='',
                                            help="Number of concordant reads (R2F1) to use for estimating the insert size distribution. Default 100000",
                                            default=100000)

            # Interval options

            interval.add_argument('-f', '--merge_fraction', type=float, metavar='',
                                    help="Merge intervals reciprocally overlapping by a fraction. Default 0.99",
                                    default=0.99)

            interval.add_argument('-P', '--interval_probability', type=float, metavar='',
                                    help="Skip edges of the graph with a probability below the threshold. Default: 0.01",
                                    default=0.01)
            interval.add_argument('-K', '--clustering_dist', type=int, metavar='',
                                    help="Cluster reads that are K nucleotides appart in the same node. Default: 500",
                                    default=500)

            interval.add_argument('-D', '--only_discordants', help="Use only discordant reads to build the graph",
                                    action='store_false')
            interval.add_argument('-F', '--allele_frequency', type=float, metavar='',
                                    help="Minimum allele frequency required to report the circle interval. Default (0.1)",
                                    default=0.1)
            # When to call a circle

            out_decision.add_argument('-S', '--split', type=int, metavar='',
                                    help="Number of required split reads to output a eccDNA. Default: 0",
                                    default=0)

            out_decision.add_argument('-O', '--number_of_discordants', type=int, metavar='',
                                    help="Number of required discordant reads for intervals with only discordants. Default: 3",
                                    default=3)
            out_decision.add_argument('-r', '--ratio', type=float, metavar='',
                                    help="Minimum in/out required coverage ratio. Default: 0.0",
                                    default=0.0)

            # coverage metrics

            coverage_metrics.add_argument('-N', '--no_coverage', help="Don't compute coverage statistics",
                                        action='store_true')

            coverage_metrics.add_argument('-b', '--bases', type=int, metavar='',
                                        help="Number of bases to extend for computing the coverage ratio. Default: 200",
                                        default=200)

            coverage_metrics.add_argument('-cq', '--cmapq', type=int, metavar='',
                                        help="Minimum mapping quality treshold for coverage computation. Default: 0",
                                        default=0)

            coverage_metrics.add_argument('-E', '--extension', type=int, metavar='',
                                        help="Number of bases inside the eccDNA breakpoint coordinates to compute the ratio. Default: 100",
                                        default=100)

            # run options

            running.add_argument('-t', '--threads', type=int, metavar='',
                                help="Number of threads to use.Default 1",
                                default=1)

            running.add_argument('-dir', '--directory', metavar='',
                                help="Working directory, default is the working directory",
                                default=os.getcwd())

            running.add_argument('-v', '--verbose', type=int, metavar='',
                                help='Verbose level, 1=error,2=warning, 3=message',
                                choices=[1, 2, 3], default=3)



        else:
            # output
            io_options.add_argument('-o', metavar='', help="Output filename")

            alignment_options.add_argument('-n', '--nhits', type=int, metavar='',
                                            help="Number of realignment attempts. Default: 10",
                                            default=10)

            alignment_options.add_argument('-p', '--cut_off', type=float, metavar='',
                                            help="Probability cut-off for considering a soft-clipped as realigned: Default: 0.99",
                                            default=0.99)

            alignment_options.add_argument('-m', '--min_sc', type=float, metavar='',
                                            help="Minimum soft-clipped length to attempt the realignment. Default: 8",
                                            default=8)

            alignment_options.add_argument('-g', '--gap_open', type=int, metavar='',
                                            help="Gap open penalty in the position specific scoring matrix. Default: 5",
                                            default=5)

            alignment_options.add_argument('-e', '--gap_ext', type=int, metavar='',
                                            help="Gap extension penalty in the position specific scoring matrix. Default: 1",
                                            default=1)

            alignment_options.add_argument('-q', '--mapq', type=int, metavar='',
                                            help="Minimum mapping quality allowed in the supplementary alignments. Default: 20",
                                            default=20)

            alignment_options.add_argument('-d', '--edit_distance-fraction', type=float, metavar='',
                                            help="Maximum edit distance fraction allowed in the first realignment. Default (0.05)",
                                            default=0.05)

            alignment_options.add_argument('-Q', '--split_quality', type=float, metavar='',
                                            help="Minium split score to output an interval. Default (0.0)",
                                            default=0.0)
            alignment_options.add_argument('-R', '--remap_splits', help="Remap probabilistacally bwa-mem split reads",
                                            action='store_true')

            # insert size

            i_size_estimate.add_argument('-iq', '--insert_mapq', type=int, metavar='',
                                            help="Mapq cutoff for stimating the insert size distribution. Default 60",
                                            default=60)

            i_size_estimate.add_argument('-sd', '--std', type=int, metavar='',
                                            help="Standard deviations of the insert size to extend the intervals. Default 5",
                                            default=5)

            i_size_estimate.add_argument('-s', '--sample_size', type=int, metavar='',
                                            help="Number of concordant reads (R2F1) to use for estimating the insert size distribution. Default 100000",
                                            default=100000)

            # Interval options

            interval.add_argument('-f', '--merge_fraction', type=float, metavar='',
                                    help="Merge intervals reciprocally overlapping by a fraction. Default 0.99",
                                    default=0.99)

            interval.add_argument('-P', '--interval_probability', type=float, metavar='',
                                    help="Skip edges of the graph with a probability below the threshold. Default: 0.01",
                                    default=0.01)
            interval.add_argument('-K', '--clustering_dist', type=int, metavar='',
                                    help="Cluster reads that are K nucleotides appart in the same node. Default: 500",
                                    default=500)
            interval.add_argument('-D', '--only_discordants', help="Use only discordant reads to build the graph",
                                    action='store_true')
            interval.add_argument('-F', '--allele_frequency', type=float, metavar='',
                                    help="Minimum allele frequency required to report the circle interval. Default (0.1)",
                                    default=0.1)

            # When to call a circle

            out_decision.add_argument('-S', '--split', type=int, metavar='',
                                        help="Number of required split reads to output a eccDNA. Default: 0",
                                        default=0)
            out_decision.add_argument('-O', '--number_of_discordants', type=int, metavar='',
                                        help="Number of required discordant reads for intervals with only discordants. Default: 3",
                                        default=3)

            out_decision.add_argument('-r', '--ratio', type=float, metavar='',
                                        help="Minimum in/out required coverage ratio. Default: 0.0",
                                        default=0.0)

            # coverage metrics

            coverage_metrics.add_argument('-N', '--no_coverage', help="Don't compute coverage statistics",
                                            action='store_true')

            coverage_metrics.add_argument('-b', '--bases', type=int, metavar='',
                                            help="Number of bases to extend for computing the coverage ratio. Default: 200",
                                            default=200)

            coverage_metrics.add_argument('-cq', '--cmapq', type=int, metavar='',
                                            help="Minimum mapping quality treshold for coverage computation. Default: 0",
                                            default=0)

            coverage_metrics.add_argument('-E', '--extension', type=int, metavar='',
                                            help="Number of bases inside the eccDNA breakpoint coordinates to compute the ratio. Default: 100",
                                            default=100)

            # Running options

            running.add_argument('-t', '--threads', type=int, metavar='',
                                help="Number of threads to use.Default 1",
                                default=1)

            running.add_argument('-dir', '--directory', metavar='',
                                help="Working directory, default is the working directory",
                                default=os.getcwd())

            running.add_argument('-v', '--verbose', type=int, metavar='',
                                help='Verbose level, 1=error,2=warning, 3=message',
                                choices=[1, 2, 3], default=3)

            # find out which arguments are missing

            parser.print_help()

            time.sleep(0.01)
            sys.stderr.write("\nInput does not match. Check that you provide the -i, -qbam and -fasta options"
                            "\nExiting\n")
            sys.exit(0)

        if len(sys.argv[2:]) == 0:
            parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo arguments given to Realign. Exiting\n")
            sys.exit(0)

        return (parser)


    def args_bam2bam(self):
        parser = self.realigner

        # declare the different groups for the parser
        parser._action_groups.pop()
        io_options = parser.add_argument_group('Required')
        alignment_options = parser.add_argument_group('Alignment options')
        i_size_estimate = parser.add_argument_group('Insert size estimation options')
        interval = parser.add_argument_group('Interval processing options')
        running = parser.add_argument_group('Running options')

        io_options.add_argument('-i', metavar='',
                                help="Input: bam file containing the reads extracted by ReadExtractor")
        io_options.add_argument('-qbam', metavar='', help="Input: query name sorted bam file")
        io_options.add_argument('-fasta', metavar='', help="Input: Reference genome fasta file")
        io_options.add_argument('-o', '--output', metavar='', help="Output BAM name")

        if "-i" and "-qbam" and "-fasta" and "-o" in sys.argv:
            # output



            # alignment
            alignment_options.add_argument('-n', '--nhits', type=int, metavar='',
                                            help="Number of realignment attempts. Default: 10",
                                            default=10)

            alignment_options.add_argument('-p', '--cut_off', type=float, metavar='',
                                            help="Probability cut-off for considering a soft-clipped as realigned: Default: 0.99",
                                            default=0.99)

            alignment_options.add_argument('-m', '--min_sc', type=float, metavar='',
                                            help="Minimum soft-clipped length to attempt the realignment. Default: 8",
                                            default=8)

            alignment_options.add_argument('-g', '--gap_open', type=int, metavar='',
                                            help="Gap open penalty in the position specific scoring matrix. Default: 5",
                                            default=5)

            alignment_options.add_argument('-e', '--gap_ext', type=int, metavar='',
                                            help="Gap extension penalty in the position specific scoring matrix. Default: 1",
                                            default=1)

            alignment_options.add_argument('-q', '--mapq', type=int, metavar='',
                                            help="Minimum mapping quality allowed in the supplementary alignments. Default: 20",
                                            default=20)

            alignment_options.add_argument('-d', '--edit_distance-fraction', type=float, metavar='',
                                            help="Maximum edit distance fraction allowed in the first realignment. Default (0.05)",
                                            default=0.05)

            alignment_options.add_argument('-Q', '--split_quality', type=float, metavar='',
                                            help="Minium split score to output an interval. Default (0.0)",
                                            default=0.0)

            alignment_options.add_argument('-R', '--remap_splits', help="Remap probabilistacally the split reads",
                                            action='store_true')

            # insert size

            i_size_estimate.add_argument('-iq', '--insert_mapq', type=int, metavar='',
                                            help="Mapq cutoff for stimating the insert size distribution. Default 60",
                                            default=60)

            i_size_estimate.add_argument('-sd', '--std', type=int, metavar='',
                                            help="Standard deviations of the insert size to extend the intervals. Default 5",
                                            default=4)

            i_size_estimate.add_argument('-s', '--sample_size', type=int, metavar='',
                                            help="Number of concordant reads (R2F1) to use for estimating the insert size distribution. Default 100000",
                                            default=100000)

            # Interval options



            interval.add_argument('-P', '--interval_probability', type=float, metavar='',
                                    help="Skip edges of the graph with a probability below the threshold. Default: 0.01",
                                    default=0.01)
            interval.add_argument('-K', '--clustering_dist', type=int, metavar='',
                                    help="Cluster reads that are K nucleotides appart in the same node. Default: 500",
                                    default=500)

            interval.add_argument('-D', '--only_discordants', help="Use only discordant reads to build the graph",
                                    action='store_false')



            # run options

            running.add_argument('-t', '--threads', type=int, metavar='',
                                    help="Number of threads to use.Default 1",
                                    default=1)

            running.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            running.add_argument('-v', '--verbose', type=int, metavar='',
                                    help='Verbose level, 1=error,2=warning, 3=message',
                                    choices=[1, 2, 3], default=3)



        else:
            # output
            alignment_options.add_argument('-n', '--nhits', type=int, metavar='',
                                            help="Number of realignment attempts. Default: 10",
                                            default=10)

            alignment_options.add_argument('-p', '--cut_off', type=float, metavar='',
                                            help="Probability cut-off for considering a soft-clipped as realigned: Default: 0.99",
                                            default=0.99)

            alignment_options.add_argument('-m', '--min_sc', type=float, metavar='',
                                            help="Minimum soft-clipped length to attempt the realignment. Default: 8",
                                            default=8)

            alignment_options.add_argument('-g', '--gap_open', type=int, metavar='',
                                            help="Gap open penalty in the position specific scoring matrix. Default: 5",
                                            default=5)

            alignment_options.add_argument('-e', '--gap_ext', type=int, metavar='',
                                            help="Gap extension penalty in the position specific scoring matrix. Default: 1",
                                            default=1)

            alignment_options.add_argument('-q', '--mapq', type=int, metavar='',
                                            help="Minimum mapping quality allowed in the supplementary alignments. Default: 20",
                                            default=20)

            alignment_options.add_argument('-d', '--edit_distance-fraction', type=float, metavar='',
                                            help="Maximum edit distance fraction allowed in the first realignment. Default (0.05)",
                                            default=0.05)

            alignment_options.add_argument('-Q', '--split_quality', type=float, metavar='',
                                            help="Minium split score to output an interval. Default (0.0)",
                                            default=0.0)
            alignment_options.add_argument('-R', '--remap_splits', help="Remap probabilistacally bwa-mem split reads",
                                            action='store_true')

            # insert size

            i_size_estimate.add_argument('-iq', '--insert_mapq', type=int, metavar='',
                                            help="Mapq cutoff for stimating the insert size distribution. Default 60",
                                            default=60)

            i_size_estimate.add_argument('-sd', '--std', type=int, metavar='',
                                            help="Standard deviations of the insert size to extend the intervals. Default 5",
                                            default=5)

            i_size_estimate.add_argument('-s', '--sample_size', type=int, metavar='',
                                            help="Number of concordant reads (R2F1) to use for estimating the insert size distribution. Default 100000",
                                            default=100000)

            # Interval options


            interval.add_argument('-P', '--interval_probability', type=float, metavar='',
                                    help="Skip edges of the graph with a probability below the threshold. Default: 0.01",
                                    default=0.01)
            interval.add_argument('-K', '--clustering_dist', type=int, metavar='',
                                    help="Cluster reads that are K nucleotides appart in the same node. Default: 500",
                                    default=500)
            interval.add_argument('-D', '--only_discordants', help="Use only discordant reads to build the graph",
                                    action='store_true')


            # Running options

            running.add_argument('-t', '--threads', type=int, metavar='',
                                    help="Number of threads to use.Default 1",
                                    default=1)

            running.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            running.add_argument('-v', '--verbose', type=int, metavar='',
                                    help='Verbose level, 1=error,2=warning, 3=message',
                                    choices=[1, 2, 3], default=3)

            # find out which arguments are missing

            parser.print_help()

            time.sleep(0.01)
            sys.stderr.write("\nInput does not match. Check that you provide the -i, -qbam and -fasta options"
                                "\nExiting\n")
            sys.exit(0)

        if len(sys.argv[2:]) == 0:
            parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo arguments given to bam2bam. Exiting\n")
            sys.exit(0)

        return (parser)

    def args_repeats(self):

        parser = self.repeats

        parser._action_groups.pop()
        required = parser.add_argument_group('required arguments')
        optional = parser.add_argument_group('optional arguments')
        # prefixing the argument with -- means it's optional
        # input and output

        required.add_argument('-i', metavar='', help="Input: coordinate name sorted bam file")

        if "-i" in sys.argv:

            optional.add_argument('-o', '--output', metavar='',
                                    help="Ouput: Reads indicating circular DNA structural variants from repeat regions",
                                    default="circle_repeats_%s" % sys.argv[sys.argv.index("-i") + 1])

            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            # coverage metrics
            optional.add_argument('-m', '--mismatch', metavar='',
                                    help="Number of mismatches allowed on the reads",
                                    default=2)

            optional.add_argument('-b', '--bases', type=int, metavar='',
                                    help="Number of bases to extend for computing the coverage ratio. Default: 200",
                                    default=200)

            optional.add_argument('-cq', '--cmapq', type=int, metavar='',
                                    help="Minimum mapping quality treshold for coverage computation. Default: 0",
                                    default=0)

            optional.add_argument('-E', '--extension', type=int, metavar='',
                                    help="Number of bases inside the eccDNA coordinates to compute the ratio. Default: 100",
                                    default=100)

            optional.add_argument('-r', '--ratio', type=float, metavar='',
                                    help="Minimum in/out required ratio. Default: 0.6",
                                    default=0.6)

            optional.add_argument('-f', '--fraction', type=float, metavar='',
                                    help="Required fraction to merge the intervals of the double mapped reads. Default 0.8",
                                    default=0.8)

            optional.add_argument('-n', '--read_number', metavar='',
                                    help="Minimum number of reads required to output",
                                    default=20)



        else:

            optional.add_argument('-o', '--output', metavar='',
                                    help="Ouput: Reads indicating circular DNA structural variants",
                                    )

            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())

            # coverage metrics

            optional.add_argument('-m', '--mismatch', metavar='',
                                    help="Number of mismatches allowed on the reads",
                                    default=2)

            optional.add_argument('-b', '--bases', type=int, metavar='',
                                    help="Number of bases to extend for computing the coverage ratio. Default: 200",
                                    default=200)

            optional.add_argument('-cq', '--cmapq', type=int, metavar='',
                                    help="Minimum mapping quality treshold for coverage computation. Default: 0",
                                    default=0.6)

            optional.add_argument('-E', '--extension', type=int, metavar='',
                                    help="Number of bases inside the eccDNA coordinates to compute the ratio. Default: 100",
                                    default=100)

            optional.add_argument('-r', '--ratio', type=float, metavar='',
                                    help="Minimum in/out required ratio. Default: 0.6",
                                    default=0.6)

            optional.add_argument('-f', '--fraction', type=float, metavar='',
                                    help="Required fraction to merge the intervals of the double mapped reads. Default 0.8",
                                    default=0.8)

            optional.add_argument('-n', '--read_number', metavar='',
                                    help="Minimum number of reads required to output",
                                    default=20)

            parser.print_help()

            time.sleep(0.01)
            sys.stderr.write("\nNo input input given to Repeats, be sure that you are providing the flag '-i'"
                                "\nExiting\n")
            sys.exit(0)

        # parse the commands

        if len(sys.argv[2:]) == 0:
            parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo arguments given to Repeats. Exiting\n")
            sys.exit(0)

        return (parser)

    def args_simulate(self):

        parser = self.simulate

        parser._action_groups.pop()
        required = parser.add_argument_group('required arguments')
        optional = parser.add_argument_group('optional arguments')
        # prefixing the argument with -- means it's optional
        # input and output

        if "-g" and "-N" in sys.argv:
            required.add_argument('-g', metavar='',
                                    help="Genome fasta file (Needs to be indexed with samtools faidx)")
            required.add_argument('-N', '--read-number', type=int, metavar='',
                                    help="Number of reads to simulate")
            optional.add_argument('-o', '--output', default='simulated.bed',
                                    help="Output file name")
            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())
            optional.add_argument('-b', '--base-name', metavar='', default='simulated',
                                    help="Fastq output basename")
            optional.add_argument('-s', '--skip-region', metavar='', default=None,
                                    help="Regions of the genome to skip the simulation. The input needs to be in bed format")
            optional.add_argument('-r', '--read-length', metavar='', type=int, default=150,
                                    help="Read length to simulate")
            optional.add_argument('-m', '--mean-insert-size', metavar='', type=int, default=300,
                                    help="Mean of the insert size distribution")
            optional.add_argument('-c', '--mean-coverage', metavar='', type=int, default=30,
                                    help="Mean sequencing coverage within the eccDNA coordinates")
            optional.add_argument('-p', '--processes', metavar='', type=int, default=1,
                                    help="Mean sequencing coverage within the eccDNA coordinates")

            optional.add_argument('-v', '--variants', action='store_true',
                                    help="If set to true, introduce mutations in the reference genome prior to simulating"
                                        "reads.")
            optional.add_argument('-S', '--substitution', metavar='', type=float, default=0.0001,
                                    help="Fraction of base substitutions to introduce on the genome. Default: 0.0001")

            optional.add_argument('-I', '--Indels', metavar='', type=float, default=0.001,
                                    help="Fraction of indels to introduce on the genome. Default: 0.001")
            optional.add_argument('-J', '--java_memory', metavar='', type=str, default="-Xmx16g",
                                    help="Java memory allocation, required for mutating the genome. Default: -Xmx16g")

            optional.add_argument('-e', '--error', action='store_true',
                                    help="Introduce sequencing errors ( Uses ART on the background)")

            optional.add_argument('-i', '--instrument', metavar='', type=str, default="HS25",
                                    help="Illumina sequecing instrument to simulate reads from (Default HiSeq 2500)")

            optional.add_argument('-ir', '--insRate', metavar='', type=float, default=0.00009,
                                    help="the first-read insertion rate (default: 0.00009)")
            optional.add_argument('-ir2', '--insRate2', metavar='', type=float, default=0.00015,
                                    help="the second-read insertion rate (default: 0.00015)")
            optional.add_argument('-dr', '--delRate', metavar='', type=float, default=0.00011,
                                    help="the first-read deletion rate (default:  0.00011)")
            optional.add_argument('-dr2', '--delRate2', metavar='', type=float, default=0.00023,
                                    help="the second-read deletion rate (default: 0.00023)")
        else:
            required.add_argument('-g', metavar='',
                                    help="Genome fasta file (Needs to be indexed with samtools faidx)")
            required.add_argument('-N', '--read-number', type=int, metavar='',
                                    help="Number of reads to simulate")
            optional.add_argument('-o', '--output', default='simulated.bed',
                                    help="Output file name")
            optional.add_argument('-dir', '--directory', metavar='',
                                    help="Working directory, default is the working directory",
                                    default=os.getcwd())
            optional.add_argument('-b', '--base-name', metavar='', default='simulated',
                                    help="Fastq output basename")
            optional.add_argument('-s', '--skip-region', metavar='', default=None,
                                    help="Regions of the genome to skip the simulation. The input needs to be in bed format")
            optional.add_argument('-r', '--read-length', metavar='', type=int, default=150,
                                    help="Read length to simulate")
            optional.add_argument('-m', '--mean-insert', metavar='', type=int, default=300,
                                    help="Mean of the insert size distribution")

            optional.add_argument('-c', '--mean-coverage', metavar='', type=int, default=30,
                                    help="Mean sequencing coverage within the eccDNA coordinates")

            optional.add_argument('-p', '--processes', metavar='', type=int, default=1,
                                    help="Number of parallel processes to use")

            optional.add_argument('-v', '--variants', action='store_true',
                                    help="If set to true, introduce mutations in the reference genome prior to simulating"
                                        "reads.")
            optional.add_argument('-S', '--substitution', metavar='', type=float, default=0.0001,
                                    help="Fraction of base substitutions to introduce on the genome. Default: 0.0001")

            optional.add_argument('-I', '--Indels', metavar='', type=float, default=0.001,
                                    help="Fraction of indels to introduce on the genome. Default: 0.001")
            optional.add_argument('-J', '--java_memory', metavar='', type=str, default="-Xmx16g",
                                    help="Java memory allocation, required for mutating the genome. Default: -Xmx16g")

            optional.add_argument('-e', '--error', action='store_true',
                                    help="Introduce sequencing errors ( Uses ART on the background)")

            optional.add_argument('-i', '--instrument', metavar='', type=str, default="HS25",
                                    help="Illumina sequecing instrument to simulate reads from (Default HiSeq 2500)")
            optional.add_argument('-ir', '--insRate', metavar='', type=float, default=0.00009,
                                    help="the first-read insertion rate (default: 0.00009). Requires -e")
            optional.add_argument('-ir2', '--insRate2', metavar='', type=float, default=0.00015,
                                    help="the second-read insertion rate (default: 0.00015). Requires -e")
            optional.add_argument('-dr', '--delRate', metavar='', type=float, default=0.00011,
                                    help="the first-read deletion rate (default:  0.00011). Requires -e")
            optional.add_argument('-dr2', '--delRate2', metavar='', type=float, default=0.00023,
                                    help="the second-read deletion rate (default: 0.00023). Requires -e")

            parser.print_help()

            time.sleep(0.01)
            sys.stderr.write(
                "\nNo input input given to Simulate, be sure that you are providing the flags '-g' and '-N'"
                "\nExiting\n")

            sys.exit(0)


        if len(sys.argv[2:]) == 0:
            parser.print_help()
            time.sleep(0.01)
            sys.stderr.write("\nNo arguments given to Simulate. Exiting\n")




        return (parser)

def main():
    run = circle_map()
    pid = run.__getpid__()
    # clean
    os.system("rm -rf temp_files_%s" % pid)

if __name__ == '__main__':
    main()

