tb-profiler

#! /usr/bin/env python3
import sys
import pathogenprofiler as pp
import argparse
import json
import tbprofiler as tbp
import os
import csv

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"}
    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)
    if args.external_db:
        conf = get_conf_dict(args.external_db)
    else:
        conf = get_conf_dict(sys.base_prefix+"/share/tbprofiler/%s" % args.db)

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

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

    ### 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 ###
    files_prefix = args.dir+"/"+args.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 args.min_depth<args.missing_cov_threshold:
        sys.stderr.write(
            "\nThe value of --min_depth (%s) is lower than than --missing_cov_threshold (%s). "
            "This could result in variants called also being present missing_positions results.\n" % (args.min_depth,args.missing_cov_threshold)
        )
    ### 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, min_depth=args.min_depth
    )

    ### Reformat the results to TB-Profiler style ###
    results = tbp.reformat(results, conf, reporting_af=args.reporting_af)
    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}

    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=args.reporting_af)
    if args.csv:
        tbp.write_csv(results,conf,csv_output,extra_columns)

    ### Move files to respective directories ###
    if not args.bam:
        pp.run_cmd("mv %(dir)s/%(prefix)s.bam* %(dir)s/bam/" % vars(args))
        if not args.no_trim:
            pp.run_cmd("rm -f %s" % " ".join(fastq_obj.files))
    pp.run_cmd("mv -f %(dir)s/%(prefix)s*.vcf.gz* %(dir)s/vcf/" % vars(args))
    if run_delly and results["delly"]=="success":
        pp.run_cmd("mv -f %(dir)s/%(prefix)s.delly.bcf* %(dir)s/vcf/" % vars(args))

    ### Add meta data to results
    if args.meta:
        for row in csv.DictReader(open(args.meta)):
            if row["id"]==results["id"]:
                for col in row:
                    results["meta_"+col] = row[col]
    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("python parse_db.py --seqname %s" % args.seqname)
    pp.run_cmd("tb-profiler load_library tbdb")
    os.chdir("../")
    pp.log("Sucessfully updated TBDB")

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.run_cmd("bwa index %s" % 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.run_cmd("gatk CreateSequenceDictionary -R %s" % 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)
    if args.snps:
        with open("%s.lineage.snps.txt" % (args.prefix),"w") as OUT:
            barcode_bed = pp.load_bed(conf["barcode"],[1,2,3,4,5],3)
            for snp in sorted(list(barcode_bed.values()),key=lambda x:x[3]):
                tmp = mutations[snp[0]][int(snp[2])]
                lineage_reads = tmp[snp[4]] if snp[4] in tmp else 0
                not_lineage_reads = sum(tmp.values()) - lineage_reads
                fraction = lineage_reads/(lineage_reads+not_lineage_reads) if (sum(tmp.values())>0) else  "NA"
                OUT.write("%s\t%s\t%s\t%s\t%s\n" % (snp[2],snp[3],lineage_reads,not_lineage_reads,fraction))

    results = {}
    results["barcode"] = pp.barcode(mutations,conf["barcode"])
    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,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)

    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 = pp.get_random_file()
        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))
        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))


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

if __name__=="__main__":
    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")

    parser_sub = subparsers.add_parser('profile', help='Run whole profiling pipeline', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser_sub.add_argument('--platform','-m',choices=["illumina","nanopore"],default="illumina",help='NGS Platform used to generate data')
    parser_sub.add_argument('--read1','-1',help='First read file')
    parser_sub.add_argument('--read2','-2',help='Second read file')
    parser_sub.add_argument('--bam','-a',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 for all results generated')
    parser_sub.add_argument('--no_trim',action="store_true",help="Don't trim files using trimmomatic")
    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('--mapper',default="bwa", choices=["bwa","minimap2","bowtie2"],help="Mapping tool to use. If you are using nanopore data it will default to minimap2",type=str)
    parser_sub.add_argument('--caller',default="bcftools", choices=["bcftools","gatk","freebayes"],help="Variant calling tool to use.",type=str)
    parser_sub.add_argument('--calling_params',type=str,help='Override default parameters for variant calling')
    parser_sub.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')
    parser_sub.add_argument('--af',default=0.1,type=float,help='Minimum allele frequency to call variants')
    parser_sub.add_argument('--reporting_af',default=0.1,type=float,help='Minimum allele frequency to use variants for prediction')
    parser_sub.add_argument('--coverage_fraction_threshold',default=0,type=int,help='Cutoff used to calculate fraction of region covered by <= this value')
    parser_sub.add_argument('--missing_cov_threshold',default=10,type=int,help='Cutoff used to positions/codons in genes which are missing')
    parser_sub.add_argument('--threads','-t',default=1,help='Threads to use',type=int)
    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('--pdf',action="store_true",help="Add PDF output. This requires pdflatex to be installed")
    parser_sub.add_argument('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and pdf results")
    parser_sub.add_argument('--meta',default=None,type=str,help="Add meta data from a CSV file to the results. The CSV file must contain a column labelled \"id\" with the same value as the prefix argument")
    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('--no_flagstat',action="store_true",help="Don't collect flagstats")
    parser_sub.add_argument('--no_delly',action="store_true",help="Don't run delly")
    parser_sub.add_argument('--version', action='version', version="TBProfiler version %s" % tbp._VERSION)
    parser_sub.set_defaults(func=main_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('--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('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and pdf 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)

    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',help='VCF file')
    parser_sub.add_argument('prefix',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('--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('--add_columns',default=None,type=str,help="Add additional columns found in the mutation database to the text and pdf 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('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 pdf', 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('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('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","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('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)