tb-profiler

#! /usr/bin/env python3
import sys
import pathogenprofiler as pp
import argparse
import json
import tbprofiler as tbp
import os
import csv
from datetime import datetime
from uuid import uuid4
import glob
import atexit

files_prefix = None
conf = None
@atexit.register
def cleanup():
    if "last_traceback" in vars(sys):
        if files_prefix and not args.no_cleanup:
            sys.stderr.write("Cleaning up after failed run\n")
            for f in glob.glob(files_prefix+"*"):
                os.remove(f)
        import traceback
        with open(args.prefix+".errlog", "w") as O:
            O.write("# tb-profiler error report\n\n")
            O.write("* OS: %s\n" % sys.platform)
            O.write("* Program version: %s\n" % tbp._VERSION)
            O.write("* Database version: %s\n" % json.load(open(conf["version"]))["commit"])
            O.write("* Program call:\n")
            O.write("```\n")
            O.write("%s\n" % vars(args))
            O.write("```\n")

            O.write("## Traceback:\n")
            O.write("```\n")
            traceback.print_tb(sys.last_traceback,file=O)
            O.write("```\n")

            O.write("## Value:\n")
            O.write("```\n")
            O.write("%s" % sys.last_value)
            O.write("```\n")
        sys.stderr.write("""\n
################################# ERROR #######################################

This run has failed. Please check all arguments and make sure all input files
exist. If no solution is found, please open up an issue at
https://github.com/jodyphelan/TBProfiler/issues/new and paste or attach the
contents of the error log (%s)

###############################################################################
""" % (args.prefix+".errlog"))


try:
    sys.base_prefix
except:
    sys.base_prefix = getattr(sys, 'base_prefix', getattr(sys, 'real_prefix', sys.prefix))

def get_conf_dict(library_prefix):
    files = {"gff":".gff","ref":".fasta","ann":".ann.txt","barcode":".barcode.bed","bed":".bed","json_db":".dr.json","version":".version.json"}
    global conf
    conf = {}
    for key in files:
        sys.stderr.write("Using %s file: %s\n" % (key,library_prefix+files[key]))
        conf[key] = pp.filecheck(library_prefix+files[key])
    return conf

def main_reprofile(args):
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)
    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    old_results = json.load(open(args.json))
    new_results = old_results.copy()
    variant_dump = {}
    for var in old_results["dr_variants"]:
        del var["drug"]
        var["gene_id"] = var["locus_tag"]
        var["change"] = var["_internal_change"]
        variant_dump[(var["locus_tag"],var["change"])] = var
    for var in old_results["other_variants"]:
        var["gene_id"] = var["locus_tag"]
        var["change"] = var["_internal_change"]
        variant_dump[(var["locus_tag"],var["change"])] = var
    new_results["variants"] = list(variant_dump.values())
    del new_results["other_variants"]
    del new_results["dr_variants"]
    new_results = pp.db_compare(db_file=conf["json_db"],mutations=new_results)
    tbp.reformat_annotations(new_results,conf)
    for var in new_results["dr_variants"]:
        del var["gene_id"]
    for var in new_results["other_variants"]:
        del var["gene_id"]
    new_results["db_version"] = json.load(open(conf["version"]))
    json.dump(new_results,open("%s.results.json"%args.prefix,"w"))

def main_profile(args):
    #### Setup conf dictionary ###
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)
    global conf
    if args.external_db:
        conf = tbp.get_conf_dict_with_path(args.external_db)
    else:
        conf = tbp.get_conf_dict(args.db)

    ### Create folders for results if they don't exist ###
    if pp.nofolder(args.dir):
        os.mkdir(args.dir)



    ### Set up platform dependant parameters ###
    if args.platform=="nanopore":
        args.mapper = "minimap2"
        args.caller = "bcftools"
        args.no_trim=True
        run_delly = False
    else:
        if args.no_delly:
            run_delly = False
        else:
            run_delly = True

    ### Setup prefix for files ###
    global files_prefix
    args.tmp_prefix = str(uuid4())
    files_prefix = args.dir+"/"+args.tmp_prefix

    ### Create bam file if fastq has been supplied ###
    if args.bam==None:
        if args.read1 and args.read2 and args.no_trim:
            # Paired + no trimming
            fastq_obj = pp.fastq(args.read1,args.read2)
        elif args.read1 and args.read2 and not args.no_trim:
            # Paired + trimming
            untrimmed_fastq_obj = pp.fastq(args.read1,args.read2)
            fastq_obj = untrimmed_fastq_obj.trim(files_prefix,threads=args.threads)
        elif args.read1 and not args.read2 and args.no_trim:
            # Unpaired + trimming
            fastq_obj = pp.fastq(args.read1,args.read2)
        elif args.read1 and not args.read2 and not args.no_trim:
            # Unpaired + trimming
            untrimmed_fastq_obj = pp.fastq(args.read1)
            fastq_obj = untrimmed_fastq_obj.trim(files_prefix,threads=args.threads)
        else:
            exit("\nPlease provide a bam file or a fastq file(s)...Exiting!\n")
        bam_obj = fastq_obj.map_to_ref(
            ref_file=conf["ref"], prefix=files_prefix,sample_name=args.prefix,
            aligner=args.mapper, platform=args.platform, threads=args.threads
        )
        bam_file = bam_obj.bam_file
    else:
        bam_file = args.bam

    if not args.missing_cov_threshold:
        args.missing_cov_threshold = args.min_depth
    else:
        sys.stderr.write("\nWARNING: The --missing_cov_threshold argument is deprecated and will be removed in future releases. This parameter can now be set with --min_depth.\n")
    ### Run profiling module from pathogen-profiler ###
    results = pp.bam_profiler(
        conf=conf, bam_file=bam_file, prefix=files_prefix, platform=args.platform,
        caller=args.caller, threads=args.threads, no_flagstat=args.no_flagstat,
        run_delly = run_delly, calling_params=args.calling_params,
        coverage_fraction_threshold=args.coverage_fraction_threshold,
        missing_cov_threshold=args.missing_cov_threshold, samclip=True,
        min_depth=args.min_depth,delly_bcf_file=args.delly_bcf,call_wg=args.call_whole_genome
    )

    ### Reformat the results to TB-Profiler style ###
    results = tbp.reformat(results, conf, reporting_af=args.reporting_af, mutation_metadata=args.add_mutation_metadata)
    results["id"] = args.prefix
    results["tbprofiler_version"] = tbp._VERSION
    results["pipeline"] = {"mapper":args.mapper if not args.bam else "N/A","variant_caller":args.caller}

    for x in ["bam","vcf","results"]:
        if pp.nofolder(args.dir+"/"+x):
            os.mkdir(args.dir+"/"+x)

    json_output = args.dir+"/results/"+args.prefix+".results.json"
    text_output = args.dir+"/results/"+args.prefix+".results.txt"
    csv_output = args.dir+"/results/"+args.prefix+".results.csv"
    pdf_output = args.dir+"/results/"+args.prefix+".results.pdf"

    results["timestamp"] = datetime.now().strftime("%d-%m-%Y %H:%M:%S")
    json.dump(results,open(json_output,"w"))
    extra_columns = [x.lower() for x in args.add_columns.split(",")] if args.add_columns else []
    if args.pdf:
        tbp.write_pdf(results,conf,pdf_output)
    if args.txt:
        tbp.write_text(results,conf,text_output,extra_columns,reporting_af=args.reporting_af)
    if args.csv:
        tbp.write_csv(results,conf,csv_output,extra_columns)

    ### Move files to respective directories ###
    result_files = {
        "%s.delly.bcf" % files_prefix: "%s/vcf/%s.delly.bcf" % (args.dir,args.prefix),
        "%s.targets.csq.vcf.gz" % files_prefix: "%s/vcf/%s.targets.csq.vcf.gz" % (args.dir,args.prefix),
        "%s.vcf.gz" % files_prefix: "%s/vcf/%s.vcf.gz" % (args.dir,args.prefix),
        "%s.bam" % files_prefix: "%s/bam/%s.bam" % (args.dir,args.prefix),
        "%s.bam.bai" % files_prefix: "%s/bam/%s.bam.bai" % (args.dir,args.prefix),
    }
    for file,target in result_files.items():
        if os.path.isfile(file):
            os.rename(file,target)

    pp.run_cmd("rm %s*" % files_prefix)

    pp.log("Profiling finished sucessfully!")

def main_update_tbdb(args):
    if pp.nofolder("tbdb"):
        pp.run_cmd("git clone https://github.com/jodyphelan/tbdb.git")
    os.chdir("tbdb")
    pp.run_cmd("git pull")
    pp.run_cmd("tb-profiler create_db --seqname %s" % args.seqname)
    pp.run_cmd("tb-profiler load_library tbdb")
    os.chdir("../")
    pp.log("Sucessfully updated TBDB")

def main_create_db(args):
    tbp.create_db(args)

def main_load_library(args):
    lib_prefix = args.prefix.split("/")[-1]
    files = {"gff":".gff","ref":".fasta","ann":".ann.txt","barcode":".barcode.bed","bed":".bed","json_db":".dr.json","version":".version.json"}
    if pp.nofolder(sys.base_prefix+"/share/tbprofiler"):
        pp.run_cmd("mkdir %s " % (sys.base_prefix+"/share/tbprofiler/"))
    for key in files:
        new_file_location = sys.base_prefix+"/share/tbprofiler/"+lib_prefix+files[key]
        pp.run_cmd("cp %s %s" % (args.prefix+files[key],new_file_location))
    pp.run_cmd("samtools faidx %s" % sys.base_prefix+"/share/tbprofiler/"+lib_prefix+".fasta")
    pp.bwa_index(sys.base_prefix+"/share/tbprofiler/"+lib_prefix+".fasta")
    if os.path.isfile("%s" % sys.base_prefix+"/share/tbprofiler/"+lib_prefix+".dict"):
        pp.run_cmd("rm %s" % sys.base_prefix+"/share/tbprofiler/"+lib_prefix+".dict")
    pp.create_seq_dict(sys.base_prefix+"/share/tbprofiler/"+lib_prefix+".fasta")
    pp.log("Sucessfully imported library")

def main_lineage(args):
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)
    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    pp.filecheck(args.bam)
    bam = pp.bam(args.bam,args.bam,conf["ref"])
    mutations = bam.get_bed_gt(conf["barcode"],conf["ref"],args.caller)

    results = {}
    snps_file = args.prefix+".lineage.snps.txt" if args.snps else None

    results["barcode"] = pp.barcode(mutations,conf["barcode"],snps_file=snps_file)
    tbp.barcode2lineage(results)
    outfile = "%s.lineage.%s" % (args.prefix,args.outfmt)
    O = open(outfile,"w")
    if args.outfmt=="json":
        json.dump(results,O)
    elif args.outfmt=="txt":
        O.write(tbp.text.lineagejson2text(results["lineage"]))
    O.close()

def main_collate(args):
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)
    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    tbp.collate_results(args.prefix,conf,sample_file=args.samples,result_dir=args.dir,reporting_af=args.reporting_af,mark_missing=args.mark_missing)

def main_version(args):
    print("\nTBProfiler version %s\n" % tbp._VERSION)

def main_reformat(args):
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)

    results = json.load(open(args.json))
    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    args.prefix = results["id"]
    tex_output = args.prefix+".results.tex"
    csv_output = args.prefix+".results.csv"
    text_output = args.prefix+".results.txt"
    # if args.pdf:
    #     tbp.write_tex(results,conf,tex_output)
    #     pp.run_cmd("pdflatex %s"%tex_output,verbose=1)
    #     pp.rm_files([tex_output, args.prefix+".results.aux",args.prefix+".results.log"])
    if args.txt:
        tbp.write_text(results,conf,text_output)
    if args.csv:
        tbp.write_csv(results,conf,csv_output)


def main_fasta_profile(args):
    for x in ["bam","vcf","results"]:
        if pp.nofolder(args.dir+"/"+x):
            os.mkdir(args.dir+"/"+x)
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)

    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    results = pp.fasta_profiler(conf, args.prefix, args.fasta)
    results = tbp.reformat(results,conf,reporting_af=0)
    results["id"] = args.prefix
    results["tbprofiler_version"] = tbp._VERSION
    results["pipeline"] = {"mapper":"N/A","variant_caller":"paftools.js"}

    json_output = args.dir+"/results/"+args.prefix+".results.json"
    tex_output = args.dir+"/results/"+args.prefix+".results.tex"
    text_output = args.dir+"/results/"+args.prefix+".results.txt"
    csv_output = args.dir+"/results/"+args.prefix+".results.csv"

    json.dump(results,open(json_output,"w"))
    extra_columns = [x.lower() for x in args.add_columns.split(",")] if args.add_columns else []
    # if args.pdf:
    #     tbp.write_tex(results,conf,tex_output,extra_columns)
    #     pp.run_cmd("pdflatex %s"%tex_output,verbose=1)
    #     pp.rm_files([tex_output, args.dir+"/"+args.prefix+".results.aux",args.dir+"/"+args.prefix+".results.log"])
    if args.txt:
        tbp.write_text(results,conf,text_output,extra_columns,reporting_af=0.7)
    if args.csv:
        tbp.write_csv(results,conf,csv_output,extra_columns)

    pp.run_cmd("mv -f %(dir)s/%(prefix)s*.vcf.gz* %(dir)s/vcf/" % vars(args))



def main_vcf_profile(args):
    for x in ["bam","vcf","results"]:
        if pp.nofolder(args.dir+"/"+x):
            os.mkdir(args.dir+"/"+x)

    if args.lofreq_sample_name:
        import re
        modified_vcf = "%s.lofreq_modified.vcf" % pp.vcf(args.vcf).prefix
        with open(modified_vcf,"w") as O:
            for l in pp.cmd_out("bcftools view %(vcf)s" % vars(args)):
                row = l.strip().split()
                if "INFO=<ID=DP" in l:
                    O.write("##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")
                    O.write("##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"Allelic depths (high-quality bases)\">\n")

                    O.write(l+"\n")
                elif l[:2]=="##":
                    O.write(l+"\n")
                elif row[0]=="#CHROM":
                    row = row + ["FORMAT", args.lofreq_sample_name]
                    O.write("%s\n" % "\t".join(row))
                else:
                    re_obj = re.search("DP4=(\d+),(\d+),(\d+),(\d+)",l)
                    dp4_1 = int(re_obj.group(1))
                    dp4_2 = int(re_obj.group(2))
                    dp4_3 = int(re_obj.group(3))
                    dp4_4 = int(re_obj.group(4))
                    ad = "%s,%s" % ( (dp4_1+dp4_2), (dp4_3+dp4_4))
                    row = row + ["GT:AD", "0/1:%s" %ad]
                    O.write("%s\n" % "\t".join(row))
        pp.run_cmd("bgzip -f %s" % modified_vcf)
        args.vcf = modified_vcf+".gz"

    vcf_obj = pp.vcf(args.vcf)
    if args.db=="tbdb" and not args.external_db and pp.nofile(sys.base_prefix+"/share/tbprofiler/tbdb.fasta"):
        pp.log("Can't find the tbdb file at %s. Please run 'tb-profiler update_tbdb' to load the default library or specify another using the '--external_db' flag" % sys.base_prefix,ext=True)

    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)
    for sample_name in vcf_obj.samples:
        args.sample_name = sample_name
        files_prefix = args.dir+"/"+sample_name
        args.tmp_vcf = "%s.vcf.gz" % uuid4()
        if args.reporting_af>0:
            pp.run_cmd("bcftools view -s %(sample_name)s -c 1 %(vcf)s | bcftools view -i 'sum(AD)>0' | bcftools +fixploidy -Oz -o %(tmp_vcf)s " % vars(args))
        else:
            pp.run_cmd("bcftools view -s %(sample_name)s -c 1 %(vcf)s | bcftools +fixploidy -Oz -o %(tmp_vcf)s " % vars(args))
        results = pp.vcf_profiler(conf, files_prefix, sample_name, args.tmp_vcf)
        results = tbp.reformat(results,conf,reporting_af=args.reporting_af)
        results["id"] = sample_name
        results["tbprofiler_version"] = tbp._VERSION
        results["pipeline"] = {"mapper": "N/A","variant_caller":"N/A"}

        json_output = args.dir+"/results/"+sample_name+".results.json"
        tex_output = args.dir+"/results/"+sample_name+".results.tex"
        text_output = args.dir+"/results/"+sample_name+".results.txt"
        csv_output = args.dir+"/results/"+sample_name+".results.csv"

        json.dump(results,open(json_output,"w"))
        extra_columns = [x.lower() for x in args.add_columns.split(",")] if args.add_columns else []
        # if args.pdf:
        #     tbp.write_tex(results,conf,tex_output,extra_columns)
        #     pp.run_cmd("pdflatex %s"%tex_output,verbose=1)
        #     pp.rm_files([tex_output, args.dir+"/"+args.prefix+".results.aux",args.dir+"/"+args.prefix+".results.log"])
        if args.txt:
            tbp.write_text(results,conf,text_output,extra_columns)
        if args.csv:
            tbp.write_csv(results,conf,csv_output,extra_columns)
        pp.run_cmd("rm %(tmp_vcf)s*" % vars(args))

    if args.lofreq_sample_name:
        pp.run_cmd("rm %(vcf)s*" % vars(args))


def main_test(args):
    pp.run_cmd("tb-profiler profile -1 %s" % (sys.base_prefix+"/share/tbprofiler/tbprofiler.test.fq.gz"),verbose=2)

#### Argument Parsing ####

parser = argparse.ArgumentParser(description='TBProfiler pipeline',formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
subparsers = parser.add_subparsers(help="Task to perform")

# Profile #
parser_sub = subparsers.add_parser('profile', help='Run whole profiling pipeline', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
input=parser_sub.add_argument_group("Input options")
group = input.add_mutually_exclusive_group(required=True)
group.add_argument('--read1','-1',help='First read file')
input.add_argument('--read2','-2',help='Second read file')
group.add_argument('--bam','-a',help='BAM file. Make sure it has been generated using the H37Rv genome (GCA_000195955.2)')
input.add_argument('--platform','-m',choices=["illumina","nanopore"],default="illumina",help='NGS Platform used to generate data')
input.add_argument('--db',default='tbdb',help='Mutation panel name')
input.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')

output=parser_sub.add_argument_group("Output options")
output.add_argument('--prefix','-p',default="tbprofiler",help='Sample prefix for all results generated')
output.add_argument('--dir','-d',default=".",help='Storage directory')
output.add_argument('--csv',action="store_true",help="Add CSV output")
output.add_argument('--txt',action="store_true",help="Add text output")
output.add_argument('--pdf',action="store_true",help="Add PDF output. This requires pdflatex to be installed")
output.add_argument('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and csv results")
output.add_argument('--add_mutation_metadata',action="store_true",help="Add mutation data from the tbdr database (not currently functional)")
output.add_argument('--call_whole_genome',action="store_true",help="Call whole genome")

algorithm=parser_sub.add_argument_group("Algorithm options")
algorithm.add_argument('--mapper',default="bwa", choices=["bwa","minimap2","bowtie2","bwa-mem2"],help="Mapping tool to use. If you are using nanopore data it will default to minimap2",type=str)
algorithm.add_argument('--caller',default="bcftools", choices=["bcftools","gatk","freebayes"],help="Variant calling tool to use.",type=str)
algorithm.add_argument('--calling_params',type=str,help='Override default parameters for variant calling')
algorithm.add_argument('--min_depth',default=10,type=int,help='Minimum depth required to call variant. Bases with depth below this cutoff will be marked as missing')
algorithm.add_argument('--af',default=0.1,type=float,help='Minimum allele frequency to call variants')
algorithm.add_argument('--reporting_af',default=0.1,type=float,help='Minimum allele frequency to use variants for prediction')
algorithm.add_argument('--coverage_fraction_threshold',default=0,type=int,help='Cutoff used to calculate fraction of region covered by <= this value')
algorithm.add_argument('--missing_cov_threshold',type=int,help='Cutoff used to positions/codons in genes which are missing (this argument has now been merged with --min_depth argument and will be deprecated in future releases)')
algorithm.add_argument('--no_trim',action="store_true",help="Don't trim files using trimmomatic")
algorithm.add_argument('--no_flagstat',action="store_true",help="Don't collect flagstats")
algorithm.add_argument('--no_delly',action="store_true",help="Don't run delly")
algorithm.add_argument('--threads','-t',default=1,help='Threads to use',type=int)

other=parser_sub.add_argument_group("Other options")
other.add_argument('--verbose','-v',default=0, choices=[0,1,2],help="Verbosity increases from 0 to 2",type=int)
other.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
other.add_argument('--no_cleanup',action="store_true",help="Don't remove temporary files on error")
other.add_argument('--delly_bcf',help=argparse.SUPPRESS)
parser_sub.set_defaults(func=main_profile)


# VCF profile #
parser_sub = subparsers.add_parser('vcf_profile', help='Run profiling pipeline on VCF file. Warning: this assumes that you have good coverage across the genome', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('vcf',help='VCF file')
parser_sub.add_argument('--db',default='tbdb',help='Mutation panel name')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--dir','-d',default=".",help='Storage directory')
parser_sub.add_argument('--reporting_af',default=0.1,type=float,help='Minimum allele frequency to call variants')
parser_sub.add_argument('--lofreq_sample_name',help='Sample name to use if VCF has been generated by lofreq')
parser_sub.add_argument('--txt',action="store_true",help="Add text output")
parser_sub.add_argument('--csv',action="store_true",help="Add CSV output")
parser_sub.add_argument('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and csv results")
parser_sub.add_argument('--verbose','-v',default=0, choices=[0,1,2],help="Verbosity increases from 0 to 2",type=int)
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_vcf_profile)


# Fasta profile #
parser_sub = subparsers.add_parser('fasta_profile', help='(BETA) Run profiling pipeline on Fasta file. Warning: this assumes that this is a good quality assembly which coveres all drug resistance loci', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('--fasta','-f',help='VCF file',required=True)
parser_sub.add_argument('--prefix','-p',help='VCF file',required=True)
parser_sub.add_argument('--db',default='tbdb',help='Mutation panel name')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--dir','-d',default=".",help='Storage directory')
parser_sub.add_argument('--txt',action="store_true",help="Add text output")
parser_sub.add_argument('--csv',action="store_true",help="Add CSV output")
parser_sub.add_argument('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and csv results")
parser_sub.add_argument('--verbose','-v',default=0, choices=[0,1,2],help="Verbosity increases from 0 to 2",type=int)
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_fasta_profile)


parser_sub = subparsers.add_parser('lineage', help='Profile only lineage', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('--bam','-a',required=True, help='BAM file. Make sure it has been generated using the H37Rv genome (GCA_000195955.2)')
parser_sub.add_argument('--prefix','-p',default="tbprofiler",help='Sample prefix')
parser_sub.add_argument('--outfmt',default='json',choices=["json","txt"],type=str,help="Output format")
parser_sub.add_argument('--snps',action="store_true",help='Sample prefix')
parser_sub.add_argument('--caller',default='bcftools',choices=["bcftools","freebayes","gatk"],type=str,help="Variant caller")
parser_sub.add_argument('--db',default='tbdb',help='Mutation panel name')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_lineage)

parser_sub = subparsers.add_parser('collate', help='Collate results form multiple samples together', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('--prefix','-p',default="tbprofiler",help='Sample prefix')
parser_sub.add_argument('--samples',help='File with samples (one per line)')
parser_sub.add_argument('--dir','-d',default="results",help='Storage directory')
parser_sub.add_argument('--full',action="store_true",help='Output mutations in main result file')
parser_sub.add_argument('--all_variants',action="store_true",help='Output all variants in variant matrix')
parser_sub.add_argument('--mark_missing',action="store_true",help='An asteriks will be use to mark predictions which are affected by missing data at a drug resistance position')
parser_sub.add_argument('--reporting_af',default=0.1,type=float,help='Minimum allele frequency to call variants')
parser_sub.add_argument('--db',default='tbdb',help='Full path to mutation database json file to use')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_collate)



parser_sub = subparsers.add_parser('reprofile', help='Reprofile previous results using a new library. The new library must have same targets and the old one.', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('json',help='JSON output file')
parser_sub.add_argument('--prefix','-p',default="tbprofiler",help='Sample prefix')
parser_sub.add_argument('--db',default='tbdb',help='Mutation panel name')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_reprofile)


parser_sub = subparsers.add_parser('reformat', help='Reformat json results into text or csv', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('json',default="tbprofiler",help='Sample prefix')
parser_sub.add_argument('--txt',action="store_true",help="Add text output")
parser_sub.add_argument('--csv',action="store_true",help="Add CSV output")
# parser_sub.add_argument('--pdf',action="store_true",help="Add PDF output. This requires pdflatex to be installed")
parser_sub.add_argument('--db',default='tbdb',help='Mutation panel name')
parser_sub.add_argument('--external_db',type=str,help='Path to db files prefix (overrides "--db" parameter)')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_reformat)

parser_sub = subparsers.add_parser('create_db', help='Generate the files required to run TBProfiler', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('--prefix','-p',default="tbdb",type=str,help='The input CSV file containing the mutations')
parser_sub.add_argument('--csv','-c',default="tbdb.csv",type=str,help='The prefix for all output files')
parser_sub.add_argument('--watchlist','-w',default="tbdb.watchlist.csv",type=str,help='A csv file containing genes to profile but without any specific associated mutations')
parser_sub.add_argument('--seqname','-s',default="Chromosome",type=str,help='The prefix for all output files')
parser_sub.add_argument('--confidence',default="tbdb.confidence.csv",type=str,help="A CSV containing gene, mutation, drug and confidence columns")
parser_sub.add_argument('--custom',action="store_true",help='Tells the script this is a custom database, this is used to alter the generation of the version file')
parser_sub.add_argument('--db-name',help='Overrides the name of the database in the version file')
parser_sub.add_argument('--db-commit',help='Overrides the commit string of the database in the version file')
parser_sub.add_argument('--db-author',help='Overrides the author of the database in the version file')
parser_sub.add_argument('--db-date',help='Overrides the date of the database in the version file')
parser_sub.set_defaults(func=main_create_db)

parser_sub = subparsers.add_parser('load_library', help='Load new library', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('prefix',type=str,help='Prefix to the library files')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_load_library)

parser_sub = subparsers.add_parser('update_tbdb', help='Pull the latest tbdb library and load', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.add_argument('--seqname',default='Chromosome',help='Mutation panel name')
parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
parser_sub.set_defaults(func=main_update_tbdb)

parser_sub = subparsers.add_parser('version', help='Output program version and exit', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser_sub.set_defaults(func=main_version)
args = parser.parse_args()
if vars(args)=={}:
    parser.print_help(sys.stderr)
else:
    args.func(args)