Combine VCF files from the same sample into a single output file.

    Handles cases where we split files into SNPs/Indels for processing then
    need to merge back into a final file.
426110
def combine_variant_files(orig_files, out_file, ref_file, config,
                          quiet_out=True, region=None):
    
    in_pipeline = False
    if isinstance(orig_files, dict):
        file_key = config["file_key"]
        in_pipeline = True
        orig_files = orig_files[file_key]
    if not utils.file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            exist_files = [x for x in orig_files if os.path.exists(x)]
            ready_files = run_multicore(p_bgzip_and_index, [[x, config] for x in exist_files], config)
            dict_file = "%s.dict" % utils.splitext_plus(ref_file)[0]
            cores = dd.get_num_cores({"config": config})
            memscale = {"magnitude": 0.9 * cores, "direction": "increase"} if cores > 1 else None
            cmd = ["picard"] + broad.get_picard_opts(config, memscale) + \
                  ["MergeVcfs", "D=%s" % dict_file, "O=%s" % tx_out_file] + \
                  ["I=%s" % f for f in ready_files]
            cmd = "%s && %s" % (utils.get_java_clprep(), " ".join(cmd))
            do.run(cmd, "Combine variant files")
    if out_file.endswith(".gz"):
        bgzip_and_index(out_file, config)
    if in_pipeline:
        return [{file_key: out_file, "region": region, "sam_ref": ref_file, "config": config}]
    else:
        return out_file
124492
def _summarize_combined(samples, vkey):
    
    validate_dir = utils.safe_makedir(os.path.join(samples[0]["dirs"]["work"], vkey))
    combined, _ = _group_validate_samples(samples, vkey, [["metadata", "validate_combine"]])
    for vname, vitems in combined.items():
        if vname:
            cur_combined = collections.defaultdict(int)
            for data in sorted(vitems, key=lambda x: x.get("lane", dd.get_sample_name(x))):
                validations = [variant.get(vkey) for variant in data.get("variants", [])]
                validations = [v for v in validations if v]
                if len(validations) == 0 and vkey in data:
                    validations = [data.get(vkey)]
                for validate in validations:
                    with open(validate["summary"]) as in_handle:
                        reader = csv.reader(in_handle)
                        next(reader)  # header
                        for _, caller, vtype, metric, value in reader:
                            cur_combined[(caller, vtype, metric)] += int(value)
            out_csv = os.path.join(validate_dir, "grading-summary-%s.csv" % vname)
            with open(out_csv, "w") as out_handle:
                writer = csv.writer(out_handle)
                header = ["sample", "caller", "vtype", "metric", "value"]
                writer.writerow(header)
                for (caller, variant_type, category), val in cur_combined.items():
                    writer.writerow(["combined-%s" % vname, caller, variant_type, category, val])
            plots = validateplot.classifyplot_from_valfile(out_csv)
1030031
def merge(args):
    
    from jcvi.formats.gff import merge as gmerge

    p = OptionParser(merge.__doc__)
    p.set_home("maker")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    outdir, outputgff = args
    fsnames, suffix = get_fsnames(outdir)
    nfs = len(fsnames)
    cmd = op.join(opts.maker_home, "bin/gff3_merge")

    outfile = "merge.sh"
    write_file(outfile, mergesh.format(suffix, cmd))

    # Generate per split directory
    # Note that gff3_merge write to /tmp, so I limit processes here to avoid
    # filling up disk space
    sh("parallel -j 8 merge.sh {} ::: " + " ".join(fsnames))

    # One final output
    gffnames = glob("*.all.gff")
    assert len(gffnames) == nfs

    # Again, DO NOT USE gff3_merge to merge with a smallish /tmp/ area
    gfflist = "gfflist"
    fw = open(gfflist, "w")
    print("\n".join(gffnames), file=fw)
    fw.close()

    nlines = sum(1 for x in open(gfflist))
    assert nlines == nfs  # Be extra, extra careful to include all results
    gmerge([gfflist, "-o", outputgff])
    logging.debug("Merged GFF file written to `{0}`".format(outputgff))
1239405
def _merge_summary(in_files, out_file, data):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                for i, in_file in enumerate(in_files):
                    with open(in_file) as in_handle:
                        for j, line in enumerate(in_handle):
                            if j == 0:
                                if i == 0:
                                    out_handle.write(line)
                            else:
                                out_handle.write(line)
    return out_file
620558
def merge(self, workdir, pot_files, out_dvdb, delete_source=True):
        
        # We work with absolute paths.
        pot_files = [os.path.abspath(s) for s in list_strings(pot_files)]
        if not os.path.isabs(out_dvdb):
            out_dvdb = os.path.join(os.path.abspath(workdir), os.path.basename(out_dvdb))

        if self.verbose:
            print("Will merge %d files into output DVDB %s" % (len(pot_files), out_dvdb))
            for i, f in enumerate(pot_files):
                print(" [%d] %s" % (i, f))

        # Handle the case of a single file since mrgddb uses 1 to denote GS files!
        if len(pot_files) == 1:
            with open(pot_files[0], "r") as inh, open(out_dvdb, "w") as out:
                for line in inh:
                    out.write(line)
            return out_dvdb

        self.stdin_fname, self.stdout_fname, self.stderr_fname = \
            map(os.path.join, 3 * [os.path.abspath(workdir)], ["mrgdvdb.stdin", "mrgdvdb.stdout", "mrgdvdb.stderr"])

        inp = StringIO()
        inp.write(out_dvdb + "\n")             # Name of the output file.
        inp.write(str(len(pot_files)) + "\n")  # Number of input POT files.

        # Names of the POT files.
        for fname in pot_files:
            inp.write(fname + "\n")

        self.stdin_data = [s for s in inp.getvalue()]

        with open(self.stdin_fname, "wt") as fh:
            fh.writelines(self.stdin_data)
            # Force OS to write data to disk.
            fh.flush()
            os.fsync(fh.fileno())

        retcode = self.execute(workdir)
        if retcode == 0 and delete_source:
            # Remove pot files.
            for f in pot_files:
                try:
                    os.remove(f)
                except IOError:
                    pass

        return out_dvdb
1196898
def combine_files(joinfile, sourcefiles, client):
    
    from mediasync.conf import msettings

    joinfile = joinfile.strip('/')

    if joinfile.endswith('.css'):
        dirname = msettings['CSS_PATH'].strip('/')
        separator = '\n'
    elif joinfile.endswith('.js'):
        dirname = msettings['JS_PATH'].strip('/')
        separator = ';\n'
    else:
        # By-pass this file since we only join CSS and JS.
        return None

    buffer = cStringIO.StringIO()

    for sourcefile in sourcefiles:
        sourcepath = os.path.join(client.media_root, dirname, sourcefile)
        if os.path.isfile(sourcepath):
            f = open(sourcepath)
            buffer.write(f.read())
            f.close()
            buffer.write(separator)

    filedata = buffer.getvalue()
    buffer.close()
    return (filedata, dirname)
844774
def subset_vcf(in_file, region, out_file, config):
    
    work_file = vcfutils.bgzip_and_index(in_file, config)
    if not file_exists(out_file):
        with file_transaction(config, out_file) as tx_out_file:
            bcftools = config_utils.get_program("bcftools", config)
            region_str = bamprep.region_to_gatk(region)
            cmd = "{bcftools} view -r {region_str} {work_file} > {tx_out_file}"
            do.run(cmd.format(**locals()), "subset %s: %s" % (os.path.basename(work_file), region_str))
    return out_file
409381
def generate_simple_vcf(filename, variant_collection):
    
    contigs = []
    positions = []
    refs = []
    alts = []
    for variant in variant_collection:
        contigs.append("chr" + variant.contig)
        positions.append(variant.start)
        refs.append(variant.ref)
        alts.append(variant.alt)
    df = pd.DataFrame()
    df["contig"] = contigs
    df["position"] = positions
    df["id"] = ["."] * len(variant_collection)
    df["ref"] = refs
    df["alt"] = alts
    df["qual"] = ["."] * len(variant_collection)
    df["filter"] = ["."] * len(variant_collection)
    df["info"] = ["."] * len(variant_collection)
    df["format"] = ["GT:AD:DP"] * len(variant_collection)
    normal_ref_depths = [randint(1, 10) for v in variant_collection]
    normal_alt_depths = [randint(1, 10) for v in variant_collection]
    df["n1"] = ["0:%d,%d:%d" % (normal_ref_depths[i], normal_alt_depths[i],
                                normal_ref_depths[i] + normal_alt_depths[i])
                for i in range(len(variant_collection))]
    tumor_ref_depths = [randint(1, 10) for v in variant_collection]
    tumor_alt_depths = [randint(1, 10) for v in variant_collection]
    df["t1"] = ["0/1:%d,%d:%d" % (tumor_ref_depths[i], tumor_alt_depths[i], tumor_ref_depths[i] + tumor_alt_depths[i])
                for i in range(len(variant_collection))]

    with open(filename, "w") as f:
        f.write("##fileformat=VCFv4.1\n")
        f.write("##reference=file:///projects/ngs/resources/gatk/2.3/ucsc.hg19.fasta\n")

    with open(filename, "a") as f:
        df.to_csv(f, sep="\t", index=None, header=None)
1219693
def get_multisample_vcf(fnames, name, caller, data):
    
    unique_fnames = []
    for f in fnames:
        if f not in unique_fnames:
            unique_fnames.append(f)
    out_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "gemini"))
    if len(unique_fnames) > 1:
        gemini_vcf = os.path.join(out_dir, "%s-%s.vcf.gz" % (name, caller))
        vrn_file_batch = None
        for variant in data.get("variants", []):
            if variant["variantcaller"] == caller and variant.get("vrn_file_batch"):
                vrn_file_batch = variant["vrn_file_batch"]
        if vrn_file_batch:
            utils.symlink_plus(vrn_file_batch, gemini_vcf)
            return gemini_vcf
        else:
            return vcfutils.merge_variant_files(unique_fnames, gemini_vcf, dd.get_ref_file(data),
                                                data["config"])
    else:
        gemini_vcf = os.path.join(out_dir, "%s-%s%s" % (name, caller, utils.splitext_plus(unique_fnames[0])[1]))
        utils.symlink_plus(unique_fnames[0], gemini_vcf)
        return gemini_vcf
1099835
def run_gvcfgenotyper(data, orig_region, vrn_files, out_file):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            regions = _find_gvcf_blocks(vrn_files[0], bamprep.region_to_gatk(orig_region),
                                        os.path.dirname(tx_out_file))
            if len(regions) == 1:
                _run_gvcfgenotyper(data, regions[0], vrn_files, tx_out_file)
            else:
                split_outs = [_run_gvcfgenotyper(data, r, vrn_files,
                                                 "%s-%s.vcf.gz" % (utils.splitext_plus(out_file)[0],
                                                                   r.replace(":", "_").replace("-", "_")))
                              for r in regions]
                vcfutils.concat_variant_files(split_outs, tx_out_file, regions,
                                              dd.get_ref_file(data), data["config"])
    return vcfutils.bgzip_and_index(out_file, data["config"])
207636
def summarize_vc(items):
    
    items = [utils.to_single_data(x) for x in utils.flatten(items)]
    items = [_normalize_vc_input(x) for x in items]
    items = validate.summarize_grading(items)
    items = [utils.to_single_data(x) for x in items]
    out = {"validate": validate.combine_validations(items),
           "variants": {"calls": [], "gvcf": [], "samples": []}}
    added = set([])
    variants_by_sample = collections.defaultdict(list)
    sample_order = []
    for data in items:
        batch_samples = data.get("batch_samples", [dd.get_sample_name(data)])
        for s in batch_samples:
            if s not in sample_order:
                sample_order.append(s)
        if data.get("vrn_file"):
            # Only get batches if we're actually doing variantcalling in bcbio
            # otherwise we'll be using the original files
            names = dd.get_batches(data) if dd.get_variantcaller(data) else None
            if not names:
                names = [dd.get_sample_name(data)]
            batch_name = names[0]
            if data.get("vrn_file_joint") is not None:
                to_add = [("vrn_file", "gvcf", dd.get_sample_name(data)),
                          ("vrn_file_joint", "calls", batch_name)]
            else:
                to_add = [("vrn_file", "calls", batch_name)]
            for vrn_key, out_key, name in to_add:
                cur_name = "%s-%s" % (name, dd.get_variantcaller(data))
                out_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data),
                                                                        "variants", out_key)),
                                        "%s.vcf.gz" % cur_name)
                for s in batch_samples:
                    variants_by_sample[s].append(out_file)
                if cur_name not in added:
                    added.add(cur_name)
                    # Ideally could symlink here but doesn't appear to work with
                    # Docker container runs on Toil where PATHs don't get remapped
                    utils.copy_plus(os.path.realpath(data[vrn_key]), out_file)
                    vcfutils.bgzip_and_index(out_file, data["config"])
                    out["variants"][out_key].append(out_file)
    for sample in sample_order:
        out["variants"]["samples"].append(variants_by_sample[sample])
    return [out]
525553
def make_simple_merged_vcf_with_no_combinations(self, ref_seq):
        ''''''
        if len(self) <= 1:
            return

        merged_vcf_record = self.vcf_records[0]

        for i in range(1, len(self.vcf_records), 1):
            if self.vcf_records[i].intersects(merged_vcf_record):
                return
            else:
                merged_vcf_record = merged_vcf_record.merge(self.vcf_records[i], ref_seq)

        self.vcf_records = [merged_vcf_record]
1032462
def cli(variant_file, vep, split, outfile, verbose, silent, check_info,
        allele_symbol, logfile, loglevel):
    
    from vcf_parser import logger, init_log

    if not loglevel:
        if verbose:
            loglevel = 'INFO'

    init_log(logger, logfile, loglevel)
    nr_of_variants = 0
    start = datetime.now()
    
    # with open(variant_file, 'r', encoding="utf-8") as f:
    #     for line in f:
    #         if not line.startswith('#'):
    #             nr_of_variants += 1
    
    if variant_file == '-':
        logger.info("Start parsing variants from stdin")
        my_parser = VCFParser(
            fsock=sys.stdin, 
            split_variants=split,
            check_info=check_info, 
            allele_symbol=allele_symbol
        )
    else:
        logger.info("Start parsing variants from file {0}".format(variant_file))
        my_parser = VCFParser(
            infile = variant_file,
            split_variants=split, 
            check_info=check_info, 
            allele_symbol=allele_symbol
        )

    if outfile:
        f = open(outfile, 'w', encoding='utf-8')
        logger.info("Printing vcf to file {0}".format(outfile))
    
    if not silent:
        logger.info("Printing vcf to stdout")
    else:
        logger.info("Skip printing since silent is active")
    
    for line in my_parser.metadata.print_header():
        if outfile:
            f.write(line+'\n')
        else:
            if not silent:
                print(line)
    try:
        for variant in my_parser:
            variant_line = '\t'.join([variant[head] for head in my_parser.header])
            if outfile:
                f.write(variant_line + '\n')
            else:
                if not silent:
                    print(variant_line)
            nr_of_variants += 1
    except SyntaxError as e:
        print(e)

    logger.info('Number of variants: {0}'.format(nr_of_variants))
    logger.info('Time to parse file: {0}'.format(str(datetime.now() - start)))
587618
def combine_count_files(files, out_file=None, ext=".fpkm"):
    
    files = list(files)
    if not files:
        return None
    assert all([file_exists(x) for x in files]), \
        "Some count files in %s do not exist." % files
    for f in files:
        assert file_exists(f), "%s does not exist or is empty." % f
    col_names = [os.path.basename(x.replace(ext, "")) for x in files]
    if not out_file:
        out_dir = os.path.join(os.path.dirname(files[0]))
        out_file = os.path.join(out_dir, "combined.counts")

    if file_exists(out_file):
        return out_file
    logger.info("Combining count files into %s." % out_file)
    row_names = []
    col_vals = defaultdict(list)
    for i, f in enumerate(files):
        vals = []
        if i == 0:
            with open(f) as in_handle:
                for line in in_handle:
                    rname, val = line.strip().split("\t")
                    row_names.append(rname)
                    vals.append(val)
        else:
            with open(f) as in_handle:
                for line in in_handle:
                    _, val = line.strip().split("\t")
                    vals.append(val)
        col_vals[col_names[i]] = vals

    df = pd.DataFrame(col_vals, index=row_names)
    df.to_csv(out_file, sep="\t", index_label="id")
    return out_file
774878
def finalize_vcf(in_file, variantcaller, items):
    
    out_file = "%s-annotated%s" % utils.splitext_plus(in_file)
    if not utils.file_uptodate(out_file, in_file):
        header_cl = _add_vcf_header_sample_cl(in_file, items, out_file)
        contig_cl = _add_contig_cl(in_file, items, out_file)
        cls = [x for x in (contig_cl, header_cl) if x]
        if cls:
            post_cl = " | ".join(cls) + " | "
        else:
            post_cl = None
        dbsnp_file = tz.get_in(("genome_resources", "variation", "dbsnp"), items[0])
        if dbsnp_file:
            out_file = _add_dbsnp(in_file, dbsnp_file, items[0], out_file, post_cl)
    if utils.file_exists(out_file):
        return vcfutils.bgzip_and_index(out_file, items[0]["config"])
    else:
        return in_file
253677
def split_snps_indels(orig_file, ref_file, config):
    
    base, ext = utils.splitext_plus(orig_file)
    snp_file = "{base}-snp{ext}".format(base=base, ext=ext)
    indel_file = "{base}-indel{ext}".format(base=base, ext=ext)
    for out_file, select_arg in [(snp_file, "--types snps"),
                                 (indel_file, "--exclude-types snps")]:
        if not utils.file_exists(out_file):
            with file_transaction(config, out_file) as tx_out_file:
                bcftools = config_utils.get_program("bcftools", config)
                output_type = "z" if out_file.endswith(".gz") else "v"
                cmd = "{bcftools} view -O {output_type} {orig_file} {select_arg} > {tx_out_file}"
                do.run(cmd.format(**locals()), "Subset to SNPs and indels")
        if out_file.endswith(".gz"):
            bgzip_and_index(out_file, config)
    return snp_file, indel_file
400458
def combine_express(samples, combined):
    """"""
    if not combined:
        return None
    to_combine = [dd.get_express_counts(x) for x in
                  dd.sample_data_iterator(samples) if dd.get_express_counts(x)]
    gtf_file = dd.get_gtf_file(samples[0][0])
    isoform_to_gene_file = os.path.join(os.path.dirname(combined), "isoform_to_gene.txt")
    isoform_to_gene_file = express.isoform_to_gene_name(
        gtf_file, isoform_to_gene_file, next(dd.sample_data_iterator(samples)))
    if len(to_combine) > 0:
        eff_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_counts"
        eff_counts_combined = count.combine_count_files(to_combine, eff_counts_combined_file, ext=".counts")
        to_combine = [dd.get_express_tpm(x) for x in
                      dd.sample_data_iterator(samples) if dd.get_express_tpm(x)]
        tpm_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_tpm"
        tpm_counts_combined = count.combine_count_files(to_combine, tpm_counts_combined_file)
        to_combine = [dd.get_express_fpkm(x) for x in dd.sample_data_iterator(samples)
                      if dd.get_express_fpkm(x)]
        fpkm_counts_combined_file = os.path.splitext(combined)[0] + ".isoform.express_fpkm"
        fpkm_counts_combined = count.combine_count_files(to_combine, fpkm_counts_combined_file, ext=".fpkm")
        return {'counts': eff_counts_combined, 'tpm': tpm_counts_combined,
                'fpkm': fpkm_counts_combined, 'isoform_to_gene': isoform_to_gene_file}
    return {}