Run cortex_var run_calls.pl, producing a VCF variant file.
973393
def _run_cortex(fastq, indexes, params, out_base, dirs, config):
    
    print(out_base)
    fastaq_index = "{0}.fastaq_index".format(out_base)
    se_fastq_index = "{0}.se_fastq".format(out_base)
    pe_fastq_index = "{0}.pe_fastq".format(out_base)
    reffasta_index = "{0}.list_ref_fasta".format(out_base)
    with open(se_fastq_index, "w") as out_handle:
        out_handle.write(fastq + "\n")
    with open(pe_fastq_index, "w") as out_handle:
        out_handle.write("")
    with open(fastaq_index, "w") as out_handle:
        out_handle.write("{0}\t{1}\t{2}\t{2}\n".format(params["sample"], se_fastq_index,
                                                       pe_fastq_index))
    with open(reffasta_index, "w") as out_handle:
        for x in indexes["fasta"]:
            out_handle.write(x + "\n")
    os.environ["PERL5LIB"] = "{0}:{1}:{2}".format(
        os.path.join(dirs["cortex"], "scripts/calling"),
        os.path.join(dirs["cortex"], "scripts/analyse_variants/bioinf-perl/lib"),
        os.environ.get("PERL5LIB", ""))
    kmers = sorted(params["kmers"])
    kmer_info = ["--first_kmer", str(kmers[0])]
    if len(kmers) > 1:
        kmer_info += ["--last_kmer", str(kmers[-1]),
                      "--kmer_step", str(kmers[1] - kmers[0])]
    subprocess.check_call(["perl", os.path.join(dirs["cortex"], "scripts", "calling", "run_calls.pl"),
                           "--fastaq_index", fastaq_index,
                           "--auto_cleaning", "yes", "--bc", "yes", "--pd", "yes",
                           "--outdir", os.path.dirname(out_base), "--outvcf", os.path.basename(out_base),
                           "--ploidy", str(config["algorithm"].get("ploidy", 2)),
                           "--stampy_hash", indexes["stampy"],
                           "--stampy_bin", os.path.join(dirs["stampy"], "stampy.py"),
                           "--refbindir", os.path.dirname(indexes["cortex"][0]),
                           "--list_ref_fasta",  reffasta_index,
                           "--genome_size", str(params["genome_size"]),
                           "--max_read_len", "30000",
                           #"--max_var_len", "4000",
                           "--format", "FASTQ", "--qthresh", "5", "--do_union", "yes",
                           "--mem_height", "17", "--mem_width", "100",
                           "--ref", "CoordinatesAndInCalling", "--workflow", "independent",
                           "--vcftools_dir", dirs["vcftools"],
                           "--logfile", "{0}.logfile,f".format(out_base)]
                          + kmer_info)
    final = glob.glob(os.path.join(os.path.dirname(out_base), "vcfs",
                                   "{0}*FINALcombined_BC*decomp.vcf".format(os.path.basename(out_base))))
    # No calls, need to setup an empty file
    if len(final) != 1:
        print("Did not find output VCF file for {0}".format(out_base))
        return None
    else:
        return final[0]
588736
def run_filter(vrn_file, align_bam, ref_file, data, items):
    
    if not should_filter(items) or not vcfutils.vcf_has_variants(vrn_file):
        return data
    else:
        raw_file = "%s-damage.vcf" % utils.splitext_plus(vrn_file)[0]
        out_plot_files = ["%s%s" % (utils.splitext_plus(raw_file)[0], ext)
                          for ext in ["_seq_bias_simplified.pdf", "_pcr_bias_simplified.pdf"]]
        if not utils.file_uptodate(raw_file, vrn_file) and not utils.file_uptodate(raw_file + ".gz", vrn_file):
            with file_transaction(items[0], raw_file) as tx_out_file:
                # Does not apply --qcSummary plotting due to slow runtimes
                cmd = ["dkfzbiasfilter.py", "--filterCycles", "1", "--passOnly",
                       "--tempFolder", os.path.dirname(tx_out_file),
                       vrn_file, align_bam, ref_file, tx_out_file]
                do.run(cmd, "Filter low frequency variants for DNA damage and strand bias")
                for out_plot in out_plot_files:
                    tx_plot_file = os.path.join("%s_qcSummary" % utils.splitext_plus(tx_out_file)[0], "plots",
                                                os.path.basename(out_plot))
                    if utils.file_exists(tx_plot_file):
                        shutil.move(tx_plot_file, out_plot)
        raw_file = vcfutils.bgzip_and_index(raw_file, items[0]["config"])
        data["vrn_file"] = _filter_to_info(raw_file, items[0])
        out_plot_files = [x for x in out_plot_files if utils.file_exists(x)]
        data["damage_plots"] = out_plot_files
        return data
175408
def _get_variant_file(x, key, suffix="", sample=None, ignore_do_upload=False):
    
    out = []
    fname = utils.get_in(x, key)
    upload_key = list(key)
    upload_key[-1] = "do_upload"
    do_upload = tz.get_in(tuple(upload_key), x, True)
    if fname and (ignore_do_upload or do_upload):
        if fname.endswith(".vcf.gz"):
            out.append({"path": fname,
                        "type": "vcf.gz",
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
            if utils.file_exists(fname + ".tbi"):
                out.append({"path": fname + ".tbi",
                            "type": "vcf.gz.tbi",
                            "index": True,
                            "ext": "%s%s" % (x["variantcaller"], suffix),
                            "variantcaller": x["variantcaller"]})
        elif fname.endswith((".vcf", ".bed", ".bedpe", ".bedgraph", ".cnr", ".cns", ".cnn", ".txt", ".tsv")):
            ftype = utils.splitext_plus(fname)[-1][1:]
            if ftype == "txt":
                extended_ftype = fname.split("-")[-1]
                if "/" not in extended_ftype:
                    ftype = extended_ftype
            out.append({"path": fname,
                        "type": ftype,
                        "ext": "%s%s" % (x["variantcaller"], suffix),
                        "variantcaller": x["variantcaller"]})
    if sample:
        out_sample = []
        for x in out:
            x["sample"] = sample
            out_sample.append(x)
        return out_sample
    else:
        return out
1234337
def _run_vardict_caller(align_bams, items, ref_file, assoc_files,
                          region=None, out_file=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            vrs = bedutils.population_variant_regions(items)
            target = shared.subset_variant_regions(
                vrs, region, out_file, items=items, do_merge=False)
            num_bams = len(align_bams)
            sample_vcf_names = []  # for individual sample names, given batch calling may be required
            for bamfile, item in zip(align_bams, items):
                # prepare commands
                sample = dd.get_sample_name(item)
                vardict = get_vardict_command(items[0])
                opts, var2vcf_opts = _vardict_options_from_config(items, config, out_file, target)
                vcfstreamsort = config_utils.get_program("vcfstreamsort", config)
                compress_cmd = "| bgzip -c" if tx_out_file.endswith("gz") else ""
                fix_ambig_ref = vcfutils.fix_ambiguous_cl()
                fix_ambig_alt = vcfutils.fix_ambiguous_cl(5)
                remove_dup = vcfutils.remove_dup_cl()
                py_cl = os.path.join(utils.get_bcbio_bin(), "py")
                jvm_opts = _get_jvm_opts(items[0], tx_out_file)
                setup = ("%s && unset JAVA_HOME &&" % utils.get_R_exports())
                contig_cl = vcfutils.add_contig_to_header_cl(ref_file, tx_out_file)
                lowfreq_filter = _lowfreq_linear_filter(0, False)
                cmd = ("{setup}{jvm_opts}{vardict} -G {ref_file} "
                       "-N {sample} -b {bamfile} {opts} "
                       "| teststrandbias.R "
                       "| var2vcf_valid.pl -A -N {sample} -E {var2vcf_opts} "
                       "| {contig_cl} | bcftools filter -i 'QUAL >= 0' | {lowfreq_filter} "
                       "| {fix_ambig_ref} | {fix_ambig_alt} | {remove_dup} | {vcfstreamsort} {compress_cmd}")
                if num_bams > 1:
                    temp_file_prefix = out_file.replace(".gz", "").replace(".vcf", "") + item["name"][1]
                    tmp_out = temp_file_prefix + ".temp.vcf"
                    tmp_out += ".gz" if out_file.endswith("gz") else ""
                    sample_vcf_names.append(tmp_out)
                    with file_transaction(item, tmp_out) as tx_tmp_file:
                        if not _is_bed_file(target):
                            vcfutils.write_empty_vcf(tx_tmp_file, config, samples=[sample])
                        else:
                            cmd += " > {tx_tmp_file}"
                            do.run(cmd.format(**locals()), "Genotyping with VarDict: Inference", {})
                else:
                    if not _is_bed_file(target):
                        vcfutils.write_empty_vcf(tx_out_file, config, samples=[sample])
                    else:
                        cmd += " > {tx_out_file}"
                        do.run(cmd.format(**locals()), "Genotyping with VarDict: Inference", {})
            if num_bams > 1:
                # N.B. merge_variant_files wants region in 1-based end-inclusive
                # coordinates. Thus use bamprep.region_to_gatk
                vcfutils.merge_variant_files(orig_files=sample_vcf_names,
                                             out_file=tx_out_file, ref_file=ref_file,
                                             config=config, region=bamprep.region_to_gatk(region))
    return out_file
172966
def run(call_file, ref_file, vrn_files, data):
    
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if includes_missingalt(data):
        logger.info("Removing variants with missing alts from %s." % call_file)
        call_file = gatk_remove_missingalt(call_file, data)

    if "gatkcnn" in dd.get_tools_on(data):
        return _cnn_filter(call_file, vrn_files, data)
    elif config_utils.use_vqsr(algs, call_file):
        if vcfutils.is_gvcf_file(call_file):
            raise ValueError("Cannot force gVCF output with joint calling using tools_on: [gvcf] and use VQSR. "
                             "Try using cutoff-based soft filtering with tools_off: [vqsr]")
        snp_file, indel_file = vcfutils.split_snps_indels(call_file, ref_file, data["config"])
        snp_filter_file = _variant_filtration(snp_file, ref_file, vrn_files, data, "SNP",
                                              vfilter.gatk_snp_cutoff)
        indel_filter_file = _variant_filtration(indel_file, ref_file, vrn_files, data, "INDEL",
                                                vfilter.gatk_indel_cutoff)
        orig_files = [snp_filter_file, indel_filter_file]
        out_file = "%scombined.vcf.gz" % os.path.commonprefix(orig_files)
        combined_file = vcfutils.combine_variant_files(orig_files, out_file, ref_file, data["config"])
        return combined_file
    else:
        snp_filter = vfilter.gatk_snp_cutoff(call_file, data)
        indel_filter = vfilter.gatk_indel_cutoff(snp_filter, data)
        return indel_filter
1098605
def print_vcf(data):
    """"""
    id_name = "."
    qual = "."
    chrom = data['chrom']
    pos = data['pre_pos']
    nt_ref = data['nt'][1]
    nt_snp = data['nt'][0]
    flt = "PASS"
    info = "ID=%s" % data['mature']
    frmt = "GT:NR:NS"
    gntp = "%s:%s:%s" % (_genotype(data), data["counts"], data["diff"])
    print("\t".join(map(str, [chrom, pos, id_name, nt_ref, nt_snp, qual, flt, info, frmt, gntp])), file=STDOUT, end="")
1123205
def _run_genotype_gvcfs_gatk3(data, region, vrn_files, ref_file, out_file):
    
    if not utils.file_exists(out_file):
        broad_runner = broad.runner_from_config(data["config"])
        with file_transaction(data, out_file) as tx_out_file:
            assoc_files = tz.get_in(("genome_resources", "variation"), data, {})
            if not assoc_files: assoc_files = {}
            params = ["-T", "GenotypeGVCFs",
                      "-R", ref_file, "-o", tx_out_file,
                      "-L", bamprep.region_to_gatk(region),
                      "--max_alternate_alleles", "4"]
            for vrn_file in vrn_files:
                params += ["--variant", vrn_file]
            if assoc_files.get("dbsnp"):
                params += ["--dbsnp", assoc_files["dbsnp"]]
            broad_runner.new_resources("gatk-haplotype")
            cores = dd.get_cores(data)
            if cores > 1:
                # GATK performs poorly with memory usage when parallelizing
                # with a large number of cores but makes use of extra memory,
                # so we cap at 6 cores.
                # See issue #1565 for discussion
                # Recent GATK 3.x versions also have race conditions with multiple
                # threads, so limit to 1 and keep memory available
                # https://gatkforums.broadinstitute.org/wdl/discussion/8718/concurrentmodificationexception-in-gatk-3-7-genotypegvcfs
                # params += ["-nt", str(min(6, cores))]
                memscale = {"magnitude": 0.9 * cores, "direction": "increase"}
            else:
                memscale = None
            broad_runner.run_gatk(params, memscale=memscale, parallel_gc=True)
    return vcfutils.bgzip_and_index(out_file, data["config"])
95450
def postprocess_variants(items):
    
    vrn_key = "vrn_file"
    if not isinstance(items, dict):
        items = [utils.to_single_data(x) for x in items]
        if "vrn_file_joint" in items[0]:
            vrn_key = "vrn_file_joint"
    data, items = _get_batch_representative(items, vrn_key)
    items = cwlutils.unpack_tarballs(items, data)
    data = cwlutils.unpack_tarballs(data, data)
    cur_name = "%s, %s" % (dd.get_sample_name(data), get_variantcaller(data, require_bam=False))
    logger.info("Finalizing variant calls: %s" % cur_name)
    orig_vrn_file = data.get(vrn_key)
    data = _symlink_to_workdir(data, [vrn_key])
    data = _symlink_to_workdir(data, ["config", "algorithm", "variant_regions"])
    if data.get(vrn_key):
        logger.info("Calculating variation effects for %s" % cur_name)
        ann_vrn_file, vrn_stats = effects.add_to_vcf(data[vrn_key], data)
        if ann_vrn_file:
            data[vrn_key] = ann_vrn_file
        if vrn_stats:
            data["vrn_stats"] = vrn_stats
        orig_items = _get_orig_items(items)
        logger.info("Annotate VCF file: %s" % cur_name)
        data[vrn_key] = annotation.finalize_vcf(data[vrn_key], get_variantcaller(data, require_bam=False),
                                                orig_items)
        if cwlutils.is_cwl_run(data):
            logger.info("Annotate with population level variation data")
            ann_file = population.run_vcfanno(data[vrn_key], data)
            if ann_file:
                data[vrn_key] = ann_file
        logger.info("Filtering for %s" % cur_name)
        data[vrn_key] = variant_filtration(data[vrn_key], dd.get_ref_file(data),
                                           tz.get_in(("genome_resources", "variation"), data, {}),
                                           data, orig_items)
        logger.info("Prioritization for %s" % cur_name)
        prio_vrn_file = prioritize.handle_vcf_calls(data[vrn_key], data, orig_items)
        if prio_vrn_file != data[vrn_key]:
            data[vrn_key] = prio_vrn_file
            logger.info("Germline extraction for %s" % cur_name)
            data = germline.extract(data, orig_items)

        if dd.get_align_bam(data):
            data = damage.run_filter(data[vrn_key], dd.get_align_bam(data), dd.get_ref_file(data),
                                     data, orig_items)
    if orig_vrn_file and os.path.samefile(data[vrn_key], orig_vrn_file):
        data[vrn_key] = orig_vrn_file
    return [[data]]
324514
def _run_internal(self,
                    context,
                    pipeline_key,
                    root_pipeline_key,
                    caller_output):
    """"""
    self._set_values_internal(
        context, pipeline_key, root_pipeline_key, caller_output,
        _PipelineRecord.RUN)
    logging.debug('Running %s(*%s, **%s)#%s',
                  self._class_path, _short_repr(self.args),
                  _short_repr(self.kwargs), self._pipeline_key.name())
    return self.run(*self.args, **self.kwargs)
1331298
def to_vcf(in_tsv, data):
    
    call_convert = {"Amp": "DUP", "Del": "DEL"}
    out_file = "%s.vcf" % utils.splitext_plus(in_tsv)[0]
    if not utils.file_uptodate(out_file, in_tsv):
        with file_transaction(data, out_file) as tx_out_file:
            with open(in_tsv) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    out_handle.write(VCF_HEADER + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t%s\n"
                                     % (dd.get_sample_name(data)))
                    header = in_handle.readline().split("\t")
                    for cur in (dict(zip(header, l.split("\t"))) for l in in_handle):
                        if cur["Amp_Del"] in call_convert:
                            svtype = call_convert[cur["Amp_Del"]]
                            info = "SVTYPE=%s;END=%s;SVLEN=%s;FOLD_CHANGE_LOG=%s;PROBES=%s;GENE=%s" % (
                                svtype, cur["End"], int(cur["End"]) - int(cur["Start"]),
                                cur["Log2ratio"], cur["Ab_Seg"], cur["Gene"])
                            out_handle.write("\t".join([cur["Chr"], cur["Start"], ".", "N", "<%s>" % (svtype),
                                                        ".", ".", info, "GT", "1/1"]) + "\n")
    return vcfutils.sort_by_ref(out_file, data)
161651
def run_file(path):
    """"""
    if PY26:
        dirpath, name, _ = splitname(path)
        found = imp.find_module(name, [dirpath])
        module = imp.load_module("__main__", *found)
        return vars(module)
    else:
        return runpy.run_path(path, run_name="__main__")
129313
def write_vcf(path, callset, rename=None, number=None, description=None,
              fill=None, write_header=True):
    """"""

    names, callset = normalize_callset(callset)

    with open(path, 'w') as vcf_file:
        if write_header:
            write_vcf_header(vcf_file, names, callset=callset, rename=rename,
                             number=number, description=description)
        write_vcf_data(vcf_file, names, callset=callset, rename=rename, fill=fill)
732476
def run_vcfanno(vcf_file, data, decomposed=False):
    
    conf_files = dd.get_vcfanno(data)
    if conf_files:
        with_basepaths = collections.defaultdict(list)
        gemini_basepath = _back_compatible_gemini(conf_files, data)
        for f in conf_files:
            name = os.path.splitext(os.path.basename(f))[0]
            if f.endswith(".lua"):
                conf_file = None
                lua_file = f
            else:
                conf_file = f
                lua_file = "%s.lua" % utils.splitext_plus(conf_file)[0]
            if lua_file and not os.path.exists(lua_file):
                lua_file = None
            data_basepath = gemini_basepath if name == "gemini" else None
            if conf_file and os.path.exists(conf_file):
                with_basepaths[(data_basepath, name)].append(conf_file)
            if lua_file and os.path.exists(lua_file):
                with_basepaths[(data_basepath, name)].append(lua_file)
        conf_files = with_basepaths.items()
    out_file = None
    if conf_files:
        VcfannoIn = collections.namedtuple("VcfannoIn", ["conf", "lua"])
        bp_files = collections.defaultdict(list)
        for (data_basepath, name), anno_files in conf_files:
            anno_files = list(set(anno_files))
            if len(anno_files) == 1:
                cur = VcfannoIn(anno_files[0], None)
            elif len(anno_files) == 2:
                lua_files = [x for x in anno_files if x.endswith(".lua")]
                assert len(lua_files) == 1, anno_files
                lua_file = lua_files[0]
                anno_files.remove(lua_file)
                cur = VcfannoIn(anno_files[0], lua_file)
            else:
                raise ValueError("Unexpected annotation group %s" % anno_files)
            bp_files[data_basepath].append(cur)
        for data_basepath, anno_files in bp_files.items():
            ann_file = vcfanno.run(vcf_file, [x.conf for x in anno_files],
                                   [x.lua for x in anno_files], data,
                                   basepath=data_basepath,
                                   decomposed=decomposed)
            if ann_file:
                out_file = ann_file
                vcf_file = ann_file
    return out_file
62815
def variantcall_sample(data, region=None, align_bams=None, out_file=None):
    
    if out_file is None or not os.path.exists(out_file) or not os.path.lexists(out_file):
        utils.safe_makedir(os.path.dirname(out_file))
        ref_file = dd.get_ref_file(data)
        config = data["config"]
        caller_fns = get_variantcallers()
        caller_fn = caller_fns[config["algorithm"].get("variantcaller")]
        if len(align_bams) == 1:
            items = [data]
        else:
            items = multi.get_orig_items(data)
            assert len(items) == len(align_bams)
        assoc_files = tz.get_in(("genome_resources", "variation"), data, {})
        if not assoc_files: assoc_files = {}
        for bam_file in align_bams:
            bam.index(bam_file, data["config"], check_timestamp=False)
        out_file = caller_fn(align_bams, items, ref_file, assoc_files, region, out_file)
    if region:
        data["region"] = region
    data["vrn_file"] = out_file
    return [data]
160592
def variants_to_vcf(variants, variant_to_metadata, out=sys.stdout):
    

    # TODO: The variant metadata dictionary (in the `VariantCollection`)
    # contains different data depending on the original input file format (VCF,
    # MAF, CSV).  It's easy to output variants originally in VCF format, but we
    # might want to consider how fields from MAF (for example) map to those in
    # VCF.

    # TODO: The VCF file we output doesn't contain any VCF metadata headers, as
    # the original headers were thrown away when the VCF file was parsed. We
    # may want to keep some of that information and/or infer some of the
    # headers based on the variants themselves. The former is difficult because
    # merge conflicts will inevitably occur; the latter is difficult because
    # the variants themselves don't contain all the information required for
    # these metadata headers (e.g., descriptions).
    #
    # As a side note, adding headers for certain fields will make them parse
    # correctly -- into an integer instead of a string (the default), for
    # example.

    # TODO: If we maintain headers (see above TODO), what should happen if
    # variant sources use different reference genomes?
    #
    # If we don't maintain headers, what should the default reference genome
    # be? This code chose one fairly arbitrarily.

    # TODO: If we end up needing more functions to "build" VCF record fields
    # (see functions below), we may want to abstract away the individual
    # functions and create a mapping from field to format function.

    def get_metadata_field(key, variant, default='.'):
        """Retrieve field from variant metadata dictionary."""
        val = variant_to_metadata[variant].get(key)
        if val is None:
            return default
        return val

    def build_filter_field(variant):
        """Build the filter field from the given variant.

        The `filter` field, as stored in the variants metadata dictionary,
        comes in 3 flavors:

            - empty list: 1+ filters were run and none failed
            - non-empty list: 1+ filters were run and 1+ failed
            - `None`: no filters were run

        This function maps each of these internal representations to their
        corresponding VCF representations.
        """
        filter_metadata = get_metadata_field('filter', variant)
        if type(filter_metadata) == list:
            return 'PASS' if filter_metadata == [] else ';'.join(filter_metadata)
        else:
            # TODO: Can the filter field ever be something other than the 3
            # possibilities described in this function's doc comment?
            bad_value_msg = (
                'Filter metadata field took on unexpected value `{}`. Update '
                'code to account for this value.').format(str(filter_metadata))
            assert filter_metadata == '.', bad_value_msg
            return filter_metadata

    def build_info_field(variant):
        """Build the info field from the given variant.

        Format is `<key>=<val>,...;<key>=<val>,...;<key>=<val>,...`.
        """

        def build_info_pair(key, val):
            """Build key/val pair for info field."""
            # Note: Different from `val == True`, which returns True when `val == 1`.
            if val is True:
                return key

            if type(val) == list:
                val = ','.join(map(str, val))
            else:
                val = str(val)
            return '{}={}'.format(key, val)

        info_dict = get_metadata_field('info', variant, default={})
        if not info_dict:
            return '.'

        return ';'.join(build_info_pair(k, v) for (k, v) in info_dict.items())

    def build_format_field(variant):
        """Build the sample format string from the given variant.

        Each sample column follows this format for the specified variant.
        """
        sample_info = get_metadata_field('sample_info', variant, default={})
        return ':'.join(list(sample_info.values())[0].keys()) if sample_info else '.'

    def build_sample_fields(variant):
        """Build the sample fields for the given variant."""

        def build_sample_field(sample):
            """Build a specific sample's field (i.e., one sample column)."""
            sample_vals = sample_info[sample].values()
            return ':'.join(build_sample_val(val) for val in sample_vals)

        def build_sample_val(sample_val):
            """Build a specific value for a sample (i.e., one value for one column)."""
            if type(sample_val) is list:
                return ','.join(map(str, sample_val))
            elif sample_val is not None:
                return str(sample_val)
            else:
                return '.'

        sample_info = get_metadata_field('sample_info', variant, default={})
        return list(build_sample_field(sample) for sample in sample_info)

    def build_vcf_record(variant, add_sample_info):
        """Return a list of all the variant's VCF record fields."""
        record = [
            str(variant.original_contig),
            str(variant.original_start),
            get_metadata_field('id', variant),
            variant.original_ref,
            variant.original_alt,
            str(get_metadata_field('qual', variant)),
            build_filter_field(variant),
            build_info_field(variant),
        ]
        if add_sample_info:
            record.append(build_format_field(variant))
            record.extend(build_sample_fields(variant))
        return record

    def merge_duplicate_variants():
        """Merge duplicate variants (according to ID) and return *list* of merged, sorted variants.

        Multiple `Variant`s can have the same VCF id (e.g., those variants which
        originally had > 1 alternate base), but we can't output a VCF file with multiple
        records having the same id. This function merges those duplicate variants.
        """

        def construct_id2variants():
            """Construct dictionary which maps variant IDs to `Variant`s."""
            id2variants = defaultdict(list)
            for variant in variants:
                # Note: For variants without IDs (e.g., those that came from
                # MAF files), we assign the missing value.
                variant_id = get_metadata_field('id', variant)
                id2variants[variant_id].append(variant)
            return id2variants

        def merge_variant_list(duplicates):
            """Merge duplicate variant list into one."""

            # TODO: Currently assumes that only alternate bases differ, but we may
            # want or need to merge variants that have the same ID but different
            # contigs, positions, or reference bases. If merging isn't possible (which
            # I think is the case for many situations), it may be easiest to assign
            # the "missing value" to the IDs and move on.

            assert len(duplicates) > 0
            assert all(v.original_contig == duplicates[0].original_contig and
                       v.original_start == duplicates[0].original_start and
                       v.original_ref == duplicates[0].original_ref for v in duplicates)
            import copy
            merged = copy.copy(duplicates[0])
            merged.original_alt = ','.join(duplicate.original_alt for duplicate in duplicates)
            return merged

        id2variants = construct_id2variants()
        variants_no_id = id2variants.pop('.', [])  # don't want to merge variants w/ no id
        merged_variants = list(map(merge_variant_list, id2variants.values())) + variants_no_id

        # groups variants by contig; relative ordering of contigs doesn't matter
        return sorted(
            merged_variants,
            key=lambda v: (str(v.original_contig), str(v.original_start)))

    def get_sample_names():
        """Return the sample names for all variants."""

        # TODO: For now, ensures that every variant has the same samples. If they didn't,
        # we'd have to use the missing value for VCF records with no data on a given
        # sample.  In and of itself, that isn't a problem until the format field contains
        # GT (genotype), which is required for all samples if present in the format field.
        #
        # A couple of ways to handle this:
        #
        #   1) Ignore this requirement of VCF files. I'd assume this would still be
        #   compatible with varcode, but may not be with other VCF parsers.
        #   2) Remove the GT field from those records where this rule would be violated.
        #   This is sad because we lose information.
        #
        # It's also important to note that if/when we combine all samples, we'll have to
        # add a `.` (missing value) for those samples in which a variant has no data.
        # Check the VCF spec to make sure this is valid; if not, we may have to write
        # `.:.:.:.` -- one `.` for each field in the format string.
        #
        # Moreover, we'll want to note whether we care about maintaining the relative
        # ordering of the sample names in the original VCF files. This probably isn't
        # important, but the code below does not do this (it effectively alphabetizes the
        # sample names because of the use of `set`).

        sample_names = set()
        for variant in variants:
            sample_info = get_metadata_field('sample_info', variant, default={})
            addl_sample_names = set(sample_info.keys())

            # Ensures all variants have the same samples.
            if sample_names and sample_names != addl_sample_names:
                return []

            sample_names.update(addl_sample_names)

        return list(sample_names)

    headers = ['#CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']
    sample_names = get_sample_names()
    if sample_names:
        headers += ['FORMAT'] + sample_names

    unique_variants_list = merge_duplicate_variants()

    # usually we have just one reference genome for the whole variant collection
    # but if variants from multiple sources have been merged then we might
    # not be able to write out a VCF since the individual variants may be using
    # different coordinate systems
    genome_names = list({v.ensembl.reference_name for v in unique_variants_list})
    if len(genome_names) > 1:
        raise ValueError(
            "Cannot create VCF for variants with multiple reference genomes: %s" % (
                genome_names,))
    genome_name = genome_names[0]
    print('##fileformat=VCFv4.2', file=out)
    print('##reference=%s' % genome_name, file=out)
    print('\t'.join(headers), file=out)
    for variant in unique_variants_list:
        record = build_vcf_record(variant, add_sample_info=bool(sample_names))
        print('\t'.join(record), file=out)
440756
def run_region(data, region, vrn_files, out_file):
    
    broad_runner = broad.runner_from_config(data["config"])
    if broad_runner.gatk_type() == "gatk4":
        genomics_db = _run_genomicsdb_import(vrn_files, region, out_file, data)
        return _run_genotype_gvcfs_genomicsdb(genomics_db, region, out_file, data)
    else:
        vrn_files = _batch_gvcfs(data, region, vrn_files, dd.get_ref_file(data), out_file)
        return _run_genotype_gvcfs_gatk3(data, region, vrn_files, dd.get_ref_file(data), out_file)
503211
def _select_final_variants(base_vcf, out_vcf, config):
    
    min_depth = int(config["algorithm"].get("min_depth", 4))
    with file_transaction(out_vcf) as tx_out_file:
        with open(base_vcf) as in_handle:
            with open(tx_out_file, "w") as out_handle:
                for line in in_handle:
                    if line.startswith("#"):
                        passes = True
                    else:
                        passes = _passes_cortex_depth(line, min_depth)
                    if passes:
                        out_handle.write(line)
    return out_vcf