<commit_before>"""Utilities for manipulating variant files in standard VCF format.
"""

from collections import namedtuple, defaultdict
import copy
import os
import pprint
import shutil
import subprocess

import toolz as tz

import six
from six.moves import zip

from bcbio import broad, utils
from bcbio.bam import ref
from bcbio.distributed.multi import run_multicore, zeromq_aware_logging
from bcbio.distributed.split import parallel_split_combine
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import config_utils, tools
from bcbio.pipeline import datadict as dd
from bcbio.provenance import do

from bcbio.log import logger

# ## Tumor/normal paired cancer analyses

PairedData = namedtuple("PairedData", ["tumor_bam", "tumor_name",
                                       "normal_bam", "normal_name", "normal_panel",
                                       "tumor_config", "tumor_data", "normal_data"])

def is_paired_analysis(align_bams, items):
    """Determine if BAMs are from a tumor/normal paired analysis.
    """
    return get_paired_bams(align_bams, items) is not None


def somatic_batches(items):
    """Group items into somatic calling batches (tumor-only or tumor/normal).

    Returns batches, where a data item may be in pairs, and somatic and non_somatic
    (which are the original list of items).
    """
    non_somatic = []
    somatic = []
    data_by_batches = defaultdict(list)
    for data in items:
        if not get_paired_phenotype(data):
            non_somatic.append(data)
        else:
            somatic.append(data)
            batches = dd.get_batches(data)
            if batches:
                for batch in batches:
                    data_by_batches[batch].append(data)
    return data_by_batches.values(), somatic, non_somatic

def get_paired(items):
    return get_paired_bams([dd.get_align_bam(d) for d in items], items)

def get_paired_bams(align_bams, items):
    """Split aligned bams into tumor / normal pairs if this is a paired analysis.
    Allows cases with only tumor BAMs to handle callers that can work without
    normal BAMs or with normal VCF panels.
    """
    tumor_bam, tumor_name, tumor_data, normal_bam, normal_name, normal_panel, tumor_config, normal_data = (None,) * 8
    for bamfile, item in zip(align_bams, items):
        phenotype = get_paired_phenotype(item)
        if phenotype == "normal":
            normal_bam = bamfile
            normal_name = dd.get_sample_name(item)
            normal_data = item
        elif phenotype == "tumor":
            tumor_bam = bamfile
            tumor_name = dd.get_sample_name(item)
            tumor_data = item
            tumor_config = item["config"]
            normal_panel = dd.get_background_variant(item)
    if tumor_bam or tumor_name:
        return PairedData(tumor_bam, tumor_name, normal_bam,
                          normal_name, normal_panel, tumor_config,
                          tumor_data, normal_data)

def get_somatic_variantcallers(items):
    """Retrieve all variant callers for somatic calling, handling somatic/germline.
    """
    out = []
    for data in items:
        vcs = dd.get_variantcaller(data)
        if isinstance(vcs, dict) and "somatic" in vcs:
            vcs = vcs["somatic"]
        if not isinstance(vcs, (list, tuple)):
            vcs = [vcs]
        out += vcs
    return set(vcs)

def check_paired_problems(items):
    """Check for incorrectly paired tumor/normal samples in a batch.
    """
    # ensure we're in a paired batch
    if not get_paired(items):
        return
    num_tumor = len([x for x in items if dd.get_phenotype(x).lower() == "tumor"])
    if num_tumor > 1:
        raise ValueError("Unsupported configuration: found multiple tumor samples in batch %s: %s" %
                         (tz.get_in(["metadata", "batch"], items[0]),
                          [dd.get_sample_name(data) for data in items]))
    elif num_tumor == 0 and any(dd.get_phenotype(data).lower() == "normal" for data in items):
        raise ValueError("Found normal sample without tumor in batch %s: %s" %
                         (tz.get_in(["metadata", "batch"], items[0]),
                          [dd.get_sample_name(data) for data in items]))
    else:
        vcs = get_somatic_variantcallers(items)
        if "mutect" in vcs or "mutect2" in vcs or "strelka2" in vcs:
            paired = get_paired(items)
            if not (paired.normal_data or paired.normal_panel):
                raise ValueError("MuTect, MuTect2 and Strelka2 somatic calling requires normal sample or panel: %s" %
                                 [dd.get_sample_name(data) for data in items])

def get_paired_phenotype(data):
    """Retrieve the phenotype for a paired tumor/normal analysis.
    """
    allowed_names = set(["tumor", "normal"])
    p = tz.get_in(["metadata", "phenotype"], data)
    return p if p in allowed_names else None

# ## General utilities

def fix_ambiguous_cl(column=4):
    """awk command to replace non-N ambiguous REF bases with N.

    Some callers include these if present in the reference genome but GATK does
    not like them.
    """
    return r"""awk -F$'\t' -v OFS='\t' '{if ($0 !~ /^#/) gsub(/[KMRYSWBVHDXkmryswbvhdx]/, "N", $%s) } {print}'""" % column

def remove_dup_cl():
    """awk command line to remove duplicate alleles where the ref and alt are the same.
    """
    return r""" awk -F$'\t' -v OFS='\t' '$1!~/^#/ && $4 == $5 {next} {print}'"""

def get_indelcaller(d_or_c):
    """Retrieve string for indelcaller to use, or empty string if not specified.
    """<commit_msg>allow t-only without pon for mutect2 (#3489)<commit_after>"""Utilities for manipulating variant files in standard VCF format.
"""

from collections import namedtuple, defaultdict
import copy
import os
import pprint
import shutil
import subprocess

import toolz as tz

import six
from six.moves import zip

from bcbio import broad, utils
from bcbio.bam import ref
from bcbio.distributed.multi import run_multicore, zeromq_aware_logging
from bcbio.distributed.split import parallel_split_combine
from bcbio.distributed.transaction import file_transaction
from bcbio.pipeline import config_utils, tools
from bcbio.pipeline import datadict as dd
from bcbio.provenance import do

from bcbio.log import logger

# ## Tumor/normal paired cancer analyses

PairedData = namedtuple("PairedData", ["tumor_bam", "tumor_name",
                                       "normal_bam", "normal_name", "normal_panel",
                                       "tumor_config", "tumor_data", "normal_data"])

def is_paired_analysis(align_bams, items):
    """Determine if BAMs are from a tumor/normal paired analysis.
    """
    return get_paired_bams(align_bams, items) is not None


def somatic_batches(items):
    """Group items into somatic calling batches (tumor-only or tumor/normal).

    Returns batches, where a data item may be in pairs, and somatic and non_somatic
    (which are the original list of items).
    """
    non_somatic = []
    somatic = []
    data_by_batches = defaultdict(list)
    for data in items:
        if not get_paired_phenotype(data):
            non_somatic.append(data)
        else:
            somatic.append(data)
            batches = dd.get_batches(data)
            if batches:
                for batch in batches:
                    data_by_batches[batch].append(data)
    return data_by_batches.values(), somatic, non_somatic

def get_paired(items):
    return get_paired_bams([dd.get_align_bam(d) for d in items], items)

def get_paired_bams(align_bams, items):
    """Split aligned bams into tumor / normal pairs if this is a paired analysis.
    Allows cases with only tumor BAMs to handle callers that can work without
    normal BAMs or with normal VCF panels.
    """
    tumor_bam, tumor_name, tumor_data, normal_bam, normal_name, normal_panel, tumor_config, normal_data = (None,) * 8
    for bamfile, item in zip(align_bams, items):
        phenotype = get_paired_phenotype(item)
        if phenotype == "normal":
            normal_bam = bamfile
            normal_name = dd.get_sample_name(item)
            normal_data = item
        elif phenotype == "tumor":
            tumor_bam = bamfile
            tumor_name = dd.get_sample_name(item)
            tumor_data = item
            tumor_config = item["config"]
            normal_panel = dd.get_background_variant(item)
    if tumor_bam or tumor_name:
        return PairedData(tumor_bam, tumor_name, normal_bam,
                          normal_name, normal_panel, tumor_config,
                          tumor_data, normal_data)

def get_somatic_variantcallers(items):
    """Retrieve all variant callers for somatic calling, handling somatic/germline.
    """
    out = []
    for data in items:
        vcs = dd.get_variantcaller(data)
        if isinstance(vcs, dict) and "somatic" in vcs:
            vcs = vcs["somatic"]
        if not isinstance(vcs, (list, tuple)):
            vcs = [vcs]
        out += vcs
    return set(vcs)

def check_paired_problems(items):
    """Check for incorrectly paired tumor/normal samples in a batch."""
    # ensure we're in a paired batch
    if not get_paired(items):
        return
    num_tumor = len([x for x in items if dd.get_phenotype(x).lower() == "tumor"])
    if num_tumor > 1:
        raise ValueError("Unsupported configuration: found multiple tumor samples in batch %s: %s" %
                         (tz.get_in(["metadata", "batch"], items[0]),
                          [dd.get_sample_name(data) for data in items]))
    elif num_tumor == 0 and any(dd.get_phenotype(data).lower() == "normal" for data in items):
        raise ValueError("Found normal sample without tumor in batch %s: %s" %
                         (tz.get_in(["metadata", "batch"], items[0]),
                          [dd.get_sample_name(data) for data in items]))
    else:
        vcs = get_somatic_variantcallers(items)
        if "mutect" in vcs or "mutect2" in vcs or "strelka2" in vcs:
            paired = get_paired(items)
            if not (paired.normal_data or paired.normal_panel):
                # give a warning for mutect2, error out for mutect and strelka
                if "mutect2" in vcs:
                    logger.info("Using a PON or a normal sample is recommended!")
                else:
                    raise ValueError("MuTect and Strelka2 somatic calling requires normal sample or panel: %s" %
                                 [dd.get_sample_name(data) for data in items])

def get_paired_phenotype(data):
    """Retrieve the phenotype for a paired tumor/normal analysis.
    """
    allowed_names = set(["tumor", "normal"])
    p = tz.get_in(["metadata", "phenotype"], data)
    return p if p in allowed_names else None

# ## General utilities

def fix_ambiguous_cl(column=4):
    """awk command to replace non-N ambiguous REF bases with N.

    Some callers include these if present in the reference genome but GATK does
    not like them.
    """
    return r"""awk -F$'\t' -v OFS='\t' '{if ($0 !~ /^#/) gsub(/[KMRYSWBVHDXkmryswbvhdx]/, "N", $%s) } {print}'""" % column

def remove_dup_cl():
    """awk command line to remove duplicate alleles where the ref and alt are the same.
    """
    return r""" awk -F$'\t' -v OFS='\t' '$1!~/^#/ && $4 == $5 {next} {print}'"""

def get_indelcaller(d_or_c):
    """Retrieve string for indelcaller to use, or empty string if not specified.
    """
{'fraction_non_alphanumeric': 0.0972987149226331, 'fraction_numerical': 0.002972287787394003, 'mean_word_length': 2.684380032206119, 'pattern_counts': {'":': 4, '<': 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': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4909027', 'n_tokens_mistral': 3653, 'n_tokens_neox': 3406, 'n_words': 1092}