<commit_before>"""Copy number detection with CNVkit with specific support for targeted sequencing.

http://cnvkit.readthedocs.org
"""
import copy
import os
import sys
import tempfile

import pybedtools
import numpy as np
import toolz as tz

from bcbio import install, utils
from bcbio.bam import ref
from bcbio.distributed.multi import run_multicore, zeromq_aware_logging
from bcbio.distributed.transaction import file_transaction
from bcbio.heterogeneity import chromhacks
from bcbio.pipeline import datadict as dd
from bcbio.pipeline import config_utils
from bcbio.variation import bedutils, effects, vcfutils
from bcbio.provenance import do
from bcbio.structural import annotate, shared, plot

def run(items, background=None):
    """Detect copy number variations from batched set of samples using CNVkit.
    """
    if not background: background = []
    return _cnvkit_by_type(items, background)

def _sv_workdir(data):
    return utils.safe_makedir(os.path.join(data["dirs"]["work"], "structural",
                                           dd.get_sample_name(data), "cnvkit"))

def _cnvkit_by_type(items, background):
    """Dispatch to specific CNVkit functionality based on input type.
    """
    if len(items + background) == 1:
        return _run_cnvkit_single(items[0])
    elif vcfutils.get_paired_phenotype(items[0]):
        return _run_cnvkit_cancer(items, background)
    else:
        return _run_cnvkit_population(items, background)

def _associate_cnvkit_out(ckouts, items):
    """Associate cnvkit output with individual items.
    """
    assert len(ckouts) == len(items)
    out = []
    for ckout, data in zip(ckouts, items):
        ckout = copy.deepcopy(ckout)
        ckout["variantcaller"] = "cnvkit"
        if utils.file_exists(ckout["cns"]):
            ckout = _add_seg_to_output(ckout, data)
            ckout = _add_gainloss_to_output(ckout, data)
            ckout = _add_segmetrics_to_output(ckout, data)
            ckout = _add_variantcalls_to_output(ckout, data)
            # ckout = _add_coverage_bedgraph_to_output(ckout, data)
            ckout = _add_cnr_bedgraph_and_bed_to_output(ckout, data)
            if "svplots" in dd.get_tools_on(data):
                ckout = _add_plots_to_output(ckout, data)
            if "sv" not in data:
                data["sv"] = []
            data["sv"].append(ckout)
        out.append(data)
    return out

def _run_cnvkit_single(data, background=None):
    """Process a single input file with BAM or uniform background.
    """
    work_dir = _sv_workdir(data)
    test_bams = [data["align_bam"]]
    if background:
        background_bams = [x["align_bam"] for x in background]
        background_name = os.path.splitext(os.path.basename(background_bams[0]))[0]
    else:
        background_bams = []
        background_name = None
    ckouts = _run_cnvkit_shared([data], test_bams, background_bams, work_dir,
                               background_name=background_name)
    if not ckouts:
        return [data]
    else:
        assert len(ckouts) == 1
        return _associate_cnvkit_out(ckouts, [data])

def _run_cnvkit_cancer(items, background):
    """Run CNVkit on a tumor/normal pair.
    """
    paired = vcfutils.get_paired_bams([x["align_bam"] for x in items], items)
    work_dir = _sv_workdir(paired.tumor_data)
    ckouts = _run_cnvkit_shared([paired.tumor_data], [paired.tumor_bam], [paired.normal_bam],
                               work_dir, background_name=paired.normal_name)
    if not ckouts:
        return items
    assert len(ckouts) == 1
    tumor_data = _associate_cnvkit_out(ckouts, [paired.tumor_data])
    normal_data = [x for x in items if dd.get_sample_name(x) != paired.tumor_name]
    return tumor_data + normal_data

def _run_cnvkit_population(items, background):
    """Run CNVkit on a population of samples.

    Tries to calculate background based on case/controls, otherwise uses
    a flat background for each sample and calls independently.
    """
    assert not background
    inputs, background = shared.find_case_control(items)
    work_dir = _sv_workdir(inputs[0])
    ckouts = _run_cnvkit_shared(inputs, [x["align_bam"] for x in inputs],
                                [x["align_bam"] for x in background], work_dir,
                                background_name=dd.get_sample_name(background[0]) if len(background) > 0 else None)
    return _associate_cnvkit_out(ckouts, inputs) + background

def _get_cmd():
    return os.path.join(os.path.dirname(sys.executable), "cnvkit.py")

def _bam_to_outbase(bam_file, work_dir):
    """Convert an input BAM file into CNVkit expected output.
    """
    out_base = os.path.splitext(os.path.basename(bam_file))[0].split(".")[0]
    return os.path.join(work_dir, out_base)

def _run_cnvkit_shared(items, test_bams, background_bams, work_dir, background_name=None):
    """Shared functionality to run CNVkit, parallelizing over multiple BAM files.
    """
    raw_work_dir = utils.safe_makedir(os.path.join(work_dir, "raw"))

    background_cnn = os.path.join(raw_work_dir,
                                  "%s_background.cnn" % (background_name if background_name else "flat"))
    ckouts = []
    for test_bam in test_bams:
        out_base = _bam_to_outbase(test_bam, raw_work_dir)
        ckouts.append({"cnr": "%s.cns" % out_base,
                       "cns": "%s.cns" % out_base,
                       "back_cnn": background_cnn})
    if not utils.file_exists(ckouts[0]["cnr"]):
        data = items[0]
        cov_interval = dd.get_coverage_interval(data)
        raw_target_bed, access_bed = _get_target_access_files(cov_interval, data, work_dir)
        # bail out if we ended up with no regions
        if not utils.file_exists(raw_target_bed):
            return {}
        raw_target_bed = annotate.add_genes(raw_target_bed, data)
        parallel = {"type": "local", "cores": dd.get_cores(data), "progs": ["cnvkit"]}
        target_bed, antitarget_bed = _cnvkit_targets(raw_target_bed, access_bed, cov_interval, raw_work_dir, data)
        def _bam_to_itype(bam):
            return "background" if bam in background_bams else "evaluate"
        split_cnns = run_multicore(_cnvkit_coverage,
                                   [(bam, bed, _bam_to_itype(bam), raw_work_dir, data)
                                    for bam in test_bams + background_bams
                                    for bed in _split_bed(target_bed, data) + _split_bed(antitarget_bed, data)],
                                   data["config"], parallel)
        coverage_cnns = _merge_coverage(split_cnns, data)
        background_cnn = _cnvkit_background([x["file"] for x in coverage_cnns if x["itype"] == "background"],
                                            background_cnn, target_bed, antitarget_bed, data)
        fixed_cnrs = run_multicore(_cnvkit_fix,
                                   [(cnns, background_cnn, data) for cnns in
                                    tz.groupby("bam", [x for x in coverage_cnns
                                                       if x["itype"] == "evaluate"]).values()],
                                      data["config"], parallel)
        called_segs = run_multicore(_cnvkit_segment,
                                    [(cnr, cov_interval, data) for cnr in fixed_cnrs],
                                    data["config"], parallel)
    return ckouts

@utils.map_wrap
@zeromq_aware_logging
def _cnvkit_segment(cnr_file, cov_interval, data):
    """Perform segmentation and copy number calling on normalized inputs
    """
    out_file = "%s.cns" % os.path.splitext(cnr_file)[0]
    if not utils.file_uptodate(out_file, cnr_file):
        with file_transaction(data, out_file) as tx_out_file:
            local_sitelib = os.path.join(install.get_defaults().get("tooldir", "/usr/local"),
                                            "lib", "R", "site-library")
            cmd = [_get_cmd(), "segment", "-o", tx_out_file, "--rlibpath", local_sitelib, cnr_file]
            if cov_interval == "genome":
                cmd += ["--threshold", "0.00001"]<commit_msg>CNVkit: do not pass old R site-library, no longer needed with conda R libraries<commit_after>"""Copy number detection with CNVkit with specific support for targeted sequencing.

http://cnvkit.readthedocs.org
"""
import copy
import os
import sys
import tempfile

import pybedtools
import numpy as np
import toolz as tz

from bcbio import install, utils
from bcbio.bam import ref
from bcbio.distributed.multi import run_multicore, zeromq_aware_logging
from bcbio.distributed.transaction import file_transaction
from bcbio.heterogeneity import chromhacks
from bcbio.pipeline import datadict as dd
from bcbio.pipeline import config_utils
from bcbio.variation import bedutils, effects, vcfutils
from bcbio.provenance import do
from bcbio.structural import annotate, shared, plot

def run(items, background=None):
    """Detect copy number variations from batched set of samples using CNVkit.
    """
    if not background: background = []
    return _cnvkit_by_type(items, background)

def _sv_workdir(data):
    return utils.safe_makedir(os.path.join(data["dirs"]["work"], "structural",
                                           dd.get_sample_name(data), "cnvkit"))

def _cnvkit_by_type(items, background):
    """Dispatch to specific CNVkit functionality based on input type.
    """
    if len(items + background) == 1:
        return _run_cnvkit_single(items[0])
    elif vcfutils.get_paired_phenotype(items[0]):
        return _run_cnvkit_cancer(items, background)
    else:
        return _run_cnvkit_population(items, background)

def _associate_cnvkit_out(ckouts, items):
    """Associate cnvkit output with individual items.
    """
    assert len(ckouts) == len(items)
    out = []
    for ckout, data in zip(ckouts, items):
        ckout = copy.deepcopy(ckout)
        ckout["variantcaller"] = "cnvkit"
        if utils.file_exists(ckout["cns"]):
            ckout = _add_seg_to_output(ckout, data)
            ckout = _add_gainloss_to_output(ckout, data)
            ckout = _add_segmetrics_to_output(ckout, data)
            ckout = _add_variantcalls_to_output(ckout, data)
            # ckout = _add_coverage_bedgraph_to_output(ckout, data)
            ckout = _add_cnr_bedgraph_and_bed_to_output(ckout, data)
            if "svplots" in dd.get_tools_on(data):
                ckout = _add_plots_to_output(ckout, data)
            if "sv" not in data:
                data["sv"] = []
            data["sv"].append(ckout)
        out.append(data)
    return out

def _run_cnvkit_single(data, background=None):
    """Process a single input file with BAM or uniform background.
    """
    work_dir = _sv_workdir(data)
    test_bams = [data["align_bam"]]
    if background:
        background_bams = [x["align_bam"] for x in background]
        background_name = os.path.splitext(os.path.basename(background_bams[0]))[0]
    else:
        background_bams = []
        background_name = None
    ckouts = _run_cnvkit_shared([data], test_bams, background_bams, work_dir,
                               background_name=background_name)
    if not ckouts:
        return [data]
    else:
        assert len(ckouts) == 1
        return _associate_cnvkit_out(ckouts, [data])

def _run_cnvkit_cancer(items, background):
    """Run CNVkit on a tumor/normal pair.
    """
    paired = vcfutils.get_paired_bams([x["align_bam"] for x in items], items)
    work_dir = _sv_workdir(paired.tumor_data)
    ckouts = _run_cnvkit_shared([paired.tumor_data], [paired.tumor_bam], [paired.normal_bam],
                               work_dir, background_name=paired.normal_name)
    if not ckouts:
        return items
    assert len(ckouts) == 1
    tumor_data = _associate_cnvkit_out(ckouts, [paired.tumor_data])
    normal_data = [x for x in items if dd.get_sample_name(x) != paired.tumor_name]
    return tumor_data + normal_data

def _run_cnvkit_population(items, background):
    """Run CNVkit on a population of samples.

    Tries to calculate background based on case/controls, otherwise uses
    a flat background for each sample and calls independently.
    """
    assert not background
    inputs, background = shared.find_case_control(items)
    work_dir = _sv_workdir(inputs[0])
    ckouts = _run_cnvkit_shared(inputs, [x["align_bam"] for x in inputs],
                                [x["align_bam"] for x in background], work_dir,
                                background_name=dd.get_sample_name(background[0]) if len(background) > 0 else None)
    return _associate_cnvkit_out(ckouts, inputs) + background

def _get_cmd():
    return os.path.join(os.path.dirname(sys.executable), "cnvkit.py")

def _bam_to_outbase(bam_file, work_dir):
    """Convert an input BAM file into CNVkit expected output.
    """
    out_base = os.path.splitext(os.path.basename(bam_file))[0].split(".")[0]
    return os.path.join(work_dir, out_base)

def _run_cnvkit_shared(items, test_bams, background_bams, work_dir, background_name=None):
    """Shared functionality to run CNVkit, parallelizing over multiple BAM files.
    """
    raw_work_dir = utils.safe_makedir(os.path.join(work_dir, "raw"))

    background_cnn = os.path.join(raw_work_dir,
                                  "%s_background.cnn" % (background_name if background_name else "flat"))
    ckouts = []
    for test_bam in test_bams:
        out_base = _bam_to_outbase(test_bam, raw_work_dir)
        ckouts.append({"cnr": "%s.cns" % out_base,
                       "cns": "%s.cns" % out_base,
                       "back_cnn": background_cnn})
    if not utils.file_exists(ckouts[0]["cnr"]):
        data = items[0]
        cov_interval = dd.get_coverage_interval(data)
        raw_target_bed, access_bed = _get_target_access_files(cov_interval, data, work_dir)
        # bail out if we ended up with no regions
        if not utils.file_exists(raw_target_bed):
            return {}
        raw_target_bed = annotate.add_genes(raw_target_bed, data)
        parallel = {"type": "local", "cores": dd.get_cores(data), "progs": ["cnvkit"]}
        target_bed, antitarget_bed = _cnvkit_targets(raw_target_bed, access_bed, cov_interval, raw_work_dir, data)
        def _bam_to_itype(bam):
            return "background" if bam in background_bams else "evaluate"
        split_cnns = run_multicore(_cnvkit_coverage,
                                   [(bam, bed, _bam_to_itype(bam), raw_work_dir, data)
                                    for bam in test_bams + background_bams
                                    for bed in _split_bed(target_bed, data) + _split_bed(antitarget_bed, data)],
                                   data["config"], parallel)
        coverage_cnns = _merge_coverage(split_cnns, data)
        background_cnn = _cnvkit_background([x["file"] for x in coverage_cnns if x["itype"] == "background"],
                                            background_cnn, target_bed, antitarget_bed, data)
        fixed_cnrs = run_multicore(_cnvkit_fix,
                                   [(cnns, background_cnn, data) for cnns in
                                    tz.groupby("bam", [x for x in coverage_cnns
                                                       if x["itype"] == "evaluate"]).values()],
                                      data["config"], parallel)
        called_segs = run_multicore(_cnvkit_segment,
                                    [(cnr, cov_interval, data) for cnr in fixed_cnrs],
                                    data["config"], parallel)
    return ckouts

@utils.map_wrap
@zeromq_aware_logging
def _cnvkit_segment(cnr_file, cov_interval, data):
    """Perform segmentation and copy number calling on normalized inputs
    """
    out_file = "%s.cns" % os.path.splitext(cnr_file)[0]
    if not utils.file_uptodate(out_file, cnr_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = [_get_cmd(), "segment", "-o", tx_out_file, cnr_file]
            if cov_interval == "genome":
                cmd += ["--threshold", "0.00001"]
{'fraction_non_alphanumeric': 0.09912045412014223, 'fraction_numerical': 0.002619923897448693, 'mean_word_length': 2.67622105021784, 'pattern_counts': {'":': 14, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7027280', 'n_tokens_mistral': 5223, 'n_tokens_neox': 5040, 'n_words': 1166}