#!/usr/bin/env python

######################################
############# DEFINE ENV #############
######################################

# general module imports
import argparse, os
from argparse import RawTextHelpFormatter
import copy as cp
import pickle
import string
import shutil 
import random
import sys
from shutil import copyfile
import time

# get the cwd were all the scripts are 
CWD = "/".join(__file__.split("/")[0:-1]); sys.path.insert(0, CWD)

# define the module name
module_name = __file__.split("/")[-1]

# define the EnvDir where the environment is defined
EnvDir = "/".join(sys.executable.split("/")[0:-2])

# import functions
import sv_functions as fun

# import perSVade-specific modules

######################################
######################################
######################################


#################################### 
############## ARGS ################
####################################

description = """
Run parameter optimization for SV calling adapted to any input. This is based on simulations of SVs in a reference genome (--ref), which resemble (in terms of coverage, read length and insert size) an input sequencing dataset of aligned reads (provided through --sortedbam). The module outpus a .json file with the optimized parameters for the input and a report of the accuracy on the simulations. You can find more information about how this works in https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs#how-does-the-parameter-optimization-for-sv-calling-work.
"""

# mandatory args
parser = argparse.ArgumentParser(description=description, formatter_class=RawTextHelpFormatter)
parser.add_argument("-o", "--outdir", dest="outdir", action="store", required=True, help="Output directory.")
parser.add_argument("-r", "--ref", dest="ref", required=True, help="Reference genome. It has to end with .fasta.")
parser.add_argument("-sbam", "--sortedbam", dest="sortedbam", required=True, help="The path to the sorted .bam file, which should have a bam.bai file in the same dir. For example, if your bam file is called 'aligned_reads.bam', there should be an 'aligned_reads.bam.bai' as well.")
parser.add_argument("-mchr", "--mitochondrial_chromosome", dest="mitochondrial_chromosome", required=True, type=str, help="The name of the mitochondrial chromosome. If there is no mitochondria just put 'no_mitochondria'. If there is more than one mitochindrial scaffold, provide them as comma-sepparated IDs, like '--mitochondrial_chromosome chr_mito_1,chr_mito_2'.")
parser.add_argument("--repeats_file", dest="repeats_file", required=True, help="A file with the repeats of the reference genome, such as the file 'combined_repeats.tab' generated by perSVade infer_repeats. You may set '--repeats_file skip' if you don't want to consider repeats for SV filtering.")
parser.add_argument("--regions_SVsimulations", dest="regions_SVsimulations", required=True, help="A .bedpe file with the regions around which to perfomr the simulations of SVs. This may be obtained through the modules 'find_homologous_regions' or 'find_knownSVs_regions'. If you want simulations randomly distributed across the genome you can specify '--regions_SVsimulations random'.")
parser.add_argument("--simulation_ploidies", dest="simulation_ploidies", type=str, required=True, help='A comma-sepparated string of the ploidies to simulate. It can have any of "haploid", "diploid_homo", "diploid_hetero", "ref:2_var:1", "ref:3_var:1", "ref:4_var:1", "ref:5_var:1", "ref:9_var:1", "ref:19_var:1", "ref:99_var:1". For example, if you are running on ploidy 1 and 2, you can specify "--simulation_ploidies haploid" or "--simulation_ploidies diploid_hetero", respectively. If you specify more than one ploidy (i.e. "--simulation_ploidies diploid_hetero,haploid") the parameters will be optimized for both types of ploidies. If you specify "--simulation_ploidies ref:3_var:1" this module will simulate a tetraploid genome where 3 chromosomes are like the reference and 1 chromosome has SVs.')

# simulations optional args
parser.add_argument("--simulation_chromosomes", dest="simulation_chromosomes", type=str, default=None, help="A comma-sepparated set of chromosomes (i.e.: chr1,chr2,chr3) in which to perform simulations. By default it takes all the chromosomes. This can be useful to speed up the computation.")
parser.add_argument("--nvars", dest="nvars", default=50, type=int, help="Number of variants to simulate for each SVtype.")
parser.add_argument("--nsimulations", dest="nsimulations", default=2, type=int, help="The number of 'replicate' simulations that will be produced.")
parser.add_argument("--range_filtering_benchmark", dest="range_filtering_benchmark", type=str, default="theoretically_meaningful", help='The range of parameters that should be tested in the SV optimisation pipeline. It can be any of large, medium, small, theoretically_meaningful, theoretically_meaningful_NoFilterRepeats or single. ')

# general optional args
parser.add_argument("--max_median_insert_size", dest="max_median_insert_size", default=10000, type=int, help="This module calculates the median insert size of the input sorted bam. This value is expected to be low (i.e. < 1000). If it is very high it suggests that the input .bam is somehow corrupt. This argument allows you to control what is considered as a very high insert size (default 10,000 bp). For example, if you want to allow a higher threshold you may specify '--max_median_insert_size 20000'.")
parser.add_argument("--min_chromosome_len", dest="min_chromosome_len", default=100000, type=int, help="The minimum length to consider chromosomes from the provided fasta for calculating the window length (used in may steps of perSVade to parallelize across fractions of the genome).")
parser.add_argument("--replace", dest="replace", action="store_true", help="Re-run all the steps by deleting the output directory.")
parser.add_argument("--verbose", dest="verbose", action="store_true", default=False, help="Print a verbose log.")
parser.add_argument("--tmpdir", dest="tmpdir", default=None, help="A full path to a directory where to write intermediate files. This is useful if you are running on a cluster that has some directories that have higher writing speed than others.")

# resources
parser.add_argument("--fraction_available_mem", dest="fraction_available_mem", default=None, type=float, help="This pipeline calculates the available RAM for several steps, and it may not work well in some systems (i.e. HPC clusters). This parameter allows you to correct possible errors. If --fraction_available_mem is not provided (default behavior), this pipeline will calculate the available RAM by filling the memory, which may give errors. If you want to use all the available memory you should specify --fraction_available_mem 1.0. See the FAQ 'How does the --fraction_available_mem work?' from https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs for more info.")

parser.add_argument("-thr", "--threads", dest="threads", default=16, type=int, help="Number of threads, Default: 16")
parser.add_argument("--fractionRAM_to_dedicate", dest="fractionRAM_to_dedicate", type=float,  default=0.5, help="This is the fraction of the available memory that will be used by several java programs that require a heap size. By default we set this to 0.5 to not overload the system.")

opt = parser.parse_args()

####################################
####################################
####################################

#################################
########### MAIN CODE ###########
#################################

# remove outdir if replace, and set replace to False
if opt.replace is True: fun.delete_folder(opt.outdir)
opt.replace = False

# make the outdir
fun.make_folder(opt.outdir)

# exit if the final file exists
final_file = "%s/perSVade_finished_file.txt"%opt.outdir

if not fun.file_is_empty(final_file): 
    fun.print_with_runtime("WARNING: %s exists, suggesting that perSVade was already  run in this folder. Remove this file if you want this command to work. Exiting..."%final_file)
    sys.exit(0)


# define the start time
start_time = time.time()

# define the verbosity. If opt.verbose is False, none of the 'print' statements of sv_functions will have an effect
fun.printing_verbose_mode = opt.verbose

# define a file that will contain all the cmds ran
fun.log_file_all_cmds = "%s/all_cmds.txt"%opt.outdir
if fun.file_is_empty(fun.log_file_all_cmds): open(fun.log_file_all_cmds, "w").write("# These are all the cmds:\n")

# get sample name
sample_name = fun.get_sampleName_from_perSVade_outdir(opt.outdir)

####### SET RESOURCES ########

# define the fraction of RAM to dedicate
if opt.fractionRAM_to_dedicate>0.95: raise ValueError("You are using >95 pct of the systems RAM, which is dangerous")
fun.fractionRAM_to_dedicate = opt.fractionRAM_to_dedicate

# define the fraction of available mem
fun.fraction_available_mem = opt.fraction_available_mem
if opt.fraction_available_mem is None: fun.print_with_runtime("WARNING: You did not specify how much RAM should be used through --fraction_available_mem. perSVade will calculate this by filling the memory, which may be dangerous. If you want to use all the allocated memory you should specify --fraction_available_mem 1.0")

# print the available resources
real_available_threads = fun.get_available_threads(opt.outdir)
if opt.threads>real_available_threads:  fun.print_with_runtime("WARNING: There are %i available threads, and you required %i."%(real_available_threads, opt.threads))

available_Gb_RAM = fun.get_availableGbRAM(opt.outdir)
fun.print_with_runtime("Running perSVade %s into %s with %.3f Gb of RAM and %i cores"%(module_name, opt.outdir, available_Gb_RAM, opt.threads))

##############################

# check that the tmpdir exists
if opt.tmpdir is not None:
    if not os.path.isdir(opt.tmpdir): raise ValueError("The folder that you specified with --tmpdir does not exist")

# prepare the reference genome
opt.ref, reference_genome_dir = fun.prepare_reference_genome_for_perSVade(opt.ref, opt.outdir, opt.mitochondrial_chromosome, opt.simulation_chromosomes, opt.replace)
fun.window_l = fun.get_perSVade_window_l(opt.ref, opt.mitochondrial_chromosome, opt.min_chromosome_len)

# get the bam file under output
sorted_bam, index_bam = fun.get_sorted_bam_in_outdir(opt.sortedbam, opt.outdir)

# prepare the repeats file
fun.prepare_repeats_file_for_perSVade(opt.repeats_file, opt.ref)

# debug the regions_SVsimulations
real_bedpe_breakpoints = fun.get_debugged_real_bedpe_breakpoints_for_SVsimulation(opt.regions_SVsimulations, opt.outdir, opt.ref)

# get the simulation ploidies
simulation_ploidies = opt.simulation_ploidies.split(",")

# log
fun.print_with_runtime("Running parameter optimization")

# define outdirs (this is done like this to be able to use fun.clean_perSVade_outdir at the end)
SVdetection_outdir = "%s/SVdetection_output"%opt.outdir; fun.make_folder(SVdetection_outdir)
outdir_gridss_final = "%s/final_gridss_running"%SVdetection_outdir; fun.make_folder(outdir_gridss_final)

# clean reference
fun.clean_reference_genome_windows_files(opt.ref)

# calculate the insert size
median_insert_size, median_insert_size_sd  = fun.get_insert_size_distribution(sorted_bam, replace=opt.replace, threads=opt.threads)
fun.print_with_runtime("The median insert size is %i, with an absolute deviation of %i"%(median_insert_size, median_insert_size_sd))

# debug very high insert sizes
if median_insert_size>opt.max_median_insert_size: raise ValueError("The median insert size is very high (%s). Maybe there was an error with the bam file generation."%(median_insert_size))

# get the optimised parameters
parameter_optimisation_dir = "%s/parameter_optimisation"%SVdetection_outdir; fun.make_folder(parameter_optimisation_dir)

gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict, max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup, df_cross_benchmark_best = fun.get_best_parameters_for_GridssClove_run(sorted_bam, opt.ref, parameter_optimisation_dir, threads=opt.threads, replace=opt.replace, n_simulated_genomes=opt.nsimulations, mitochondrial_chromosome=opt.mitochondrial_chromosome, simulation_ploidies=simulation_ploidies, range_filtering_benchmark=opt.range_filtering_benchmark, nvars=opt.nvars, real_bedpe_breakpoints=real_bedpe_breakpoints, median_insert_size=median_insert_size, median_insert_size_sd=median_insert_size_sd, tmpdir=opt.tmpdir, simulation_chromosomes=opt.simulation_chromosomes)

# write these optimal parameters
json_file = "%s/optimized_parameters.json"%opt.outdir
fun.write_gridss_parameters_as_json(gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict, max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup, json_file, replace=opt.replace)

#################################
#################################
#################################


##################################
########## CLEAN OUTPUT ##########
##################################

# log
fun.print_with_runtime("cleaning")

# clean normal perSVade outdir
fun.clean_perSVade_outdir(opt.outdir)

# remove extra files from outdir
for f in os.listdir(opt.outdir): 
    file_path = "%s/%s"%(opt.outdir, f)
    if os.path.isfile(file_path) and f not in {fun.get_file(json_file), "all_cmds.txt"}: fun.remove_file(file_path)

# move to outir from parameter_optimisation_dir
for f in os.listdir(parameter_optimisation_dir): 
    file_path = "%s/%s"%(parameter_optimisation_dir, f)
    os.rename(file_path, "%s/%s"%(opt.outdir, f))

# remove specific dirs
for f in [reference_genome_dir, SVdetection_outdir]: fun.delete_folder(f)

# move the parameter accuracy to the outdir
df_accuracy = fun.get_tab_as_df_or_empty_df("%s/benchmarking_all_filters_for_all_genomes_and_ploidies/df_cross_benchmark_best.tab"%(opt.outdir))
df_accuracy = df_accuracy[["test_simName", "test_ploidy", "test_svtype", "precision", "recall", "Fvalue", "nevents", "TP", "FP", "FN"]].rename(columns={"test_simName":"simulation_genome", "test_ploidy":"simulation_ploidy", "test_svtype":"SV_type"})
fun.save_df_as_tab(df_accuracy, "%s/optimized_parameters_accuracy.tab"%opt.outdir)

# wite final file
fun.generate_final_file_report_one_module(final_file, start_time, time.time())

# print the message
fun.print_with_runtime("perSVade %s finished correctly"%module_name)

##################################
##################################
##################################

d61f77738c2d102281bb5a8929cecf65ec7ea33f
<mask token>

<mask token>
sys.path.insert(0, CWD)
<mask token>
parser.add_argument('-o', '--outdir', dest='outdir', action='store',
    required=True, help='Output directory.')
parser.add_argument('-r', '--ref', dest='ref', required=True, help=
    'Reference genome. It has to end with .fasta.')
parser.add_argument('-sbam', '--sortedbam', dest='sortedbam', required=True,
    help=
    "The path to the sorted .bam file, which should have a bam.bai file in the same dir. For example, if your bam file is called 'aligned_reads.bam', there should be an 'aligned_reads.bam.bai' as well."
    )
parser.add_argument('-mchr', '--mitochondrial_chromosome', dest=
    'mitochondrial_chromosome', required=True, type=str, help=
    "The name of the mitochondrial chromosome. If there is no mitochondria just put 'no_mitochondria'. If there is more than one mitochindrial scaffold, provide them as comma-sepparated IDs, like '--mitochondrial_chromosome chr_mito_1,chr_mito_2'."
    )
parser.add_argument('--repeats_file', dest='repeats_file', required=True,
    help=
    "A file with the repeats of the reference genome, such as the file 'combined_repeats.tab' generated by perSVade infer_repeats. You may set '--repeats_file skip' if you don't want to consider repeats for SV filtering."
    )
parser.add_argument('--regions_SVsimulations', dest='regions_SVsimulations',
    required=True, help=
    "A .bedpe file with the regions around which to perfomr the simulations of SVs. This may be obtained through the modules 'find_homologous_regions' or 'find_knownSVs_regions'. If you want simulations randomly distributed across the genome you can specify '--regions_SVsimulations random'."
    )
parser.add_argument('--simulation_ploidies', dest='simulation_ploidies',
    type=str, required=True, help=
    'A comma-sepparated string of the ploidies to simulate. It can have any of "haploid", "diploid_homo", "diploid_hetero", "ref:2_var:1", "ref:3_var:1", "ref:4_var:1", "ref:5_var:1", "ref:9_var:1", "ref:19_var:1", "ref:99_var:1". For example, if you are running on ploidy 1 and 2, you can specify "--simulation_ploidies haploid" or "--simulation_ploidies diploid_hetero", respectively. If you specify more than one ploidy (i.e. "--simulation_ploidies diploid_hetero,haploid") the parameters will be optimized for both types of ploidies. If you specify "--simulation_ploidies ref:3_var:1" this module will simulate a tetraploid genome where 3 chromosomes are like the reference and 1 chromosome has SVs.'
    )
parser.add_argument('--simulation_chromosomes', dest=
    'simulation_chromosomes', type=str, default=None, help=
    'A comma-sepparated set of chromosomes (i.e.: chr1,chr2,chr3) in which to perform simulations. By default it takes all the chromosomes. This can be useful to speed up the computation.'
    )
parser.add_argument('--nvars', dest='nvars', default=50, type=int, help=
    'Number of variants to simulate for each SVtype.')
parser.add_argument('--nsimulations', dest='nsimulations', default=2, type=
    int, help="The number of 'replicate' simulations that will be produced.")
parser.add_argument('--range_filtering_benchmark', dest=
    'range_filtering_benchmark', type=str, default=
    'theoretically_meaningful', help=
    'The range of parameters that should be tested in the SV optimisation pipeline. It can be any of large, medium, small, theoretically_meaningful, theoretically_meaningful_NoFilterRepeats or single. '
    )
parser.add_argument('--max_median_insert_size', dest=
    'max_median_insert_size', default=10000, type=int, help=
    "This module calculates the median insert size of the input sorted bam. This value is expected to be low (i.e. < 1000). If it is very high it suggests that the input .bam is somehow corrupt. This argument allows you to control what is considered as a very high insert size (default 10,000 bp). For example, if you want to allow a higher threshold you may specify '--max_median_insert_size 20000'."
    )
parser.add_argument('--min_chromosome_len', dest='min_chromosome_len',
    default=100000, type=int, help=
    'The minimum length to consider chromosomes from the provided fasta for calculating the window length (used in may steps of perSVade to parallelize across fractions of the genome).'
    )
parser.add_argument('--replace', dest='replace', action='store_true', help=
    'Re-run all the steps by deleting the output directory.')
parser.add_argument('--verbose', dest='verbose', action='store_true',
    default=False, help='Print a verbose log.')
parser.add_argument('--tmpdir', dest='tmpdir', default=None, help=
    'A full path to a directory where to write intermediate files. This is useful if you are running on a cluster that has some directories that have higher writing speed than others.'
    )
parser.add_argument('--fraction_available_mem', dest=
    'fraction_available_mem', default=None, type=float, help=
    "This pipeline calculates the available RAM for several steps, and it may not work well in some systems (i.e. HPC clusters). This parameter allows you to correct possible errors. If --fraction_available_mem is not provided (default behavior), this pipeline will calculate the available RAM by filling the memory, which may give errors. If you want to use all the available memory you should specify --fraction_available_mem 1.0. See the FAQ 'How does the --fraction_available_mem work?' from https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs for more info."
    )
parser.add_argument('-thr', '--threads', dest='threads', default=16, type=
    int, help='Number of threads, Default: 16')
parser.add_argument('--fractionRAM_to_dedicate', dest=
    'fractionRAM_to_dedicate', type=float, default=0.5, help=
    'This is the fraction of the available memory that will be used by several java programs that require a heap size. By default we set this to 0.5 to not overload the system.'
    )
<mask token>
if opt.replace is True:
    fun.delete_folder(opt.outdir)
<mask token>
fun.make_folder(opt.outdir)
<mask token>
if not fun.file_is_empty(final_file):
    fun.print_with_runtime(
        'WARNING: %s exists, suggesting that perSVade was already  run in this folder. Remove this file if you want this command to work. Exiting...'
         % final_file)
    sys.exit(0)
<mask token>
if fun.file_is_empty(fun.log_file_all_cmds):
    open(fun.log_file_all_cmds, 'w').write('# These are all the cmds:\n')
<mask token>
if opt.fractionRAM_to_dedicate > 0.95:
    raise ValueError(
        'You are using >95 pct of the systems RAM, which is dangerous')
<mask token>
if opt.fraction_available_mem is None:
    fun.print_with_runtime(
        'WARNING: You did not specify how much RAM should be used through --fraction_available_mem. perSVade will calculate this by filling the memory, which may be dangerous. If you want to use all the allocated memory you should specify --fraction_available_mem 1.0'
        )
<mask token>
if opt.threads > real_available_threads:
    fun.print_with_runtime(
        'WARNING: There are %i available threads, and you required %i.' % (
        real_available_threads, opt.threads))
<mask token>
fun.print_with_runtime(
    'Running perSVade %s into %s with %.3f Gb of RAM and %i cores' % (
    module_name, opt.outdir, available_Gb_RAM, opt.threads))
if opt.tmpdir is not None:
    if not os.path.isdir(opt.tmpdir):
        raise ValueError(
            'The folder that you specified with --tmpdir does not exist')
<mask token>
fun.prepare_repeats_file_for_perSVade(opt.repeats_file, opt.ref)
<mask token>
fun.print_with_runtime('Running parameter optimization')
<mask token>
fun.make_folder(SVdetection_outdir)
<mask token>
fun.make_folder(outdir_gridss_final)
fun.clean_reference_genome_windows_files(opt.ref)
<mask token>
fun.print_with_runtime(
    'The median insert size is %i, with an absolute deviation of %i' % (
    median_insert_size, median_insert_size_sd))
if median_insert_size > opt.max_median_insert_size:
    raise ValueError(
        'The median insert size is very high (%s). Maybe there was an error with the bam file generation.'
         % median_insert_size)
<mask token>
fun.make_folder(parameter_optimisation_dir)
<mask token>
fun.write_gridss_parameters_as_json(gridss_blacklisted_regions,
    gridss_maxcoverage, gridss_filters_dict,
    max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup,
    json_file, replace=opt.replace)
fun.print_with_runtime('cleaning')
fun.clean_perSVade_outdir(opt.outdir)
for f in os.listdir(opt.outdir):
    file_path = '%s/%s' % (opt.outdir, f)
    if os.path.isfile(file_path) and f not in {fun.get_file(json_file),
        'all_cmds.txt'}:
        fun.remove_file(file_path)
for f in os.listdir(parameter_optimisation_dir):
    file_path = '%s/%s' % (parameter_optimisation_dir, f)
    os.rename(file_path, '%s/%s' % (opt.outdir, f))
for f in [reference_genome_dir, SVdetection_outdir]:
    fun.delete_folder(f)
<mask token>
fun.save_df_as_tab(df_accuracy, '%s/optimized_parameters_accuracy.tab' %
    opt.outdir)
fun.generate_final_file_report_one_module(final_file, start_time, time.time())
fun.print_with_runtime('perSVade %s finished correctly' % module_name)

<mask token>
CWD = '/'.join(__file__.split('/')[0:-1])
sys.path.insert(0, CWD)
module_name = __file__.split('/')[-1]
EnvDir = '/'.join(sys.executable.split('/')[0:-2])
<mask token>
description = """
Run parameter optimization for SV calling adapted to any input. This is based on simulations of SVs in a reference genome (--ref), which resemble (in terms of coverage, read length and insert size) an input sequencing dataset of aligned reads (provided through --sortedbam). The module outpus a .json file with the optimized parameters for the input and a report of the accuracy on the simulations. You can find more information about how this works in https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs#how-does-the-parameter-optimization-for-sv-calling-work.
"""
parser = argparse.ArgumentParser(description=description, formatter_class=
    RawTextHelpFormatter)
parser.add_argument('-o', '--outdir', dest='outdir', action='store',
    required=True, help='Output directory.')
parser.add_argument('-r', '--ref', dest='ref', required=True, help=
    'Reference genome. It has to end with .fasta.')
parser.add_argument('-sbam', '--sortedbam', dest='sortedbam', required=True,
    help=
    "The path to the sorted .bam file, which should have a bam.bai file in the same dir. For example, if your bam file is called 'aligned_reads.bam', there should be an 'aligned_reads.bam.bai' as well."
    )
parser.add_argument('-mchr', '--mitochondrial_chromosome', dest=
    'mitochondrial_chromosome', required=True, type=str, help=
    "The name of the mitochondrial chromosome. If there is no mitochondria just put 'no_mitochondria'. If there is more than one mitochindrial scaffold, provide them as comma-sepparated IDs, like '--mitochondrial_chromosome chr_mito_1,chr_mito_2'."
    )
parser.add_argument('--repeats_file', dest='repeats_file', required=True,
    help=
    "A file with the repeats of the reference genome, such as the file 'combined_repeats.tab' generated by perSVade infer_repeats. You may set '--repeats_file skip' if you don't want to consider repeats for SV filtering."
    )
parser.add_argument('--regions_SVsimulations', dest='regions_SVsimulations',
    required=True, help=
    "A .bedpe file with the regions around which to perfomr the simulations of SVs. This may be obtained through the modules 'find_homologous_regions' or 'find_knownSVs_regions'. If you want simulations randomly distributed across the genome you can specify '--regions_SVsimulations random'."
    )
parser.add_argument('--simulation_ploidies', dest='simulation_ploidies',
    type=str, required=True, help=
    'A comma-sepparated string of the ploidies to simulate. It can have any of "haploid", "diploid_homo", "diploid_hetero", "ref:2_var:1", "ref:3_var:1", "ref:4_var:1", "ref:5_var:1", "ref:9_var:1", "ref:19_var:1", "ref:99_var:1". For example, if you are running on ploidy 1 and 2, you can specify "--simulation_ploidies haploid" or "--simulation_ploidies diploid_hetero", respectively. If you specify more than one ploidy (i.e. "--simulation_ploidies diploid_hetero,haploid") the parameters will be optimized for both types of ploidies. If you specify "--simulation_ploidies ref:3_var:1" this module will simulate a tetraploid genome where 3 chromosomes are like the reference and 1 chromosome has SVs.'
    )
parser.add_argument('--simulation_chromosomes', dest=
    'simulation_chromosomes', type=str, default=None, help=
    'A comma-sepparated set of chromosomes (i.e.: chr1,chr2,chr3) in which to perform simulations. By default it takes all the chromosomes. This can be useful to speed up the computation.'
    )
parser.add_argument('--nvars', dest='nvars', default=50, type=int, help=
    'Number of variants to simulate for each SVtype.')
parser.add_argument('--nsimulations', dest='nsimulations', default=2, type=
    int, help="The number of 'replicate' simulations that will be produced.")
parser.add_argument('--range_filtering_benchmark', dest=
    'range_filtering_benchmark', type=str, default=
    'theoretically_meaningful', help=
    'The range of parameters that should be tested in the SV optimisation pipeline. It can be any of large, medium, small, theoretically_meaningful, theoretically_meaningful_NoFilterRepeats or single. '
    )
parser.add_argument('--max_median_insert_size', dest=
    'max_median_insert_size', default=10000, type=int, help=
    "This module calculates the median insert size of the input sorted bam. This value is expected to be low (i.e. < 1000). If it is very high it suggests that the input .bam is somehow corrupt. This argument allows you to control what is considered as a very high insert size (default 10,000 bp). For example, if you want to allow a higher threshold you may specify '--max_median_insert_size 20000'."
    )
parser.add_argument('--min_chromosome_len', dest='min_chromosome_len',
    default=100000, type=int, help=
    'The minimum length to consider chromosomes from the provided fasta for calculating the window length (used in may steps of perSVade to parallelize across fractions of the genome).'
    )
parser.add_argument('--replace', dest='replace', action='store_true', help=
    'Re-run all the steps by deleting the output directory.')
parser.add_argument('--verbose', dest='verbose', action='store_true',
    default=False, help='Print a verbose log.')
parser.add_argument('--tmpdir', dest='tmpdir', default=None, help=
    'A full path to a directory where to write intermediate files. This is useful if you are running on a cluster that has some directories that have higher writing speed than others.'
    )
parser.add_argument('--fraction_available_mem', dest=
    'fraction_available_mem', default=None, type=float, help=
    "This pipeline calculates the available RAM for several steps, and it may not work well in some systems (i.e. HPC clusters). This parameter allows you to correct possible errors. If --fraction_available_mem is not provided (default behavior), this pipeline will calculate the available RAM by filling the memory, which may give errors. If you want to use all the available memory you should specify --fraction_available_mem 1.0. See the FAQ 'How does the --fraction_available_mem work?' from https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs for more info."
    )
parser.add_argument('-thr', '--threads', dest='threads', default=16, type=
    int, help='Number of threads, Default: 16')
parser.add_argument('--fractionRAM_to_dedicate', dest=
    'fractionRAM_to_dedicate', type=float, default=0.5, help=
    'This is the fraction of the available memory that will be used by several java programs that require a heap size. By default we set this to 0.5 to not overload the system.'
    )
opt = parser.parse_args()
if opt.replace is True:
    fun.delete_folder(opt.outdir)
opt.replace = False
fun.make_folder(opt.outdir)
final_file = '%s/perSVade_finished_file.txt' % opt.outdir
if not fun.file_is_empty(final_file):
    fun.print_with_runtime(
        'WARNING: %s exists, suggesting that perSVade was already  run in this folder. Remove this file if you want this command to work. Exiting...'
         % final_file)
    sys.exit(0)
start_time = time.time()
fun.printing_verbose_mode = opt.verbose
fun.log_file_all_cmds = '%s/all_cmds.txt' % opt.outdir
if fun.file_is_empty(fun.log_file_all_cmds):
    open(fun.log_file_all_cmds, 'w').write('# These are all the cmds:\n')
sample_name = fun.get_sampleName_from_perSVade_outdir(opt.outdir)
if opt.fractionRAM_to_dedicate > 0.95:
    raise ValueError(
        'You are using >95 pct of the systems RAM, which is dangerous')
fun.fractionRAM_to_dedicate = opt.fractionRAM_to_dedicate
fun.fraction_available_mem = opt.fraction_available_mem
if opt.fraction_available_mem is None:
    fun.print_with_runtime(
        'WARNING: You did not specify how much RAM should be used through --fraction_available_mem. perSVade will calculate this by filling the memory, which may be dangerous. If you want to use all the allocated memory you should specify --fraction_available_mem 1.0'
        )
real_available_threads = fun.get_available_threads(opt.outdir)
if opt.threads > real_available_threads:
    fun.print_with_runtime(
        'WARNING: There are %i available threads, and you required %i.' % (
        real_available_threads, opt.threads))
available_Gb_RAM = fun.get_availableGbRAM(opt.outdir)
fun.print_with_runtime(
    'Running perSVade %s into %s with %.3f Gb of RAM and %i cores' % (
    module_name, opt.outdir, available_Gb_RAM, opt.threads))
if opt.tmpdir is not None:
    if not os.path.isdir(opt.tmpdir):
        raise ValueError(
            'The folder that you specified with --tmpdir does not exist')
opt.ref, reference_genome_dir = fun.prepare_reference_genome_for_perSVade(opt
    .ref, opt.outdir, opt.mitochondrial_chromosome, opt.
    simulation_chromosomes, opt.replace)
fun.window_l = fun.get_perSVade_window_l(opt.ref, opt.
    mitochondrial_chromosome, opt.min_chromosome_len)
sorted_bam, index_bam = fun.get_sorted_bam_in_outdir(opt.sortedbam, opt.outdir)
fun.prepare_repeats_file_for_perSVade(opt.repeats_file, opt.ref)
real_bedpe_breakpoints = (fun.
    get_debugged_real_bedpe_breakpoints_for_SVsimulation(opt.
    regions_SVsimulations, opt.outdir, opt.ref))
simulation_ploidies = opt.simulation_ploidies.split(',')
fun.print_with_runtime('Running parameter optimization')
SVdetection_outdir = '%s/SVdetection_output' % opt.outdir
fun.make_folder(SVdetection_outdir)
outdir_gridss_final = '%s/final_gridss_running' % SVdetection_outdir
fun.make_folder(outdir_gridss_final)
fun.clean_reference_genome_windows_files(opt.ref)
median_insert_size, median_insert_size_sd = fun.get_insert_size_distribution(
    sorted_bam, replace=opt.replace, threads=opt.threads)
fun.print_with_runtime(
    'The median insert size is %i, with an absolute deviation of %i' % (
    median_insert_size, median_insert_size_sd))
if median_insert_size > opt.max_median_insert_size:
    raise ValueError(
        'The median insert size is very high (%s). Maybe there was an error with the bam file generation.'
         % median_insert_size)
parameter_optimisation_dir = '%s/parameter_optimisation' % SVdetection_outdir
fun.make_folder(parameter_optimisation_dir)
(gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict,
    max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup,
    df_cross_benchmark_best) = (fun.get_best_parameters_for_GridssClove_run
    (sorted_bam, opt.ref, parameter_optimisation_dir, threads=opt.threads,
    replace=opt.replace, n_simulated_genomes=opt.nsimulations,
    mitochondrial_chromosome=opt.mitochondrial_chromosome,
    simulation_ploidies=simulation_ploidies, range_filtering_benchmark=opt.
    range_filtering_benchmark, nvars=opt.nvars, real_bedpe_breakpoints=
    real_bedpe_breakpoints, median_insert_size=median_insert_size,
    median_insert_size_sd=median_insert_size_sd, tmpdir=opt.tmpdir,
    simulation_chromosomes=opt.simulation_chromosomes))
json_file = '%s/optimized_parameters.json' % opt.outdir
fun.write_gridss_parameters_as_json(gridss_blacklisted_regions,
    gridss_maxcoverage, gridss_filters_dict,
    max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup,
    json_file, replace=opt.replace)
fun.print_with_runtime('cleaning')
fun.clean_perSVade_outdir(opt.outdir)
for f in os.listdir(opt.outdir):
    file_path = '%s/%s' % (opt.outdir, f)
    if os.path.isfile(file_path) and f not in {fun.get_file(json_file),
        'all_cmds.txt'}:
        fun.remove_file(file_path)
for f in os.listdir(parameter_optimisation_dir):
    file_path = '%s/%s' % (parameter_optimisation_dir, f)
    os.rename(file_path, '%s/%s' % (opt.outdir, f))
for f in [reference_genome_dir, SVdetection_outdir]:
    fun.delete_folder(f)
df_accuracy = fun.get_tab_as_df_or_empty_df(
    '%s/benchmarking_all_filters_for_all_genomes_and_ploidies/df_cross_benchmark_best.tab'
     % opt.outdir)
df_accuracy = df_accuracy[['test_simName', 'test_ploidy', 'test_svtype',
    'precision', 'recall', 'Fvalue', 'nevents', 'TP', 'FP', 'FN']].rename(
    columns={'test_simName': 'simulation_genome', 'test_ploidy':
    'simulation_ploidy', 'test_svtype': 'SV_type'})
fun.save_df_as_tab(df_accuracy, '%s/optimized_parameters_accuracy.tab' %
    opt.outdir)
fun.generate_final_file_report_one_module(final_file, start_time, time.time())
fun.print_with_runtime('perSVade %s finished correctly' % module_name)

import argparse, os
from argparse import RawTextHelpFormatter
import copy as cp
import pickle
import string
import shutil
import random
import sys
from shutil import copyfile
import time
CWD = '/'.join(__file__.split('/')[0:-1])
sys.path.insert(0, CWD)
module_name = __file__.split('/')[-1]
EnvDir = '/'.join(sys.executable.split('/')[0:-2])
import sv_functions as fun
description = """
Run parameter optimization for SV calling adapted to any input. This is based on simulations of SVs in a reference genome (--ref), which resemble (in terms of coverage, read length and insert size) an input sequencing dataset of aligned reads (provided through --sortedbam). The module outpus a .json file with the optimized parameters for the input and a report of the accuracy on the simulations. You can find more information about how this works in https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs#how-does-the-parameter-optimization-for-sv-calling-work.
"""
parser = argparse.ArgumentParser(description=description, formatter_class=
    RawTextHelpFormatter)
parser.add_argument('-o', '--outdir', dest='outdir', action='store',
    required=True, help='Output directory.')
parser.add_argument('-r', '--ref', dest='ref', required=True, help=
    'Reference genome. It has to end with .fasta.')
parser.add_argument('-sbam', '--sortedbam', dest='sortedbam', required=True,
    help=
    "The path to the sorted .bam file, which should have a bam.bai file in the same dir. For example, if your bam file is called 'aligned_reads.bam', there should be an 'aligned_reads.bam.bai' as well."
    )
parser.add_argument('-mchr', '--mitochondrial_chromosome', dest=
    'mitochondrial_chromosome', required=True, type=str, help=
    "The name of the mitochondrial chromosome. If there is no mitochondria just put 'no_mitochondria'. If there is more than one mitochindrial scaffold, provide them as comma-sepparated IDs, like '--mitochondrial_chromosome chr_mito_1,chr_mito_2'."
    )
parser.add_argument('--repeats_file', dest='repeats_file', required=True,
    help=
    "A file with the repeats of the reference genome, such as the file 'combined_repeats.tab' generated by perSVade infer_repeats. You may set '--repeats_file skip' if you don't want to consider repeats for SV filtering."
    )
parser.add_argument('--regions_SVsimulations', dest='regions_SVsimulations',
    required=True, help=
    "A .bedpe file with the regions around which to perfomr the simulations of SVs. This may be obtained through the modules 'find_homologous_regions' or 'find_knownSVs_regions'. If you want simulations randomly distributed across the genome you can specify '--regions_SVsimulations random'."
    )
parser.add_argument('--simulation_ploidies', dest='simulation_ploidies',
    type=str, required=True, help=
    'A comma-sepparated string of the ploidies to simulate. It can have any of "haploid", "diploid_homo", "diploid_hetero", "ref:2_var:1", "ref:3_var:1", "ref:4_var:1", "ref:5_var:1", "ref:9_var:1", "ref:19_var:1", "ref:99_var:1". For example, if you are running on ploidy 1 and 2, you can specify "--simulation_ploidies haploid" or "--simulation_ploidies diploid_hetero", respectively. If you specify more than one ploidy (i.e. "--simulation_ploidies diploid_hetero,haploid") the parameters will be optimized for both types of ploidies. If you specify "--simulation_ploidies ref:3_var:1" this module will simulate a tetraploid genome where 3 chromosomes are like the reference and 1 chromosome has SVs.'
    )
parser.add_argument('--simulation_chromosomes', dest=
    'simulation_chromosomes', type=str, default=None, help=
    'A comma-sepparated set of chromosomes (i.e.: chr1,chr2,chr3) in which to perform simulations. By default it takes all the chromosomes. This can be useful to speed up the computation.'
    )
parser.add_argument('--nvars', dest='nvars', default=50, type=int, help=
    'Number of variants to simulate for each SVtype.')
parser.add_argument('--nsimulations', dest='nsimulations', default=2, type=
    int, help="The number of 'replicate' simulations that will be produced.")
parser.add_argument('--range_filtering_benchmark', dest=
    'range_filtering_benchmark', type=str, default=
    'theoretically_meaningful', help=
    'The range of parameters that should be tested in the SV optimisation pipeline. It can be any of large, medium, small, theoretically_meaningful, theoretically_meaningful_NoFilterRepeats or single. '
    )
parser.add_argument('--max_median_insert_size', dest=
    'max_median_insert_size', default=10000, type=int, help=
    "This module calculates the median insert size of the input sorted bam. This value is expected to be low (i.e. < 1000). If it is very high it suggests that the input .bam is somehow corrupt. This argument allows you to control what is considered as a very high insert size (default 10,000 bp). For example, if you want to allow a higher threshold you may specify '--max_median_insert_size 20000'."
    )
parser.add_argument('--min_chromosome_len', dest='min_chromosome_len',
    default=100000, type=int, help=
    'The minimum length to consider chromosomes from the provided fasta for calculating the window length (used in may steps of perSVade to parallelize across fractions of the genome).'
    )
parser.add_argument('--replace', dest='replace', action='store_true', help=
    'Re-run all the steps by deleting the output directory.')
parser.add_argument('--verbose', dest='verbose', action='store_true',
    default=False, help='Print a verbose log.')
parser.add_argument('--tmpdir', dest='tmpdir', default=None, help=
    'A full path to a directory where to write intermediate files. This is useful if you are running on a cluster that has some directories that have higher writing speed than others.'
    )
parser.add_argument('--fraction_available_mem', dest=
    'fraction_available_mem', default=None, type=float, help=
    "This pipeline calculates the available RAM for several steps, and it may not work well in some systems (i.e. HPC clusters). This parameter allows you to correct possible errors. If --fraction_available_mem is not provided (default behavior), this pipeline will calculate the available RAM by filling the memory, which may give errors. If you want to use all the available memory you should specify --fraction_available_mem 1.0. See the FAQ 'How does the --fraction_available_mem work?' from https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs for more info."
    )
parser.add_argument('-thr', '--threads', dest='threads', default=16, type=
    int, help='Number of threads, Default: 16')
parser.add_argument('--fractionRAM_to_dedicate', dest=
    'fractionRAM_to_dedicate', type=float, default=0.5, help=
    'This is the fraction of the available memory that will be used by several java programs that require a heap size. By default we set this to 0.5 to not overload the system.'
    )
opt = parser.parse_args()
if opt.replace is True:
    fun.delete_folder(opt.outdir)
opt.replace = False
fun.make_folder(opt.outdir)
final_file = '%s/perSVade_finished_file.txt' % opt.outdir
if not fun.file_is_empty(final_file):
    fun.print_with_runtime(
        'WARNING: %s exists, suggesting that perSVade was already  run in this folder. Remove this file if you want this command to work. Exiting...'
         % final_file)
    sys.exit(0)
start_time = time.time()
fun.printing_verbose_mode = opt.verbose
fun.log_file_all_cmds = '%s/all_cmds.txt' % opt.outdir
if fun.file_is_empty(fun.log_file_all_cmds):
    open(fun.log_file_all_cmds, 'w').write('# These are all the cmds:\n')
sample_name = fun.get_sampleName_from_perSVade_outdir(opt.outdir)
if opt.fractionRAM_to_dedicate > 0.95:
    raise ValueError(
        'You are using >95 pct of the systems RAM, which is dangerous')
fun.fractionRAM_to_dedicate = opt.fractionRAM_to_dedicate
fun.fraction_available_mem = opt.fraction_available_mem
if opt.fraction_available_mem is None:
    fun.print_with_runtime(
        'WARNING: You did not specify how much RAM should be used through --fraction_available_mem. perSVade will calculate this by filling the memory, which may be dangerous. If you want to use all the allocated memory you should specify --fraction_available_mem 1.0'
        )
real_available_threads = fun.get_available_threads(opt.outdir)
if opt.threads > real_available_threads:
    fun.print_with_runtime(
        'WARNING: There are %i available threads, and you required %i.' % (
        real_available_threads, opt.threads))
available_Gb_RAM = fun.get_availableGbRAM(opt.outdir)
fun.print_with_runtime(
    'Running perSVade %s into %s with %.3f Gb of RAM and %i cores' % (
    module_name, opt.outdir, available_Gb_RAM, opt.threads))
if opt.tmpdir is not None:
    if not os.path.isdir(opt.tmpdir):
        raise ValueError(
            'The folder that you specified with --tmpdir does not exist')
opt.ref, reference_genome_dir = fun.prepare_reference_genome_for_perSVade(opt
    .ref, opt.outdir, opt.mitochondrial_chromosome, opt.
    simulation_chromosomes, opt.replace)
fun.window_l = fun.get_perSVade_window_l(opt.ref, opt.
    mitochondrial_chromosome, opt.min_chromosome_len)
sorted_bam, index_bam = fun.get_sorted_bam_in_outdir(opt.sortedbam, opt.outdir)
fun.prepare_repeats_file_for_perSVade(opt.repeats_file, opt.ref)
real_bedpe_breakpoints = (fun.
    get_debugged_real_bedpe_breakpoints_for_SVsimulation(opt.
    regions_SVsimulations, opt.outdir, opt.ref))
simulation_ploidies = opt.simulation_ploidies.split(',')
fun.print_with_runtime('Running parameter optimization')
SVdetection_outdir = '%s/SVdetection_output' % opt.outdir
fun.make_folder(SVdetection_outdir)
outdir_gridss_final = '%s/final_gridss_running' % SVdetection_outdir
fun.make_folder(outdir_gridss_final)
fun.clean_reference_genome_windows_files(opt.ref)
median_insert_size, median_insert_size_sd = fun.get_insert_size_distribution(
    sorted_bam, replace=opt.replace, threads=opt.threads)
fun.print_with_runtime(
    'The median insert size is %i, with an absolute deviation of %i' % (
    median_insert_size, median_insert_size_sd))
if median_insert_size > opt.max_median_insert_size:
    raise ValueError(
        'The median insert size is very high (%s). Maybe there was an error with the bam file generation.'
         % median_insert_size)
parameter_optimisation_dir = '%s/parameter_optimisation' % SVdetection_outdir
fun.make_folder(parameter_optimisation_dir)
(gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict,
    max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup,
    df_cross_benchmark_best) = (fun.get_best_parameters_for_GridssClove_run
    (sorted_bam, opt.ref, parameter_optimisation_dir, threads=opt.threads,
    replace=opt.replace, n_simulated_genomes=opt.nsimulations,
    mitochondrial_chromosome=opt.mitochondrial_chromosome,
    simulation_ploidies=simulation_ploidies, range_filtering_benchmark=opt.
    range_filtering_benchmark, nvars=opt.nvars, real_bedpe_breakpoints=
    real_bedpe_breakpoints, median_insert_size=median_insert_size,
    median_insert_size_sd=median_insert_size_sd, tmpdir=opt.tmpdir,
    simulation_chromosomes=opt.simulation_chromosomes))
json_file = '%s/optimized_parameters.json' % opt.outdir
fun.write_gridss_parameters_as_json(gridss_blacklisted_regions,
    gridss_maxcoverage, gridss_filters_dict,
    max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup,
    json_file, replace=opt.replace)
fun.print_with_runtime('cleaning')
fun.clean_perSVade_outdir(opt.outdir)
for f in os.listdir(opt.outdir):
    file_path = '%s/%s' % (opt.outdir, f)
    if os.path.isfile(file_path) and f not in {fun.get_file(json_file),
        'all_cmds.txt'}:
        fun.remove_file(file_path)
for f in os.listdir(parameter_optimisation_dir):
    file_path = '%s/%s' % (parameter_optimisation_dir, f)
    os.rename(file_path, '%s/%s' % (opt.outdir, f))
for f in [reference_genome_dir, SVdetection_outdir]:
    fun.delete_folder(f)
df_accuracy = fun.get_tab_as_df_or_empty_df(
    '%s/benchmarking_all_filters_for_all_genomes_and_ploidies/df_cross_benchmark_best.tab'
     % opt.outdir)
df_accuracy = df_accuracy[['test_simName', 'test_ploidy', 'test_svtype',
    'precision', 'recall', 'Fvalue', 'nevents', 'TP', 'FP', 'FN']].rename(
    columns={'test_simName': 'simulation_genome', 'test_ploidy':
    'simulation_ploidy', 'test_svtype': 'SV_type'})
fun.save_df_as_tab(df_accuracy, '%s/optimized_parameters_accuracy.tab' %
    opt.outdir)
fun.generate_final_file_report_one_module(final_file, start_time, time.time())
fun.print_with_runtime('perSVade %s finished correctly' % module_name)

#!/usr/bin/env python

######################################
############# DEFINE ENV #############
######################################

# general module imports
import argparse, os
from argparse import RawTextHelpFormatter
import copy as cp
import pickle
import string
import shutil 
import random
import sys
from shutil import copyfile
import time

# get the cwd were all the scripts are 
CWD = "/".join(__file__.split("/")[0:-1]); sys.path.insert(0, CWD)

# define the module name
module_name = __file__.split("/")[-1]

# define the EnvDir where the environment is defined
EnvDir = "/".join(sys.executable.split("/")[0:-2])

# import functions
import sv_functions as fun

# import perSVade-specific modules

######################################
######################################
######################################


#################################### 
############## ARGS ################
####################################

description = """
Run parameter optimization for SV calling adapted to any input. This is based on simulations of SVs in a reference genome (--ref), which resemble (in terms of coverage, read length and insert size) an input sequencing dataset of aligned reads (provided through --sortedbam). The module outpus a .json file with the optimized parameters for the input and a report of the accuracy on the simulations. You can find more information about how this works in https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs#how-does-the-parameter-optimization-for-sv-calling-work.
"""

# mandatory args
parser = argparse.ArgumentParser(description=description, formatter_class=RawTextHelpFormatter)
parser.add_argument("-o", "--outdir", dest="outdir", action="store", required=True, help="Output directory.")
parser.add_argument("-r", "--ref", dest="ref", required=True, help="Reference genome. It has to end with .fasta.")
parser.add_argument("-sbam", "--sortedbam", dest="sortedbam", required=True, help="The path to the sorted .bam file, which should have a bam.bai file in the same dir. For example, if your bam file is called 'aligned_reads.bam', there should be an 'aligned_reads.bam.bai' as well.")
parser.add_argument("-mchr", "--mitochondrial_chromosome", dest="mitochondrial_chromosome", required=True, type=str, help="The name of the mitochondrial chromosome. If there is no mitochondria just put 'no_mitochondria'. If there is more than one mitochindrial scaffold, provide them as comma-sepparated IDs, like '--mitochondrial_chromosome chr_mito_1,chr_mito_2'.")
parser.add_argument("--repeats_file", dest="repeats_file", required=True, help="A file with the repeats of the reference genome, such as the file 'combined_repeats.tab' generated by perSVade infer_repeats. You may set '--repeats_file skip' if you don't want to consider repeats for SV filtering.")
parser.add_argument("--regions_SVsimulations", dest="regions_SVsimulations", required=True, help="A .bedpe file with the regions around which to perfomr the simulations of SVs. This may be obtained through the modules 'find_homologous_regions' or 'find_knownSVs_regions'. If you want simulations randomly distributed across the genome you can specify '--regions_SVsimulations random'.")
parser.add_argument("--simulation_ploidies", dest="simulation_ploidies", type=str, required=True, help='A comma-sepparated string of the ploidies to simulate. It can have any of "haploid", "diploid_homo", "diploid_hetero", "ref:2_var:1", "ref:3_var:1", "ref:4_var:1", "ref:5_var:1", "ref:9_var:1", "ref:19_var:1", "ref:99_var:1". For example, if you are running on ploidy 1 and 2, you can specify "--simulation_ploidies haploid" or "--simulation_ploidies diploid_hetero", respectively. If you specify more than one ploidy (i.e. "--simulation_ploidies diploid_hetero,haploid") the parameters will be optimized for both types of ploidies. If you specify "--simulation_ploidies ref:3_var:1" this module will simulate a tetraploid genome where 3 chromosomes are like the reference and 1 chromosome has SVs.')

# simulations optional args
parser.add_argument("--simulation_chromosomes", dest="simulation_chromosomes", type=str, default=None, help="A comma-sepparated set of chromosomes (i.e.: chr1,chr2,chr3) in which to perform simulations. By default it takes all the chromosomes. This can be useful to speed up the computation.")
parser.add_argument("--nvars", dest="nvars", default=50, type=int, help="Number of variants to simulate for each SVtype.")
parser.add_argument("--nsimulations", dest="nsimulations", default=2, type=int, help="The number of 'replicate' simulations that will be produced.")
parser.add_argument("--range_filtering_benchmark", dest="range_filtering_benchmark", type=str, default="theoretically_meaningful", help='The range of parameters that should be tested in the SV optimisation pipeline. It can be any of large, medium, small, theoretically_meaningful, theoretically_meaningful_NoFilterRepeats or single. ')

# general optional args
parser.add_argument("--max_median_insert_size", dest="max_median_insert_size", default=10000, type=int, help="This module calculates the median insert size of the input sorted bam. This value is expected to be low (i.e. < 1000). If it is very high it suggests that the input .bam is somehow corrupt. This argument allows you to control what is considered as a very high insert size (default 10,000 bp). For example, if you want to allow a higher threshold you may specify '--max_median_insert_size 20000'.")
parser.add_argument("--min_chromosome_len", dest="min_chromosome_len", default=100000, type=int, help="The minimum length to consider chromosomes from the provided fasta for calculating the window length (used in may steps of perSVade to parallelize across fractions of the genome).")
parser.add_argument("--replace", dest="replace", action="store_true", help="Re-run all the steps by deleting the output directory.")
parser.add_argument("--verbose", dest="verbose", action="store_true", default=False, help="Print a verbose log.")
parser.add_argument("--tmpdir", dest="tmpdir", default=None, help="A full path to a directory where to write intermediate files. This is useful if you are running on a cluster that has some directories that have higher writing speed than others.")

# resources
parser.add_argument("--fraction_available_mem", dest="fraction_available_mem", default=None, type=float, help="This pipeline calculates the available RAM for several steps, and it may not work well in some systems (i.e. HPC clusters). This parameter allows you to correct possible errors. If --fraction_available_mem is not provided (default behavior), this pipeline will calculate the available RAM by filling the memory, which may give errors. If you want to use all the available memory you should specify --fraction_available_mem 1.0. See the FAQ 'How does the --fraction_available_mem work?' from https://github.com/Gabaldonlab/perSVade/wiki/8.-FAQs for more info.")

parser.add_argument("-thr", "--threads", dest="threads", default=16, type=int, help="Number of threads, Default: 16")
parser.add_argument("--fractionRAM_to_dedicate", dest="fractionRAM_to_dedicate", type=float,  default=0.5, help="This is the fraction of the available memory that will be used by several java programs that require a heap size. By default we set this to 0.5 to not overload the system.")

opt = parser.parse_args()

####################################
####################################
####################################

#################################
########### MAIN CODE ###########
#################################

# remove outdir if replace, and set replace to False
if opt.replace is True: fun.delete_folder(opt.outdir)
opt.replace = False

# make the outdir
fun.make_folder(opt.outdir)

# exit if the final file exists
final_file = "%s/perSVade_finished_file.txt"%opt.outdir

if not fun.file_is_empty(final_file): 
    fun.print_with_runtime("WARNING: %s exists, suggesting that perSVade was already  run in this folder. Remove this file if you want this command to work. Exiting..."%final_file)
    sys.exit(0)


# define the start time
start_time = time.time()

# define the verbosity. If opt.verbose is False, none of the 'print' statements of sv_functions will have an effect
fun.printing_verbose_mode = opt.verbose

# define a file that will contain all the cmds ran
fun.log_file_all_cmds = "%s/all_cmds.txt"%opt.outdir
if fun.file_is_empty(fun.log_file_all_cmds): open(fun.log_file_all_cmds, "w").write("# These are all the cmds:\n")

# get sample name
sample_name = fun.get_sampleName_from_perSVade_outdir(opt.outdir)

####### SET RESOURCES ########

# define the fraction of RAM to dedicate
if opt.fractionRAM_to_dedicate>0.95: raise ValueError("You are using >95 pct of the systems RAM, which is dangerous")
fun.fractionRAM_to_dedicate = opt.fractionRAM_to_dedicate

# define the fraction of available mem
fun.fraction_available_mem = opt.fraction_available_mem
if opt.fraction_available_mem is None: fun.print_with_runtime("WARNING: You did not specify how much RAM should be used through --fraction_available_mem. perSVade will calculate this by filling the memory, which may be dangerous. If you want to use all the allocated memory you should specify --fraction_available_mem 1.0")

# print the available resources
real_available_threads = fun.get_available_threads(opt.outdir)
if opt.threads>real_available_threads:  fun.print_with_runtime("WARNING: There are %i available threads, and you required %i."%(real_available_threads, opt.threads))

available_Gb_RAM = fun.get_availableGbRAM(opt.outdir)
fun.print_with_runtime("Running perSVade %s into %s with %.3f Gb of RAM and %i cores"%(module_name, opt.outdir, available_Gb_RAM, opt.threads))

##############################

# check that the tmpdir exists
if opt.tmpdir is not None:
    if not os.path.isdir(opt.tmpdir): raise ValueError("The folder that you specified with --tmpdir does not exist")

# prepare the reference genome
opt.ref, reference_genome_dir = fun.prepare_reference_genome_for_perSVade(opt.ref, opt.outdir, opt.mitochondrial_chromosome, opt.simulation_chromosomes, opt.replace)
fun.window_l = fun.get_perSVade_window_l(opt.ref, opt.mitochondrial_chromosome, opt.min_chromosome_len)

# get the bam file under output
sorted_bam, index_bam = fun.get_sorted_bam_in_outdir(opt.sortedbam, opt.outdir)

# prepare the repeats file
fun.prepare_repeats_file_for_perSVade(opt.repeats_file, opt.ref)

# debug the regions_SVsimulations
real_bedpe_breakpoints = fun.get_debugged_real_bedpe_breakpoints_for_SVsimulation(opt.regions_SVsimulations, opt.outdir, opt.ref)

# get the simulation ploidies
simulation_ploidies = opt.simulation_ploidies.split(",")

# log
fun.print_with_runtime("Running parameter optimization")

# define outdirs (this is done like this to be able to use fun.clean_perSVade_outdir at the end)
SVdetection_outdir = "%s/SVdetection_output"%opt.outdir; fun.make_folder(SVdetection_outdir)
outdir_gridss_final = "%s/final_gridss_running"%SVdetection_outdir; fun.make_folder(outdir_gridss_final)

# clean reference
fun.clean_reference_genome_windows_files(opt.ref)

# calculate the insert size
median_insert_size, median_insert_size_sd  = fun.get_insert_size_distribution(sorted_bam, replace=opt.replace, threads=opt.threads)
fun.print_with_runtime("The median insert size is %i, with an absolute deviation of %i"%(median_insert_size, median_insert_size_sd))

# debug very high insert sizes
if median_insert_size>opt.max_median_insert_size: raise ValueError("The median insert size is very high (%s). Maybe there was an error with the bam file generation."%(median_insert_size))

# get the optimised parameters
parameter_optimisation_dir = "%s/parameter_optimisation"%SVdetection_outdir; fun.make_folder(parameter_optimisation_dir)

gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict, max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup, df_cross_benchmark_best = fun.get_best_parameters_for_GridssClove_run(sorted_bam, opt.ref, parameter_optimisation_dir, threads=opt.threads, replace=opt.replace, n_simulated_genomes=opt.nsimulations, mitochondrial_chromosome=opt.mitochondrial_chromosome, simulation_ploidies=simulation_ploidies, range_filtering_benchmark=opt.range_filtering_benchmark, nvars=opt.nvars, real_bedpe_breakpoints=real_bedpe_breakpoints, median_insert_size=median_insert_size, median_insert_size_sd=median_insert_size_sd, tmpdir=opt.tmpdir, simulation_chromosomes=opt.simulation_chromosomes)

# write these optimal parameters
json_file = "%s/optimized_parameters.json"%opt.outdir
fun.write_gridss_parameters_as_json(gridss_blacklisted_regions, gridss_maxcoverage, gridss_filters_dict, max_rel_coverage_to_consider_del, min_rel_coverage_to_consider_dup, json_file, replace=opt.replace)

#################################
#################################
#################################


##################################
########## CLEAN OUTPUT ##########
##################################

# log
fun.print_with_runtime("cleaning")

# clean normal perSVade outdir
fun.clean_perSVade_outdir(opt.outdir)

# remove extra files from outdir
for f in os.listdir(opt.outdir): 
    file_path = "%s/%s"%(opt.outdir, f)
    if os.path.isfile(file_path) and f not in {fun.get_file(json_file), "all_cmds.txt"}: fun.remove_file(file_path)

# move to outir from parameter_optimisation_dir
for f in os.listdir(parameter_optimisation_dir): 
    file_path = "%s/%s"%(parameter_optimisation_dir, f)
    os.rename(file_path, "%s/%s"%(opt.outdir, f))

# remove specific dirs
for f in [reference_genome_dir, SVdetection_outdir]: fun.delete_folder(f)

# move the parameter accuracy to the outdir
df_accuracy = fun.get_tab_as_df_or_empty_df("%s/benchmarking_all_filters_for_all_genomes_and_ploidies/df_cross_benchmark_best.tab"%(opt.outdir))
df_accuracy = df_accuracy[["test_simName", "test_ploidy", "test_svtype", "precision", "recall", "Fvalue", "nevents", "TP", "FP", "FN"]].rename(columns={"test_simName":"simulation_genome", "test_ploidy":"simulation_ploidy", "test_svtype":"SV_type"})
fun.save_df_as_tab(df_accuracy, "%s/optimized_parameters_accuracy.tab"%opt.outdir)

# wite final file
fun.generate_final_file_report_one_module(final_file, start_time, time.time())

# print the message
fun.print_with_runtime("perSVade %s finished correctly"%module_name)

##################################
##################################
##################################
