<commit_before>#!/usr/bin/env python2.7
"""
vg_calleval.py: Compare vcfs with vcfeval.  Option to make freebayes calls to use as baseline.  Can
run on vg_mapeval.py output. 
"""
from __future__ import print_function
import argparse, sys, os, os.path, errno, random, subprocess, shutil, itertools, glob, tarfile
import doctest, re, json, collections, time, timeit
import logging, logging.handlers, SocketServer, struct, socket, threading
import string, math
import urlparse
import getpass
import pdb
import gzip
import logging
import copy

from math import ceil
from subprocess import Popen, PIPE

try:
    import numpy as np
    from sklearn.metrics import roc_auc_score, average_precision_score, r2_score, roc_curve
    have_sklearn = True
except:
    have_sklearn = False

import tsv
import vcf

from toil.common import Toil
from toil.job import Job
from toil.realtimeLogger import RealtimeLogger
from toil_vg.vg_common import *
from toil_vg.vg_call import chunked_call_parse_args, run_all_calling, run_merge_vcf
from toil_vg.vg_vcfeval import vcfeval_parse_args, run_extract_sample_truth_vcf, run_vcfeval, run_vcfeval_roc_plot, run_happy, run_sv_eval
from toil_vg.context import Context, run_write_info_to_outstore
from toil_vg.vg_construct import run_unzip_fasta
from toil_vg.vg_surject import run_surjecting

logger = logging.getLogger(__name__)

def calleval_subparser(parser):
    """
    Create a subparser for calleval.  Should pass in results of subparsers.add_parser()
    """

    # Add the Toil options so the job store is the first argument
    Job.Runner.addToilOptions(parser)
    
    # Add the out_store
    # TODO: do this at a higher level?
    # Or roll into Context?
    parser.add_argument('out_store',
                        help='output store.  All output written here. Path specified using same syntax as toil jobStore')

    parser.add_argument("--chroms", nargs='+', required=True,
                        help="Name(s) of reference path in graph(s) (separated by space).")
    # todo: move to chunked_call_parse_args and share with toil-vg run
    parser.add_argument("--gams", nargs='+', type=make_url,
                        help="GAMs to call.  Each GAM treated as separate input (and must contain all chroms)")
    parser.add_argument("--sample_name", type=str, required=True,
                        help="sample name (ex NA12878)")
    parser.add_argument("--gam_index_cores", type=int,
                        help="number of threads used for gam indexing")

    # Add common options shared with everybody
    add_common_vg_parse_args(parser)

    # Add common call options shared with toil_vg pipeline
    chunked_call_parse_args(parser)
    
    # Add common vcfeval options shared with toil_vg pipeline
    vcfeval_parse_args(parser)

    # Add common calleval options shared with toil_vg pipeline
    calleval_parse_args(parser)

    # Add common docker options shared with toil_vg pipeline
    add_container_tool_parse_args(parser)
    
def calleval_parse_args(parser):
    """
    Add the calleval options to the given argparse parser.
    """
    parser.add_argument('--gam_names', nargs='+',
                        help='names of vg runs (corresponds to gams and xg_paths)')
    parser.add_argument('--xg_paths', nargs='+', type=make_url,
                        help='xg indexes for the different graphs')
    parser.add_argument('--freebayes', action='store_true',
                        help='run freebayes as a baseline')
    parser.add_argument("--caller_fasta", type=make_url,
                        help="Use this FASTA instead of the vcfeval fasta for Freebayes. Maybe be gzipped")
    parser.add_argument('--platypus', action='store_true',
                        help='run platypus as a baseline')
    parser.add_argument('--bam_names', nargs='+', default=[],
                        help='names of bwa runs (corresponds to bams)')
    parser.add_argument('--bams', nargs='+', type=make_url, default=[],
                        help='bam inputs for freebayes')
    parser.add_argument('--filter_opts_gt',
                        help='override filter-opts for genotype only')
    parser.add_argument('--clip_only', action='store_true',
                        help='only compute accuracy clipped to --vcfeval_bed_regions')
    parser.add_argument('--plot_sets', nargs='+', default=[],
                        help='comma-separated lists of condition-tagged BAM/GAM names (primary-mp-pe, etc.) with colon-separated title prefixes')
    parser.add_argument('--call', action='store_true',
                        help='run vg call (even if --genotype used)')    
    parser.add_argument("--surject", action="store_true",
                        help="surject GAMs to BAMs, adding the latter to the comparison")
    parser.add_argument("--interleaved", action="store_true", default=False,
                        help="assume GAM files are interleaved when surjecting")
        
def validate_calleval_options(options):
    """
    Throw an error if an invalid combination of options has been selected.
    """
    if options.gams or options.gam_names or options.xg_paths:        
        require(options.gams and options.gam_names and options.xg_paths and
                len(options.gam_names) == len(options.xg_paths) == len(options.gams),
                '--gam_names, --xg_paths, --gams must all contain same number of elements')
        require(options.call or options.genotype or options.surject,
                '--call and/or --genotype and/or --surject required with --gams')
    if options.freebayes or options.platypus:
        require(options.bams or options.surject, '--bams and/or --surject needed with --freebayes or --platypus')
    if options.bams or options.bam_names:
        require(options.bams and options.bam_names and len(options.bams) == len(options.bam_names),
                '--bams and --bam_names must be same length')
    require(options.vcfeval_baseline, '--vcfeval_baseline required')
    require(options.vcfeval_fasta, '--vcfeval_fasta required')
    if options.call or options.genotype or options.surject:
        require(options.gams, '--gams must be given with --call, --genotype, and --surject')
    require(options.vcfeval_bed_regions is not None or not options.clip_only,
            '--vcfeval_bed_regions must be given with --clip_only')
    if options.surject or options.bams:
        require(options.freebayes or options.platypus,
                '--freebayes and/or --platypus must be used with --bams and --surject')

def run_bam_index(job, context, bam_file_id, bam_name):
    """
    sort and index a bam.  return the sorted bam and its idx
    """
    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # download the input
    bam_path = os.path.join(work_dir, bam_name + '.bam')
    job.fileStore.readGlobalFile(bam_file_id, bam_path)

    sort_bam_path = os.path.join(work_dir, 'sort.bam')
    sort_cmd = ['samtools', 'sort', os.path.basename(bam_path), '-o',
                os.path.basename(sort_bam_path), '-O', 'BAM', '-@', max(0, job.cores - 1)]
    context.runner.call(job, sort_cmd, work_dir=work_dir)
    bam_index_cmd = ['samtools', 'index', os.path.basename(sort_bam_path)]
    context.runner.call(job, bam_index_cmd, work_dir=work_dir)

    out_bam_id = context.write_intermediate_file(job, sort_bam_path)
    out_idx_id = context.write_intermediate_file(job, sort_bam_path + '.bai')
    return out_bam_id, out_idx_id
    
def run_all_bam_caller(job, context, fasta_file_id, bam_file_id, bam_idx_id,
                       sample_name, chroms, offsets, out_name, bam_caller,
                       bam_caller_opts = []):
    """
    run freebayes or platypus on a set of chromosomal regions.  this is done by sending each region to a 
    child job and farming off the entire input to each (ie not splitting the input)
    """
    # to encapsulate everything under this job
    child_job = Job()
    job.addChild(child_job)

    fb_vcf_ids = []
    fb_tbi_ids = []
    fb_timers = []
    assert chroms
    if not offsets:
        offsets = [None] * len(chroms)
    for chrom, offset in zip(chroms, offsets):
        fb_job = child_job.addChildJobFn(run_bam_caller, context, fasta_file_id, bam_file_id, bam_idx_id,
                                         sample_name, chrom, offset, out_name, bam_caller, bam_caller_opts,
                                         cores=context.config.calling_cores,
                                         memory=context.config.calling_mem,
                                         disk=context.config.calling_disk)
        fb_vcf_ids.append(fb_job.rv(0))
        fb_tbi_ids.append(fb_job.rv(1))
        fb_timers.append([fb_job.rv(2)])

    merge_vcf_job = child_job.addFollowOnJobFn(run_merge_vcf, context, out_name, zip(fb_vcf_ids, fb_tbi_ids), fb_timers)
    return merge_vcf_job.rv()
    
def run_bam_caller(job, context, fasta_file_id, bam_file_id, bam_idx_id,
                   sample_name, chrom, offset, out_name, bam_caller, bam_caller_opts):
    """
    run freebayes or platypus to make a vcf
    """

    assert bam_caller in ['freebayes', 'platypus']
    
    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # download the input
    fasta_path = os.path.join(work_dir, 'ref.fa')
    bam_path = os.path.join(work_dir, 'alignment.bam')
    bam_idx_path = bam_path + '.bai'
    job.fileStore.readGlobalFile(fasta_file_id, fasta_path)
    job.fileStore.readGlobalFile(bam_file_id, bam_path)
    job.fileStore.readGlobalFile(bam_idx_id, bam_idx_path)

    # output
    vcf_path = os.path.join(work_dir, '{}-raw.vcf'.format(out_name))
    
    if bam_caller == 'freebayes':
        fb_cmd = ['freebayes', '-f', os.path.basename(fasta_path), os.path.basename(bam_path)]
        if bam_caller_opts:
            fb_cmd += bam_caller_opts
        if chrom:
            fb_cmd += ['-r', chrom]
        timer = TimeTracker('freebayes')
        with open(vcf_path, 'w') as out_vcf:
            context.runner.call(job, fb_cmd, work_dir=work_dir, outfile=out_vcf)
        timer.stop()
            
    elif bam_caller == 'platypus':
        context.runner.call(job, ['samtools', 'faidx', 'ref.fa'], work_dir=work_dir)
        plat_cmd = ['Platypus.py', 'callVariants', '--refFile', os.path.basename(fasta_path),
                    '--bamFiles', os.path.basename(bam_path), '-o', os.path.basename(vcf_path)]
        if bam_caller_opts:
            plat_cmd += bam_caller_opts
        if chrom:
            plat_cmd += ['--regions', chrom]
        timer = TimeTracker('Platypus')
        context.runner.call(job, plat_cmd, work_dir=work_dir)
        timer.stop()
    else:
        assert False

    context.write_intermediate_file(job, vcf_path)

    vcf_fix_path = os.path.join(work_dir, '{}.vcf'.format(out_name))
    
    # apply offset and sample name
    vcf_reader = vcf.Reader(open(vcf_path))
    if sample_name:
        # Freebayes always outputs "unknown" for the sample if the BAM
        # doesn't specify a name, and can't be convinced to do otherwise:
        # https://github.com/ekg/freebayes/issues/471
        # So we hack the VCFReader to think the sample names are what we want them to be
        assert(len(vcf_reader.samples) == 1)
        # TODO: asserting that sample 0 is named unknown can fail. Why?
        RealtimeLogger.info('Correcting Freebayes samples {} to [{}]'.format(vcf_reader.samples, sample_name))
        vcf_reader.samples = [sample_name]
        # Rebuild the secret sample index
        vcf_reader._sample_indexes = {sample_name: 0}
    vcf_writer = vcf.Writer(open(vcf_fix_path, 'w'), vcf_reader)
    
    # We insist that Freebayes produce at least one record, or we will abort and dump files.
    have_records = False
    for record in vcf_reader:
        have_records = True
        if offset:
            record.POS += int(offset)
        vcf_writer.write_record(record)
    vcf_writer.flush()
    vcf_writer.close()
    
    if not have_records:
        # This is a problem. Stop here for debugging.
        RealtimeLogger.error('{} produced no calls. Dumping files...'.format(bam_caller))
        for dump_path in [fasta_path, bam_path, bam_idx_path]:
            if dump_path and os.path.isfile(dump_path):
                context.write_output_file(job, dump_path)
        raise RuntimeError('{} produced no calls'.format(bam_caller))

    context.runner.call(job, ['bgzip', os.path.basename(vcf_fix_path)], work_dir = work_dir)
    context.runner.call(job, ['tabix', '-p', 'vcf', os.path.basename(vcf_fix_path) + '.gz'], work_dir = work_dir)

    return (context.write_output_file(job, vcf_fix_path + '.gz'),
            context.write_output_file(job, vcf_fix_path + '.gz.tbi'),
            timer)

def run_calleval_plots(job, context, names, eval_results_dict, plot_sets):
    """
    
    Make and output calleval ROC plots.
    
    Takes a "names" list of all conditions. Condition names in the list (or in
    plot_sets) do not include an "-unclipped" tag; both clipped and unclipped
    plots are made if the data is available.
    
    Takes a nested dict by condition name, then clipping status ("clipped",
    "unclipped"), and then variant type ("snp", "non_snp", "weighted").
    Eventual entries are to ROC data file ids (.tsv.gz).
    
    plot_sets is a data structure of collections of conditions to plot against
    each other, as produced by parse_plot_sets.
    
    Returns a list of created plot file IDs.
    
    """
    
    # make some roc plots
    roc_plot_ids = []
    for roc_type in ['snp', 'non_snp', 'weighted']:
        # For each type of ROC to plot
    
        for mode in ['unclipped', 'clipped']:
            # For each clipping mode

            # What kind of ROC plot is it?
            # It should be this unless there is a clipped mode for this ROC and we are the unclipped one.
            roc_kind = roc_type
            if mode == 'unclipped' and True in [eval_results_dict.has_key('clipped') for eval_result in eval_results_dict.itervalues()]:
                # We are the unclipped mode and there will be a clipped mode
                roc_kind += '-unclipped'
            
            # Get all the eval results for this mode, by condition name
            mode_results = {name: eval_result.get(mode) for name, eval_result in eval_results_dict.iteritems()}
            
            if None in mode_results.itervalues():
                # We can't do this mode since it wasn't run
                continue
                
            # Extract out all the stats file IDs for this ROC type, by
            # condition name. This dict is implicitly for the clipped or
            # unclipped mode, depending on which mode we are doing. The
            # condition name keys don't have clipping tags.
            roc_table_ids = {name: result.get(roc_type) for name, result in mode_results.iteritems()}
            
            for subset_number, plot_set in enumerate(plot_sets):
                # For each plot set specifier
                
                # Unpack plot_set
                plot_title, plot_conditions = plot_set
                
                if plot_conditions is None:
                    # Use all the conditions instead
                    plot_conditions = names
                    
                # Add the ROC plot kind (snp, clipped, etc.) to the plot title
                if plot_title is None:
                    plot_title = roc_kind
                else:
                    plot_title = '{} ({})'.format(plot_title, roc_kind)
                
                for name in plot_conditions:
                    if name not in names:
                        # Complain if any of the names in the subset aren't conditions that actually ran
                        message = 'Condition {} not found in list of available conditions {}'.format(name, names)
                        RealtimeLogger.error(message)
                        raise RuntimeError(message)
                    if not roc_table_ids.has_key(name):
                        # Complain if any of the names in the subset lacks ROC data
                        message = 'Condition {} has no ROC data; data only available for {}'.format(name, roc_table_ids.keys())
                        RealtimeLogger.error(message)
                        raise RuntimeError(message)

                # Make a list of roc table IDs for the subset, in plot_conditions order
                subset_ids = [roc_table_ids[name] for name in plot_conditions]
                
                # Make the plot
                roc_plot_ids.append(job.addChildJobFn(run_vcfeval_roc_plot, context, subset_ids, names=plot_conditions,
                                                      kind=roc_kind, number=subset_number, title=plot_title).rv())
                                                      
    return roc_plot_ids
    

def run_calleval_results(job, context, names, vcf_tbi_pairs, eval_results_dict, happy_results_dict, sveval_results_dict,
                         timing_results, plot_sets):
    """
    
    output the calleval results
    
    Requires that, if any result in eval_results has clipped results, all
    results have clipped results, and similarly for unclipped results.

    plot_sets is a data structure of collections of conditions to plot against
    each other, as produced by parse_plot_sets.
    
    """

    RealtimeLogger.info('Handling results for conditions {} in sets {}'.format(names, plot_sets))

    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # make a simple tsv
    stats_path = os.path.join(work_dir, 'calleval_stats.tsv')
    with open(stats_path, 'w') as stats_file:
        for name in names:
            eval_result = eval_results_dict[name]
            happy_result = happy_results_dict[name]
            # Find the best result (clipped if present, unclipped if not).
            best_result = eval_result.get("clipped", eval_result.get("unclipped", None))

            # Same for the happy results
            best_happy_result = happy_result.get("clipped", happy_result.get("unclipped", None))
            if best_happy_result is not None:
                happy_snp_f1 = best_happy_result['parsed_summary']['SNP']['METRIC.F1_Score']
                happy_indel_f1 =  best_happy_result['parsed_summary']['INDEL']['METRIC.F1_Score']
            else:
                happy_snp_f1, happy_indel_f1 = -1, -1

            # Same for the sveval results
            sveval_result = sveval_results_dict[name]
            best_sveval_result = sveval_result.get("clipped", sveval_result.get("unclipped", None))
            sveval_f1 = best_sveval_result['F1'] if best_sveval_result is not None else -1
                
            # Output the F1 scores
            stats_file.write('{}\t{}\t{}\t{}\t{}\n'.format(name, best_result['f1'], happy_snp_f1, happy_indel_f1, sveval_f1))

    # Make the roc plots
    roc_plot_job = job.addChildJobFn(run_calleval_plots, context, names, eval_results_dict, plot_sets)
    roc_plot_ids = roc_plot_job.rv()

    # write some times
    times_path = os.path.join(work_dir, 'call_times.tsv')
    
    # organize our expected labels a bit
    calling_labels = ['call', 'genotype', 'freebayes', 'platypus']
    augmenting_labels = ['call-filter-augment', 'call-augment', 'call-filter']
    all_labels = set()
    for timer in timing_results:
        for label in timer.names():
            all_labels.add(label)
    other_labels = all_labels - set(calling_labels + augmenting_labels)
    with open(times_path, 'w') as times_file:
        # write the header<commit_msg>Remove sveval_f1 from output again because the tests prohibit it<commit_after>#!/usr/bin/env python2.7
"""
vg_calleval.py: Compare vcfs with vcfeval.  Option to make freebayes calls to use as baseline.  Can
run on vg_mapeval.py output. 
"""
from __future__ import print_function
import argparse, sys, os, os.path, errno, random, subprocess, shutil, itertools, glob, tarfile
import doctest, re, json, collections, time, timeit
import logging, logging.handlers, SocketServer, struct, socket, threading
import string, math
import urlparse
import getpass
import pdb
import gzip
import logging
import copy

from math import ceil
from subprocess import Popen, PIPE

try:
    import numpy as np
    from sklearn.metrics import roc_auc_score, average_precision_score, r2_score, roc_curve
    have_sklearn = True
except:
    have_sklearn = False

import tsv
import vcf

from toil.common import Toil
from toil.job import Job
from toil.realtimeLogger import RealtimeLogger
from toil_vg.vg_common import *
from toil_vg.vg_call import chunked_call_parse_args, run_all_calling, run_merge_vcf
from toil_vg.vg_vcfeval import vcfeval_parse_args, run_extract_sample_truth_vcf, run_vcfeval, run_vcfeval_roc_plot, run_happy, run_sv_eval
from toil_vg.context import Context, run_write_info_to_outstore
from toil_vg.vg_construct import run_unzip_fasta
from toil_vg.vg_surject import run_surjecting

logger = logging.getLogger(__name__)

def calleval_subparser(parser):
    """
    Create a subparser for calleval.  Should pass in results of subparsers.add_parser()
    """

    # Add the Toil options so the job store is the first argument
    Job.Runner.addToilOptions(parser)
    
    # Add the out_store
    # TODO: do this at a higher level?
    # Or roll into Context?
    parser.add_argument('out_store',
                        help='output store.  All output written here. Path specified using same syntax as toil jobStore')

    parser.add_argument("--chroms", nargs='+', required=True,
                        help="Name(s) of reference path in graph(s) (separated by space).")
    # todo: move to chunked_call_parse_args and share with toil-vg run
    parser.add_argument("--gams", nargs='+', type=make_url,
                        help="GAMs to call.  Each GAM treated as separate input (and must contain all chroms)")
    parser.add_argument("--sample_name", type=str, required=True,
                        help="sample name (ex NA12878)")
    parser.add_argument("--gam_index_cores", type=int,
                        help="number of threads used for gam indexing")

    # Add common options shared with everybody
    add_common_vg_parse_args(parser)

    # Add common call options shared with toil_vg pipeline
    chunked_call_parse_args(parser)
    
    # Add common vcfeval options shared with toil_vg pipeline
    vcfeval_parse_args(parser)

    # Add common calleval options shared with toil_vg pipeline
    calleval_parse_args(parser)

    # Add common docker options shared with toil_vg pipeline
    add_container_tool_parse_args(parser)
    
def calleval_parse_args(parser):
    """
    Add the calleval options to the given argparse parser.
    """
    parser.add_argument('--gam_names', nargs='+',
                        help='names of vg runs (corresponds to gams and xg_paths)')
    parser.add_argument('--xg_paths', nargs='+', type=make_url,
                        help='xg indexes for the different graphs')
    parser.add_argument('--freebayes', action='store_true',
                        help='run freebayes as a baseline')
    parser.add_argument("--caller_fasta", type=make_url,
                        help="Use this FASTA instead of the vcfeval fasta for Freebayes. Maybe be gzipped")
    parser.add_argument('--platypus', action='store_true',
                        help='run platypus as a baseline')
    parser.add_argument('--bam_names', nargs='+', default=[],
                        help='names of bwa runs (corresponds to bams)')
    parser.add_argument('--bams', nargs='+', type=make_url, default=[],
                        help='bam inputs for freebayes')
    parser.add_argument('--filter_opts_gt',
                        help='override filter-opts for genotype only')
    parser.add_argument('--clip_only', action='store_true',
                        help='only compute accuracy clipped to --vcfeval_bed_regions')
    parser.add_argument('--plot_sets', nargs='+', default=[],
                        help='comma-separated lists of condition-tagged BAM/GAM names (primary-mp-pe, etc.) with colon-separated title prefixes')
    parser.add_argument('--call', action='store_true',
                        help='run vg call (even if --genotype used)')    
    parser.add_argument("--surject", action="store_true",
                        help="surject GAMs to BAMs, adding the latter to the comparison")
    parser.add_argument("--interleaved", action="store_true", default=False,
                        help="assume GAM files are interleaved when surjecting")
        
def validate_calleval_options(options):
    """
    Throw an error if an invalid combination of options has been selected.
    """
    if options.gams or options.gam_names or options.xg_paths:        
        require(options.gams and options.gam_names and options.xg_paths and
                len(options.gam_names) == len(options.xg_paths) == len(options.gams),
                '--gam_names, --xg_paths, --gams must all contain same number of elements')
        require(options.call or options.genotype or options.surject,
                '--call and/or --genotype and/or --surject required with --gams')
    if options.freebayes or options.platypus:
        require(options.bams or options.surject, '--bams and/or --surject needed with --freebayes or --platypus')
    if options.bams or options.bam_names:
        require(options.bams and options.bam_names and len(options.bams) == len(options.bam_names),
                '--bams and --bam_names must be same length')
    require(options.vcfeval_baseline, '--vcfeval_baseline required')
    require(options.vcfeval_fasta, '--vcfeval_fasta required')
    if options.call or options.genotype or options.surject:
        require(options.gams, '--gams must be given with --call, --genotype, and --surject')
    require(options.vcfeval_bed_regions is not None or not options.clip_only,
            '--vcfeval_bed_regions must be given with --clip_only')
    if options.surject or options.bams:
        require(options.freebayes or options.platypus,
                '--freebayes and/or --platypus must be used with --bams and --surject')

def run_bam_index(job, context, bam_file_id, bam_name):
    """
    sort and index a bam.  return the sorted bam and its idx
    """
    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # download the input
    bam_path = os.path.join(work_dir, bam_name + '.bam')
    job.fileStore.readGlobalFile(bam_file_id, bam_path)

    sort_bam_path = os.path.join(work_dir, 'sort.bam')
    sort_cmd = ['samtools', 'sort', os.path.basename(bam_path), '-o',
                os.path.basename(sort_bam_path), '-O', 'BAM', '-@', max(0, job.cores - 1)]
    context.runner.call(job, sort_cmd, work_dir=work_dir)
    bam_index_cmd = ['samtools', 'index', os.path.basename(sort_bam_path)]
    context.runner.call(job, bam_index_cmd, work_dir=work_dir)

    out_bam_id = context.write_intermediate_file(job, sort_bam_path)
    out_idx_id = context.write_intermediate_file(job, sort_bam_path + '.bai')
    return out_bam_id, out_idx_id
    
def run_all_bam_caller(job, context, fasta_file_id, bam_file_id, bam_idx_id,
                       sample_name, chroms, offsets, out_name, bam_caller,
                       bam_caller_opts = []):
    """
    run freebayes or platypus on a set of chromosomal regions.  this is done by sending each region to a 
    child job and farming off the entire input to each (ie not splitting the input)
    """
    # to encapsulate everything under this job
    child_job = Job()
    job.addChild(child_job)

    fb_vcf_ids = []
    fb_tbi_ids = []
    fb_timers = []
    assert chroms
    if not offsets:
        offsets = [None] * len(chroms)
    for chrom, offset in zip(chroms, offsets):
        fb_job = child_job.addChildJobFn(run_bam_caller, context, fasta_file_id, bam_file_id, bam_idx_id,
                                         sample_name, chrom, offset, out_name, bam_caller, bam_caller_opts,
                                         cores=context.config.calling_cores,
                                         memory=context.config.calling_mem,
                                         disk=context.config.calling_disk)
        fb_vcf_ids.append(fb_job.rv(0))
        fb_tbi_ids.append(fb_job.rv(1))
        fb_timers.append([fb_job.rv(2)])

    merge_vcf_job = child_job.addFollowOnJobFn(run_merge_vcf, context, out_name, zip(fb_vcf_ids, fb_tbi_ids), fb_timers)
    return merge_vcf_job.rv()
    
def run_bam_caller(job, context, fasta_file_id, bam_file_id, bam_idx_id,
                   sample_name, chrom, offset, out_name, bam_caller, bam_caller_opts):
    """
    run freebayes or platypus to make a vcf
    """

    assert bam_caller in ['freebayes', 'platypus']
    
    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # download the input
    fasta_path = os.path.join(work_dir, 'ref.fa')
    bam_path = os.path.join(work_dir, 'alignment.bam')
    bam_idx_path = bam_path + '.bai'
    job.fileStore.readGlobalFile(fasta_file_id, fasta_path)
    job.fileStore.readGlobalFile(bam_file_id, bam_path)
    job.fileStore.readGlobalFile(bam_idx_id, bam_idx_path)

    # output
    vcf_path = os.path.join(work_dir, '{}-raw.vcf'.format(out_name))
    
    if bam_caller == 'freebayes':
        fb_cmd = ['freebayes', '-f', os.path.basename(fasta_path), os.path.basename(bam_path)]
        if bam_caller_opts:
            fb_cmd += bam_caller_opts
        if chrom:
            fb_cmd += ['-r', chrom]
        timer = TimeTracker('freebayes')
        with open(vcf_path, 'w') as out_vcf:
            context.runner.call(job, fb_cmd, work_dir=work_dir, outfile=out_vcf)
        timer.stop()
            
    elif bam_caller == 'platypus':
        context.runner.call(job, ['samtools', 'faidx', 'ref.fa'], work_dir=work_dir)
        plat_cmd = ['Platypus.py', 'callVariants', '--refFile', os.path.basename(fasta_path),
                    '--bamFiles', os.path.basename(bam_path), '-o', os.path.basename(vcf_path)]
        if bam_caller_opts:
            plat_cmd += bam_caller_opts
        if chrom:
            plat_cmd += ['--regions', chrom]
        timer = TimeTracker('Platypus')
        context.runner.call(job, plat_cmd, work_dir=work_dir)
        timer.stop()
    else:
        assert False

    context.write_intermediate_file(job, vcf_path)

    vcf_fix_path = os.path.join(work_dir, '{}.vcf'.format(out_name))
    
    # apply offset and sample name
    vcf_reader = vcf.Reader(open(vcf_path))
    if sample_name:
        # Freebayes always outputs "unknown" for the sample if the BAM
        # doesn't specify a name, and can't be convinced to do otherwise:
        # https://github.com/ekg/freebayes/issues/471
        # So we hack the VCFReader to think the sample names are what we want them to be
        assert(len(vcf_reader.samples) == 1)
        # TODO: asserting that sample 0 is named unknown can fail. Why?
        RealtimeLogger.info('Correcting Freebayes samples {} to [{}]'.format(vcf_reader.samples, sample_name))
        vcf_reader.samples = [sample_name]
        # Rebuild the secret sample index
        vcf_reader._sample_indexes = {sample_name: 0}
    vcf_writer = vcf.Writer(open(vcf_fix_path, 'w'), vcf_reader)
    
    # We insist that Freebayes produce at least one record, or we will abort and dump files.
    have_records = False
    for record in vcf_reader:
        have_records = True
        if offset:
            record.POS += int(offset)
        vcf_writer.write_record(record)
    vcf_writer.flush()
    vcf_writer.close()
    
    if not have_records:
        # This is a problem. Stop here for debugging.
        RealtimeLogger.error('{} produced no calls. Dumping files...'.format(bam_caller))
        for dump_path in [fasta_path, bam_path, bam_idx_path]:
            if dump_path and os.path.isfile(dump_path):
                context.write_output_file(job, dump_path)
        raise RuntimeError('{} produced no calls'.format(bam_caller))

    context.runner.call(job, ['bgzip', os.path.basename(vcf_fix_path)], work_dir = work_dir)
    context.runner.call(job, ['tabix', '-p', 'vcf', os.path.basename(vcf_fix_path) + '.gz'], work_dir = work_dir)

    return (context.write_output_file(job, vcf_fix_path + '.gz'),
            context.write_output_file(job, vcf_fix_path + '.gz.tbi'),
            timer)

def run_calleval_plots(job, context, names, eval_results_dict, plot_sets):
    """
    
    Make and output calleval ROC plots.
    
    Takes a "names" list of all conditions. Condition names in the list (or in
    plot_sets) do not include an "-unclipped" tag; both clipped and unclipped
    plots are made if the data is available.
    
    Takes a nested dict by condition name, then clipping status ("clipped",
    "unclipped"), and then variant type ("snp", "non_snp", "weighted").
    Eventual entries are to ROC data file ids (.tsv.gz).
    
    plot_sets is a data structure of collections of conditions to plot against
    each other, as produced by parse_plot_sets.
    
    Returns a list of created plot file IDs.
    
    """
    
    # make some roc plots
    roc_plot_ids = []
    for roc_type in ['snp', 'non_snp', 'weighted']:
        # For each type of ROC to plot
    
        for mode in ['unclipped', 'clipped']:
            # For each clipping mode

            # What kind of ROC plot is it?
            # It should be this unless there is a clipped mode for this ROC and we are the unclipped one.
            roc_kind = roc_type
            if mode == 'unclipped' and True in [eval_results_dict.has_key('clipped') for eval_result in eval_results_dict.itervalues()]:
                # We are the unclipped mode and there will be a clipped mode
                roc_kind += '-unclipped'
            
            # Get all the eval results for this mode, by condition name
            mode_results = {name: eval_result.get(mode) for name, eval_result in eval_results_dict.iteritems()}
            
            if None in mode_results.itervalues():
                # We can't do this mode since it wasn't run
                continue
                
            # Extract out all the stats file IDs for this ROC type, by
            # condition name. This dict is implicitly for the clipped or
            # unclipped mode, depending on which mode we are doing. The
            # condition name keys don't have clipping tags.
            roc_table_ids = {name: result.get(roc_type) for name, result in mode_results.iteritems()}
            
            for subset_number, plot_set in enumerate(plot_sets):
                # For each plot set specifier
                
                # Unpack plot_set
                plot_title, plot_conditions = plot_set
                
                if plot_conditions is None:
                    # Use all the conditions instead
                    plot_conditions = names
                    
                # Add the ROC plot kind (snp, clipped, etc.) to the plot title
                if plot_title is None:
                    plot_title = roc_kind
                else:
                    plot_title = '{} ({})'.format(plot_title, roc_kind)
                
                for name in plot_conditions:
                    if name not in names:
                        # Complain if any of the names in the subset aren't conditions that actually ran
                        message = 'Condition {} not found in list of available conditions {}'.format(name, names)
                        RealtimeLogger.error(message)
                        raise RuntimeError(message)
                    if not roc_table_ids.has_key(name):
                        # Complain if any of the names in the subset lacks ROC data
                        message = 'Condition {} has no ROC data; data only available for {}'.format(name, roc_table_ids.keys())
                        RealtimeLogger.error(message)
                        raise RuntimeError(message)

                # Make a list of roc table IDs for the subset, in plot_conditions order
                subset_ids = [roc_table_ids[name] for name in plot_conditions]
                
                # Make the plot
                roc_plot_ids.append(job.addChildJobFn(run_vcfeval_roc_plot, context, subset_ids, names=plot_conditions,
                                                      kind=roc_kind, number=subset_number, title=plot_title).rv())
                                                      
    return roc_plot_ids
    

def run_calleval_results(job, context, names, vcf_tbi_pairs, eval_results_dict, happy_results_dict, sveval_results_dict,
                         timing_results, plot_sets):
    """
    
    output the calleval results
    
    Requires that, if any result in eval_results has clipped results, all
    results have clipped results, and similarly for unclipped results.

    plot_sets is a data structure of collections of conditions to plot against
    each other, as produced by parse_plot_sets.
    
    """

    RealtimeLogger.info('Handling results for conditions {} in sets {}'.format(names, plot_sets))

    # make a local work directory
    work_dir = job.fileStore.getLocalTempDir()

    # make a simple tsv
    stats_path = os.path.join(work_dir, 'calleval_stats.tsv')
    with open(stats_path, 'w') as stats_file:
        for name in names:
            eval_result = eval_results_dict[name]
            happy_result = happy_results_dict[name]
            # Find the best result (clipped if present, unclipped if not).
            best_result = eval_result.get("clipped", eval_result.get("unclipped", None))

            # Same for the happy results
            best_happy_result = happy_result.get("clipped", happy_result.get("unclipped", None))
            if best_happy_result is not None:
                happy_snp_f1 = best_happy_result['parsed_summary']['SNP']['METRIC.F1_Score']
                happy_indel_f1 =  best_happy_result['parsed_summary']['INDEL']['METRIC.F1_Score']
            else:
                happy_snp_f1, happy_indel_f1 = -1, -1

            # Same for the sveval results
            sveval_result = sveval_results_dict[name]
            best_sveval_result = sveval_result.get("clipped", sveval_result.get("unclipped", None))
            sveval_f1 = best_sveval_result['F1'] if best_sveval_result is not None else -1
                
            # Output the F1 scores
            # TODO: We amy want sveval_f1 in here, but the tests prohibit tis presence
            stats_file.write('{}\t{}\t{}\t{}\n'.format(name, best_result['f1'], happy_snp_f1, happy_indel_f1))

    # Make the roc plots
    roc_plot_job = job.addChildJobFn(run_calleval_plots, context, names, eval_results_dict, plot_sets)
    roc_plot_ids = roc_plot_job.rv()

    # write some times
    times_path = os.path.join(work_dir, 'call_times.tsv')
    
    # organize our expected labels a bit
    calling_labels = ['call', 'genotype', 'freebayes', 'platypus']
    augmenting_labels = ['call-filter-augment', 'call-augment', 'call-filter']
    all_labels = set()
    for timer in timing_results:
        for label in timer.names():
            all_labels.add(label)
    other_labels = all_labels - set(calling_labels + augmenting_labels)
    with open(times_path, 'w') as times_file:
        # write the header
{'fraction_non_alphanumeric': 0.09225428839283452, 'fraction_numerical': 0.001798971292472192, 'mean_word_length': 2.5120128136679125, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5494957', 'n_tokens_mistral': 12434, 'n_tokens_neox': 11958, 'n_words': 3294}