New upstream version 1.1.9

## 1.1.8

## 1.1.9 (in progress)

- Fix for get VEP cache.

- Support Picard's new syntax for ReorderSam (REFERENCE -> SEQUENCE_DICTIONARY).

- Remove mitochondrial reads from ChIP/ATAC-seq calling.

- Add documentation describing ATAC-seq outputs.

- Add ENCODE library complexity metrics for ATAC/ChIP-seq to MultiQC report 

  (see https://www.encodeproject.org/data-standards/terms/#library for a description of the metrics)

- Add STAR sample-specific 2-pass. This helps assign a moderate number of reads per genes. Thanks

  to @naumenko-sa for the intial implementation and push to get this going.

- Index transcriptomes only once for pseudo/quasi aligner tools. This fixes race conditions that

  can happen.

- Add --buildversion option, for tracking which version of a gene build was used. This is used

  during `bcbio_setup_genome.py`. Suggested formats are source_version, so Ensembl_94, 

  EnsemblMetazoa_25, FlyBase_26, etc.

- Sort MACS2 bedgraph files before compressing. Thanks to @LMannarino for the suggestion.

- Check for the reserved field `sample` in RNA-seq metadata and quit with a useful error message. 

  Thanks to @marypiper for suggesting this.

- Split ATAC-seq BAM files into nucleosome-free and mono/di/tri nucleosome files, so we can call 

  peaks on them separately.

- Call peaks on NF/MN/DN/TN regions separately for each caller during ATAC-seq.

- Allow viral contamination to be assasyed on non tumor/normal samples.

- Ensure EBV coverage is calculated when run on genomes with it included as a contig.

## 1.1.8 (28 October 2019)

- Add `antibody` configuration option. Setting a specific antibody for ChIP-seq will use appropriate

  settings for that antibody. See the documentation for supported antibodies.

- Add `use_lowfreq_filter` for forcing vardict to report variants with low allelic frequency,

    return out_bam

def sort(in_bam, config, order="coordinate", out_dir=None):

def sort(in_bam, config, order="coordinate", out_dir=None, force=False):

    """Sort a BAM file, skipping if already present.

    """

    assert is_bam(in_bam), "%s in not a BAM file" % in_bam

    if bam_already_sorted(in_bam, config, order):

    if not force and bam_already_sorted(in_bam, config, order):

        return in_bam

    sort_stem = _get_sort_stem(in_bam, order, out_dir)

        do.run(cmd, message)

    index(out_bam, dd.get_config(data))

    return out_bam

def count_alignments(in_bam, data, filter=None):

    """

    count alignments in a BAM file passing a given filter. valid filter

    strings are available in the sambamba documentation:

    https://github.com/biod/sambamba/wiki/%5Bsambamba-view%5D-Filter-expression-syntax

    """

    sambamba = config_utils.get_program("sambamba", dd.get_config(data))

    num_cores = dd.get_num_cores(data)

    if not filter:

        filter_string = ""

        message = f"Counting alignments in {in_bam}."

    else:

        filter_string = "--filter {filter}"

        message = f"Counting alignments in {in_bam} matching {filter}."

    cmd = f"{sambamba} view -c --nthreads {num_cores} -f bam {filter_string} {in_bam}"

    logger.info(message)

    result = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE,

                            stderr=subprocess.PIPE)

    return int(result.stdout.decode().strip())

    if not file_exists(out_file):

        with tx_tmpdir(picard._config) as tmp_dir:

            with file_transaction(picard._config, out_file) as tx_out_file:

                dict_file = "%s.dict" % os.path.splitext(ref_file)[0]

                opts = [("INPUT", in_bam),

                        ("OUTPUT", tx_out_file),

                        ("REFERENCE", ref_file),

                        ("SEQUENCE_DICTIONARY", dict_file),

                        ("ALLOW_INCOMPLETE_DICT_CONCORDANCE", "true"),

                        ("TMP_DIR", tmp_dir)]

                picard.run("ReorderSam", opts)

import os

import shutil

import subprocess

import sys

import toolz as tz

from bcbio.distributed.transaction import file_transaction

from bcbio.provenance import do

from bcbio.log import logger

from bcbio.heterogeneity.chromhacks import get_mitochondrial_chroms

from bcbio.heterogeneity.chromhacks import (get_mitochondrial_chroms,

                                            get_nonmitochondrial_chroms)

from bcbio.chipseq import atac

def clean_chipseq_alignment(data):

    # lcr_bed = utils.get_in(data, ("genome_resources", "variation", "lcr"))

    method = dd.get_chip_method(data)

    if method == "atac":

        data = shift_ATAC(data)

    work_bam = dd.get_work_bam(data)

    work_bam = bam.sort(work_bam, dd.get_config(data))

    bam.index(work_bam, dd.get_config(data))

    clean_bam = remove_nonassembled_chrom(work_bam, data)

    clean_bam = remove_mitochondrial_reads(clean_bam, data)

    data = atac.calculate_complexity_metrics(clean_bam, data)

    if not dd.get_keep_multimapped(data):

        clean_bam = remove_multimappers(clean_bam, data)

    if not dd.get_keep_duplicates(data):

        clean_bam = bam.remove_duplicates(clean_bam, data)

    data["work_bam"] = clean_bam

    # for ATAC-seq, brewak alignments into NF, mono/di/tri nucleosome BAM files

    if method == "atac":

        data  = atac.split_ATAC(data)

    encode_bed = tz.get_in(["genome_resources", "variation", "encode_blacklist"], data)

    if encode_bed:

        data["work_bam"] = remove_blacklist_regions(dd.get_work_bam(data), encode_bed, data['config'])

                                       dd.get_work_bam(data), data)

    return [[data]]

def remove_mitochondrial_reads(bam_file, data):

    mito = get_mitochondrial_chroms(data)

    if not mito:

        logger.info(f"Mitochondrial chromosome not identified, skipping removal of "

                    "mitochondrial reads from {bam_file}.")

        return bam_file

    nonmito = get_nonmitochondrial_chroms(data)

    mito_bam = os.path.splitext(bam_file)[0] + "-noMito.bam"

    if utils.file_exists(mito_bam):

        return mito_bam

    samtools = config_utils.get_program("samtools", dd.get_config(data))

    nonmito_flag = " ".join(nonmito)

    num_cores = dd.get_num_cores(data)

    with file_transaction(mito_bam) as tx_out_bam:

        cmd = (f"{samtools} view -bh -@ {num_cores} {bam_file} {nonmito_flag} "

               f"> {tx_out_bam}")

        message = f"Removing mitochondrial reads on {','.join(mito)} from {bam_file}."

        do.run(cmd, message)

    return mito_bam

def remove_multimappers(bam_file, data):

    aligner = dd.get_aligner(data)

    if aligner:

        unique_bam = bam_file

        logger.warn("When a BAM file is given as input, bcbio skips removal of "

                    "multimappers.")

    logger.warn("ChIP/ATAC-seq usually requires duplicate marking, but it was disabled.")

    return unique_bam

def remove_nonassembled_chrom(bam_file, data):

    """Remove non-assembled contigs from the BAM file"""

    ref_file =  dd.get_ref_file(data)

    config = dd.get_config(data)

    bam.index(bam_file, config)

    fai = "%s.fai" % ref_file

    chrom = []

    with open(fai) as inh:

def _compute_deeptools_matrix(data):

    pass

def extract_NF_regions(data):

def shift_ATAC(data):

    """

    extract the nucleosome free regions from the work_bam. These regions will

    be < 100 bases

    shift the ATAC-seq alignments

    """

    MAX_FRAG_LENGTH = 100

    sieve = config_utils.get_program("alignmentSieve", data)

    work_bam = dd.get_work_bam(data)

    num_cores = dd.get_num_cores(data)

    out_file = os.path.splitext(work_bam)[0] + "-NF.bam"

    log_file = os.path.splitext(work_bam)[0] + "-NF.log"

    if file_exists(out_file):

        data["NF_bam"] = out_file

    out_file = os.path.splitext(work_bam)[0] + "-shifted.bam"

    log_file = os.path.splitext(work_bam)[0] + "-shifted.log"

    if utils.file_exists(out_file):

        data["work_bam"] = out_file

        return data

    unsorted_bam = os.path.splitext(out_file)[0] + ".unsorted.bam"

    if not utils.file_exists(out_file):

        with file_transaction(out_file) as tx_out_file, \

            file_transaction(log_file) as tx_log_file:

        tx_unsorted_bam = tx_out_file + ".unsorted"

            tx_unsorted_file = os.path.splitext(tx_out_file)[0] + ".tmp.bam"

            cmd = (

            f"{sieve} --bam ${work_bam} --outFile {tx_unsorted_bam} --ATACshift "

            f"--numberOfProcessors {num_cores} --maxFragmentLength {MAX_FRAG_LENGTH} "

                f"{sieve} --verbose --bam {work_bam} --outFile {tx_unsorted_file} --ATACshift "

                f"--numberOfProcessors {num_cores} --maxFragmentLength 0 "

                f"--minFragmentLength 0 "

                f"--minMappingQuality 10 "

                f"--filterMetrics {tx_log_file} ")

        do.run(cmd, "Extract NF regions from {work_bam} to {tx_unsorted_bam}.")

        tx_out_file = bam.sort(tx_unsorted_bam)

    data["NF_bam"] = out_file

            do.run(cmd, f"Shifting ATAC-seq alignments in {work_bam} to {tx_unsorted_file}.")

            # shifting can cause the file to become unsorted

            sorted_file = bam.sort(tx_unsorted_file, dd.get_config(data), force=True)

            shutil.move(sorted_file, tx_out_file)

    bam.index(out_file, dd.get_config(data))

    data["work_bam"] = out_file

    return data

import os

import toolz as tz

from collections import namedtuple

from bcbio import utils

from bcbio.pipeline import datadict as dd

from bcbio.pipeline import config_utils

from bcbio.log import logger

from bcbio.distributed.transaction import file_transaction

from bcbio.provenance import do

from bcbio import bam

# ranges taken from Buenrostro, Nat. Methods 10, 1213–1218 (2013).

ATACRange = namedtuple('ATACRange', ['label', 'min', 'max'])

ATACRanges = {"NF": ATACRange("NF", 0, 100),

              "MN": ATACRange("MN", 180, 247),

              "DN": ATACRange("DN", 315, 473),

              "TN": ATACRange("TN", 558, 615)}

def calculate_complexity_metrics(work_bam, data):

    """

    the work_bam should have duplicates marked but not removed

    mitochondrial reads should be removed 

    """

    bedtools = config_utils.get_program("bedtools", dd.get_config(data))

    work_dir = dd.get_work_dir(data)

    metrics_dir = os.path.join(work_dir, "metrics", "atac")

    utils.safe_makedir(metrics_dir)

    metrics_file = os.path.join(metrics_dir,

                                f"{dd.get_sample_name(data)}-atac-metrics.csv")

    if utils.file_exists(metrics_file):

        data = tz.assoc_in(data, ['atac', 'complexity_metrics_file'], metrics_file)

        return data

    # BAM file must be sorted by read name

    work_bam = bam.sort(work_bam, dd.get_config(data), order="queryname")

    with file_transaction(metrics_file) as tx_metrics_file:

        with open(tx_metrics_file, "w") as out_handle:

            out_handle.write("mt,m0,m1,m2\n")

        cmd = (f"{bedtools} bamtobed -bedpe -i {work_bam} | "

               "awk 'BEGIN{OFS=\"\\t\"}{print $1,$2,$4,$6,$9,$10}' | "

               "sort | "

               "uniq -c | "

               "awk 'BEGIN{mt=0;m0=0;m1=0;m2=0}($1==1){m1=m1+1} "

               "($1==2){m2=m2+1}{m0=m0+1}{mt=mt+$1}END{printf \"%d,%d,%d,%d\\n\", mt,m0,m1,m2}' >> "

               f"{tx_metrics_file}")

        message = f"Calculating ATAC-seq complexity metrics on {work_bam}, saving as {metrics_file}."

        do.run(cmd, message)

    data = tz.assoc_in(data, ['atac', 'complexity_metrics_file'], metrics_file)

    return data

def calculate_encode_complexity_metrics(data):

    metrics_file = tz.get_in(['atac', 'complexity_metrics_file'], data, None)

    if not metrics_file:

        return {}

    else:

        with open(metrics_file) as in_handle:

            header = next(in_handle).strip().split(",")

            values = next(in_handle).strip().split(",")

    raw_metrics = {h: int(v) for h, v in zip(header, values)}

    metrics = {"PBC1": raw_metrics["m1"] / raw_metrics["m0"],

               "NRF": raw_metrics["m0"] / raw_metrics["mt"]}

    if raw_metrics["m2"] == 0:

        PBC2 = 0

    else:

        PBC2 = raw_metrics["m1"] / raw_metrics["m2"]

    metrics["PBC2"] = PBC2

    return(metrics)

def split_ATAC(data, bam_file=None):

    """

    splits a BAM into nucleosome-free (NF) and mono/di/tri nucleosome BAMs based

    on the estimated insert sizes

    uses the current working BAM file if no BAM file is supplied

    """

    sambamba = config_utils.get_program("sambamba", data)

    num_cores = dd.get_num_cores(data)

    base_cmd = f'{sambamba} view --format bam --nthreads {num_cores} '

    bam_file = bam_file if bam_file else dd.get_work_bam(data)

    out_stem = os.path.splitext(bam_file)[0]

    split_files = {}

    for arange in ATACRanges.values():

        out_file = f"{out_stem}-{arange.label}.bam"

        if not utils.file_exists(out_file):

            with file_transaction(out_file) as tx_out_file:

                cmd = base_cmd +\

                    f'-F "template_length > {arange.min} and template_length < {arange.max}" ' +\

                    f'{bam_file} > {tx_out_file}'

                message = f'Splitting {arange.label} regions from {bam_file}.'

                do.run(cmd, message)

        split_files[arange.label] = out_file

    data = tz.assoc_in(data, ['atac', 'align'], split_files)

    return data