<commit_before>"""Structural variation detection for split and paired reads using lumpy.

Uses lumpyexpress for lumpy integration and samblaster for read preparation:
https://github.com/GregoryFaust/samblaster
https://github.com/arq5x/lumpy-sv
"""
import contextlib
import os
import sys
import shutil

import vcf

from bcbio import utils
from bcbio.distributed.transaction import file_transaction, tx_tmpdir
from bcbio.pipeline import datadict as dd
from bcbio.provenance import do
from bcbio.structural import shared as sshared
from bcbio.variation import effects, vcfutils, vfilter

# ## Lumpy main

def _run_lumpy(full_bams, sr_bams, disc_bams, work_dir, items):
    """Run lumpy-sv, using speedseq pipeline.
    """
    batch = sshared.get_cur_batch(items)
    ext = "-%s-svs" % batch if batch else "-svs"
    out_file = os.path.join(work_dir, "%s%s.vcf"
                            % (os.path.splitext(os.path.basename(items[0]["align_bam"]))[0], ext))
    sv_exclude_bed = sshared.prepare_exclude_file(items, out_file)
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            with tx_tmpdir(items[0]) as tmpdir:
                full_bams = ",".join(full_bams)
                sr_bams = ",".join(sr_bams)
                disc_bams = ",".join(disc_bams)
                exclude = "-x %s" % sv_exclude_bed if utils.file_exists(sv_exclude_bed) else ""
                ref_file = dd.get_ref_file(items[0])
                # use our bcbio python for runs within lumpyexpress
                curpython_dir = os.path.dirname(sys.executable)
                cmd = ("export PATH={curpython_dir}:$PATH && "
                       "lumpyexpress -v -B {full_bams} -S {sr_bams} -D {disc_bams} "
                       "{exclude} -T {tmpdir} -o {tx_out_file}")
                do.run(cmd.format(**locals()), "lumpyexpress", items[0])
    return vcfutils.sort_by_ref(out_file, items[0]), sv_exclude_bed

def _filter_by_support(in_file, data):
    """Filter call file based on supporting evidence, adding FILTER annotations to VCF.

    Filters based on the following criteria:
      - Minimum read support for the call (SU = total support)
      - Large calls need split read evidence.
    """
    rc_filter = ("FORMAT/SU < 4 || "
                 "(FORMAT/SR == 0 && ABS(SVLEN)>50000)")
    return vfilter.hard_w_expression(in_file, rc_filter, data, name="ReadCountSupport",
                                     limit_regions=None)

def _filter_by_background(base_samples, back_samples, gt_vcfs, data):
    """Filter base samples, marking any also present in the background.
    """
    filtname = "InBackground"
    filtdoc = "Variant also present in background samples with same genotype"
    for base_name in base_samples:
        orig_vcf = gt_vcfs[base_name]
        out_file = "%s-backfilter.vcf" % (utils.splitext_plus(orig_vcf)[0])
        if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
            with file_transaction(data, out_file) as tx_out_file:
                with utils.open_gzipsafe(orig_vcf) as in_handle:
                    with _vcf_readers([gt_vcfs[n] for n in back_samples]) as back_readers:
                        inp = vcf.Reader(in_handle, orig_vcf)
                        inp.filters[filtname] = vcf.parser._Filter(filtname, filtdoc)
                        with open(tx_out_file, "w") as out_handle:
                            outp = vcf.Writer(out_handle, inp)
                            for rec in inp:
                                back_recs = [r.next() for r in back_readers]
                                if _genotype_in_background(rec, back_recs):
                                    rec.add_filter(filtname)
                                outp.write_record(rec)
        gt_vcfs[base_name] = vcfutils.bgzip_and_index(out_file, data["config"])
    return gt_vcfs

def _genotype_in_background(rec, back_recs):
    """Check if the genotype in the record of interest is present in the background records.
    """
    def passes(rec):
        return not rec.FILTER or len(rec.FILTER) == 0
    return any([passes(brec) and passes(rec) and rec.samples[0].gt_alleles == brec.samples[0].gt_alleles
                for brec in back_recs])
<commit_msg>Support lumpy re-runs with pre-bgzipped background files<commit_after>"""Structural variation detection for split and paired reads using lumpy.

Uses lumpyexpress for lumpy integration and samblaster for read preparation:
https://github.com/GregoryFaust/samblaster
https://github.com/arq5x/lumpy-sv
"""
import contextlib
import os
import sys
import shutil

import vcf

from bcbio import utils
from bcbio.distributed.transaction import file_transaction, tx_tmpdir
from bcbio.pipeline import datadict as dd
from bcbio.provenance import do
from bcbio.structural import shared as sshared
from bcbio.variation import effects, vcfutils, vfilter

# ## Lumpy main

def _run_lumpy(full_bams, sr_bams, disc_bams, work_dir, items):
    """Run lumpy-sv, using speedseq pipeline.
    """
    batch = sshared.get_cur_batch(items)
    ext = "-%s-svs" % batch if batch else "-svs"
    out_file = os.path.join(work_dir, "%s%s.vcf"
                            % (os.path.splitext(os.path.basename(items[0]["align_bam"]))[0], ext))
    sv_exclude_bed = sshared.prepare_exclude_file(items, out_file)
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            with tx_tmpdir(items[0]) as tmpdir:
                full_bams = ",".join(full_bams)
                sr_bams = ",".join(sr_bams)
                disc_bams = ",".join(disc_bams)
                exclude = "-x %s" % sv_exclude_bed if utils.file_exists(sv_exclude_bed) else ""
                ref_file = dd.get_ref_file(items[0])
                # use our bcbio python for runs within lumpyexpress
                curpython_dir = os.path.dirname(sys.executable)
                cmd = ("export PATH={curpython_dir}:$PATH && "
                       "lumpyexpress -v -B {full_bams} -S {sr_bams} -D {disc_bams} "
                       "{exclude} -T {tmpdir} -o {tx_out_file}")
                do.run(cmd.format(**locals()), "lumpyexpress", items[0])
    return vcfutils.sort_by_ref(out_file, items[0]), sv_exclude_bed

def _filter_by_support(in_file, data):
    """Filter call file based on supporting evidence, adding FILTER annotations to VCF.

    Filters based on the following criteria:
      - Minimum read support for the call (SU = total support)
      - Large calls need split read evidence.
    """
    rc_filter = ("FORMAT/SU < 4 || "
                 "(FORMAT/SR == 0 && ABS(SVLEN)>50000)")
    return vfilter.hard_w_expression(in_file, rc_filter, data, name="ReadCountSupport",
                                     limit_regions=None)

def _filter_by_background(base_samples, back_samples, gt_vcfs, data):
    """Filter base samples, marking any also present in the background.
    """
    filtname = "InBackground"
    filtdoc = "Variant also present in background samples with same genotype"
    for base_name in base_samples:
        orig_vcf = gt_vcfs[base_name]
        out_file = "%s-backfilter.vcf" % (utils.splitext_plus(orig_vcf)[0])
        if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
            with file_transaction(data, out_file) as tx_out_file:
                with utils.open_gzipsafe(orig_vcf) as in_handle:
                    with _vcf_readers([gt_vcfs[n] for n in back_samples]) as back_readers:
                        inp = vcf.Reader(in_handle, orig_vcf)
                        inp.filters[filtname] = vcf.parser._Filter(filtname, filtdoc)
                        with open(tx_out_file, "w") as out_handle:
                            outp = vcf.Writer(out_handle, inp)
                            for rec in inp:
                                back_recs = [r.next() for r in back_readers]
                                if _genotype_in_background(rec, back_recs):
                                    rec.add_filter(filtname)
                                outp.write_record(rec)
        if utils.file_exists(out_file + ".gz"):
            out_file = out_file + ".gz"
        gt_vcfs[base_name] = vcfutils.bgzip_and_index(out_file, data["config"])
    return gt_vcfs

def _genotype_in_background(rec, back_recs):
    """Check if the genotype in the record of interest is present in the background records.
    """
    def passes(rec):
        return not rec.FILTER or len(rec.FILTER) == 0
    return any([passes(brec) and passes(rec) and rec.samples[0].gt_alleles == brec.samples[0].gt_alleles
                for brec in back_recs])

{'fraction_non_alphanumeric': 0.09973480917790845, 'fraction_numerical': 0.004381413582382105, 'mean_word_length': 2.4654414702357172, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6579908', 'n_tokens_mistral': 2812, 'n_tokens_neox': 2721, 'n_words': 676}