<commit_before>"""Examine callable regions following genome mapping of short reads.

Identifies callable analysis regions surrounded by larger regions lacking
aligned bases. This allows parallelization of smaller chromosome chunks
through post-processing and variant calling, with each sub-section
mapping handled separately.

Regions are split to try to maintain relative uniformity across the
genome and avoid extremes of large blocks or large numbers of
small blocks.
"""
import contextlib
import copy
import operator
import os
import shutil

import pybedtools
import pysam
from py_descriptive_statistics import Enum as Stats

from bcbio import bam, broad, utils
from bcbio.log import logger, setup_local_logging
from bcbio.distributed.messaging import parallel_runner
from bcbio.distributed.split import parallel_split_combine
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import shared

def parallel_callable_loci(in_bam, ref_file, config):
    num_cores = config["algorithm"].get("num_cores", 1)
    config = copy.deepcopy(config)
    config["algorithm"]["memory_adjust"] = {"direction": "decrease", "magnitude": 2}
    data = {"work_bam": in_bam, "sam_ref": ref_file, "config": config}
    parallel = {"type": "local", "cores": num_cores, "module": "bcbio.distributed"}
    runner = parallel_runner(parallel, {}, config)
    split_fn = shared.process_bam_by_chromosome("-callable.bed", "work_bam")
    out = parallel_split_combine([[data]], split_fn, runner,
                                 "calc_callable_loci", "combine_bed",
                                 "callable_bed", ["config"])[0]
    return out[0]["callable_bed"]

def combine_bed(in_files, out_file, config):
    """Combine multiple BED files into a single output.
    """
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                for in_file in in_files:
                    with open(in_file) as in_handle:
                        shutil.copyfileobj(in_handle, out_handle)
    return out_file

def calc_callable_loci(data, region=None, out_file=None):
    """Determine callable bases for input BAM using Broad's CallableLoci walker.

    http://www.broadinstitute.org/gatk/gatkdocs/
    org_broadinstitute_sting_gatk_walkers_coverage_CallableLoci.html
    """
    if data["config"].get("parallel", {}).get("log_queue"):
        handler = setup_local_logging(data["config"], data["config"]["parallel"])
    else:
        handler = None
    broad_runner = broad.runner_from_config(data["config"])
    if out_file is None:
        out_file = "%s-callable.bed" % os.path.splitext(data["work_bam"])[0]
    out_summary = "%s-callable-summary.txt" % os.path.splitext(data["work_bam"])[0]
    variant_regions = data["config"]["algorithm"].get("variant_regions", None)
    # set a maximum depth to avoid calling in repetitive regions with excessive coverage
    max_depth = int(1e6 if data["config"]["algorithm"].get("coverage_depth", "").lower() == "super-high"
                    else 2.5e4)
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            bam.index(data["work_bam"], data["config"])
            params = ["-T", "CallableLoci",
                      "-R", data["sam_ref"],
                      "-I", data["work_bam"],
                      "--minDepth", "0",
                      "--downsample_to_coverage", str(max_depth + 1000),
                      "--minMappingQuality", "0",
                      "--maxFractionOfReadsWithLowMAPQ", "1.1",
                      "--maxDepth", str(max_depth),
                      "--out", tx_out_file,
                      "--summary", out_summary]
            ready_region = shared.subset_variant_regions(variant_regions, region, tx_out_file)
            if ready_region:
                params += ["-L", ready_region]
            if ((variant_regions and ready_region and os.path.isfile(ready_region))
                 or not variant_regions or not region):
                broad_runner.run_gatk(params)
            else:
                with open(out_file, "w") as out_handle:
                    for tregion in get_ref_bedtool(data["sam_ref"], data["config"]):
                        if tregion.chrom == region:
                            out_handle.write("%s\t%s\t%s\tNO_COVERAGE\n" %
                                             (tregion.chrom, tregion.start, tregion.stop))
    if handler and hasattr(handler, "close"):
        handler.close()
    return [{"callable_bed": out_file, "config": data["config"], "work_bam": data["work_bam"]}]

def sample_callable_bed(bam_file, ref_file, config):
    """Retrieve callable regions for a sample subset by defined analysis regions.
    """
    out_file = "%s-callable_sample.bed" % os.path.splitext(bam_file)[0]
    callable_bed = parallel_callable_loci(bam_file, ref_file, config)
    input_regions_bed = config["algorithm"].get("variant_regions", None)
    if not utils.file_uptodate(out_file, callable_bed):
        with file_transaction(out_file) as tx_out_file:
            callable_regions = pybedtools.BedTool(callable_bed)
            filter_regions = callable_regions.filter(lambda x: x.name == "CALLABLE")
            if input_regions_bed:
                if not utils.file_uptodate(out_file, input_regions_bed):
                    input_regions = pybedtools.BedTool(input_regions_bed)
                    filter_regions.intersect(input_regions).saveas(tx_out_file)
            else:<commit_msg>Provide ability to recover from java memory issues when running callability assessment. It can reliably fail on quick running processes in memory limited environments<commit_after>"""Examine callable regions following genome mapping of short reads.

Identifies callable analysis regions surrounded by larger regions lacking
aligned bases. This allows parallelization of smaller chromosome chunks
through post-processing and variant calling, with each sub-section
mapping handled separately.

Regions are split to try to maintain relative uniformity across the
genome and avoid extremes of large blocks or large numbers of
small blocks.
"""
import contextlib
import copy
import operator
import os
import time
import shutil
import subprocess

import pybedtools
import pysam
from py_descriptive_statistics import Enum as Stats

from bcbio import bam, broad, utils
from bcbio.log import logger, setup_local_logging
from bcbio.distributed.messaging import parallel_runner
from bcbio.distributed.split import parallel_split_combine
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import shared

def parallel_callable_loci(in_bam, ref_file, config):
    num_cores = config["algorithm"].get("num_cores", 1)
    config = copy.deepcopy(config)
    config["algorithm"]["memory_adjust"] = {"direction": "decrease", "magnitude": 2}
    data = {"work_bam": in_bam, "sam_ref": ref_file, "config": config}
    parallel = {"type": "local", "cores": num_cores, "module": "bcbio.distributed"}
    runner = parallel_runner(parallel, {}, config)
    split_fn = shared.process_bam_by_chromosome("-callable.bed", "work_bam")
    out = parallel_split_combine([[data]], split_fn, runner,
                                 "calc_callable_loci", "combine_bed",
                                 "callable_bed", ["config"])[0]
    return out[0]["callable_bed"]

def combine_bed(in_files, out_file, config):
    """Combine multiple BED files into a single output.
    """
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                for in_file in in_files:
                    with open(in_file) as in_handle:
                        shutil.copyfileobj(in_handle, out_handle)
    return out_file

def calc_callable_loci(data, region=None, out_file=None):
    """Determine callable bases for input BAM using Broad's CallableLoci walker.

    http://www.broadinstitute.org/gatk/gatkdocs/
    org_broadinstitute_sting_gatk_walkers_coverage_CallableLoci.html
    """
    if data["config"].get("parallel", {}).get("log_queue"):
        handler = setup_local_logging(data["config"], data["config"]["parallel"])
    else:
        handler = None
    broad_runner = broad.runner_from_config(data["config"])
    if out_file is None:
        out_file = "%s-callable.bed" % os.path.splitext(data["work_bam"])[0]
    out_summary = "%s-callable-summary.txt" % os.path.splitext(data["work_bam"])[0]
    variant_regions = data["config"]["algorithm"].get("variant_regions", None)
    # set a maximum depth to avoid calling in repetitive regions with excessive coverage
    max_depth = int(1e6 if data["config"]["algorithm"].get("coverage_depth", "").lower() == "super-high"
                    else 2.5e4)
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out_file:
            bam.index(data["work_bam"], data["config"])
            params = ["-T", "CallableLoci",
                      "-R", data["sam_ref"],
                      "-I", data["work_bam"],
                      "--minDepth", "0",
                      "--downsample_to_coverage", str(max_depth + 1000),
                      "--minMappingQuality", "0",
                      "--maxFractionOfReadsWithLowMAPQ", "1.1",
                      "--maxDepth", str(max_depth),
                      "--out", tx_out_file,
                      "--summary", out_summary]
            ready_region = shared.subset_variant_regions(variant_regions, region, tx_out_file)
            if ready_region:
                params += ["-L", ready_region]
            if ((variant_regions and ready_region and os.path.isfile(ready_region))
                 or not variant_regions or not region):
                _memory_safe_run(broad_runner, params, data, region)
            else:
                with open(out_file, "w") as out_handle:
                    for tregion in get_ref_bedtool(data["sam_ref"], data["config"]):
                        if tregion.chrom == region:
                            out_handle.write("%s\t%s\t%s\tNO_COVERAGE\n" %
                                             (tregion.chrom, tregion.start, tregion.stop))
    if handler and hasattr(handler, "close"):
        handler.close()
    return [{"callable_bed": out_file, "config": data["config"], "work_bam": data["work_bam"]}]

def _memory_safe_run(runner, params, data, region):
    """Detect java memory errors which occur during large scale processing.

    When running on multiple simultaneous projects or with small regions, this
    can error out with temporary memory issues. This detects and recovers from those.
    """
    max_runs = 5
    num_runs = 0
    while 1:
        try:
            runner.run_gatk(params, log_error=False)
            break
        except subprocess.CalledProcessError, msg:
            if num_runs < max_runs and ("insufficient memory" in str(msg) or
                                        "did not provide enough memory" in str(msg) or
                                        "A fatal error has been detected" in str(msg)):
                logger.info("Memory issue with callability assessment on %s, %s. Retrying."
                            % (region, data["work_bam"]))
                time.sleep(30)
                num_runs += 1
            else:
                logger.exception()
                raise

def sample_callable_bed(bam_file, ref_file, config):
    """Retrieve callable regions for a sample subset by defined analysis regions.
    """
    out_file = "%s-callable_sample.bed" % os.path.splitext(bam_file)[0]
    callable_bed = parallel_callable_loci(bam_file, ref_file, config)
    input_regions_bed = config["algorithm"].get("variant_regions", None)
    if not utils.file_uptodate(out_file, callable_bed):
        with file_transaction(out_file) as tx_out_file:
            callable_regions = pybedtools.BedTool(callable_bed)
            filter_regions = callable_regions.filter(lambda x: x.name == "CALLABLE")
            if input_regions_bed:
                if not utils.file_uptodate(out_file, input_regions_bed):
                    input_regions = pybedtools.BedTool(input_regions_bed)
                    filter_regions.intersect(input_regions).saveas(tx_out_file)
            else:
{'fraction_non_alphanumeric': 0.10493976877678066, 'fraction_numerical': 0.003718974856496079, 'mean_word_length': 2.7587359465208143, 'pattern_counts': {'":': 22, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1163842', 'n_tokens_mistral': 3664, 'n_tokens_neox': 3536, 'n_words': 891}