gotcha_beta_2

#!/usr/bin/env python

#####################################################################################
#                                                                                   #
#                    gotha catch'em all (envronimental DNA baits)                   #
#                                                                                   #
#####################################################################################

#notes Kevin;
        # Taxon filtering in combination with clustering sequences
# bug notes
        # python3 ../gotcha_git/gotcha/gotcha.3.0.1.py -bp ./bold-cli -t "Apis,Eucera,Ceratina" -m COI-5P -c 5 -o Bees_beta_3.1 -ft species -e -v -cl 0.99
        # size selection bug

# new
        # tree output with branch length.
        # Fixed gap problem with selecting baits (no bait should have a gap)
        # Fixed problem with node selection. There is a problem that sequence close to the root are placed on tips.
        # Fixed, for each tip sequence look back into the tree to the closest ancestral node.

import re
import os
import re
import sys
import glob
import shutil
import os.path
import datetime
import argparse
import operator
import subprocess
from os import path
import pandas as pd
from tqdm import tqdm
from Bio import SeqIO
from Bio import AlignIO
from Bio.Seq import Seq
from Bio.SeqIO import FastaIO
from Bio.SeqRecord import SeqRecord
from Bio.Align.Applications import MafftCommandline
import ete3 #new v2.2.2
from ete3 import Tree  #new v2.2.2
from statistics import mean #new v2.2.2
import statistics #new v2.3.2
from collections import Counter #new v2.2.2
from numpy import random # new v3
import matplotlib.pyplot as plt # new v3
import numpy as np
##################################################################################################### parsing arguments and setting variables

class Range(object):
    def __init__(self, start, end):
        self.start = start
        self.end = end
    def __eq__(self, other):
        return self.start <= other <= self.end

def check_positive(value):
    try:
        value = int(value)
        if value <= 0:
            raise argparse.ArgumentTypeError("{} is not a positive integer".format(value))
    except ValueError:
        raise Exception("{} is not an integer".format(value))
    return value

parser = argparse.ArgumentParser(prog='gotcha', description='gotcha version 1.0.0 built Feb 2023\ndeveloped by Kevin Nota & Giobbe Forni\n',add_help=False)

parser.formatter_class = lambda prog: argparse.RawTextHelpFormatter(prog, max_help_position=150)
formatter_class=argparse.RawTextHelpFormatter

man = parser.add_argument_group('MARKER')
man.add_argument('-m', '--marker', metavar='\b', choices=['COI-5P', 'COI-3P', 'rbcL', 'matK', 'ITS', 'coding' , 'noncoding'], help='can be COI/rbcL/matK/ITS or coding/noncoding using custom .fna/.aln')
man.add_argument('-c', '--code', metavar='\b', help='genetic code - e.g. 1 for plastid and nuclear, 5 for mitochondrial invertebrate')

bld = parser.add_argument_group('BOLD DOWNLOAD')
bld.add_argument('-t', '--taxa', metavar='\b', help='input taxa')
bld.add_argument('-g', '--geo', metavar='\b', default="", help='geographic location of samples')
bld.add_argument('-ft', '--filter_tax', metavar='\b', default="", choices=['species', 'genus', 'family'], help='can be species/genus/family - defeault is none')
bld.add_argument("-bp", "--boldcli_path", metavar='\b', default="bold-cli", help="path to bold-cli")

ctm = parser.add_argument_group('CUSTOM INPUTS')
ctm.add_argument('-cf', '--custom_fna', metavar='\b', default="", help='custom nucleotide sequences file (.fasta format)')
ctm.add_argument('-ca', '--custom_aln', metavar='\b', default="", help='custom nucleotide alignment file (.fasta format) - gotcha will jump to size selection')
ctm.add_argument('-cn', '--custom_nwk', metavar='\b', default="", help='custom newick file - tree can be multifurcating and can lack taxa')

par = parser.add_argument_group('BAITS SEARCH PARAMETERS')
par.add_argument('-fl', '--filter_len', metavar='\b', default=100, help='minimum amminoacids / nucleotides length - default is 100', type=check_positive)
par.add_argument('-sc', '--seq_collapse', metavar='\b', type=float, default=1, help='percent identity to collapse sequences - default is 1')
par.add_argument('-bc', '--bts_collapse', metavar='\b', type=float, default=1, help='percent identity to collapse baits - defeault is 1')
par.add_argument('-bl', "--baitlength", metavar='\b', type=int, default=80, help='lenght of the bait sequences - default is 80')
par.add_argument('-tl', "--tiling", metavar='\b', default=30, type=int, help='kmer tiling for ancestral node selection - default is 30') #changed the default from 10 to 30 (think most cases ~3x coverage is alright and for a 80mer, that would be ~30)
par.add_argument('-ds', "--distance", metavar='\b', type=float, default=0.09, help='maximal distance to the ancestral node - default is 0.09')
par.add_argument('-so', "--trim_seqoverlap", metavar='\b', default="", help='trimal param between 0 and 100 - if not specified will skip')
par.add_argument('-ro', "--trim_resoverlap", metavar='\b', default="", help='trimal param between 1 and 0 - if not specified will skip')
par.add_argument('-dt', "--dust_threshold", metavar='\b', default=10, help='score threshold for low complexity region masking - default is 10')
par.add_argument('-mf', "--max_farg_lenght", metavar='\b', default=3000, type=int, help='this value can restrain the selected fragment size by choosing the optimal fragment below this value, \
                                                                                                                                                                                \ndefault will take the highest value if number of sequences multiplied by fragment size') #new in version 2.3.2
par.add_argument('-GC', "--GC_threshold", metavar='\b', default="0.30-0.70", help='GC range for bait filtering - default is 0.30-0.70')
par.add_argument('-ks', "--keep_short_sequences",  help='keep short sequences in size selection step', action='store_true')
par.add_argument('-NS', "--no_size_selection",  help='preforms size selection', action='store_false', default=True)
par.add_argument('-fs', '--fast', action='store_true', help='performs the fast baits search w/out ancestral state reconstruction')
par.add_argument('-mc', '--manual_check', action='store_true', help='stops prior to tree inference for a manual check of the alignment')

fal = parser.add_argument_group('VALIDATION')
fal.add_argument('-VR', "--run_validation", help='Run bait validation with simulated data', action='store_true')
fal.add_argument('-VO', "--run_validation_only",  help='Run only the bait validation with simulated data', action='store_true')
fal.add_argument('-ns', "--n_simulated",metavar='\b',  default=10000, type=int, help='number of reads to simulate')
fal.add_argument('-dp', "--distribution_par",metavar='\b',  default="50:50:35",  help='distribution for reads simulation min_read_length:mean:sd of normal distribution')
fal.add_argument('-mo', "--min_overlap", metavar='\b', default=20, type=int, help='min overlap read to bait to be considered in calculation')

otr = parser.add_argument_group('OTHER')
otr.add_argument("-h", "--help", action="help", help="show this help message and exit")
otr.add_argument("-path", default=".", help=argparse.SUPPRESS)
otr.add_argument('-o', '--out', metavar='\b', default="probes", help='basename of output file and folders - default is probes')
otr.add_argument('-v', '--verbose', action='store_false', help='keeps temporary folder and files')
otr.add_argument('-e', '--erase', action='store_true', help='erases and rewrites a pre existing output folder')
otr.add_argument('-th', '--threads', metavar='\b', type=check_positive, default=1, help='number of threads used for tree inference - default is 1')
otr.add_argument('--version', action='store_true', help="show program's version number")

args=parser.parse_args()

# fix output path so that it can handle "../"
args.out = re.sub("/$","", args.out)

def check_if_all_is_installed():

    if args.version == 1:
        print('\ngotcha version 5.0.0 built March 2023\ndeveloped by Kevin Nota & Giobbe Forni\n')
        exit()

    if args.boldcli_path != "bold-cli" :
        args.boldcli_path = args.boldcli_path
    rc = subprocess.call(['which', args.boldcli_path],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0 and not args.custom_fna:
        print('bold-cli\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'baseml'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('baseml\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'transeq'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('transeq\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'cd-hit'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('cdhit\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'iqtree'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('iqtree\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'trimal'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('trimal\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'translatorx'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('translatorx\tmissing in path')
        exit()

    rc = subprocess.call(['which', 'muscle'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)
    if rc != 0:
        print('muscle\tmissing in path')
        exit()


    ##################################################################################################### errors !!!

    if not args.custom_fna and not args.custom_aln and (args.taxa is None):
        print("\n WARNING! Either a custom alignment (--custom_fna) or a lineage to download from BOLD (--taxa) is required! \n")
        quit()

    if  args.custom_fna and args.custom_aln:
        print("\n WARNING! both custom sequences and alignment have been specified! \n")
        quit()

    if not args.custom_fna and (args.marker is None):
        print("\n WARNING! When not using custom alignment --marker flags is required! \n")
        quit()

    if  args.custom_nwk and (args.custom_fna is None):
        print("\n WARNING! When using a custom tree (--custom_nwk) either custom sequences (--custom_fna) or alignment (--custom_aln) has to be specified! \n")
        quit()

    if args.erase == False and path.exists(args.out):
        print("\n WARNING! An output folder with the same name already exists! Use --erase to overwrite \n")
        quit()
    elif args.erase == True and path.exists(args.out):
        shutil.rmtree(args.out)
        os.makedirs(args.out)
    else:
        os.makedirs(args.out)

    if (args.code != None):
        if int(args.code) not in range (1, 33):
            print("\n WARNING! Unknown geneitc code has been specified - please refer to NCBI! \n")
            quit()

    if not 0.01 <= args.distance <= 0.2:
        print("\n WARNING! A value between 0.01 and 0.2 should be specified for parameter distance! \n")
        quit()

    if not 1 >= args.seq_collapse >= 0.7:
        print("\n WARNING! A value between 1 and 0.7 should be specified for parameter seq_collapse! \n")
        quit()    
    
    if not 1 >= args.bts_collapse >= 0.7:
        print("\n WARNING! A value between 1 and 0.7 should be specified for parameter bts_collapse! \n")
        quit()  
    
    gc_thresholds=args.GC_threshold.split("-")
    if float(min(gc_thresholds))<=0 and float(min(gc_thresholds))>=1 :
        print("\n WARNING! GC threshold should be a float between 0-1")

    if args.trim_seqoverlap != "" and args.trim_resoverlap == "":
        print("\n WARNING! You have to specify both seqoverlap and resoverlap parameters for the trimal step \n")
        quit()

    if args.trim_seqoverlap == "" and args.trim_resoverlap != "":
        print("\n WARNING! You have to specify both seqoverlap and resoverlap parameters for the trimal step \n")
        quit()

    if args.trim_seqoverlap != "":
        if not 100 >= float(args.trim_seqoverlap) >= 1:
            print("\n WARNING! A value between 100 and 1 should be specified for parameter seqoverlap! \n")
            quit()

    if args.trim_resoverlap != "":
        if not 1 >= float(args.trim_resoverlap) >= 0:
            print("\n WARNING! A value between 1 and 0 should be specified for parameter resoverlap! \n")
            quit()

    os.makedirs(args.out + "/tmp")
    os.chdir(args.out + "/tmp")

    ##################################################################################################### define markers
    marker_list=[]
    if ( args.marker == "coding" ):
        coding = True
    if ( args.marker == "non-coding" ):
        coding = False

    if ( args.marker == "COI-5P" ):
        marker_list=[ "COI-5P" ]
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    if ( args.marker == "COI-3P" ):
        marker_list=[ "COI-3P" ]
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    elif ( args.marker == "rbcL" ):
        marker_list=[ "rbcL", "rbcl" , "RBCL" , "Rbcl" , "rbcla" , "rbcl-a" ]
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    elif ( args.marker == "matK" ):
        marker_list=[ "matK" , "matk" , "MATK" ]
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    elif ( args.marker == "ITS" ):
        marker_list=[ "ITS" , "its" ]
        coding = False
        if ( args.code != None ):
            print("\n WARNING! a genetic code has been specified when using a non-coding marker! \n")
            quit()

    if ( args.custom_fna != "" ) and ( args.marker == "coding" ):
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    elif ( args.custom_fna != "" ) and ( args.marker == "noncoding" ):
        coding = False
        if ( args.code != None ):
            print("\n WARNING! a genetic code has been specified for a non-coding marker! \n")
            quit()

    if ( args.custom_aln != "" ) and ( args.marker == "coding" ):
        coding = True
        if ( args.code == None ):
            print("\n WARNING! a genetic code has to be specified when using a coding marker! \n")
            quit()
    elif ( args.custom_aln != "" ) and ( args.marker == "noncoding" ):
        coding = False
        if ( args.code != None ):
            print("\n WARNING! a genetic code has been specified for a non-coding marker! \n")
            quit()
    print("\nanalysis started on" , datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S"), "\n")
    
    return coding, marker_list

def write_log_file_01():
    gc_thresholds=args.GC_threshold.split("-")
    with open("log.txt", "a+") as log:

        log.write("# analysis started on " + datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + "\n")

    # INPUTS

        log.write("\n\n\t INPUTS:" + "\n\n")

        if ( args.custom_fna != "" ):
            log.write("\t BOLD TAXONOMY \t\tN" + "\n")
            log.write("\t CUSTOM SEQUENCES \tY" + "\n")
            log.write("\t CUSTOM ALIGNMENT \tN" + "\n")
        elif ( args.custom_aln != "" ):
            log.write("\t BOLD TAXONOMY \t\tN" + "\n")
            log.write("\t CUSTOM SEQUENCES \tN" + "\n")
            log.write("\t CUSTOM ALIGNMENT \tY" + "\n")
        else:
            if ( args.geo == "" ):
                log.write("\t BOLD TAXONOMY \t\t" + str(args.taxa) + "\n")
                log.write("\t GEOGRAPHY \t\tN" + "\n")
                log.write("\t CUSTOM SEQUENCES \tN" + "\n")
                log.write("\t CUSTOM ALIGNMENT \tN" + "\n")
            else:
                log.write("\t BOLD TAXONOMY \t\t" + str(args.taxa) + "\n")
                log.write("\t GEOGRAPHY \t\t" + args.geo + "\n")
                log.write("\t CUSTOM SEQUENCES \tN" + "\n")
                log.write("\t CUSTOM ALIGNMENT \tN" + "\n")

        if ( args.custom_nwk == "" ):
            log.write("\t CUSTOM TREE \t\tN" + "\n")
        else:
            log.write("\t CUSTOM TREE \t\tY" + "\n")

        log.write("\t MARKER \t\t" + str(args.marker) + "\n")

        if (coding == True):
            log.write("\t GEN CODE \t\t" + str(args.code) + "\n")
        else:
            log.write("\t GEN CODE \t\t" + "-" + "\n")

        if (args.fast == False):
            log.write("\t FAST MODE \t\t" + "N" + "\n")
        else:
            log.write("\t FAST MODE \t\t" + "Y" + "\n")

    # PARAMETER

        log.write("\n\n\t PARAMETERS:" + "\n\n")
        if ( args.filter_tax == "" ):
            log.write("\t TAXONOMIC FILTER \tN" + "\n")
        else:
            log.write("\t TAXONOMIC FILTER \t" + str(args.filter_tax) + "\n")
        log.write("\t MIN MARKER LENGTH \t" + str(args.filter_len) + "\n")
        log.write("\t MAX MARKER IDENTITY \t" + str(args.seq_collapse) + "\n")
        log.write("\t BAITS LENGTH \t\t" + str(args.baitlength) + "\n")
        log.write("\t BAITS TILING \t\t" + str(args.tiling) + "\n")
        log.write("\t GC THRESHOLD \t\t" + min(gc_thresholds) + "-" + max(gc_thresholds) + "\n")
        if ( args.trim_seqoverlap != "" ) and ( args.trim_resoverlap != "" ):
            log.write("\t SEQ OVERLAP \t\t" + str(args.trim_seqoverlap) + "\n")
            log.write("\t RES OVERLAP \t\t" + str(args.trim_resoverlap) + "\n")
        else:
            log.write("\t TRIMAL \t\tN" + "\n")
        log.write("\t DUST TRESHOLD \t\t" + str(args.dust_threshold) + "\n")

def bold_download():

    if args.boldcli_path != "bold-cli" :
        args.boldcli_path = "../../" + args.boldcli_path

    if ( args.geo == "" ):
        print("downloading" , args.marker , "sequences for taxa" , args.taxa)
        subprocess.run([args.boldcli_path , "-taxon" , args.taxa , "-marker" , args.marker , "-output" , "tmp.bold"] , stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)
        if os.path.getsize("tmp.bold") == 0 :
            print("\n WARNING! The BOLD search for marker" , args.marker , "and taxa" , args.taxa , "returned nothing! \n")
            quit()
    else:
        print("downloading" , args.marker , "sequences for taxa" , args.taxa, "from" , args.geo)
        subprocess.run([args.boldcli_path , "-taxon" , args.taxa , "-marker" , args.marker , "-geo" , args.geo, "-output" , "tmp.bold"] , stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)
        if os.path.getsize("tmp.bold") == 0 :
            print("\n WARNING! The BOLD search for marker" , args.marker , "and taxa" , args.taxa , "in" , args.geo , "returned nothing! \n")
            quit()

    num0 = sum(1 for line in open('tmp.bold', errors='ignore'))
    print("#        downloaded: ", num0)

    with open("log.txt", "a+") as log:
        log.write("\n\n\t FILTERING:" + "\n\n")
        log.write("\t " + str(num0) + "\t starting sequences" + "\n")
    log.close()

    with open("tmp.bold", 'r', errors='ignore') as file:
        if 'Fatal Error' in file.read():
            print("\n WARNING! The taxonomy specified appears to be not correct! \n")
            quit()

##################################################################################################### filter marker and taxonomy

    tmp_1_fasta=[]

    if (args.filter_tax != ""):
        print("filtering sequences with" , args.filter_tax , "identification")
    else:
        print("skipping taxonomic filtering")

    with open("tmp.bold", 'r', errors='ignore') as file:
        header = file.readline()
        for l in file :
            sl = l.split('\t')
            if (args.filter_tax == "family"):
                if sl[15]:
                    header = (">" + sl[0])
                    seq = sl[71].replace('-','')
                    if sl[69] in marker_list:
                        tmp_1_line=[header, seq]
                        tmp_1_fasta.append(tmp_1_line)
            elif (args.filter_tax == "genus"):
                if sl[19]:
                    header = (">" + sl[0])
                    seq = sl[71].replace('-','')
                    if sl[69] in marker_list:
                        tmp_1_line=[header, seq]
                        tmp_1_fasta.append(tmp_1_line)
            elif (args.filter_tax == "species"):
                if sl[21]:
                    header = (">" + sl[0])
                    seq = sl[71].replace('-','')
                    if sl[69] in marker_list:
                        tmp_1_line=[header, seq]
                        tmp_1_fasta.append(tmp_1_line)
            else:
                header = (">" + sl[0])
                seq = sl[71].replace('-','')
                if sl[69] in marker_list:
                    tmp_1_line=[header, seq]
                    tmp_1_fasta.append(tmp_1_line)

    with open('tmp1.fna', 'w') as tmp_1:
        for line in tmp_1_fasta:
            for element in line:
                tmp_1.write(str(element) + '\n')

    if os.path.getsize("tmp1.fna") == 0 :
        print("\n WARNING! No sequence passed the taxonomic filter! \n")

    num1 = len([1 for line in open("tmp1.fna") if line.startswith(">")])                                                    # count tmp1.fna
    print("#        kept: ", num1)
    
    with open("log.txt", "a+") as log:
        log.write("\t " + str(num1) + "\t after taxonomic filtering" + "\n")
    log.close()

def seq_quality_filters():
    if (args.custom_fna != ""):
        print("custom sequences detected\n")
        custom_fna =  "../../" + args.custom_fna
        shutil.copy(custom_fna,"tmp1.fna")

    if (args.custom_aln == ""):
        print("collapsing identical sequences with a percent identitiy of", args.seq_collapse)
        subprocess.run(["cd-hit","-i","tmp1.fna","-o","tmp2a.fna","-c", str(args.seq_collapse)] , stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)

        num2 = len([1 for line in open("tmp2a.fna") if line.startswith(">")])                                                   # count tmp2.fna
        print("#        kept: ", num2)
        with open("log.txt", "a+") as log:
            log.write("\t " + str(num2) + "\t after identity collapse" + "\n")
        log.close()
    ##################################################################################################### remove sequences with non ATGC nucleotides

        print("removing sequences with non ATGC nucleotides")

        selected_seqs_ns = list()

        tmp2b_fna=open("tmp2b.fna",'w')

        for record in SeqIO.parse("tmp2a.fna", "fasta"):
    #               if record.seq.count('N') == 0:
            matches = ['B','D','E','F','H','I','J','K','L','M','N','O','P','Q','R','S','U','V','W','X','Y','Z','b','d','e','f','h','i','j','k','l','m','n','o','p','q','r','s','u','v','w','x','y','z']
            if not any(x in record.seq.upper() for x in matches):
                selected_seqs_ns.append(record)

        nseq=(len(selected_seqs_ns))
        SeqIO.write(selected_seqs_ns, tmp2b_fna , "fasta")
        tmp2b_fna.close()
        print("#        kept: ", nseq)

        with open("log.txt", "a+") as log:
            log.write("\t " + str(nseq) + "\t removing sequences with non ATGC nucleotides" + "\n")
        log.close()
    ##################################################################################################### length cutoff

        print("filtering sequences with a minimum length of", args.filter_len)

        selected_seqs_nt = list()
        tmp2_fna=open("tmp2.fna",'w')

        for record in SeqIO.parse("tmp2b.fna", "fasta"):
            if len(record.seq) >= int(args.filter_len):
                selected_seqs_nt.append(record)
        nseq=(len(selected_seqs_nt))
        SeqIO.write(selected_seqs_nt, tmp2_fna , "fasta")

        print("#        kept: ", nseq)
        tmp2_fna.close()

        with open("log.txt", "a+") as log:
            log.write("\t " + str(nseq) + "\t after length filter" + "\n")
        log.close()
    ##################################################################################################### too few sequences to proceed?

        nseq = len([1 for line in open("tmp2.fna") if line.startswith(">")])
        if nseq < 4:
            print("\n WARNING! Less than 4 sequences passed the filtering steps! \n" )
            quit()

    ##################################################################################################### if the marker is coding

    if (coding == True and args.custom_aln == ""):
        #print("aligning sequences")

        # adjust direction with mafft
        with open('tmp3.fna', 'w') as tmp3_fna, open('tmp2.fna', 'r') as tmp2_fna:
            if 'MAFFT_BINARIES' in os.environ:
                os.environ.pop('MAFFT_BINARIES')
            subprocess.call(["mafft" , "--thread", str(args.threads), "--adjustdirection" , "--reorder" , "tmp2.fna"], stdout=tmp3_fna, stderr=subprocess.DEVNULL)

        # remove "N" and "-"
        with open("tmp4.fna", "w") as tmp4_fna:
            for record in SeqIO.parse("tmp3.fna", "fasta"):
                record.seq = record.seq.replace("-", "")
                record.seq = record.seq.replace("n", "")
                SeqIO.write(record, tmp4_fna, "fasta")

        # protein-guided nucleotide alignment with translatorx
        with open('tmp4.fna', 'r') as tmp4_fna, open('translatorx.log', 'w') as translatorx_log:
            subprocess.call(["translatorx" , "-i" , "tmp4.fna", "-t", "T", "-o", "tmp", "-c", args.code], stdout=translatorx_log, stderr=translatorx_log)
        os.rename("tmp.nt_ali.fasta", "tmp5.fna")

        # translate sequences to find stop codons
        with open('tmp5.fna', 'r') as tmp5_fna:
            subprocess.call(["transeq" , "-sequence" , "tmp5.fna", "-outseq", "translated4stop.fna", "-table", args.code], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

        # find sequences with stopcodons
        identifiers = set()
        with open('translated4stop.fna', 'r') as translated4stop_fna:
            for seq_record in SeqIO.parse(translated4stop_fna, "fasta"):
                if "*" in (seq_record.seq[0:(len(seq_record.seq)-2)]):
                    identifiers.add(seq_record.id[:-2])

        print("removing sequences with stop codons")

        # find sequences with stopcodons
        with open('tmp5.fna', 'r') as tmp5_fna, open('tmp6.fna', 'w') as tmp6_fna:
            records = SeqIO.parse(tmp5_fna, 'fasta')
            for record in records:
                if record.id not in identifiers:
                    SeqIO.write(record, tmp6_fna, 'fasta')

        num6 = len([1 for line in open("tmp6.fna") if line.startswith(">")])                                                    # count tmp6.fna
        print("#        kept: ", num6)

        with open("log.txt", "a+") as log:
            log.write("\t " + str(num6) + "\t the removal of sequences with stop codons" + "\n")
        log.close()

        os.rename("tmp6.fna", "tmp.aln")

    ##################################################################################################### if the marker is not coding

    if (coding == False and args.custom_aln == ""):

        print("aligning non-coding sequences")

        with open('tmp.aln', 'w') as tmp_aln, open('tmp2.fna', 'r') as tmp2_fna:
            if 'MAFFT_BINARIES' in os.environ:
                os.environ.pop('MAFFT_BINARIES')
            subprocess.call(["mafft" , "--thread", str(args.threads) , "--adjustdirection" , "tmp2.fna"], stdout=tmp_aln, stderr=subprocess.DEVNULL)

    ##################################################################################################### reformat after mafft adjustdirection

    os.rename('tmp.aln', 'tmptmp.aln')

    tmp_aln=open("tmp.aln",'w')

    for record in SeqIO.parse("tmptmp.aln", "fasta"):
        record.description = (record.id).replace('_R_', '')
        record.id=record.description
        SeqIO.write(record, tmp_aln , "fasta")

    tmp_aln.close()

    os.remove('tmptmp.aln')

    ##################################################################################################### can enter custom .aln

    if (args.custom_aln != ""):
        print("custom alignment detected")
        custom_aln =  "../../" + args.custom_aln
        shutil.copy(custom_aln,"tmp.aln")
        num = len([1 for line in open("tmp.aln") if line.startswith(">")])

def size_selection():

    input_sequences=list(SeqIO.parse("tmp.aln", "fasta")) # path to the fasta file
    fragment_selections= pd.DataFrame()
    #counting the number of gaps in the alignment file
    gaps = pd.DataFrame()

    threshold_list=[]
    number_of_sequences_list=[]
    seq_length_list=[]
    number_X_length_list=[]
    start_position_list=[]
    end_position_list=[]
    position_list=[]
    proportion_gaps=[]
    print("\nselecting optimal fragment size based on", len(input_sequences),"sequences :")# "sequences with max length of",len(input_sequences[0].seq), "\n")

    #print("\nnumber of input sequences =", len(input_sequences))
    #print("max sequences lenght =", len(input_sequences[0].seq))

    #print("calculating missing data")
    #progress=tqdm(total=None, desc="Alignment scanner ")
    for p in range(len(input_sequences[0].seq)):
        gap_positions = 0
        for i in range(len(input_sequences)):
            if (input_sequences[i].seq[p] == "-"):
                gap_positions += 1
            #data = pd.Series({'position':p, 'proportion gaps':gap_positions/len(input_sequences)})

        position_list.append(p)
        proportion_gaps.append(gap_positions/len(input_sequences))
        #gaps=gaps.append(pd.DataFrame(data), ignore_index=True)

        #gaps= pd.concat([gaps, data.to_frame().T], ignore_index=True)
        #progress.update()
    #progress.close()
    gaps=pd.DataFrame({'position':position_list, 'proportion gaps':proportion_gaps})
    #progress=tqdm(total=None, desc="Calculating distance = ")
    for treshold in range(0, 100):
        threshold_found=False
        i=0

        # find the first and last position where thresholds values are met
        # tressholds that are tested are 0-1 with 0.01 increments (101 loops)
        
        while not (threshold_found==True) :
            if (gaps.loc[i,"proportion gaps"]<=float(treshold/100)):
                threshold_found=True
                start_position=int(gaps.loc[i,"position"])
            i=i+1
            if(i==len(gaps)):
                threshold_found=True
                start_position=0

        threshold_found=False
        i=len(gaps)-1
        while not (threshold_found==True) :
            if (gaps.loc[i,"proportion gaps"]<=float(treshold/100)):
                threshold_found=True
                finish_position=int(gaps.loc[i,"position"])
            i=i-1
            if(i==0):
                threshold_found=True
                finish_position=0
        count=0

        lenght_no_gaps=[]
        # this counts the number of sequences within the threshold with no gaps (-)
        for i in range(len(input_sequences)):
            if(len(re.findall("^-|-$|^N|N$", str(input_sequences[i].seq[start_position:finish_position])))==0):
                count+=1
                lenght_no_gaps.append(len(re.sub("-", "", str(input_sequences[i].seq[start_position:finish_position]))))
                
        if len(lenght_no_gaps) == 0 :
            continue
        
        threshold_list.append(str(treshold)+"%")
        number_of_sequences_list.append(count)
        seq_length_list.append(int(re.sub(",.*", "", re.sub("\(", "", str(Counter(lenght_no_gaps).most_common()[0])))))
        number_X_length_list.append(int(re.sub(",.*", "", re.sub("\(", "", str(Counter(lenght_no_gaps).most_common()[0]))))*count)
        start_position_list.append(start_position)
        end_position_list.append(finish_position)
        
        #data = pd.Series({'threshold':str(treshold)+"%",
        #    'number of sequences':count,
        #    'seq_length':int(re.sub(",.*", "", re.sub("\(", "", str(Counter(lenght_no_gaps).most_common()[0])))), #statistics.mode(lenght_no_gaps), 
        #    'number * lenght':count*len(input_sequences[i].seq[start_position:finish_position]),
        #    'start position':start_position,
        #    'finish_position':finish_position})
        #fragment_selections = fragment_selections.append(pd.DataFrame(data), ignore_index=True)
        #fragment_selections = fragment_selections.append(pd.DataFrame(data), ignore_index=True)
        #fragment_selections = pd.concat([fragment_selections, data.to_frame().T], ignore_index=True)
        #progress.update()
    fragment_selections = pd.DataFrame({'threshold':threshold_list, 
                                        'number of sequences':number_of_sequences_list,
                                        'seq_length':seq_length_list,
                                        'number * lenght':number_X_length_list,
                                        'start position':start_position_list,
                                        'finish_position':end_position_list})

    fragment_selections.to_csv("size_selection_summary.df", sep="\t") # new version 2.3.2
    #progress.close()    

    # loop to find the most optimal fragment lenght
    find_optimal_fragment=False
    i=0

    fragment_selections = fragment_selections.loc[fragment_selections['seq_length'] <= args.max_farg_lenght] # new version 2.3.2 (gaps is a problem here)
    
    while not (find_optimal_fragment==True):
        if((fragment_selections.loc[i,"number * lenght"]==fragment_selections['number * lenght'].max())==True):
            #optimal_threshold=fragment_selections.loc[i, 'threshold']
            find_optimal_fragment=True
            start_position=fragment_selections.loc[i,'start position']
            finish_position=fragment_selections.loc[i,'finish_position']
            selected_size=fragment_selections.loc[i,'seq_length']
        seq_count=fragment_selections.loc[i, 'number of sequences']
        i=i+1

    print(selected_size)
    while not (((finish_position-start_position)/3.0).is_integer()):
        finish_position-=1
        lenght=((finish_position-start_position)/3.0)

    #correcting size for the triplet cutoff (prevent sequences being lost due to gaps)
    #count_eng_gap=0
    frag_count=0
    
    while not (seq_count<=(frag_count)):
        #count_eng_gap=0
        frag_count=0
        for i in range(len(input_sequences)):
            #if(len(re.findall("-$", (str(input_sequences[i].seq[start_position:finish_position]))))!=0):
            #    count_eng_gap+=1
            # ^-|-$|^N|N$
            if(len(re.findall("$-|$N", (str(input_sequences[i].seq[start_position:finish_position]))))==0):
                frag_count+=1
        if seq_count>=frag_count :
            finish_position=finish_position-3
    
    frag_count=0
    while not (seq_count<=(frag_count)):
        #count_eng_gap=0
        frag_count=0
        for i in range(len(input_sequences)):
            #if(len(re.findall("-$", (str(input_sequences[i].seq[start_position:finish_position]))))!=0):
            #    count_eng_gap+=1
            # ^-|-$|^N|N$
            if(len(re.findall("^-|^N", (str(input_sequences[i].seq[start_position:finish_position]))))==0):
                frag_count+=1
        if seq_count>=frag_count :
            start_position=start_position+3

    #Creating new fasta file with original labels but new sequence
    new_fasta=list()
    frag_count=0
    for i in range(len(input_sequences)):
        if args.keep_short_sequences == True : 
            new_seq=SeqRecord(Seq(str(input_sequences[i].seq[start_position:finish_position])),
            id=input_sequences[i].id,
            description="")
            new_fasta.append(new_seq)
            frag_count+=1
        if args.keep_short_sequences != True : 
            if(len(re.findall("^-|-$|^N|N$", str(input_sequences[i].seq[start_position:finish_position])))==0): #removes all sequences that start with w gap, and sequences that contain N's
                new_seq=SeqRecord(Seq(str(input_sequences[i].seq[start_position:finish_position])),
                id=input_sequences[i].id,
                description="")
                new_fasta.append(new_seq)
                frag_count+=1

    #Writing the new fasta to tmp folder
    SeqIO.write(new_fasta, "size_selected.fasta", "fasta") #path to the output file
    os.rename("tmp.aln", "no_size_selection.fasta")
    shutil.copy('size_selected.fasta', 'tmp.aln')
    #os.rename("size_selected.fasta", "tmp.aln")

    #print("optimal missing data threshold = "+str(optimal_threshold)) #removed v2.3.2
    print("#\tkept:\t\t"+str(frag_count))
    print("#\tomitted:\t "+str(len(input_sequences)-frag_count))
    print("#\tfragment size:\t "+str(selected_size)+"\n")
    #print("start position = "+str(start_position)) #removed v2.3.2
    #print("end position = "+str(finish_position)) #removed v2.3.2

    num7=str(frag_count)
    with open("log.txt", "a+") as log:
        log.write("\t " + str(num7) + "\t after size selection" + "\n")
    log.close()

def extra_trimal_fiter():
    print("Misaligned sequences removal step with seqoverlap", args.trim_seqoverlap, "and resoverlap", args.trim_resoverlap)
    with open('trimal.aln', 'w') as trimal_aln, open('tmp.aln', 'r') as tmp_aln, open('trimal.log', 'w') as trimal_log:
        subprocess.call(["trimal" , "-in" , "tmp.aln", "-out", "trimal.aln", "-noallgaps", "-seqoverlap",args.trim_seqoverlap, "-resoverlap", args.trim_resoverlap], stdout=trimal_log, stderr=trimal_log)
    os.rename('trimal.aln' , 'tmp.aln')

    if os.stat('tmp.aln').st_size == 0:
        print("\nThe misaligned sequence removal step eliminated all sequences - use different sequoverlap and resoverlap parameters!")
        exit()

    num8 = len([1 for line in open("tmp.aln") if line.startswith(">")])                                                     # count 8
    print("#        kept: ", num8)
        
    with open("log.txt", "a+") as log:
        log.write("\t " + str(num8) + "\t after the removal of misaligned sequences" + "\n")
    log.close()

def standard_trimal_filter():
    with open('trimal.aln', 'w') as trimal_aln, open('tmp.aln', 'r') as tmp_aln, open('trimal.log', 'w') as trimal_log:
        subprocess.call(["trimal" , "-in" , "tmp.aln", "-out", "trimal.aln", "-noallgaps"], stdout=trimal_log, stderr=trimal_log)
    os.rename('trimal.aln' , 'tmp.aln')

def tree_inverence():

    if (args.custom_nwk != ""):

        custom_tre = '../../' + args.custom_nwk
        shutil.copy(custom_tre, './custom.tre')
        custom_tre = './custom.tre'

        brlen = re.findall(":[0-9]*.[0-9]*", custom_tre)

        with open(custom_tre, 'r') as tre_file:

            if brlen != "":
                content = tre_file.read()
                tre_file = re.sub(":[0-9]*.[0-9]*", '', content )

            nwk_sp_list=tre_file.replace(',', ' ').replace('(', '').replace(')', '').replace(';', '').split()

        aln_sp_list = []
        for record in SeqIO.parse("tmp.aln", "fasta"):
            aln_sp_list.append(record.id)

        check = all(sp in aln_sp_list for sp in nwk_sp_list)

        if  check == False:
            print("custom nwk detected with more species in the tree than sequences - these tips will be pruned.")
            toprune_list = list(set(nwk_sp_list) - set(aln_sp_list))
            with open('toprune.lst', 'w') as f:
                for line in toprune_list:
                    f.write(line)
                    f.write('\n')

            tmp_rscript_pruning=[]
            tmp_rscript_pruning.append("library(phytools)")
            tmp_rscript_pruning.append("tree <- read.newick(\"custom.tre\")")
            tmp_rscript_pruning.append("spec <- read.table(\"toprune.lst\")")
            tmp_rscript_pruning.append("ptre<-drop.tip(tree,tree$tip.label[match(spec$V1, tree$tip.label)])")
            tmp_rscript_pruning.append("write.tree(ptre, file=\"custom.tre\")")

            with open('tmp_rscript_pruning.R', 'w') as tmp_rscript_pruning_file:
                for line in tmp_rscript_pruning:
                    tmp_rscript_pruning_file.write(line + '\n')

            with open('tmp_rscript_pruning.R_log', 'w') as tmp_rscript_pruning_log:
                subprocess.run(["Rscript" , "tmp_rscript_pruning.R"], stdout=tmp_rscript_pruning_log, stderr=tmp_rscript_pruning_log)

        elif len(nwk_sp_list) == len(aln_sp_list) and check == True:
            print("custom nwk detected with 1:1 species match to sequences - skipping topology inference.")
        elif len(nwk_sp_list) < len(aln_sp_list)  and check == True:
            print("custom nwk detected with fewer species than sequences - using constrained topology.")

##################################################################################################### if the custom .nwk is coding

        if (coding == True):

            tmp_partitions=[]
            filename = "tmp.aln"
            format = "fasta"
            tmp_aln = AlignIO.read("tmp.aln", "fasta")
            tmp_partitions.append("DNA, st = 1-%i\\3" % tmp_aln.get_alignment_length())
            tmp_partitions.append("DNA, nd = 2-%i\\3" % tmp_aln.get_alignment_length())
            tmp_partitions.append("DNA, rd = 3-%i\\3" % tmp_aln.get_alignment_length())

            with open('tmp.partitions', 'w') as tmp_partitions_file:
                for line in tmp_partitions:
                    tmp_partitions_file.write(line + '\n')

            subprocess.run(["iqtree" , "-s" , "tmp.aln" , "-nt" , str(args.threads) , "-spp" , "tmp.partitions" , "-g" , "custom.tre", "-fast"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

            def_tre_file = args.out + ".nwk"
            os.rename('tmp.partitions.treefile' , def_tre_file)
            def_aln_file = args.out + ".aln"
            os.rename('tmp.aln' , def_aln_file)

##################################################################################################### if the custom .nwk is not coding

        if (coding == False):

            with open('tmp.nwk', 'w') as tmp_nwk :
                subprocess.run(["iqtree" , "-s" , "tmp.aln" , "-nt" , str(args.threads) , "-g" , custom_tre, "-fast"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

            def_tre_file = args.out + ".nwk"
            os.rename('tmp.aln.treefile' , def_tre_file)

            def_aln_file = args.out + ".aln"
            os.rename('tmp.aln' , def_aln_file)

##################################################################################################### inferring a coding marker tree

    elif (args.custom_nwk == ""):

        print("\ninferring tree\n")

        if (coding == True):

            tmp_partitions=[]
            filename = "tmp.aln"
            format = "fasta"
            tmp_aln = AlignIO.read("tmp.aln", "fasta")
            tmp_partitions.append("DNA, st = 1-%i\\3" % tmp_aln.get_alignment_length())
            tmp_partitions.append("DNA, nd = 2-%i\\3" % tmp_aln.get_alignment_length())
            tmp_partitions.append("DNA, rd = 3-%i\\3" % tmp_aln.get_alignment_length())

            with open('tmp.partitions', 'w') as tmp_partitions_file:
                for line in tmp_partitions:
                    tmp_partitions_file.write(line + '\n')

            with open('tmp.nwk', 'w') as tmp_nwk :
                subprocess.run(["iqtree" , "-s" , "tmp.aln" , "-nt" , str(args.threads) , "-spp" , "tmp.partitions", "-fast"], stdout=tmp_nwk, stderr=subprocess.DEVNULL)

            def_tre_file = args.out + ".nwk"
            os.rename('tmp.partitions.treefile' , def_tre_file)

            def_aln_file = args.out + ".aln"
            os.rename('tmp.aln' , def_aln_file)

##################################################################################################### inferring a not coding marker tree

        if (coding == False):

            with open('tmp.nwk', 'w') as tmp_nwk :
                subprocess.run(["iqtree" , "-s" , "tmp.aln" , "-nt" , str(args.threads), "-fast"], stdout=tmp_nwk, stderr=subprocess.DEVNULL)

            def_tre_file = args.out + ".nwk"
            os.rename('tmp.aln.treefile' , def_tre_file)

            def_aln_file = args.out + ".aln"
            os.rename('tmp.aln' , def_aln_file)

##################################################################################################### ancestral sequences inference

    if (coding == True):
        print("inferring ancestral sequences using genetic code" , args.code , "\n")
    if (coding == False):
        print("inferring ancestral sequences\n")

    tmp_ctl=[]

    aln_line="seqfile = " + def_aln_file
    tmp_ctl.append(aln_line)
    tmp_ctl.append("outfile = tmp_baseml.out")
    tre_line="treefile = " + def_tre_file
    tmp_ctl.append(tre_line)
    tmp_ctl.append("noisy = 3")
    tmp_ctl.append("verbose = 1")
    tmp_ctl.append("runmode = 0")
    tmp_ctl.append("model = 7")
    tmp_ctl.append("Mgene = 0")
    tmp_ctl.append("clock = 0")
    tmp_ctl.append("fix_kappa = 0")
    tmp_ctl.append("kappa = 2.5")
    tmp_ctl.append("fix_alpha = 1")
    tmp_ctl.append("alpha = 0.")
    tmp_ctl.append("Malpha = 0")
    tmp_ctl.append("ncatG = 5")
    tmp_ctl.append("fix_rho = 1")
    tmp_ctl.append("rho = 0.")
    tmp_ctl.append("nparK = 0")
    tmp_ctl.append("nhomo = 0")
    tmp_ctl.append("getSE = 0")
    tmp_ctl.append("RateAncestor = 1")
    tmp_ctl.append("Small_Diff = 1e-6")
    tmp_ctl.append("cleandata = 0")
    if coding == True :
        icode = int(args.code) - 1
    #       icode_line="icode = " + str(icode)
        icode_line="icode = 11"
        tmp_ctl.append(icode_line)
    tmp_ctl.append("fix_blength = 2")
    tmp_ctl.append("method = 0")

    with open('tmp.ctl', 'w') as tmp_ctl_file:
        for line in tmp_ctl:
            tmp_ctl_file.write(line + '\n')

    with open('tmp.nwk', 'w') as tmp_nwk, open('tmp_baseml.log', 'w') as tmp_baseml_log:
        subprocess.run(["baseml" , "tmp.ctl"], stdout=tmp_baseml_log, stderr=subprocess.DEVNULL)

        tmp_baseml_log = './tmp_baseml.log'
        baseml_error = re.findall("stop codon TAA", tmp_baseml_log)
        if(len(baseml_error) >= 1):
            print("\n A VERY RARE EVENT HAPPENED! baseml crashed because it inferred an ancestral stop codon ... quite funny right? \n")
            quit()

##################################################################################################### extract baseml trees

    with open("rst", 'r') as rst, open("brlen.phylo.txt", 'w') as brlen_phylo, open("nodes.phylo.txt", 'w') as nodes_phylo:
        brlen = [10]
        nodes = [17]
        for position, line in enumerate(rst):
            if position in brlen:
                brlen_phylo.write(line.replace(" ", ""))
            if position in nodes:
                nodes_phylo.write(line.replace(" ", ""))

    return def_tre_file, def_aln_file

def bait_finder():
    gc_thresholds=args.GC_threshold.split("-")
    print("finding baits of", args.baitlength, "nt with max distance of", args.distance,"% and tiling",args.tiling)

    input_seq = list(SeqIO.parse(args.out + ".aln", "fasta"))
    max_n_miss = int(args.distance*args.baitlength)
    anc_seq=[]
    fasta_dict = {}
    tip_list=[]
    for seq in input_seq:
        fasta_dict[str(seq.name)]=str(seq.seq)
        tip_list.append(seq.id)
                        
    file = open("rst")
    patern=r" {3,}"
    node_seq_dict = {}
    c=0
    for row in file :
        if "List of extant and reconstructed sequences" in row:
            c=c+1
        if(c==2):
            if "node #" in row :

                node_name=re.split(patern, row)[0]
                node_seq=re.sub("\n", "", re.sub(" ", "", re.split(patern, row)[1]))
                node_seq_dict[str(node_name)] = str(node_seq)
                anc_seq.append(">"+node_name+"\n"+node_seq+"\n")
    node_seq_dict.update(fasta_dict)

    fasta_file = open("Anc_seqs.fasta", "w")
    fasta_file.writelines(anc_seq)
    fasta_file.close()
    
    t=Tree("nodes.phylo.txt",1)
    node_content = t.get_cached_content(store_attr='name')

    pos= 0  
    seq_len = len(seq)
    use_node="NO"
    bait_fasta=[]
    summary_stats=["Pos\tMean\tmin\tmax\tuniq\tgc\tmin_gc\tmax_gc\n"]
    progress=tqdm(total=None, desc="Tile number ")
    iii=0

    def get_seq(dict, key):
        seq_out=str()
        for value in ungap_pos:
            seq_out=seq_out+dict[key][value]
        return(seq_out)

    def find_ungaped_seq(start, key):
        gap_counter=0
        tip_seq=fasta_dict[key][start:(args.baitlength+start)]
        while (tip_seq.count("-") != 0) and (len(re.sub("-", "",tip_seq))!=args.baitlength):
            gap_counter = gap_counter+(args.baitlength-len(re.sub("-", "",tip_seq)))
            tip_seq=fasta_dict[key][start:(args.baitlength+start+gap_counter)]
            if (args.baitlength+start+gap_counter) >= seq_len:
                break
        return(tip_seq, gap_counter)  
    
    def update_progress_bar (bar_nam, lenght):
        if lenght > 200 :
            bar_nam.update()

    while (pos+args.baitlength) < seq_len:

        used_baits_dict={}
        
        tips_included = []
        tips_included_2 = []
        miss_list_sum_bait = []
        baitmer_list = []
        used_nodes_it=[]
        c_baits=0
        gap_counter=0
        for node in t.traverse():
            bv=0
            if not node.is_leaf():
                node_name ="node #"+str(node.name)
                level="NODE"
            else:
                node_name =re.sub("^\d*_", "", str(node.name))
                node_name_tip =str(node.name)
                level="TIP"

            anc_str= node_seq_dict[node_name][pos:(args.baitlength+pos)]
            clade=node_content[node]
            miss_list=[]
            gap_counter_all=[]
            for seq in list(clade):
                gap_counter=0
                tip_seq, gap_counter = find_ungaped_seq(pos, re.sub("^\d*_", "", seq))
            
                gap_counter_all.append(gap_counter)
            
            if(min(gap_counter_all))!=max(gap_counter_all):
                continue

            ungap_pos=[]
            for nt in range(pos, args.baitlength+pos+gap_counter):
                if args.baitlength+pos+gap_counter > seq_len:
                    continue
                gaps=[]
                
                for seq in list(clade):
                    gaps.append(fasta_dict[re.sub("^\d*_", "", seq)][nt])
                if gaps.count("-") == 0:
                    ungap_pos.append(nt)
                    if len(ungap_pos) != args.baitlength:
                        continue
            if len(ungap_pos) != args.baitlength:
                continue

            anc_str = get_seq(node_seq_dict, node_name)

            for seq in list(clade):
                tip_seq=get_seq(fasta_dict, re.sub("^\d*_", "", seq))
                miss=0
                for nt in range(len(tip_seq)):
                    if tip_seq[nt].upper()!= anc_str[nt].upper() and tip_seq[nt] != "-" :
                        miss=miss+1
                miss_list.append(miss)

            gc=((anc_str.upper().count("G")+anc_str.upper().count("C"))/args.baitlength)
            if(max(miss_list)<=max_n_miss) and len(Counter(gap_counter_all).keys())==1 and gc>=float(min(gc_thresholds)) and gc<=float(max(gc_thresholds)) : #gc
                if level == "NODE" :
                    clade_tips=node.get_leaf_names()
                elif level == "TIP":
                    continue
                    clade_tips=[node_name_tip]
                for tip in clade_tips:
                    if tip in tips_included :
                        continue
                    else:
                        tips_included.append(tip)
                        tips_included_2.append(re.sub("[0-9].*_","", tip))
                    if level == "NODE":
                        use_node="YES"

            if(use_node=="YES"):
                bait_fasta.append(">"+node_name+ " "+ str(ungap_pos[0])+":"+str(ungap_pos[args.baitlength-1])  +"\n" +anc_str+"\n")
                baitmer_list.append(anc_str)
                use_node="NO"
                miss_list_sum_bait=miss_list_sum_bait+miss_list
                c_baits=c_baits+1
                used_nodes_it.append(node_name)
                used_baits_dict.update({node_name+ " "+ str(ungap_pos[0])+":"+str(ungap_pos[args.baitlength-1]):anc_str})

        tip_no_node = []
        tip_no_node_name = []
        gap_cor_tip_dict={}
        for tip in tip_list:
            if tip in tips_included_2:
                continue
            #tip_no_node.append(get_seq(fasta_dict, tip))

            gap_cor_tip = re.sub("-", "", fasta_dict[tip][pos:pos+args.baitlength])
            
            end_of_seq=False
            gap_counter=0
            #print("seq  = " + gap_cor_tip)
            if(gap_cor_tip != args.baitlength):
                while len (gap_cor_tip) < (args.baitlength-1) :
                    gap_counter+=1
                    gap_cor_tip = re.sub("-", "", fasta_dict[tip][pos:(pos+args.baitlength+gap_counter)])
                    #print(gap_cor_tip)
                    #print(len (gap_cor_tip))
                    #print("\n"+str(seq_len))
                    #print(pos+args.baitlength+gap_counter)
                    if pos+args.baitlength+gap_counter >= seq_len:
                        end_of_seq=True
                        break
                    #print("\ntip seq : " + gap_cor_tip)
                    #print("seq len " + str(len(gap_cor_tip)))

            if gap_cor_tip.upper() in tip_no_node:
                continue
            if end_of_seq == True :
                continue
            
            ungap_pos=[]
            for nt in range(pos, (pos+args.baitlength+gap_counter)):
                if fasta_dict[tip][nt] != "-" :
                    ungap_pos.append(nt)
            
            #print(ungap_pos)
            if len(ungap_pos) != args.baitlength:# hack; not sure why some areas behave like this, looks like it has somehting to do with ending on a gap
                ungap_pos.append(pos+args.baitlength+gap_counter+1)

            #print(len(ungap_pos))
            #print("\n")
            #print(ungap_pos)
            #print(len(ungap_pos[0:(args.baitlength)-1]))
            #print("size = "+str(len(ungap_pos)))

            tip_no_node.append(gap_cor_tip.upper())
            tip_no_node_name.append(str(tip))
            gap_cor_tip_dict.update({tip:ungap_pos})
            # here the tips sequences without node need to become ungapped based on the starting position.
        #print(gap_cor_tip_dict)
        progress_int_2=[]
        if len(tip_no_node_name) > 200 :
            progress_int_2=tqdm(total=len(tip_no_node_name)+1, desc=str("Finding nearest ancestral node ; pos = " + str(pos)))
        
        # loop through the tips with no node yet to find a good pairing among the nodes already used.
        if len(tip_no_node_name)!=0:
            I=0
            k=0
            for seq in tip_no_node:
                tip_name=tip_no_node_name[I]
                I+=1
                miss_list=[]
                for node in used_baits_dict.keys():
                    miss=0
                    if len(seq)!=args.baitlength:
                        anc_str=used_baits_dict[node]
                        for nt in range(len(seq)):
                            if seq[nt].upper()!=anc_str[nt].upper() and seq[nt] != "-":
                                miss=miss+1
                    else:
                        anc_str=used_baits_dict[node]
                        if len(anc_str) != args.baitlength:
                            continue
                        for nt in range(len(seq)):
                            if seq[nt].upper()!=anc_str[nt].upper():
                                miss=miss+1
                    miss_list.append(miss)
            
                if len(miss_list)==0:
                    miss_list=[999]
            
                if min(miss_list)<= max_n_miss:
                    update_progress_bar(progress_int_2, len(tip_no_node_name))
                    continue
                
                #If no node is present, this loop will continue to look at all ancestral sequences to find the right node.

                miss_list=[]
                miss_dic={}
                seqs_already_done=[]
                if not seq in seqs_already_done:
                    seqs_already_done.append(seq)
 
                else :
                    update_progress_bar(progress_int_2, len(tip_no_node_name))
                    continue
                
                for key in node_seq_dict.keys():
                    if not "node" in key :
                        continue
                    
                    miss=0
                    if len(seq) == 0 :
                        continue

                    if len(seq)!=args.baitlength:
                        if ungap_pos[args.baitlength-1] >= seq_len:
                            continue                        
                        anc_str=get_seq(node_seq_dict, key)
                        
                        gc=((anc_str.upper().count("G")+anc_str.upper().count("C"))/args.baitlength)
                        print(gc)
                        if gc>=float(min(gc_thresholds)) and gc<=float(max(gc_thresholds)): 
                            for nt in range(len(seq)):
                                if seq[nt].upper()!=anc_str[nt].upper() and seq[nt] != "-":
                                    miss=miss+1
                        else:
                            continue
                    else:
                        if len(anc_str) != args.baitlength:
                            continue
                        ungap_pos=gap_cor_tip_dict[tip_name]
                        anc_str=get_seq(node_seq_dict, key)
                        gc=((anc_str.upper().count("G")+anc_str.upper().count("C"))/args.baitlength)
                        if gc>=float(min(gc_thresholds)) and gc<=float(max(gc_thresholds)): 
                            for nt in range(len(seq)):
                                if seq[nt].upper()!=anc_str[nt].upper():
                                    miss=miss+1
                        else:
                            continue

                    miss_dic.update({key:miss})
                if(len(miss_dic)!=0):
                    if min(miss_dic.values())>max_n_miss:
                        ungap_pos = gap_cor_tip_dict[tip_name]
                        if ungap_pos[args.baitlength-1] >= seq_len:
                            continue
                        #print(ungap_pos)
                        #print(len(ungap_pos))
                        #print(get_seq(fasta_dict, tip_name).upper())
                        bait_fasta.append(">"+tip_name+ " "+ str(ungap_pos[0])+":"+str(ungap_pos[args.baitlength-1])+"\n" +get_seq(fasta_dict, tip_name).upper()+"\n")
                        baitmer_list.append(anc_str)
                        c_baits=c_baits+1
                    else:  
                        ungap_pos = gap_cor_tip_dict[tip_name]
                        if ungap_pos[args.baitlength-1] >= seq_len:
                            continue                        
                        needed_key = [key for key in miss_dic.keys() if miss_dic[key] == min(miss_dic.values())]
                        bait_fasta.append(">"+str(sorted(needed_key)[0])+ " "+ str(ungap_pos[0])+":"+str(ungap_pos[args.baitlength-1]) + "\n" +get_seq(node_seq_dict, sorted(needed_key)[0])+"\n")
                        baitmer_list.append(anc_str)
                        c_baits=c_baits+1
                        used_nodes_it.append(sorted(needed_key)[0])

                if len(tip_no_node_name) > 200 :
                   progress_int_2.update()
        if len(tip_no_node_name) > 200 :
            progress_int_2.close()

        bait_gc=[]
        bait_gc_uniq=[]
        for bait in baitmer_list :
            if not bait in bait_gc_uniq:
                bait_gc.append((bait.upper().count("G")+bait.upper().count("C"))/args.baitlength)
                bait_gc_uniq.append(bait)
        if len(bait_gc)==0:
            print("pos in alignment no suitible baits = "+ str(pos)+":"+str(args.baitlength+pos+gap_counter))
            bait_gc.append(float(0))
            miss_list_sum_bait.append(float(0))

        if len(miss_list_sum_bait) != 0 and len(bait_gc) !=0:
            summary_stats.append(str(pos)+":"+str(args.baitlength+pos)+"\t" + \
                str(mean(miss_list_sum_bait))+ "\t" + \
                str(min(miss_list_sum_bait))+ "\t" + \
                str(max(miss_list_sum_bait))+ "\t" + \
                str(len(Counter(baitmer_list).keys()))+"\t" + \
                str(mean(bait_gc))+ "\t" + \
                str(min(bait_gc))+ "\t" + \
                str(max(bait_gc))+ "\n")
        else:
            summary_stats.append(str(pos)+":"+str(args.baitlength+pos)+"\t"+\
                str(0)+ "\t" + \
                str(0)+ "\t" + \
                str(0)+ "\t" + \
                str(0)+ "\t" + \
                str(0)+ "\t" + \
                str(0)+ "\t" + \
                str(0)+ "\n")         
        
        pos=pos+int(args.tiling)
        progress.update()
    progress.close()    
    fasta_file = open("bait_tmp.fasta", "w")
    fasta_file.writelines(bait_fasta)
    fasta_file.close()
    
    summary_stats_file = open("node_selection_stats.tsv", "w")
    summary_stats_file.writelines(summary_stats)
    summary_stats_file.close()

def bait_gc_complexity_flagger():
   ## do GC filteing
    gc_thresholds=args.GC_threshold.split("-")
    baits=[]
    for record in (SeqIO.parse("bait_tmp.fasta", "fasta")):
        bait=record.seq.upper()
        gc=((bait.count("G")+bait.count("C"))/args.baitlength)
        if(gc<float(min(gc_thresholds))):
            flag="GC="+str(gc)+"; "+"GC_flag=low;"
        elif(gc>float(max(gc_thresholds))):
            flag="GC="+str(gc)+"; "+"GC_flag=high;"
        else:
            flag="GC="+str(gc)+"; "+"GC_flag=good;"

        baits.append(SeqRecord(Seq(str(bait)), id=record.id, description=record.description+" "+flag))
    SeqIO.write(baits, "tmp_f.fasta", "fasta")

    subprocess.run(["dustmasker", "-in", "tmp_f.fasta", "-out", "tmp.dust.fasta" ,"-outfmt", "fasta", "-window", str(args.baitlength), "-level", "15"], stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)

    flagged_baits=[]
    file=open("tmp.dust.fasta")
    for record in (SeqIO.parse(file, "fasta")):
        seq=record.seq
        if not str(seq).islower() and not str(seq).isupper():
            count=0
            for nucl in seq:
                if(nucl.islower()):
                    count=count+1
            if(count>=int(args.dust_threshold)):
                record.description=record.description+"; dust_flag=failed; low_complex_lenght="+str(count)+";"
            elif(count<int(args.dust_threshold)):
                record.description=record.description+"; dust_flag=pass; low_complex_lenght="+str(count)+";"
        else:
            record.description=record.description+"; dust_flag=pass;"
        flagged_baits.append(record)

    SeqIO.write(flagged_baits, "bait.fasta", "fasta")

def summary_info_r_script():
    gc_thresholds=args.GC_threshold.split("-")
    tmp_rscript=[]
    tmp_rscript.append("library(ggplot2)")
    tmp_rscript.append("library(cowplot)")
    tmp_rscript.append("library(ggtree)")
    tmp_rscript.append("library(treeio)")
    tmp_rscript.append("library(Biostrings)")
    tmp_rscript.append("library(data.table)")
    tmp_rscript.append("library(stringr)")
    tmp_rscript.append("")
    tmp_rscript.append("fasta <- readDNAStringSet(\""+ "final_clusterd_baits.fasta" +"\")") # read the alignment file
    tmp_rscript.append("node_sum<- data.table()")
    tmp_rscript.append("node_sum<- node_sum[,.(label=gsub(\" [0-9]{1,}:[0-9]{1,}.*\", \"\", fasta@ranges@NAMES),")
    tmp_rscript.append("                       position=str_split_fixed(gsub(\"node #\", \"node_#\", fasta@ranges@NAMES), \" \", 3)[,2])]")
    tmp_rscript.append("")
    tmp_rscript.append("node_sum.dt<- node_sum[,.N,by=(label)]")
    tmp_rscript.append("node_sum.dt[, label:=(gsub(\"^[a-z]{1,} #\", \"\", label))]")
    tmp_rscript.append("")
    tmp_rscript.append("tree <- read.tree(\"" + "nodes.phylo.txt" + "\")") # read the tree with node lavels nodes.phylo.txt
    tmp_rscript.append("tree.df <-  as_tibble(tree)")
    tmp_rscript.append("tree.df$label <- gsub(\"[0-9]{1,}_\", \"\", tree.df$label)")
    tmp_rscript.append("")
    tmp_rscript.append("brlen.df <- as_tibble(read.tree(\"brlen.phylo.txt\"))")
    tmp_rscript.append("tree.df$branch.length <- brlen.df$branch.length")
    tmp_rscript.append("")
    tmp_rscript.append("node_sum.dt$label <- gsub(\"[0-9]{1,}_\", \"\", node_sum.dt$label)")
    tmp_rscript.append("tree.df$N <- NA")
    tmp_rscript.append("for (node in node_sum.dt$label){")
    tmp_rscript.append("tree.df$N[tree.df$label==node] <- node_sum.dt$N[node_sum.dt$label==node]")
    tmp_rscript.append("}")
    tmp_rscript.append("")
    tmp_rscript.append("if(sum(list.files()==\"tmp.bold\")==1){")
    tmp_rscript.append("")
    tmp_rscript.append("bold<- read.delim(\"tmp.bold\")")
    tmp_rscript.append("for (id in tree.df$label){")
    tmp_rscript.append("  tree.df$label[tree.df$label==id] <- paste(id, \" species = \", bold$species_name[bold$processid==id][1])")
    tmp_rscript.append("}")
    tmp_rscript.append("}")
    tmp_rscript.append("tree.df$label <- gsub(\"  species =  NA\" , \"\", tree.df$label)")
    tmp_rscript.append("")
    tmp_rscript.append("write.tree(as.phylo(tree.df), \"node_tree_branch_lenght.nwk\")")
    tmp_rscript.append("tree_plot <- as.treedata(tree.df)")
    tmp_rscript.append("e <- ggtree(tree_plot, layout=\"fan\")+")
    tmp_rscript.append("  geom_tippoint(aes(size=N), col=\"red\")+")
    tmp_rscript.append("  geom_tiplab(align=T, size=(4/log(length(tree_plot@phylo$tip.label))+1))+")
    tmp_rscript.append("  geom_nodepoint(aes(size=N),col=\"black\")+ggtitle(\"Tree with node counts\")+")
    tmp_rscript.append("  geom_nodelab(colour=\"blue\", size=2)")
    tmp_rscript.append("")
    tmp_rscript.append("n.df <- read.delim(\""+ "node_selection_stats.tsv" + "\")") # read summary table from node selection .py
    tmp_rscript.append("    ")
    tmp_rscript.append("plot_theme <- theme(panel.background = element_rect(fill = \"white\", colour = \"black\"),")
    tmp_rscript.append("                    strip.background = element_blank(),")
    tmp_rscript.append("                    panel.grid.minor = element_blank(),")
    tmp_rscript.append("                    panel.grid.major = element_blank(),")
    tmp_rscript.append("                    strip.text = element_text(face = \"bold\", size = 11, angle = 0, hjust = 0.5),")
    tmp_rscript.append("                    strip.background.x = element_blank(),")
    tmp_rscript.append("                    text = element_text(size = 12),")
    tmp_rscript.append("                    axis.text.x = element_text(size = , angle = 0),")
    tmp_rscript.append("                    axis.title = element_text(size=12),")
    tmp_rscript.append("                    axis.text.y = element_text(size=12),")
    tmp_rscript.append("                    plot.margin = margin(0.5, 0.5, 0.5, 0.5, \"cm\"),")
    tmp_rscript.append("                    legend.position = \"None\",")
    tmp_rscript.append("                    legend.key = element_blank(),")
    tmp_rscript.append("                    legend.text = element_text(size=12),")
    tmp_rscript.append("                    plot.title = element_text(hjust = 0.5))")
    tmp_rscript.append("")
    tmp_rscript.append("n.df$Pos <- factor(n.df$Pos, levels = n.df$Pos)")
    tmp_rscript.append("")
    tmp_rscript.append("a <- ggplot()+plot_theme+ggtitle(\"Range and mean number missmatches between baits and tip\")+")
    tmp_rscript.append("  theme(axis.text.x = element_blank())+")
    tmp_rscript.append("        #plot.margin = margin(t=0.5, r=0.5, b=-1, l=0.5, \"cm\"))+")
    tmp_rscript.append("  geom_errorbar(data=n.df, aes(Pos, ymin=min, ymax=max), alpha=0.3)+")
    tmp_rscript.append("  geom_point(data=n.df, aes(Pos, Mean), size=3)+")
    tmp_rscript.append("  scale_y_continuous(limits = c(0,7))+")
    tmp_rscript.append("  xlab(\"\")+ylab(\"Mean distance bait (in nt.)\")")
    tmp_rscript.append("  ")
    tmp_rscript.append("b <- ggplot()+plot_theme+ggtitle(\"Number of uniq baits\")+")
    tmp_rscript.append("  theme(axis.text.x = element_blank())+")
    tmp_rscript.append("        #plot.margin = margin(t=-1, r=0.5, b=-1, l=0.5, \"cm\"))+")
    tmp_rscript.append("  geom_point(data=n.df, aes(Pos, uniq), size=3)+")
    tmp_rscript.append("  xlab(\"\")+")
    tmp_rscript.append("  ylab(\"Number of uniq baits\")")
    tmp_rscript.append("")
    tmp_rscript.append("c <- ggplot()+plot_theme+ggtitle(\"Range and mean GC-content baits\")+")
    tmp_rscript.append("  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12))+")
    tmp_rscript.append("        #plot.margin = margin(t=-1, r=0.5, b=0.5, l=0.5, \"cm\"))+")
    tmp_rscript.append("  geom_errorbar(data=n.df, aes(Pos, ymin=min_gc, ymax=max_gc), alpha=0.3)+")
    tmp_rscript.append("  geom_point(data=n.df, aes(Pos, gc), size=3)+")
    tmp_rscript.append("  scale_y_continuous(limits = c(0,1))+")
    tmp_rscript.append("  geom_hline(yintercept = c("+str(min(gc_thresholds))+", "+str(max(gc_thresholds))+"), col=\"red\", linetype=\"dashed\")+")
    tmp_rscript.append("  ylab(\"Mean GC\")+xlab(\"Bait position in the selected fragment\")")
    tmp_rscript.append("")
    tmp_rscript.append("pdf(\"summary_stats.pdf\", width = 8.27, height = 11.69)")
    tmp_rscript.append("plot_grid(a,b,c, ncol = 1, align = \"vh\", axis = \"rltb\")")
    tmp_rscript.append("plot_grid(e)")
    tmp_rscript.append("dev.off()")
    tmp_rscript.append("svg(\"node_tree.svg\",  height = 20, width = 20)")
    tmp_rscript.append("plot(e)")
    tmp_rscript.append("dev.off()")


    with open('tmp.R', 'w') as tmp_rscript_file:
        for line in tmp_rscript:
            tmp_rscript_file.write(line + '\n')
    with open('tmp.R', 'w') as tmp_rscript_file:
        for line in tmp_rscript:
            tmp_rscript_file.write(line + '\n')

    with open('tmp.R_log', 'w') as tmp_rscript_log:
        subprocess.run(["Rscript" , "tmp.R"], stdout=tmp_rscript_log, stderr=subprocess.DEVNULL)

def fast_mode():
    print("\n\n###finding baits with the fast approach###")

    rc = subprocess.call(['which', 'dustmasker'],stderr=subprocess.DEVNULL, stdout=subprocess.DEVNULL)

    baits=[]
    baitsize=int(args.baitlength)

    gc_thresholds=args.GC_threshold.split("-")

    run="TRUE"
    if(run=="TRUE"):
        for record in (SeqIO.parse("tmp.aln", "fasta")):
            seq=record.seq.replace("-", "")
            for t in range(0, len(seq), int(args.tiling)):
                if(len(seq[t:(t+baitsize)])!=baitsize):
                    break
                    print(done)
                bait=seq[t:(t+baitsize)].upper()
                gc=((bait.count("G")+bait.count("C"))/baitsize)
                if(gc<float(min(gc_thresholds))):
                    flag="GC="+str(gc)+"; "+"GC_flag=low;"
                elif(gc>float(max(gc_thresholds))):
                    flag="GC="+str(gc)+"; "+"GC_flag=high;"
                else:
                    flag="GC="+str(gc)+"; "+"GC_flag=good;"

                baits.append(SeqRecord(Seq(str(bait)), id=record.id+" "+str(t)+":"+str(t+baitsize), description=flag))

        SeqIO.write(baits, "tmp_f.fasta", "fasta")

    subprocess.run(["dustmasker", "-in", "tmp_f.fasta", "-out", "tmp.dust.fasta" ,"-outfmt", "fasta", "-window", str(baitsize), "-level", "15"], stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)

    flagged_baits=[]
    file=open("tmp.dust.fasta")
    for record in (SeqIO.parse(file, "fasta")):
        seq=record.seq
        if not str(seq).islower() and not str(seq).isupper():
            count=0
            for nucl in seq:
                if(nucl.islower()):
                    count=count+1
            if(count>=int(args.dust_threshold)):
                record.description=record.description+"; dust_flag=failed; low_complex_lenght="+str(count)+";"
            elif(count<int(args.dust_threshold)):
                record.description=record.description+"; dust_flag=pass; low_complex_lenght="+str(count)+";"
        else:
            record.description=record.description+"; dust_flag=pass;"
        flagged_baits.append(record)
    SeqIO.write(flagged_baits, "bait.fasta", "fasta")

def bait_gc_complexity_filter():
    filtered_baits = []

    for record in (SeqIO.parse("bait.fasta", "fasta")):
        flags = record.description.split(" ")
        if "GC_flag=good;" and "dust_flag=pass;" in flags:
            filtered_baits.append(record)

    SeqIO.write(filtered_baits, "filtered_baits.fasta", "fasta")

    subprocess.run(["cd-hit","-i","filtered_baits.fasta","-o","filtere_uniq_baits.fasta","-c", "1"], stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)

def final_clustering():

    subprocess.run(["cd-hit","-i","filtere_uniq_baits.fasta","-o","final_clusterd_baits.fasta","-c", str(args.bts_collapse)], stdout=subprocess.DEVNULL , stderr=subprocess.DEVNULL)

def write_log_file_02():

    with open("log.txt", "a+") as log:
        log.write("\n\n\t BAITS:" + "\n\n")

        baits = len([1 for line in open("bait.fasta") if line.startswith(">")])
        log.write("\t " + str(baits) + "\t baits"+ "\n")
        baits_filterd = len([1 for line in open("filtered_baits.fasta") if line.startswith(">")])
        log.write("\t " + str(baits_filterd) + "\t baits passing gc+complexity filter"+ "\n")
        baits_colapse = len([1 for line in open("filtere_uniq_baits.fasta") if line.startswith(">")])
        log.write("\t " + str(baits_colapse) + "\t unique baits"+ "\n\n") 
        baits_colapse = len([1 for line in open("final_clusterd_baits.fasta") if line.startswith(">")])
        log.write("\t " + str(baits_colapse) + "\t unique clustered baits"+ "\n\n") 
        log.write("# analysis finished on " + datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + "\n\n")
        print("\n\nGotcha run successfully\n\na totall of " + str(baits_colapse)+ " unique and filtered baits are produced\n")

def clean_up():
    #def_prb_file = args.out + "raw_.fasta"
    os.rename('bait.fasta' , args.out + "_raw_baits.fasta")
    os.rename('filtered_baits.fasta' , args.out + "_filtered_baits.fasta")
    os.rename('filtere_uniq_baits.fasta' , args.out + "_filtered_uniq_baits.fasta")
    os.rename('final_clusterd_baits.fasta' , args.out + "_filtered_uniq_clustered_baits.fasta")

    #shutil.copy(def_prb_file, '..')
    shutil.copy("log.txt", '..')
    shutil.copy(args.out + "_raw_baits.fasta", '..')
    shutil.copy(args.out + "_filtered_baits.fasta", '..')
    shutil.copy(args.out + "_filtered_uniq_baits.fasta", '..')
    shutil.copy(args.out + "_filtered_uniq_clustered_baits.fasta", '..')
    os.rename("tmp1.fna", "original.fasta")
    shutil.copy("original.fasta", '..')
    
    def_aln_file = args.out + ".aln"
    if args.fast != True :
        shutil.copy('summary_stats.pdf', '..')
        #shutil.copy(def_tre_file, '..')
        shutil.copy("node_tree_branch_lenght.nwk", '..')
        os.rename("node_tree_branch_lenght.nwk", args.out + "_node_tree_branch_lenght.nwk")
        shutil.copy(def_aln_file, '..')
        shutil.copy("node_tree.svg", '..')

    if args.fast == True :    
        def_aln_file = args.out + ".aln"
        os.rename("tmp.aln", def_aln_file)
        shutil.copy(def_aln_file, '..')
        
    os.chdir('..')

    if args.verbose == False :
        pass
    else :
        shutil.rmtree("tmp/", ignore_errors=True) # this is not sovling the problem with folder not being removed
        #files = glob.glob("tmp/*")
        #for f in files:
        #       os.remove(f)
        #os.rmdir("tmp")

    print("\nanalysis finished on" , datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")+"\n")


###################################################
### run validator
###################################################

def hamming_dist(bait, read): 
    distance = 0
    for i in range(len(bait)):
        if len(read) < len(bait):
            break
        if bait[i].upper() != read[i].upper():
            distance += 1
    return(1-(distance/len(read)))

def make_fasta_dict():
    fasta = open("original.fasta")
    #fasta = open(def_aln_file)
    i=0
    fasta_dict={}
    for row in fasta:
        if ">" in row and i == 0:
            header=re.sub("\n|>", "", row)
            seq=str()
        elif ">" in row :
            fasta_dict.update({header:seq.upper()})
            seq=str()
            header=re.sub("\n|>", "", row)
        else:
            seq=seq+re.sub("\n", "", row)
        
        i+=1
    return(fasta_dict)

def make_aln_dict():
    alignment = open(def_aln_file)
    i=0
    alignment_dict={}
    for row in alignment:
        if ">" in row and i == 0:
            header=re.sub("\n|>", "", row)
            seq=str()
        elif ">" in row :
            alignment_dict.update({header:seq.upper()})
            seq=str()
            header=re.sub("\n|>", "", row)
        else:
            seq=seq+re.sub("\n", "", row)
        
        i+=1
    return(alignment_dict)

def make_simulated_dict_costume():
    progress = tqdm(total=args.n_simulated, desc="reads simulated ")
    #simulated_dict={}
    simulated_dict = manager.dict(simulated_dict)
    count=1

    while len(simulated_dict) != args.n_simulated :
        ran_seq=random.choice(list(fasta_dict.keys()))
        
        read_length=int(random.normal(int(args.distribution_par.split(":")[0]), int(args.distribution_par.split(":")[2])))
        pos = random.randint(0, len(fasta_dict[ran_seq]))
        if read_length <= int(args.distribution_par.split(":")[1]):
            continue
        if pos+read_length >= len(fasta_dict[ran_seq]) :
            continue
        
        simulated_dict.update({ran_seq+" sim_read_"+str(count)+"_start="+str(pos):fasta_dict[ran_seq][pos:(pos)+read_length].upper()})
        count+=1
        progress.update()
    progress.close()

    manager = multiprocessing.Manager()
    simulated_dict = manager.dict()

    return(simulated_dict)       

def make_bait_dict(): # old script
    bait_dict={}
    baits = open(args.out + "_filtered_uniq_clustered_baits.fasta")
    #baits = open(args.out + "_filtered_uniq_clustered_95_baits.fasta")
    #baits = open(args.out + "_raw_baits.fasta")
    for row in baits :
        if ">" in row:
            header = re.sub("\n|>", "", row)
            seq=str()
        else:
            seq=seq+re.sub("\n|>", "", row)

        if len(seq) == args.baitlength:
            bait_dict.update({header:seq})
    return(bait_dict)

def find_position(key): # old script
    seq_name = key.split(" ")[0]
    if seq_name in alignment_dict.keys():
        fasta_list = [alignment_dict[seq_name].upper()]
    else :
        fasta_list=[]
        for keys in alignment_dict.keys():
              fasta_list.append(alignment_dict[keys].upper())

    simulated_read = simulated_dict[key]
    print(simulated_read)
    i=0
    for fasta in fasta_list :
        for look in ["for", "rev"]:
            for pos in range(len(fasta)):
                
                seq=re.sub("-", "", fasta[pos:(pos+len(simulated_read))])

                gap_counter=0
                while(len(seq) != len(simulated_read)):
                    
                    gap_counter= fasta[pos:(pos+(len(simulated_read)+gap_counter))].count("-")
                    seq=re.sub("-", "", fasta[pos:(pos+(len(simulated_read)+gap_counter))])
                    
                    if pos+(len(simulated_read)+gap_counter) > len(fasta):
                        return(-1,-1)
                #print(gap_counter)
                #print(seq)
                if len(seq) < len(simulated_read):
                    continue

                print(fasta_list)
                if len(fasta_list)==1:
                    #print(seq)
                    #print(simulated_read)
                    if hamming_dist(seq, simulated_read) == 1 and look == "for":
                        print(seq)
                        print(simulated_read)
                        print(pos+1)
                        print(len(simulated_read))
                        return(pos+1, len(simulated_read)) 
                
                    elif hamming_dist(seq, Seq(simulated_read).reverse_complement()) == 1  and look == "rev" :

                        return(pos+1, len(simulated_read))                
                elif len(fasta_list)!=1:    
                    if hamming_dist(seq, simulated_read) >= 0.85 and look == "for":
                        return(pos+1, len(simulated_read))
                    elif hamming_dist(seq, Seq(simulated_read).reverse_complement()) >= 0.85 and look == "rev":
                        return(pos+1, len(simulated_read))
                if pos == len(fasta)-1:
                    return(-1, -1)
            i+=1

def find_pos_2(): # New much better, made with chatGPD
    progress = tqdm(total=len(simulated_dict), desc="reads parsed ")
    summary_dict={"read_name":[],
                  "bait_name":[],
                  "read_length":[],
                  "closest_dist":[],
                  "align_len":[],
                  "align_gc":[],
                  "seq":[]}
    
    for key in simulated_dict.keys():
        ham_dist_list=[]
        the_read = simulated_dict[key]
        #print(the_read)
        #print(the_read)
        max_hamming_distance = int(len(the_read)*0.2)
        #print(max_hamming_distance)
        for bait in bait_dict.keys():
            bait_seq = bait_dict[bait]
        for i in range(len(the_read)):
            for j in range(i+1, min(i+81, len(the_read)+1)):
                substring = the_read[i:j]
                if len(substring)<=args.min_overlap:
                    continue
                #print(substring)
                for k in range(len(bait_seq)):
                    bait_seq_subseq=bait_seq[k:len(substring)]
                # check if the substring matches the compare sequence within the maximum Hamming distance
                    if sum([a != b for a, b in zip(substring, bait_seq_subseq)]) <= max_hamming_distance:
                        ham_dist_list.append((len(substring),
                                              1-(sum([a != b for a, b in zip(substring, bait_seq_subseq)])/len(substring)), 
                                              bait,
                                              substring))
                        
        ham_dist_list.sort(key=lambda x: (-x[0], -x[1]))
        summary_dict["read_name"].append(key)
        summary_dict["bait_name"].append(ham_dist_list[0][2])
        summary_dict["read_length"].append(len(the_read))
        summary_dict["closest_dist"].append(ham_dist_list[0][1])
        summary_dict["align_len"].append(ham_dist_list[0][0])
        gc = (ham_dist_list[0][3].upper().count("G") + ham_dist_list[0][3].upper().count("C"))/len(ham_dist_list[0][3])
        summary_dict["align_gc"].append(gc)
        summary_dict["seq"].append(ham_dist_list[0][3].upper())
        progress.update()
    return(summary_dict)

import multiprocessing
from tqdm import tqdm

def process_key(key, simulated_dict, summary_dict, bait_dict, args):
    ham_dist_list=[]
    the_read = simulated_dict[key]
    weight=0.5
    for bait in bait_dict.keys():
        bait_seq = bait_dict[bait]
        for i in range(len(the_read)):
            for j in range(i+1, min(i+81, len(the_read)+1)):
                substring = the_read[i:j].upper()
                max_hamming_distance = int(len(substring)*0.2)
                if len(substring)<=args.min_overlap:
                    continue
                for k in range(len(bait_seq)):
                    bait_seq_subseq=bait_seq[k:len(substring)].upper()
                    if len(bait_seq_subseq) <= args.min_overlap:
                        continue
                    hamming_count=sum([a != b for a, b in zip(substring, bait_seq_subseq)])
                    if hamming_count <= max_hamming_distance:
                        value = 1-(hamming_count/len(substring)) * weight + len(bait_seq_subseq) * (1 - weight)
                        ham_dist_list.append((len(bait_seq_subseq),
                                            1-(sum([a != b for a, b in zip(substring, bait_seq_subseq)])/len(substring)), 
                                            value, #bait,
                                            substring,
                                            bait_seq_subseq))

    ham_dist_list.sort(key=lambda x: (-x[2]))

    summary_dict["read_name"].append(key)
    summary_dict["bait_name"].append(ham_dist_list[0][2])
    summary_dict["read_length"].append(len(the_read))
    summary_dict["closest_dist"].append(ham_dist_list[0][1])
    summary_dict["align_len"].append(ham_dist_list[0][0])
    gc = (ham_dist_list[0][3].upper().count("G") + ham_dist_list[0][3].upper().count("C"))/len(ham_dist_list[0][3])
    summary_dict["align_gc"].append(gc)
    summary_dict["seq"].append(ham_dist_list[0][3].upper())
    summary_dict["bait_seq_subseq"].append(ham_dist_list[0][4].upper())

def multi_process():
    fasta_dict = make_fasta_dict()
    aln_dict = make_aln_dict()
    simulated_dict = make_simulated_dict_costume(fasta_dict, aln_dict)

    p = multiprocessing.Pool(args.threads)
    results = []
    for simulated_read, simulated_seq in tqdm(simulated_dict.items(), desc="mapping simulated reads"):
        results.append(p.apply_async(map_read_to_reference, args=(simulated_read, simulated_seq, aln_dict,)))
    p.close()
    p.join()

    summary_dict = {}
    for r in results:
        summary_dict.update(r.get())

    return summary_dict


def calculate_mismatches(): # old code
    i=0
    dist_neg_count=0
    with open("validation_stats.tsv", "w") as file:
        ham_dist_list_dict={}
        progress = tqdm(total=len(simulated_dict), desc="reads parsed ")
        file.write("read_name"+ "\t"+ \
                   "bait_match" + "\t" + \
                    "read_length" + "\t" + \
                               "closest_dist" + "\t"+
                                  "align_len" + "\t" +
                                    "align_gc" + "\t" + 
                                    "seq"+"\n")
        
        for key in simulated_dict.keys():
            dist_value="-1"
            ham_dist_list=[]
            bait_list=[]
            bait_name="NA"
            pos, read_length=find_position(key)
            print("___>____>")
            print(pos)
            print(read_length)
            print("___<____<")
            if pos == -1:
                file.write(key.split(" ")[0]+"\t"+ \
                            "NA" + "\t" + \
                                 str(pos) + "\t" + \
                                         "NA" + "\t" + \
                                         "NA" + "\t" + \
                                         "NA" + "\t" + \
                                             "NA" +"\n")
                progress.update()  
                continue
            

            print(bait_dict.keys())
            miss_min = 0
            for bait in bait_dict.keys():
                the_read = simulated_dict[key]
                tile_start = int([x for x in bait.split(" ") if ":" in x][0].split(":")[0])
                bait_read_pos = len(the_read)-(tile_start-pos)
                
                if (bait_read_pos<0) or ((bait_read_pos>len(the_read) + (args.baitlength-args.min_overlap))):
                    continue
                if (bait_read_pos<args.min_overlap):
                    continue
                tile_correction = (tile_start-pos)+1

                print("\n")
                print("pos = " + str(pos))
                print("tile = " + str(tile_start))
                print('the_read = ' + the_read)
                print('bait_seq = ' + bait_dict[bait][0:(len(the_read)-tile_correction)] )

                if tile_correction>0:
                    bait_seq=bait_dict[bait][0:(len(the_read)-tile_correction)] 
                    the_read=the_read[tile_correction:len(the_read)]
                elif pos - tile_start == 1:
                    bait_seq=bait_seq
                    the_read=the_read
                else :
                    if pos+read_length > tile_start+args.baitlength:
                        if tile_correction < 0:
                            bait_seq=bait_dict[bait][(tile_correction*-1):(tile_start+args.baitlength)]
                        else: 
                            bait_seq=bait_dict[bait][(pos):(tile_start+args.baitlength)]
                    else:
                        if tile_correction < 0:
                            bait_seq=bait_dict[bait][(tile_correction*-1):(tile_correction*-1+read_length)]
                        else:
                            bait_seq=bait_dict[bait][(pos):(pos+read_length)]
                    if len(bait_seq) < len(the_read):
                        the_read=the_read[0:len(bait_seq)]
                    elif len(bait_seq) > len(the_read):
                        bait_seq=the_read[0:len(the_read)] 
                
                ham_dist_list.append(str(len(bait_seq))+":"+str('%.3f'%(hamming_dist(bait_seq, the_read))))
                bait_list.append(bait)
                miss_min = max(hamming_dist(bait_seq, the_read), miss_min)
                ham_dist_list_dict.update({str(len(bait_seq))+":"+str('%.3f'%(hamming_dist(bait_seq, the_read))):[the_read, bait,bait_seq]})
                print('cor read = ' + the_read)
                print('cor bait = ' + bait_seq)

            print("TRUE")
            print(sorted(ham_dist_list,reverse=True))
            print(ham_dist_list_dict)
            quit()    
            if len(ham_dist_list) != 0:
                dist_value = sorted(ham_dist_list,reverse=True)[0]
                for I in range(len(bait_list)):
                    if ham_dist_list[I]==dist_value:
                        if "node" in bait_list[I]:
                            bait_name=bait_list[I]
                        elif bait_name == "NA":
                            bait_name=bait_list[I]
            if dist_value == "-1" :
                progress.update()  
                continue

            if float(dist_value.split(":")[1])< 0.60 :
                ali_lenght = int(sorted(ham_dist_list,reverse=True)[0].split(":")[0])

                for i in range(1, len(sorted(ham_dist_list,reverse=True))):
                    if int(sorted(ham_dist_list,reverse=True)[i].split(":")[0]) != ali_lenght:
                        ali_lenght = int(sorted(ham_dist_list,reverse=True)[i].split(":")[0])
                    if float(sorted(ham_dist_list,reverse=True)[i].split(":")[1]) > 0.8 :
                        dist_value=sorted(ham_dist_list,reverse=True)[i]
                        break
                
            progress.update()                   
            if dist_value=="-1":
                dist_neg_count+=1
                continue
            
            file.write(key.split(" ")[0]+"\t"+ \
                        re.sub(" ", "_", bait_name) + "\t" + \
                             str(read_length) + "\t" + \
                                    str(dist_value.split(":")[1]) + "\t" + \
                                            str(dist_value.split(":")[0]) + "\t" + \
                                               str(sum(the_read.count(x) for x in ("C", "G"))/len(the_read)) + "\t" +\
                                                   the_read + "\n")            
            i+=1    
        progress.close()

def print_pandas():
    stats_pf = pd.read_csv("validation_stats.tsv",  delimiter= "\t")
    gc_thresholds=args.GC_threshold.split("-")
    plot=1
    plt.rcParams['figure.figsize'] = [15,10]
    for i in range(50, 170, 10):

        sub_bin = stats_pf.loc[(stats_pf['read_length']>=i) &(stats_pf['read_length']<i+10)]
        sub_bin_gc = sub_bin.loc[(sub_bin['align_gc']<float(min(gc_thresholds))) | (sub_bin['align_gc']>float(max(gc_thresholds)))]
        sub_bin_gc_right = sub_bin.loc[(sub_bin['align_gc']>float(min(gc_thresholds))) & (sub_bin['align_gc']<float(max(gc_thresholds)))]
        #print(len(sub_bin))
        #print(len(sub_bin_gc_right))
        plt.subplot(3,4,plot)
        plt.plot(sub_bin['closest_dist'], 
                 sub_bin['align_len'], 'ko', markersize=4, alpha=0.05)
        plt.plot(sub_bin_gc['closest_dist'],
                 sub_bin_gc['align_len'], "ro", markersize=2, alpha=0.25)
        plt.ylim(ymin = 0, ymax = 81)
        plt.xlim(xmin = 0.80, xmax = 1.005)
        plt.gca().spines.top.set_visible(False)
        plt.gca().spines.right.set_visible(False)
        plt.xticks(np.arange(0.80, 1.05, 0.05))
        if len(sub_bin) == 0 :
            continue
        plt.text(0.81, 81.5, "Read bin = "+str(i)+":"+str(i+10))
        plt.text(0.81, 20.0, "N = "+str(len(sub_bin))+\
                                   " ("+str(len(sub_bin_gc_right))+")")
        plt.text(0.81, 15.0, "91% = "+str(round(sum(sub_bin['closest_dist']>=0.90)/len(sub_bin),2))+" ("+\
                                      str(round(sum(sub_bin_gc_right['closest_dist']>=0.90)/len(sub_bin_gc_right),2))+")")
        plt.text(0.81, 10.0, "50nt, 91% = "+str(round(len(sub_bin.loc[(sub_bin['closest_dist']>=0.90) & (sub_bin['align_len']>=50)])/len(sub_bin), 2)) +\
                                           " ("+str(round(len(sub_bin_gc_right.loc[(sub_bin_gc_right['closest_dist']>=0.90) & (sub_bin_gc_right['align_len']>=50)])/len(sub_bin_gc_right), 2))+")")
        plt.text(0.81, 5.0, "GC "+str(int(float(min(gc_thresholds))*100))+"-"+ \
                 str(int(float(max(gc_thresholds))*100))+"% = "+\
                    str(round(len(sub_bin.loc[(sub_bin['align_gc']>float(min(gc_thresholds))) & 
                                              (sub_bin['align_gc']<float(max(gc_thresholds)))])/len(sub_bin),2))+\
                    "( "+ str(round(len(sub_bin_gc_right.loc[(sub_bin_gc_right['align_gc']>float(min(gc_thresholds))) & 
                                              (sub_bin_gc_right['align_gc']<float(max(gc_thresholds)))])/len(sub_bin_gc_right),2))+")")

        plot+=1
    plt.savefig("validation_plot.svg")
if __name__ == '__main__':
    if (args.run_validation_only != True):
        print("Run Gotcha Main")
        coding, marker_list = check_if_all_is_installed()

        #####################################################################
        #                         starts the log file                       #
        #####################################################################
        
        write_log_file_01()

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


        #####################################################################
        #  downloads sequences from bold if no costume alingent or fasta is given
        #####################################################################
        
        if (args.custom_fna == "") and (args.custom_aln == ""):
            bold_download()

        #####################################################################
        #                 set of quality filters fo any input               #
        #####################################################################
    
        seq_quality_filters()

        #####################################################################
        #            size selection of the reference sequences              #
        #####################################################################
        if args.no_size_selection == True:
            size_selection()
        
        #####################################################################


        #####################################################################
        # an aditonal optional timal filter, otherwise standard filer is applied
        #####################################################################
    
        if args.trim_seqoverlap != "" and args.trim_seqoverlap != "" :
            extra_trimal_fiter()
        else:
            standard_trimal_filter()
        

        #####################################################################
        #           stop the gotch for a manual quality check               #
        #####################################################################
    
        if (args.manual_check == True):
            shutil.copy('tmp.aln', '..')
            print("\nstopping prior to tree inference for a manual check of the alignment\n")
            quit()


        #####################################################################
        #               Phylogentic backbone bait finding                   #
        #####################################################################
        
        if (args.fast == False):
            def_tre_file, def_aln_file = tree_inverence()
            bait_finder()
            bait_gc_complexity_flagger()


        #####################################################################
        #                           fast mode                               #
        #####################################################################

        if (args.fast == True):
            fast_mode()

        #####################################################################
        
        
        #####################################################################
                            # final filtering and clean up
        #####################################################################
        bait_gc_complexity_filter()
        final_clustering()
        write_log_file_02()
        if (args.fast == False):
            summary_info_r_script()
        clean_up()
        #####################################################################
    
    if (args.run_validation == True) or (args.run_validation_only == True):
        def_aln_file = args.out + ".aln"
        if(args.run_validation_only == True):
            os.chdir(args.out)
        print("Run Gotcha Validator")
        fasta_dict = make_fasta_dict()
        simulated_dict = make_simulated_dict_costume()
        bait_dict = make_bait_dict()
        alignment_dict = make_aln_dict()

        progress = tqdm(total=len(simulated_dict), desc="reads parsed ")
        summary_dict={"read_name":[],
                    "bait_name":[],
                    "read_length":[],
                    "closest_dist":[],
                    "align_len":[],
                    "align_gc":[],
                    "seq":[],
                    "bait_seq_subseq":[]}
        #summary_dict=find_pos_2()
        summary_dict = multi_process()
        df = pd.DataFrame(summary_dict)
        df.to_csv("validation_stats.tsv", sep='\t', index=False)

        #calculate_mismatches()
        print_pandas()