diff --git a/Dockerfile b/Dockerfile index 4518bf7..1df81e6 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ -FROM debian:13.4-slim AS downloader +FROM debian:13.7-slim AS downloader # Tool versions and pinned SHA-256 checksums. Bumping a version is a # two-line edit (version + hash) in this block. @@ -17,8 +17,8 @@ ARG VG_VERSION=1.73.0 ARG VG_SHA256=d5752977237e801d971f0d044f43fde4f3005da35f0451b4f452d1c1c9b8414b # Sentieon's KMC fork adds the `-fa /dev/stdin` patch used by the pangenome pipeline. ARG KMC_VERSION=3.2.4 -ARG KMC_RELEASE_TAG=v3.2.4-pipe2 -ARG KMC_SHA256=d8bdf8edcd0577dba32e86e8f194b2eb04eb168c75a3e6d46721d5bef515ac96 +ARG KMC_RELEASE_TAG=v3.2.4-pipe3 +ARG KMC_SHA256=6ae37ac28c215919a67bd77edf004f9e323013bba5700bc44ad23c0b59da9a91 ARG EXPANSIONHUNTER_VERSION=5.0.0 ARG EXPANSIONHUNTER_SHA256=ebf3ec0ace6e6e3bbce12c26463da5d9f8e16374eff1ad10f0f1a9123050fa86 # T1K and segdup-caller are built from git; pin the tag and verify the @@ -155,7 +155,7 @@ RUN /opt/poetry-venv/bin/pip install --no-cache-dir \ # Build the Python environment in an isolated stage. The venv is # later COPY'd over as-is. -FROM debian:13.4-slim AS python-builder +FROM debian:13.7-slim AS python-builder ARG SENTIEON_VERSION COPY --from=downloader /opt/sentieon-cli/dist /opt/sentieon-cli/dist @@ -185,7 +185,7 @@ RUN apt-get update && \ -depth -type d -name __pycache__ -exec rm -rf {} + # Build the container -FROM debian:13.4-slim +FROM debian:13.7-slim ARG SENTIEON_VERSION ARG BEDTOOLS_VERSION=2.30.0 ARG MOSDEPTH_VERSION=0.3.9 @@ -205,7 +205,7 @@ LABEL org.opencontainers.image.title="sentieon-cli" \ org.opencontainers.image.url="https://www.sentieon.com/" \ org.opencontainers.image.source="https://github.com/Sentieon/sentieon-cli" \ org.opencontainers.image.licenses="BSD-2-Clause" \ - org.opencontainers.image.base.name="debian:13.4-slim" \ + org.opencontainers.image.base.name="debian:13.7-slim" \ org.opencontainers.image.version="${SENTIEON_VERSION}" \ org.opencontainers.image.revision="${VCS_REF}" \ org.opencontainers.image.created="${BUILD_DATE}" diff --git a/README.md b/README.md index c5e56e8..36179bf 100644 --- a/README.md +++ b/README.md @@ -15,9 +15,9 @@ pip install sentieon_cli Install the sentieon-cli through a container image: ```sh -docker pull sentieon/sentieon_cli:latest +docker pull sentieon/sentieon-cli:latest docker run --rm -v "$PWD:/data" -w /data \ - sentieon/sentieon_cli:latest \ + sentieon/sentieon-cli:latest \ sentieon-cli dnascope-pangenome --help ``` diff --git a/pyproject.toml b/pyproject.toml index 0b2a465..30f9473 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "sentieon_cli" -version = "2.0.0" +version = "2.0.1" description = "Pipeline implementations for the Sentieon software" authors = [ {name = "Don Freed", email = "don.freed@sentieon.com"}, diff --git a/sentieon_cli/_indel2cnv_utils.py b/sentieon_cli/_indel2cnv_utils.py deleted file mode 100644 index f366cbd..0000000 --- a/sentieon_cli/_indel2cnv_utils.py +++ /dev/null @@ -1,113 +0,0 @@ -"""Pure-Python replacements for ``Levenshtein.ratio`` and -``scipy.signal.find_peaks`` used by ``scripts/indel2cnv.py``. - -Kept in the package (not ``scripts/``) so pytest can import them — -``scripts/`` is excluded from collection by ``conftest.py``. -""" - -from __future__ import annotations - -from collections import Counter -from typing import Sequence - - -def ratio(a: str, b: str, threshold: float = 0.0) -> float: - """Return ``2 * LCS(a, b) / (len(a) + len(b))``. - - Equivalent to ``Levenshtein.ratio(a, b)`` (which uses insert/delete - cost 1 and substitution cost 2, making the underlying distance - equal to ``|a| + |b| - 2 * LCS(a, b)``). - - When ``threshold > 0`` and a cheap upper bound proves the ratio - cannot reach ``threshold``, returns ``0.0`` as a sentinel without - running the full LCS. Otherwise returns the exact ratio. - """ - la, lb = len(a), len(b) - if la == 0 and lb == 0: - return 1.0 - if la == 0 or lb == 0: - return 0.0 - if a == b: - return 1.0 - - tot = la + lb - - # Length-based upper bound: 2*LCS <= 2*min(|a|,|b|). - if threshold > 0.0 and 2 * min(la, lb) < threshold * tot: - return 0.0 - - # Character-multiset intersection is a true upper bound on 2*LCS. - ub_matches = sum((Counter(a) & Counter(b)).values()) - if ub_matches == 0: - return 0.0 - if threshold > 0.0 and (2.0 * ub_matches) < threshold * tot: - return 0.0 - - return 2.0 * _bit_parallel_lcs(a, b) / tot - - -def _bit_parallel_lcs(a: str, b: str) -> int: - """LCS length via the Allison-Dix bit-parallel algorithm. - - Each iteration is O(len(a) / word_size) big-int operations; Python - int ops exploit this implicitly. Puts the shorter string in ``a`` - to minimize bitmap width. - """ - if len(a) > len(b): - a, b = b, a - m = len(a) - - char_mask: dict[str, int] = {} - for i, ch in enumerate(a): - char_mask[ch] = char_mask.get(ch, 0) | (1 << i) - - full = (1 << m) - 1 - s = full - for ch in b: - t = char_mask.get(ch) - if t is None: - continue - u = s & t - if u == 0: - continue - s = ((s + u) | (s - u)) & full - return m - bin(s).count("1") - - -def find_peaks( - values: Sequence[float], *, distance: int, height: float -) -> list[int]: - """Return indices of local maxima >= ``height`` with no two kept - peaks closer than ``distance`` samples. - - Matches ``scipy.signal.find_peaks(x, distance=d, height=h)`` for the - one call site that uses it. Plateaus map to their midpoint index - (``(start + end) // 2``), matching scipy's behavior. Boundary - values are never peaks. - """ - n = len(values) - candidates: list[int] = [] - i = 1 - while i < n - 1: - if values[i] <= values[i - 1] or values[i] < height: - i += 1 - continue - j = i - while j + 1 < n and values[j + 1] == values[i]: - j += 1 - if j + 1 < n and values[j + 1] < values[i]: - candidates.append((i + j) // 2) - i = j + 1 - - if distance <= 1 or len(candidates) <= 1: - return candidates - - # Keep highest peaks first; drop any within `distance` of one - # already kept. Tie-break on index so earlier wins. - ordered = sorted(candidates, key=lambda k: (-values[k], k)) - kept: list[int] = [] - for p in ordered: - if all(abs(p - q) >= distance for q in kept): - kept.append(p) - kept.sort() - return kept diff --git a/sentieon_cli/base_pangenome.py b/sentieon_cli/base_pangenome.py index 5b2f328..4e29f87 100644 --- a/sentieon_cli/base_pangenome.py +++ b/sentieon_cli/base_pangenome.py @@ -244,6 +244,8 @@ def __init__(self) -> None: self.t1k_kir_coord: Optional[pathlib.Path] = None self.t1k_kir_locus: str = DEFAULT_T1K_KIR_LOCUS self.has_cnv_model = False + # The bundle has a `combine_cnv.json` preset for `combine_cnv.py` + self.has_combine_preset = False # The pangenome reference, resolved from the graph's metadata by # `resolve_pangenome_reference`. The two arguments are hidden # overrides, so they stay `None` unless the user supplies them. @@ -789,10 +791,22 @@ def add_pangenome_cnv( sv_vcf=sv_vcf, sample_sex=self.sample_sex, par_bed=self.cnv_par_bed, + combine_preset=self.cnv_combine_preset(), interval=interval, replace_rg=replace_rg, ).add_to(dag) + def cnv_combine_preset(self) -> Optional[pathlib.Path]: + """The `combine_cnv.py` preset for the bundle's CNVscope model. + + The bundle's `combine_cnv.json` when it has one. Without it, + `combine_cnv.py` uses its defaults, the paired-end values. + """ + if not self.has_combine_preset: + return None + bundle = self.required(self.model_bundle, "model_bundle") + return bundle.joinpath("combine_cnv.json") + def build_kmc_job( self, kmer_prefix: pathlib.Path, job_threads: int ) -> Job: diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index c208463..deb6517 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -48,6 +48,26 @@ def cmd_bedtools_subtract( return Pipeline(Command(*cmd), file_output=unphased_bed) +def cmd_bedtools_intersect( + out_bed: pathlib.Path, + in_bed: pathlib.Path, + restrict_bed: pathlib.Path, + ref_fai: pathlib.Path, +) -> Pipeline: + """The parts of `in_bed` inside `restrict_bed`, sorted and merged""" + intersect_cmd = Command( + "bedtools", + "intersect", + "-a", + str(in_bed), + "-b", + str(restrict_bed), + ) + sort_cmd = Command("bedtools", "sort", "-faidx", str(ref_fai), "-i", "-") + merge_cmd = Command("bedtools", "merge") + return Pipeline(intersect_cmd, sort_cmd, merge_cmd, file_output=out_bed) + + def cmd_bedtools_slop( in_bed: pathlib.Path, out_bed: pathlib.Path, @@ -284,7 +304,14 @@ def cmd_pyexec_hybrid_anno( bed: pathlib.Path, hybrid_anno: pathlib.Path, threads: int, + stream: bool = False, ) -> Pipeline: + """Annotate the hybrid calls with their long-read haplotype count + + With `stream`, the script writes uncompressed VCF that + `sentieon util vcfconvert` compresses and indexes, rather than + compressing and indexing every record in its own parent process. + """ cmd = [ sys.executable, str(hybrid_anno), @@ -294,9 +321,49 @@ def cmd_pyexec_hybrid_anno( str(bed), "-t", str(threads), - str(out_vcf), ] - return Pipeline(Command(*cmd)) + if not stream: + return Pipeline(Command(*cmd, str(out_vcf))) + return Pipeline( + Command(*cmd, "-"), + Command("sentieon", "util", "vcfconvert", "-", str(out_vcf)), + ) + + +def cmd_pyexec_gvcf_clip( + out_gvcf: pathlib.Path, + gvcf: pathlib.Path, + regions_bed: pathlib.Path, + reference: pathlib.Path, + gvcf_clip: pathlib.Path, + threads: int, + exclude: bool = False, +) -> Pipeline: + """Restrict a gVCF to the BED regions, or with `exclude` to the + positions outside them, clipping reference blocks at the boundaries + + The script writes uncompressed VCF that `sentieon util vcfconvert` + compresses and indexes, rather than compressing and indexing every + record in its own parent process. + """ + cmd = [ + sys.executable, + str(gvcf_clip), + "-v", + str(gvcf), + "-b", + str(regions_bed), + "-r", + str(reference), + "-t", + str(threads), + ] + if exclude: + cmd.append("--exclude") + return Pipeline( + Command(*cmd, "-"), + Command("sentieon", "util", "vcfconvert", "-", str(out_gvcf)), + ) def cmd_pyexec_sad_lad_update( @@ -401,8 +468,18 @@ def cmd_pyexec_combine_cnv( cnv_vcf: pathlib.Path, converted_vcf: pathlib.Path, combine_script: pathlib.Path, + *, + preset: Optional[pathlib.Path] = None, + sex: Optional[str] = None, + par: Optional[pathlib.Path] = None, + blocks: Optional[pathlib.Path] = None, ) -> Pipeline: - """Combine CNVscope and converted SV calls""" + """Combine CNVscope and converted SV calls. + + `preset` is the JSON file of rule overrides; without it the script + uses its defaults. `sex` and `par`, or `blocks`, are the ploidy + CNVscope was called with. + """ cmd = [ sys.executable, str(combine_script), @@ -410,9 +487,16 @@ def cmd_pyexec_combine_cnv( str(cnv_vcf), "--converted", str(converted_vcf), - "-o", - str(out_vcf), ] + if preset: + cmd.extend(["--preset", str(preset)]) + if sex: + cmd.extend(["--sex", sex]) + if par: + cmd.extend(["--par", str(par)]) + if blocks: + cmd.extend(["--blocks", str(blocks)]) + cmd.extend(["-o", str(out_vcf)]) return Pipeline(Command(*cmd)) @@ -563,6 +647,27 @@ def bcftools_subset( return Pipeline(Command("sh", "-c", shell_cmd)) +def bcftools_remove_annotations( + out_vcf: pathlib.Path, + in_vcf: pathlib.Path, + annotations: List[str], +) -> Pipeline: + """Remove annotations, such as `FORMAT/LAD`, with bcftools annotate""" + return Pipeline( + Command( + "bcftools", + "annotate", + "-x", + ",".join(annotations), + "-Oz", + "-W=tbi", + "-o", + str(out_vcf), + str(in_vcf), + ) + ) + + def bcftools_concat( out_vcf: pathlib.Path, in_vcfs: List[pathlib.Path], @@ -585,14 +690,9 @@ def filter_norm( out_vcf: pathlib.Path, in_vcf: pathlib.Path, reference: pathlib.Path, - exclude_homref: bool = True, ) -> Pipeline: """Trim and normalize""" - view_args = ["-a"] - if exclude_homref: - view_args.extend(["-e", 'GT="0/0"']) - view_args.append(str(in_vcf)) - view_cmd = Command("bcftools", "view", *view_args) + view_cmd = Command("bcftools", "view", "-a", "-e", 'GT="0/0"', str(in_vcf)) norm_cmd = Command("bcftools", "norm", "-f", str(reference)) convert_cmd = Command("sentieon", "util", "vcfconvert", "-", str(out_vcf)) return Pipeline(view_cmd, norm_cmd, convert_cmd) @@ -987,6 +1087,10 @@ def cmd_mosdepth( return Pipeline(Command(*cmd)) +# KMC fails when run with more than 128 threads +KMC_MAX_THREADS = 128 + + def cmd_kmc( output_prefix: pathlib.Path, file_list: pathlib.Path, @@ -1000,7 +1104,7 @@ def cmd_kmc( f"-k{k}", f"-m{memory}", "-okff", - f"-t{threads}", + f"-t{min(threads, KMC_MAX_THREADS)}", "-hp", f"@{file_list}", str(output_prefix), @@ -1274,8 +1378,15 @@ def cmd_bwa_extract( threads=1, unzip: str = "gzip", bwa_k: str = "2000000", + numa: Optional[str] = None, + split: Optional[str] = None, ) -> Pipeline: - """BWA alignment with pgutil extract""" + """BWA alignment with pgutil extract + + With ``split`` ("j/n") the job aligns one share of the read pairs, read + interleaved through `fqidx extract` with the same -K, so the batches and + the alignments match an unsplit run; ``numa`` pins it with `taskset`. + """ fq1_cmd = Pipeline( Command( unzip, @@ -1291,7 +1402,30 @@ def cmd_bwa_extract( ) ) + pin = ["taskset", "-c", numa] if numa else [] + if split: + reads = [ + InputProcSub( + Pipeline( + Command( + "sentieon", + "fqidx", + "extract", + "-F", + split, + "-K", + bwa_k, + InputProcSub(fq1_cmd), + InputProcSub(fq2_cmd), + ) + ) + ) + ] + else: + reads = [InputProcSub(fq1_cmd), InputProcSub(fq2_cmd)] + bwa_cmd = Command( + *pin, "sentieon", "bwa", "mem", @@ -1303,12 +1437,13 @@ def cmd_bwa_extract( str(bwa_model), "-K", bwa_k, + *(["-p"] if split else []), str(reference), - InputProcSub(fq1_cmd), - InputProcSub(fq2_cmd), + *reads, ) extract_cmd = Command( + *pin, "sentieon", "pgutil", "extract", @@ -1390,9 +1525,12 @@ def cmd_minimap2_lift( threads=1, mm2_xargs: List[str] = [], lift_prefix: str = "GRCh38#0#", + numa: Optional[str] = None, ) -> Pipeline: - """Minimap2 alignment with pgutil lift""" + """Minimap2 alignment with pgutil lift; ``numa`` pins it with `taskset`""" + pin = ["taskset", "-c", numa] if numa else [] mm2_cmd = Command( + *pin, "sentieon", "minimap2", "-t", @@ -1411,6 +1549,7 @@ def cmd_minimap2_lift( ) lift_cmd = Command( + *pin, "sentieon", "pgutil", "lift", @@ -1425,6 +1564,7 @@ def cmd_minimap2_lift( ) sort_cmd = Command( + *pin, "sentieon", "util", "sort", @@ -1480,7 +1620,7 @@ def _kmc_stdin_cmd( f"-k{k}", f"-m{memory}", "-okff", - f"-t{threads}", + f"-t{min(threads, KMC_MAX_THREADS)}", "-fa", "/dev/stdin", str(output_prefix), @@ -1618,27 +1758,70 @@ def cmd_samtools_faidx( return Pipeline(Command("samtools", "faidx", str(fasta))) +# An awk program that drops every `##contig` header line except the one +# for the `contig` variable, whether its ID comes first or not +KEEP_CONTIG_AWK = ( + "/^##contig=/ { " + 'if (index($0, "") ' + '|| index($0, ",ID=" contig ",") || index($0, ",ID=" contig ">")) ' + "print; next } 1" +) + + def cmd_bcftools_merge_trim( shard_vcf: pathlib.Path, raw_vcf: pathlib.Path, pop_vcf: pathlib.Path, trim_script: pathlib.Path, subset_bed: pathlib.Path, + contig: str, merge_rules: str, merge_xargs: List[str] = [], view_xargs: List[str] = [], ): - """Merge and trim variants""" + """Merge and trim the variants of one shard, on `contig` + + `bcftools merge` sets QUAL to the highest QUAL of the merged records + and, under its default `--filter-logic +`, combines their FILTERs. The + population VCF's QUAL is cleared before the merge (the merge ignores + a missing QUAL), so the output keeps the input VCF's, and FILTER is + cleared from both inputs. Each input is restricted to `subset_bed` as + it is read, since `-R` cannot be combined with `--no-index`. + + `--no-index` refuses inputs whose headers list different contigs, or + the same contigs in a different order, as the input VCF (the + reference's contigs) and a population VCF generally do. Both inputs + keep only the `##contig` line of the shard's contig, so the shard VCF + lists that contig alone; concatenating the shards in reference order + lists the contigs in that order again. + """ + + def subset(vcf: pathlib.Path, remove: str) -> Pipeline: + return Pipeline( + Command( + "bcftools", + "annotate", + "--no-version", + "--regions-file", + str(subset_bed), + "--regions-overlap", + "pos", + "-x", + remove, + str(vcf), + ), + Command("awk", "-v", f"contig={contig}", KEEP_CONTIG_AWK), + ) + merge_cmd = Command( "bcftools", "merge", - "--regions-file", - str(subset_bed), + "--no-index", *merge_xargs, "-i", merge_rules, - str(raw_vcf), - str(pop_vcf), + InputProcSub(subset(raw_vcf, "FILTER")), + InputProcSub(subset(pop_vcf, "QUAL,FILTER")), ) trim_cmd = Command( @@ -1658,28 +1841,28 @@ def cmd_bcftools_merge_trim( return Pipeline(merge_cmd, trim_cmd, view_cmd) -def cmd_bcftools_view_regions( +def cmd_bcftools_subset_remove( out_vcf: pathlib.Path, in_vcf: pathlib.Path, - regions: Optional[str] = None, - regions_file: Optional[pathlib.Path] = None, -): - xargs = [] - if regions: - xargs.extend(["--regions", regions]) - if regions_file: - xargs.extend(["--regions-file", str(regions_file)]) + regions_file: pathlib.Path, + remove: str, +) -> Pipeline: + """Restrict a VCF to `regions_file` and remove the `remove` + annotations, such as `FILTER`""" return Pipeline( Command( "bcftools", - "view", + "annotate", "--no-version", + "--regions-file", + str(regions_file), + "-x", + remove, "-W=tbi", "-O", "z", "-o", str(out_vcf), - *xargs, str(in_vcf), ) ) diff --git a/sentieon_cli/dnascope.py b/sentieon_cli/dnascope.py index c0c6006..87af165 100644 --- a/sentieon_cli/dnascope.py +++ b/sentieon_cli/dnascope.py @@ -5,14 +5,16 @@ import argparse import copy import itertools +import json import os import pathlib import shutil import sys -from typing import List, Optional, Set +from typing import Dict, List, Optional, Set import packaging.version +from .archive import ar_load from .dag import DAG from .driver import ( BaseAlgo, @@ -21,6 +23,7 @@ ) from .job import Job from .pipeline import BasePipeline +from .shard import Shard, determine_shards_from_fai, parse_fai, vcf_contigs from .stages.alignment import ( BWA_FASTQ_MIN_VERSIONS, BWA_REALIGN_MIN_VERSIONS, @@ -39,7 +42,9 @@ DNAscopeStage, GVCFtyperStage, TransferApplyStage, + TransferSpec, ) +from .stages.transfer import TransferConfig from .util import ( library_preloaded, parse_rg_line, @@ -48,12 +53,17 @@ set_bwt_max_mem, split_alignment, total_memory, + vcf_id, ) VARIANTS_MIN_VERSIONS = { "sentieon driver": packaging.version.Version("202308"), } +TRANSFER_MIN_VERSIONS = { + "bcftools": packaging.version.Version("1.22"), +} + class DNAscopePipeline(BasePipeline): """The DNAscope pipeline""" @@ -163,6 +173,13 @@ class DNAscopePipeline(BasePipeline): "help": "Use arguments for PCR-free data processing", "action": "store_true", }, + "pop_vcf": { + "flags": ["--pop_vcf"], + "help": ( + "A VCF containing annotations for use with DNAModelApply." + ), + "type": path_arg(exists=True, is_file=True), + }, "skip_metrics": { "help": "Skip all metrics collection and multiQC", "action": "store_true", @@ -202,6 +219,10 @@ class DNAscopePipeline(BasePipeline): "help": argparse.SUPPRESS, "action": "store_true", }, + "skip_pop_vcf_id_check": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, "util_sort_args": { # "help": "Extra arguments for sentieon util sort", "help": argparse.SUPPRESS, @@ -221,6 +242,7 @@ def __init__(self) -> None: self.bed: Optional[pathlib.Path] = None self.interval_padding: Optional[int] = None self.pcr_free = False + self.pop_vcf: Optional[pathlib.Path] = None self.gvcf = False self.duplicate_marking = "markdup" self.assay = "WGS" @@ -239,6 +261,15 @@ def __init__(self) -> None: self.bwt_max_mem: Optional[str] = None self.no_ramdisk = False self.no_split_alignment = False + self.skip_pop_vcf_id_check = False + # The annotation transfer's inputs, set by `validate` when a + # `--pop_vcf` is supplied + self.fai_data: Dict[str, Dict[str, int]] = {} + self.pop_vcf_contigs: Dict[str, Optional[int]] = {} + self.shards: List[Shard] = [] + # The bundle's sequencing platform, upper-cased, set by + # `validate_bundle` from its `platform` or `shortReadPlatform` + self.tech = "" def validate(self) -> None: self.required(self.output_vcf, "output_vcf") @@ -288,6 +319,61 @@ def validate(self) -> None: if self.r1_fastq or self.align or self.collate_align: self.validate_bwa_index() + if self.pop_vcf: + self.fai_data = parse_fai( + pathlib.Path(str(self.reference) + ".fai") + ) + self.pop_vcf_contigs = vcf_contigs(self.pop_vcf) + self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) + self.shards = determine_shards_from_fai( + self.fai_data, 10 * 1000 * 1000 + ) + self.validate_bundle() + + def validate_bundle(self) -> None: + """Read the bundle's platform and check the `--pop_vcf` against it""" + bundle = self.required(self.model_bundle, "model_bundle") + bundle_info_bytes = ar_load(str(bundle) + "/bundle_info.json") + if isinstance(bundle_info_bytes, list): + bundle_info_bytes = b"{}" + bundle_info = json.loads(bundle_info_bytes.decode()) + # Bundles name the platform `platform` or, like the hybrid and + # pangenome bundles, `shortReadPlatform` + self.tech = ( + bundle_info.get("platform") + or bundle_info.get("shortReadPlatform") + or "" + ).upper() + bundle_vcf_id = bundle_info.get("SentieonVcfID") + + if not bundle_vcf_id: + if self.pop_vcf: + self.logger.error( + "The model bundle does not require a population VCF. " + "Please use this model bundle without the `--pop_vcf` " + "argument." + ) + sys.exit(2) + return + + if not self.pop_vcf: + self.logger.error( + "The model bundle requires a population VCF. Please supply " + "the `--pop_vcf` argument." + ) + sys.exit(2) + if self.skip_pop_vcf_id_check: + return + pop_vcf_id = vcf_id(self.pop_vcf) + if pop_vcf_id != bundle_vcf_id: + self.logger.error( + "The population VCF provided does not match the population " + "VCF required by the model bundle. Expected: %s, Found: %s", + bundle_vcf_id, + pop_vcf_id, + ) + sys.exit(2) + def configure(self) -> None: self.configure_alignment() @@ -474,6 +560,10 @@ def add_small_variant_calling( bundle = self.required(self.model_bundle, "model_bundle") require_versions(VARIANTS_MIN_VERSIONS, skip=self.skip_version_check) + if self.pop_vcf: + require_versions( + TRANSFER_MIN_VERSIONS, skip=self.skip_version_check + ) out_gvcf = pathlib.Path( str(ctx.output_vcf).replace(".vcf.gz", ".g.vcf.gz") @@ -484,6 +574,9 @@ def add_small_variant_calling( out_svs = pathlib.Path( str(ctx.output_vcf).replace(".vcf.gz", "_svs.vcf.gz") ) + read_filters: List[str] = [] + if self.tech == "ULTIMA": + read_filters.append("UltimaReadFilter") pcr_indel_model = "NONE" if self.pcr_free else "CONSERVATIVE" model = bundle.joinpath("dnascope.model") @@ -523,12 +616,26 @@ def add_small_variant_calling( inputs=deduped, interval=self.bed, interval_padding=self.interval_padding, + read_filter=read_filters, ).add_to(dag, upstream) - # Genotyping and filtering with DNAModelApply + # Transfer annotations from the pop_vcf, then genotype and filter + # with DNAModelApply + transfer: Optional[TransferSpec] = None + if self.pop_vcf: + transfer_vcf = ctx.tmp_dir.joinpath( + "sample-dnascope_transfer.g.vcf.gz" + if self.gvcf + else "sample-dnascope_transfer.vcf.gz" + ) + transfer = TransferSpec( + config=TransferConfig.from_pipeline(self), + out_vcf=transfer_vcf, + ) transfer_apply = TransferApplyStage( ctx=ctx, raw_vcf=tmp_vcf, + transfer=transfer, apply=ApplySpec(model=model, output=ds_out), ).add_to(dag, call.terminal) @@ -545,6 +652,7 @@ def add_small_variant_calling( gvcf=out_gvcf, output=ctx.output_vcf, interval=self.bed, + dbsnp=self.dbsnp, ).add_to(dag, transfer_apply.terminal) # Call SVs diff --git a/sentieon_cli/dnascope_hybrid.py b/sentieon_cli/dnascope_hybrid.py index 8ad702b..0d671f9 100644 --- a/sentieon_cli/dnascope_hybrid.py +++ b/sentieon_cli/dnascope_hybrid.py @@ -18,6 +18,7 @@ from . import command_strings as cmds from .dag import DAG from .driver import ( + BaseAlgo, Driver, DNAscope, HybridStage1, @@ -46,6 +47,7 @@ from .stages.small_variants import ( ApplySpec, DNAscopeStage, + GVCFtyperStage, TransferApplyStage, TransferSpec, ) @@ -140,6 +142,17 @@ def __init__( ) +# Per-allele FORMAT fields of the hybrid model that GVCFtyper copies +# unchanged from the gVCF, still sized for and any allele it +# drops. They are removed from the VCF until GVCFtyper subsets them. +GVCFTYPER_UNSUBSET_FIELDS = [ + "FORMAT/LAD", + "FORMAT/LPL", + "FORMAT/SAD", + "FORMAT/SPL", +] + + class DNAscopeHybridPipeline(BasePipeline): """The DNAscope Hybrid pipeline""" @@ -872,21 +885,42 @@ def add_variant_calling( hybrid_model = bundle.joinpath("hybrid.model") - # First pass - combined variant calling + # With --gvcf, both calling passes emit gVCFs, which are merged, + # annotated and filtered as gVCFs into the output gVCF. GVCFtyper + # then genotypes the output gVCF to produce the VCF. vcf_suffix = ".g.vcf.gz" if self.gvcf else ".vcf.gz" - combined_vcf = self.tmp_dir.joinpath("initial" + vcf_suffix) emit_mode = "gvcf" if self.gvcf else None - call = DNAscopeStage( - ctx=ctx, - algos=[ + out_gvcf = pathlib.Path( + str(ctx.output_vcf).replace(".vcf.gz", ".g.vcf.gz") + ) + small_variants_out = out_gvcf if self.gvcf else ctx.output_vcf + + # First pass - combined variant calling. Region selection reads a + # VCF-mode call set, so with --gvcf the same driver run also writes + # one with a second DNAscope algo. + combined_vcf = self.tmp_dir.joinpath("initial" + vcf_suffix) + select_vcf = self.tmp_dir.joinpath("initial.vcf.gz") + pass1_algos: List[BaseAlgo] = [ + DNAscope( + combined_vcf, + dbsnp=self.dbsnp, + model=hybrid_model, + pcr_indel_model="none", + emit_mode=emit_mode, + ) + ] + if self.gvcf: + pass1_algos.append( DNAscope( - combined_vcf, + select_vcf, dbsnp=self.dbsnp, model=hybrid_model, pcr_indel_model="none", - emit_mode=emit_mode, ) - ], + ) + call = DNAscopeStage( + ctx=ctx, + algos=pass1_algos, inputs=lr_aln + sr_aln, interval=self.bed, read_filter=rg_info.ultima_read_filter @@ -904,7 +938,7 @@ def add_variant_calling( select_job = Job( cmds.cmd_pyexec_hybrid_select( out_bed=selected_bed, - vcf=combined_vcf, + vcf=select_vcf, ref_fai=ref_fai, hybrid_select=hybrid_select, threads=self.cores, @@ -1119,8 +1153,28 @@ def add_variant_calling( {third_stage_job}, ) + # The stage-2 BED covers the whole genome, so with --bed the second + # pass and the merge use its intersection with --bed. Stage 3 and + # the realigned short-read alignment keep the whole stage-2 BED. + recall_bed = stage2_bed + call2_upstream: Set[Job] = {third_stage_job} + merge_upstream: Set[Job] = {second_stage_job} + if self.bed: + recall_bed = self.tmp_dir.joinpath("hybrid_recall.bed") + intersect_job = Job( + cmds.cmd_bedtools_intersect( + recall_bed, stage2_bed, self.bed, ref_fai + ), + "intersect-stage2-bed", + 0, + task_name="region-selection", + ) + dag.add_job(intersect_job, {second_stage_job}) + call2_upstream.add(intersect_job) + merge_upstream = {intersect_job} + # pass 2 of variant calling - pass2_vcf = self.tmp_dir.joinpath("hybrid_pass2.vcf.gz") + pass2_vcf = self.tmp_dir.joinpath("hybrid_pass2" + vcf_suffix) call2 = DNAscopeStage( ctx=ctx, algos=[ @@ -1133,39 +1187,94 @@ def add_variant_calling( ) ], inputs=lr_aln + [stage3_aln], - interval=stage2_bed, + interval=recall_bed, read_filter=rg_info.ultima_read_filter + rg_info.lr_rg_read_filter, replace_rg=rg_info.replace_rg_args[0], name="dnascope-2", - ).add_to(dag, {third_stage_job}) - - # Merge and normalize the VCFs - subset_vcf = self.tmp_dir.joinpath("mix_subset.vcf.gz") - combined_tmp_vcf = self.tmp_dir.joinpath("combined_tmp.vcf.gz") - subset_job = Job( - cmds.bcftools_subset( - subset_vcf, + ).add_to(dag, call2_upstream) + + # Merge the calls, taking the second pass within the realigned + # regions and the first pass elsewhere + subset_vcf = self.tmp_dir.joinpath("mix_subset" + vcf_suffix) + combined_tmp_vcf = self.tmp_dir.joinpath("combined_tmp" + vcf_suffix) + if self.gvcf: + # Reference blocks from either pass can cross the region + # boundaries, so both gVCFs are clipped at the boundaries + gvcf_clip = pathlib.Path( + str(files("sentieon_cli.scripts").joinpath("gvcf_clip.py")) + ).resolve() + subset_job = Job( + cmds.cmd_pyexec_gvcf_clip( + subset_vcf, + combined_vcf, + recall_bed, + ctx.reference, + gvcf_clip, + self.cores, + exclude=True, + ), + "subset-calls", + 0, + task_name="vcf-merge", + ) + dag.add_job(subset_job, merge_upstream) + pass2_subset_vcf = self.tmp_dir.joinpath( + "hybrid_pass2_subset.g.vcf.gz" + ) + pass2_subset_job = Job( + cmds.cmd_pyexec_gvcf_clip( + pass2_subset_vcf, + pass2_vcf, + recall_bed, + ctx.reference, + gvcf_clip, + self.cores, + ), + "subset-pass2-calls", + 0, + task_name="vcf-merge", + ) + dag.add_job(pass2_subset_job, call2.terminal) + concat_inputs = [subset_vcf, pass2_subset_vcf] + concat_upstream = {subset_job, pass2_subset_job} + merge_tmp = [ combined_vcf, - stage2_bed, - ), - "subset-calls", - 0, - task_name="vcf-merge", - ) - dag.add_job(subset_job, {second_stage_job}) + select_vcf, + subset_vcf, + pass2_vcf, + pass2_subset_vcf, + ] + else: + subset_job = Job( + cmds.bcftools_subset( + subset_vcf, + combined_vcf, + recall_bed, + ), + "subset-calls", + 0, + task_name="vcf-merge", + ) + dag.add_job(subset_job, merge_upstream) + concat_inputs = [subset_vcf, pass2_vcf] + concat_upstream = {subset_job} | call2.terminal + merge_tmp = [combined_vcf, subset_vcf, pass2_vcf] + # The gVCFs are large enough for bcftools to use threads + concat_args = ["-aD"] + if self.gvcf: + concat_args += ["--threads", str(self.cores)] concat_job = Job( cmds.bcftools_concat( - combined_tmp_vcf, - [subset_vcf, pass2_vcf], + combined_tmp_vcf, concat_inputs, xargs=concat_args ), "concat-calls", 0, task_name="vcf-merge", ) - dag.add_job(concat_job, {subset_job} | call2.terminal) + dag.add_job(concat_job, concat_upstream) self.add_cleanup( dag, - [combined_vcf, subset_vcf, pass2_vcf], + merge_tmp, "rm-tmp5", {concat_job}, ) @@ -1174,9 +1283,9 @@ def add_variant_calling( hybrid_anno = pathlib.Path( str(files("sentieon_cli.scripts").joinpath("hybrid_anno.py")) ) - anno_target = self.tmp_dir.joinpath("combined_tmp_anno.vcf.gz") + anno_target = self.tmp_dir.joinpath("combined_tmp_anno" + vcf_suffix) if self.skip_model_apply and not self.pop_vcf: - anno_target = ctx.output_vcf + anno_target = small_variants_out anno_job = Job( cmds.cmd_pyexec_hybrid_anno( @@ -1185,54 +1294,95 @@ def add_variant_calling( stage1_hap_bed, hybrid_anno, self.cores, + stream=self.gvcf, ), "anno-calls", 0, task_name="annotation", ) dag.add_job(anno_job, {concat_job}) - - if not self.pop_vcf and self.skip_model_apply: - # `anno_target` is already the output VCF - return + last_jobs: Set[Job] = {anno_job} # Transfer annotations and apply the model - transfer_spec: Optional[TransferSpec] = None - if self.pop_vcf: - transfer_target = self.tmp_dir.joinpath( - "combined_tmp_transfer.vcf.gz" - ) - if self.skip_model_apply: - transfer_target = ctx.output_vcf - transfer_spec = TransferSpec( - TransferConfig.from_pipeline(self), transfer_target - ) + apply_vcf = ( + out_gvcf + if self.gvcf + else self.tmp_dir.joinpath("combined_apply.vcf.gz") + ) + if self.pop_vcf or not self.skip_model_apply: + transfer_spec: Optional[TransferSpec] = None + if self.pop_vcf: + transfer_target = self.tmp_dir.joinpath( + "combined_tmp_transfer" + vcf_suffix + ) + if self.skip_model_apply: + transfer_target = small_variants_out + transfer_spec = TransferSpec( + TransferConfig.from_pipeline(self), transfer_target + ) - apply_spec: Optional[ApplySpec] = None - apply_vcf = self.tmp_dir.joinpath("combined_apply.vcf.gz") - if not self.skip_model_apply: - apply_spec = ApplySpec(hybrid_model, apply_vcf) + apply_spec: Optional[ApplySpec] = None + if not self.skip_model_apply: + apply_spec = ApplySpec(hybrid_model, apply_vcf) - transfer_apply = TransferApplyStage( - ctx=ctx, - raw_vcf=anno_target, - transfer=transfer_spec, - apply=apply_spec, - ).add_to(dag, {anno_job}) + last_jobs = ( + TransferApplyStage( + ctx=ctx, + raw_vcf=anno_target, + transfer=transfer_spec, + apply=apply_spec, + ) + .add_to(dag, last_jobs) + .terminal + ) + + # Genotype the output gVCF to produce the VCF. GVCFtyper copies the + # hybrid per-allele fields unchanged from the gVCF, so they are + # removed from its output. That job writes the input of the final + # normalization, or under --skip_model_apply, which skips that + # step, the output VCF itself. + norm_input = apply_vcf + if self.gvcf: + gvcftyper_vcf = self.tmp_dir.joinpath("combined_gvcftyper.vcf.gz") + typer = GVCFtyperStage( + ctx=ctx, + gvcf=out_gvcf, + output=gvcftyper_vcf, + interval=self.bed, + dbsnp=self.dbsnp, + ).add_to(dag, last_jobs) + + norm_input = ( + ctx.output_vcf + if self.skip_model_apply + else self.tmp_dir.joinpath("combined_typed.vcf.gz") + ) + drop_job = Job( + cmds.bcftools_remove_annotations( + norm_input, + gvcftyper_vcf, + GVCFTYPER_UNSUBSET_FIELDS, + ), + "drop-allele-fields", + 0, + task_name="gvcftyper", + ) + dag.add_job(drop_job, typer.terminal) + last_jobs = {drop_job} if self.skip_model_apply: + # The output VCF is written without the final normalization return # Final normalize norm_job = Job( cmds.filter_norm( ctx.output_vcf, - apply_vcf, + norm_input, ctx.reference, - exclude_homref=not self.gvcf, ), "final-norm", 0, task_name="vcf-norm", ) - dag.add_job(norm_job, transfer_apply.terminal) + dag.add_job(norm_job, last_jobs) diff --git a/sentieon_cli/driver.py b/sentieon_cli/driver.py index 53b6131..1ddd13a 100644 --- a/sentieon_cli/driver.py +++ b/sentieon_cli/driver.py @@ -239,9 +239,11 @@ def __init__( self, output: pathlib.Path, vcf: pathlib.Path, + dbsnp: Optional[pathlib.Path] = None, ): self.output = output self.vcf = vcf + self.dbsnp = dbsnp class SVSolver(BaseAlgo): @@ -500,12 +502,14 @@ def __init__( *, sex: Optional[str] = None, par: Optional[pathlib.Path] = None, + blocks: Optional[pathlib.Path] = None, dump_probes: Optional[pathlib.Path] = None, ): self.output = output self.model = model self.sex = sex self.par = par + self.blocks = blocks self.dump_probes = dump_probes @@ -519,10 +523,13 @@ def __init__( output: Union[pathlib.Path, str], model: pathlib.Path, vcf: pathlib.Path, + *, + emit_conf: Optional[float] = None, ): self.output = output self.model = model self.vcf = vcf + self.emit_conf = emit_conf class PangenomeSV(BaseAlgo): diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index d3ae092..0968ac4 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -457,11 +457,12 @@ def validate_bundle(self) -> None: self.extract_model_name = extract_candidate required_members = { - "dnascope.model", "longreadsv.model", "minimap2.model", self.extract_model_name, } + if not self.skip_small_variants: + required_members.add("dnascope.model") if self.r1_fastq: required_members.add("bwa.model") if self.lr_align_input: @@ -479,6 +480,7 @@ def validate_bundle(self) -> None: sys.exit(2) self.has_cnv_model = "cnv.model" in bundle_members + self.has_combine_preset = "combine_cnv.json" in bundle_members bundle_vcf_id = bundle_info.get("SentieonVcfID") if bundle_vcf_id and not self.skip_pop_vcf_id_check: diff --git a/sentieon_cli/scripts/combine_cnv.py b/sentieon_cli/scripts/combine_cnv.py index eca4372..263c1de 100644 --- a/sentieon_cli/scripts/combine_cnv.py +++ b/sentieon_cli/scripts/combine_cnv.py @@ -1,164 +1,725 @@ #!/usr/bin/env python3 """ -Combine CNVscope calls with indel2cnv-converted SV calls (both in CNV space). +Combine CNVscope calls with SV calls converted to CNV by indel2cnv.py into one CNV VCF. -The converted CNV calls have REFSTART/REFSTOP defining the repeat array extent. -CNVscope loss calls overlapping these ranges are marked FILTER=SVdup (converted -call is preferred — it has precise breakpoints from the SV caller). -Gain calls are never marked. - -Both call sets use CN for evaluation: CNVscope has CN natively, converted calls -get CN = neutral_cn + CNDIFF. - -Usage: - python3 combine_cnv.py \\ - --cnv test/pred.HG002_30x_hs38.vcf.gz \\ - --converted test/hg002-svcaller.cnv.vcf.gz \\ - -o test/combined_cnv.vcf.gz +The rule values are DEFAULTS, those for paired-end CNVscope models, updated by the JSON file given with --preset +(presets/SE.json for single-end models), a file or a member of a Sentieon model bundle (MODEL.bundle/combine_cnv.json). +- Converted calls are written with SOURCE=SV_CNV, CN = CN_NEUTRAL + CNDIFF (floored at 0; on haploid sequence a + homozygous call counts one copy) and FILTER PASS. Dropped: converted losses longer than max_converted_loss and, + on haploid sequence, non-homozygous losses and gains longer than max_converted_gain_haploid. +- A converted gain whose array is at least nodepth_size long, whose array-average CN shift reaches + nodepth_min_shift, and over which CNVscope's raw segmentation shows depth but almost no gain, gets FILTER NoDepth. + The segments are the records of --raw (CNVscope's output) or, without --raw, every record of --cnv, which + CNVModelApply then wrote with --emit_conf 0. The gain evidence is the segments' HMM state (HMMCN) where they carry + it, else the PASS calls of --cnv. +- A CNVscope PASS call that repeats a kept converted call gets FILTER SVdup (INFO DEDUP names the array). +- With loss_support_size set (SE), a shorter CNVscope PASS loss that no converted loss overlaps gets FILTER NoSV. +- Other CNVscope PASS calls are written with SOURCE=CNV. Non-PASS CNVscope records (FILTER MLrejected: the segments + the model rejected) are written as CNVModelApply wrote them, with SOURCE=CNV appended to INFO. +The final call set is the PASS records. CN_NEUTRAL follows CNVscope's --sex/--par or --blocks ploidy. """ import argparse +import bisect +import datetime +import io +import json +import logging +import math +import os +import re import sys from collections import defaultdict import vcflib +VERSION = '2.3.0' # combine_sv_cnv release this file belongs to + +# The rule values for paired-end CNVscope models; a --preset JSON file overrides any of them (presets/SE.json) +DEFAULTS = { + 'max_converted_loss': 1500, # converted losses longer than this are dropped (0 = keep all) + 'loss_support_size': 0, # NoSV: CNVscope losses shorter than this need a converted loss (0 = off) + 'max_converted_gain_haploid': 3000, # converted gains on haploid sequence longer than this are dropped + 'loss_min_overlap': 0.5, # SVdup: share of a CNVscope loss inside one converted loss array + 'gain_min_overlap': 0.8, # SVdup: share of a CNVscope gain inside one converted gain array, + 'gain_slack': 600, # or at most this many bp of it outside the array, + 'gain_min_shift': 0.4, # and the CN shift the converted gains must explain over it + 'nodepth_size': 10000, # NoDepth: converted gains with arrays at least this long are checked, + 'nodepth_min_shift': 0.5, # if their array-average CN shift is visible to depth, + 'nodepth_min_segmented': 0.5, # and raw segments of any state cover at least this share of the array; + 'nodepth_min_cov': 0.2, # they fail when raw gain segments cover less than this share +} +FRACTIONS = ('loss_min_overlap', 'gain_min_overlap', 'nodepth_min_segmented', 'nodepth_min_cov') # rules in [0, 1] +AR_MAGIC = b'!\n' # a Sentieon model bundle is an ar archive + +AUTOSOME_CN = 2 +SEX_CHROMS = {'chrX': 'X', 'X': 'X', 'chrY': 'Y', 'Y': 'Y'} +log = logging.getLogger(__name__) + +USAGE = ('combine_cnv.py --cnv VCF --converted VCF -o VCF [--raw VCF] [--preset FILE]\n' + ' [--sex {M,F} --par BED | --blocks BED] [-v]') +EXAMPLE = '''Production use, one sample: + sentieon driver -r REF.fa --algo CNVModelApply --model MODEL --vcf S.cnv.raw.vcf.gz --emit_conf 0 S.cnv.vcf.gz + python3 indel2cnv.py REF.fa S.sv.vcf.gz S.sv.cnv.vcf.gz -t 8 + python3 combine_cnv.py --preset presets/SE.json --cnv S.cnv.vcf.gz --converted S.sv.cnv.vcf.gz \\ + --sex M --par hs38.PAR.bed -o S.combined.cnv.vcf.gz +The final call set is the PASS records: bcftools view -f PASS S.combined.cnv.vcf.gz''' + + +def tool_version(): + """COMBINE_SV_CNV_VERSION when a driver set it, else combine_sv_cnv-v, followed by + ' (sentieon-cli-)' when this file sits in the sentieon_cli package (sentieon_cli/scripts/).""" + if os.environ.get('COMBINE_SV_CNV_VERSION'): + return os.environ['COMBINE_SV_CNV_VERSION'] + ver = 'combine_sv_cnv-v' + VERSION + if os.path.basename(os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) == 'sentieon_cli': + try: + from sentieon_cli import __version__ as scli_version + except Exception: + try: + from importlib.metadata import version + scli_version = version('sentieon_cli') + except Exception: + scli_version = 'unknown' + ver += ' (sentieon-cli-%s)' % scli_version + return ver + + +def provenance_line(tool='combine_cnv'): + """The ##CommandLine header line, in the form of Sentieon's ##SentieonCommandLine.""" + date = datetime.datetime.now(datetime.timezone.utc).strftime('%Y-%m-%dT%H:%M:%SZ') + cmd = ' '.join(sys.argv).replace('\\', '\\\\').replace('"', '\\"') + return '##CommandLine.%s=' % (tool, tool, tool_version(), date, cmd) + + +def fmt_value(v): + return '%g' % v if isinstance(v, float) else str(v) + + +def describe(rules): + return ', '.join(f'{k}={fmt_value(v)}' for k, v in rules.items()) + -def eprint(*args, **kwargs): - print(*args, file=sys.stderr, **kwargs) +def ar_member(archive, name): + """The bytes of the member called name in an ar archive (a model bundle), reading only the member headers and + that member. Member names as sentieon_cli.archive.ar_load reads them: GNU/SysV "name/" and "/" into the + "//" table, BSD "#1/" and plain names; the symbol tables are skipped. Raises ValueError naming the + archive.""" + try: + f = open(archive, 'rb') + except OSError as e: + raise ValueError(f'{archive}: cannot read ({e.strerror})') from None + with f: + if f.read(len(AR_MAGIC)) != AR_MAGIC: + raise ValueError(f'{archive}: not an ar archive') + file_size, want = os.fstat(f.fileno()).st_size, os.fsencode(name) + names, table = [], None + while True: + hdr = f.read(60) + if not hdr: + break + if len(hdr) < 60: + raise ValueError(f'{archive}: truncated ar archive') + field, size = hdr[:16].rstrip(b' '), hdr[48:58].strip() + if hdr[58:] != b'`\n' or not size.isdigit(): + raise ValueError(f'{archive}: corrupted ar archive') + start, size = f.tell(), int(size) + if start + size > file_size: + raise ValueError(f'{archive}: truncated ar archive') + member = None + if field in (b'/', b'/SYM64/'): # symbol tables + pass + elif field == b'//': # GNU long-name table + table = f.read(size) + elif field[:1] == b'/': # GNU long name: / into the table + off = field[1:] + stop = table.find(b'/\n', int(off)) if table is not None and off.isdigit() else -1 + if stop < 0: + raise ValueError(f'{archive}: corrupted ar archive') + member = table[int(off):stop] + elif field[-1:] == b'/': # GNU/SysV short name + member = field[:-1] + elif field[:3] == b'#1/': # BSD long name: the first bytes of the data + n = field[3:] + if not n.isdigit() or int(n) > size: + raise ValueError(f'{archive}: corrupted ar archive') + member = f.read(int(n)).rstrip(b'\0') + else: # BSD short name + member = field + if member == want: + return f.read(start + size - f.tell()) + if member is not None: + names.append(member.decode('utf-8', 'replace')) + f.seek(start + size + size % 2) + raise ValueError(f'{archive}: no member "{name}"; members: {", ".join(names) or "none"}') -def load_converted_intervals(converted_vcf): - """Load REFSTART/REFSTOP intervals from converted calls, keyed by (chrom, direction).""" +def read_preset(path): + """The bytes of a preset: the file path, or member M of the ar archive (model bundle) A when path is A/M.""" + archive, member = os.path.split(path) + if os.path.isfile(archive): + return ar_member(archive, member) + try: + with open(path, 'rb') as f: + head = f.read(len(AR_MAGIC)) + if head == AR_MAGIC: + raise ValueError(f'{path}: an ar archive (model bundle); name a member, e.g. {path}/combine_cnv.json') + return head + f.read() + except (FileNotFoundError, NotADirectoryError): # NotADirectoryError: the path runs through a file + raise ValueError(f'{path}: no such file') from None + except OSError as e: + raise ValueError(f'{path}: cannot read ({e.strerror})') from None + + +def load_preset(path): + """DEFAULTS updated by a preset: one JSON object of DEFAULTS keys and an optional "description" string, in a file + or a member of a model bundle (MODEL.bundle/combine_cnv.json). Raises ValueError naming the path and the key.""" + def unique(pairs): + keys = [k for k, _ in pairs] + for k in keys: + if keys.count(k) > 1: + raise ValueError(f'{path}: duplicate key "{k}"') + return dict(pairs) + + text = io.TextIOWrapper(io.BytesIO(read_preset(path)), encoding='utf-8') # decoded as open(path) would + try: + data = json.load(text, object_pairs_hook=unique) + except UnicodeDecodeError: + raise ValueError(f'{path}: invalid JSON (not UTF-8 text)') from None + except json.JSONDecodeError as e: + raise ValueError(f'{path}: invalid JSON ({e})') from None + if not isinstance(data, dict): + raise ValueError(f'{path}: the top-level value is not a JSON object') + rules = dict(DEFAULTS) + for key, value in data.items(): + shown = json.dumps(value) + if key == 'description': + if not isinstance(value, str): + raise ValueError(f'{path}: "description" is {shown}, not a string') + continue + if key not in DEFAULTS: + raise ValueError(f'{path}: unknown key "{key}"; valid keys: description, {", ".join(DEFAULTS)}') + is_int = isinstance(DEFAULTS[key], int) + if isinstance(value, bool): + raise ValueError(f'{path}: "{key}" is {shown}, a boolean, not a number') + if is_int and not isinstance(value, int): + raise ValueError(f'{path}: "{key}" is {shown}, not an integer') + if not isinstance(value, (int, float)) or isinstance(value, float) and not math.isfinite(value): + raise ValueError(f'{path}: "{key}" is {shown}, not a number') + if key in FRACTIONS and not 0 <= value <= 1: + raise ValueError(f'{path}: "{key}" is {shown}, outside [0, 1]') + if value < 0: + raise ValueError(f'{path}: "{key}" is {shown}, negative') + rules[key] = value if is_int else float(value) + return rules + + +def preset_rules(parser, path): + """The rule values in effect and where they come from: DEFAULTS, updated by the --preset file when given.""" + if path is None: + return dict(DEFAULTS), "the code's defaults (the PE values)" + try: + return load_preset(path), f'--preset {path}' + except ValueError as e: + parser.error(f'--preset {e}') + + +def parse_par_bed(path): + """chrom -> [(start, end)] 0-based half-open, same as cnv_eval.py.""" + intervals = {} + with open(path) as f: + for line in f: + flds = line.rstrip().split('\t') + if not line.strip() or line.startswith('#') or len(flds) < 3: + continue + intervals.setdefault(flds[0], []).append((int(flds[1]), int(flds[2]))) + return intervals + + +def parse_blocks_bed(path): + """CNVscope --blocks BED: chrom -> [(start, end, cn_neutral)] plus the sex= header hint.""" + blocks, sex = {}, None + with open(path) as f: + for line in f: + if line.startswith('#'): + m = re.search(r'\bsex=([MF])\b', line) + if m: + sex = m.group(1) + continue + flds = line.rstrip().split('\t') + if not line.strip() or len(flds) < 4: + continue + blocks.setdefault(flds[0], []).append((int(flds[1]), int(flds[2]), int(flds[3]))) + return blocks, sex + + +class Ploidy: + """Neutral CN lookup mirroring CNVscope: the --blocks block covering pos, else the --sex/--par preset.""" + def __init__(self, sex, par_intervals, blocks=None): + self.sex = sex + self.par = par_intervals or {} + self.blocks = {} + for chrom, ivs in (blocks or {}).items(): + ivs = sorted(ivs) + self.blocks[chrom] = ([s for s, _, _ in ivs], ivs) + + def cn_neutral(self, chrom, pos): + """-1 = uncovered (female chrY), else the local neutral CN.""" + if chrom in self.blocks: + starts, ivs = self.blocks[chrom] + i = bisect.bisect_right(starts, pos) - 1 + if i >= 0 and pos < ivs[i][1]: + return ivs[i][2] + sc = SEX_CHROMS.get(chrom) + if not self.sex or sc is None: + return AUTOSOME_CN + if self.sex == 'F': + return -1 if sc == 'Y' else AUTOSOME_CN + if sc == 'X' and any(s <= pos < e for s, e in self.par.get(chrom, ())): + return AUTOSOME_CN + return 1 + + +def header_sex(vcf): + for h in vcf.headers: + if h.startswith('##SampleSex='): + return h.split('=', 1)[1].strip() + return None + + +def gt_alleles(v): + """The sample's GT alleles as strings ('1|1' -> ['1', '1'], '1' -> ['1']).""" + gt = v.samples[0].get('GT', '.') if v.samples else '.' + return re.split(r'[/|]', str(gt)) + + +def is_hom_alt(v): + """Every allele non-reference and called: 1/1, 1|1, 1 (haploid), 1|2.""" + return all(a not in ('0', '.', '') for a in gt_alleles(v)) + + +def over_loss_cap(v, max_loss): + return v.info.get('CNDIFF', 0) < 0 and bool(max_loss) and v.end - v.pos > max_loss + + +def haploid_drop(v, cn0, max_gain_haploid): + """Why a converted call on haploid sequence (CN_NEUTRAL 1) is dropped, or None.""" + if cn0 != 1: + return None + cndiff = v.info.get('CNDIFF', 0) + # the SV caller genotypes male X/Y as diploid: a het deletion on single-copy sequence cannot be real + if cndiff < 0 and not is_hom_alt(v): + return 'haploid_het_loss' + if cndiff > 0 and max_gain_haploid and v.end - v.pos > max_gain_haploid: + return 'haploid_gain_oversize' + return None + + +def converted_cn(v, cn0): + """CN_NEUTRAL + CNDIFF floored at 0, and whether it was floored. On haploid sequence a hom call's CNDIFF counts + the one copy twice (the converter sums both diploid haplotypes), so it is halved there.""" + cndiff = v.info.get('CNDIFF', 0) + if cn0 == 1 and is_hom_alt(v) and len(gt_alleles(v)) > 1: + half = (abs(cndiff) + 1) // 2 # an odd count (merged units) rounds away from zero + cndiff = half if cndiff > 0 else -half + return max(0, cn0 + cndiff), cn0 + cndiff < 0 + + +def load_converted_intervals(converted_vcf, max_loss=0, ploidy=None, max_gain_haploid=0): + """(chrom, direction) -> sorted (refstart, refstop, cndiff, svlen) of the converted calls; with ploidy, only the + calls the size cap and the haploid rules keep.""" intervals = defaultdict(list) - n = 0 for chrom in converted_vcf.contigs: for v in converted_vcf.range(chrom): cndiff = v.info.get('CNDIFF', 0) refstart = v.info.get('REFSTART', 0) refstop = v.info.get('REFSTOP', 0) - if cndiff == 0 or refstart >= refstop: + if cndiff == 0 or refstart >= refstop or over_loss_cap(v, max_loss): continue - d = -1 if cndiff < 0 else 1 - intervals[(chrom, d)].append((refstart, refstop)) - n += 1 + if ploidy is not None and haploid_drop(v, ploidy.cn_neutral(chrom, v.pos), max_gain_haploid): + continue + intervals[(chrom, -1 if cndiff < 0 else 1)].append((refstart, refstop, cndiff, v.info.get('SVLEN', 0))) for key in intervals: intervals[key].sort() - return intervals, n + return intervals -def overlaps_any(intervals, pos, end, min_frac): - """Check if [pos, end) overlaps any interval by >= min_frac of the call.""" +def overlaps_any(intervals, pos, end, min_frac, slack=0): + """The first interval holding >= min_frac of [pos, end), or leaving at most `slack` bp of it outside.""" length = end - pos if length <= 0: return None - for rs, re_ in intervals: + for rs, re_, _, _ in intervals: if rs >= end: break if re_ <= pos: continue ovl = min(end, re_) - max(pos, rs) - if ovl / length >= min_frac: + if ovl / length >= min_frac or length - ovl <= slack: return (rs, re_) return None -def main(): +def explained_shift(intervals, pos, end): + """Expected CN shift over [pos, end) from converted calls (CNDIFF*SVLEN spread over each array), and the + largest contribution.""" + length = end - pos + shift = 0.0 + best = None + for rs, re_, cndiff, svlen in intervals: + if rs >= end: + break + if re_ <= pos or re_ <= rs: + continue + part = cndiff * svlen * (min(end, re_) - max(pos, rs)) / ((re_ - rs) * length) + shift += part + if best is None or abs(part) > abs(best[2]): + best = (rs, re_, part) + return shift, best + + +def loss_duplicate(intervals, pos, end, rules): + """The converted loss array holding loss_min_overlap of the CNVscope loss, or None.""" + return overlaps_any(intervals, pos, end, rules['loss_min_overlap']) + + +def gain_duplicate(intervals, pos, end, rules): + """The converted gain array holding the CNVscope gain (gain_min_overlap, or gain_slack bp outside at most), + when the converted gains also explain gain_min_shift of CN over it; else None.""" + hit = overlaps_any(intervals, pos, end, rules['gain_min_overlap'], rules['gain_slack']) + if hit is None: + return None + shift, best = explained_shift(intervals, pos, end) + return hit if best is not None and abs(shift) >= rules['gain_min_shift'] else None + + +def loss_supported(intervals, pos, end): + """Does any converted loss array overlap [pos, end)?""" + if end <= pos: + return False + for rs, re_, _, _ in intervals: + if rs >= end: + break + if re_ > pos: + return True + return False + + +def load_raw_segments(records, chrom, ploidy): + """CNVscope's raw segments among the records on chrom: -1 losses, 1 gains, 0 every segment (where CNVscope had + depth to judge).""" + segs = {-1: [], 0: [], 1: []} + for v in records: + cn0 = v.info.get('CN_NEUTRAL', ploidy.cn_neutral(chrom, v.pos)) + cn = v.info.get('HMMCN', cn0) + if cn0 < 0: + continue + segs[0].append((v.pos, v.end, cn)) + if cn != cn0: + segs[-1 if cn < cn0 else 1].append((v.pos, v.end, cn)) + for d in segs: + segs[d].sort() + return segs + + +def load_pass_calls(records, chrom, ploidy): + """CNVscope's PASS calls among the records on chrom: -1 losses, 1 gains (the gain evidence for raw segments without + HMM states).""" + calls = {-1: [], 1: []} + for v in records: + if v.filter and set(v.filter).difference(('PASS',)): + continue + cn0 = v.info.get('CN_NEUTRAL', ploidy.cn_neutral(chrom, v.pos)) + cn = v.info.get('CN', cn0) + if cn0 >= 0 and cn != cn0: + calls[-1 if cn < cn0 else 1].append((v.pos, v.end, cn)) + for d in calls: + calls[d].sort() + return calls + + +def with_source(line): + """A record line as it was, with SOURCE=CNV appended to its INFO.""" + f = line.split('\t') + f[7] = 'SOURCE=CNV' if f[7] == '.' else f[7] + ';SOURCE=CNV' + return '\t'.join(f) + + +def has_line(vcf, prefix): + return any(h.startswith(prefix) for h in vcf.headers) + + +def emit_conf(vcf): + """The --emit_conf value on the last CNVModelApply command line in the header, or None.""" + lines = [h for h in vcf.headers if h.startswith('##SentieonCommandLine.CNVModelApply=')] + values = re.findall(r'--emit_conf[=\s]+([^\s">]+)', lines[-1]) if lines else [] + try: + return float(values[-1]) if values else None + except ValueError: + return None + + +def has_hmm_states(raw_vcf): + """Whether the raw segments carry their HMM state (HMMCN).""" + return has_line(raw_vcf, '##INFO== end: + break + if ce <= start: + continue + lo, hi = max(cs, reach), min(ce, end) + if hi > lo: + covered += hi - lo + reach = hi + return covered / length + + +def no_depth(cndiff, svlen, rs, re_, raw, rules): + """A long converted gain over which CNVscope had depth but saw (almost) no gain.""" + if raw is None or cndiff <= 0 or re_ - rs < rules['nodepth_size']: + return False + if abs(cndiff * svlen) / (re_ - rs) < rules['nodepth_min_shift']: + return False + return (covered_fraction(raw[0], rs, re_) >= rules['nodepth_min_segmented'] + and covered_fraction(raw[1], rs, re_) < rules['nodepth_min_cov']) + + +def build_parser(): parser = argparse.ArgumentParser( - description='Combine CNVscope and converted SV calls (both CNV space)') - parser.add_argument('--cnv', required=True, help='CNVscope call VCF') - parser.add_argument('--converted', required=True, - help='indel2cnv.py output VCF (converted SV calls)') - parser.add_argument('-o', '--output', required=True, help='Output VCF (.vcf.gz)') - parser.add_argument('--min_overlap', type=float, default=0.5, - help='Min overlap fraction to mark SVdup (default: 0.5)') - parser.add_argument('--neutral_cn', type=int, default=2, - help='Neutral copy number (default: 2)') - args = parser.parse_args() + prog='combine_cnv.py', usage=USAGE, epilog=EXAMPLE, formatter_class=argparse.RawDescriptionHelpFormatter, + description='Combine CNVscope calls with SV calls converted to CNV (indel2cnv.py) into one CNV VCF.') + req = parser.add_argument_group('required') + req.add_argument('--cnv', required=True, metavar='VCF', help='CNVscope calls (CNVModelApply output)') + req.add_argument('--converted', required=True, metavar='VCF', help='SV calls converted by indel2cnv.py') + req.add_argument('-o', '--output', required=True, metavar='VCF', help='combined CNV VCF (.vcf.gz)') + parser.add_argument('--raw', metavar='VCF', + help="CNVscope's own output (the input of CNVModelApply); default: the segments in --cnv, " + 'which CNVModelApply must have written with --emit_conf 0') + parser.add_argument('--preset', metavar='FILE', + help='JSON file of rule values, or a member of a model bundle, e.g. ' + 'MODEL.bundle/combine_cnv.json; without it the PE values apply (paired-end models). ' + 'presets/SE.json is for single-end models (cnv.se.model: Roche SBX, Ultima)') + parser.add_argument('--sex', choices=['M', 'F'], help='sample sex (default: ##SampleSex in --cnv)') + parser.add_argument('--par', metavar='BED', help='PAR BED, required for male samples, as for CNVscope') + parser.add_argument('--blocks', metavar='BED', help='instead of --sex/--par: the CNVscope --blocks file') + parser.add_argument('-v', '--verbose', action='store_true', help='log progress and counts') + return parser - conv_vcf = vcflib.VCF(args.converted, 'r') - cnv_vcf = vcflib.VCF(args.cnv, 'r') - # Build dedup intervals from converted calls - eprint('Loading converted intervals...') - conv_intervals, n_conv = load_converted_intervals(conv_vcf) - eprint(f' {n_conv} intervals') +def setup(parser, args): + """Check the inputs and resolve the ploidy as CNVscope does; returns (ploidy, sex).""" + cnv_vcf = vcflib.VCF(args.cnv, 'r') + raw_vcf = vcflib.VCF(args.raw, 'r') if args.raw is not None else None + try: + if not has_line(cnv_vcf, '##SentieonCommandLine.CNVModelApply='): + parser.error(f'--cnv {args.cnv} is not CNVModelApply output') + if raw_vcf is None: + conf = emit_conf(cnv_vcf) + if conf is None or conf > 0: + parser.error(f'--cnv {args.cnv} was not written with CNVModelApply --emit_conf 0 and so lacks most ' + "rejected segments: pass --raw (CNVscope's output) or rerun CNVModelApply with " + '--emit_conf 0') + elif (not has_line(raw_vcf, '##SentieonCommandLine.CNVscope=') + or has_line(raw_vcf, '##SentieonCommandLine.CNVModelApply=')): + parser.error(f'--raw {args.raw} is not CNVscope\'s own output (the input of CNVModelApply)') + if args.blocks and (args.sex or args.par): + parser.error('--sex/--par and --blocks are mutually exclusive (as for CNVscope)') + blocks, blocks_sex = parse_blocks_bed(args.blocks) if args.blocks else ({}, None) + hsex = header_sex(cnv_vcf) + if args.blocks and hsex: + log.warning('--blocks given; ignoring ##SampleSex=%s in %s', hsex, args.cnv) + sex = blocks_sex if args.blocks else (args.sex or hsex) + if args.sex and hsex and args.sex != hsex: + log.warning('--sex %s differs from ##SampleSex=%s in %s', args.sex, hsex, args.cnv) + if not sex and not args.blocks: + # read as diploid, a male's chrX/chrY would turn CNVscope's single-copy calls into losses + parser.error(f'--cnv {args.cnv} names no sex (no ##SampleSex): give --sex (with --par for a male sample) ' + 'or --blocks') + if sex == 'M' and not args.par and not args.blocks: + parser.error('--sex M requires --par ') + par_intervals = parse_par_bed(args.par) if args.par else {} + for chrom in par_intervals: + if chrom not in cnv_vcf.contigs: + parser.error(f'PAR BED contig {chrom!r} not in {args.cnv}') + for chrom, ivs in blocks.items(): + if chrom not in cnv_vcf.contigs: + parser.error(f'blocks BED contig {chrom!r} not in {args.cnv}') + if any(cn < 1 or cn > 4 for _, _, cn in ivs): + parser.error(f'--blocks cn_neutral must be in [1,4] ({chrom})') + finally: + cnv_vcf.close() + if raw_vcf is not None: + raw_vcf.close() + if args.blocks: + log.info('Ploidy: blocks=%s (%d blocks; uncovered sex chroms use preset %s)', args.blocks, + sum(len(b) for b in blocks.values()), sex or 'diploid') + else: + log.info('Ploidy: sex=%s par=%s', sex or 'unset (diploid)', args.par or 'none') + return Ploidy(sex, par_intervals, blocks), sex - # Reopen for full iteration - conv_vcf.close() - conv_vcf = vcflib.VCF(args.converted, 'r') - # Merged header: CNV as base, add converted-specific headers - conv_extra = [h for h in conv_vcf.headers if h.startswith('##')] - new_hdrs = ( +def output_headers(sex, tool='combine_cnv'): + hdrs = [ + '##INFO=', '##INFO=', '##INFO=', - '##FILTER=', - ) - out_vcf = vcflib.VCF(args.output, 'wb') - out_vcf.copy_header(cnv_vcf, update=tuple(conv_extra) + new_hdrs) - out_vcf.emit_header() + '##FILTER=', + '##FILTER=', + '##FILTER=', + provenance_line(tool), + ] + if sex: + hdrs.append(f'##SampleSex={sex}') + return hdrs - all_contigs = list(cnv_vcf.contigs.keys()) - for c in conv_vcf.contigs: - if c not in cnv_vcf.contigs: - all_contigs.append(c) - conv_stats = {'total': 0} - cnv_stats = {'kept': 0, 'marked': 0, 'neutral': 0, 'filtered': 0} +def combine(cnv_path, conv_path, raw_path, out_path, rules, ploidy, headers, on_converted=None): + """Write the combined VCF and return the counts. Without raw_path, cnv_path (an --emit_conf 0 file) gives the + segments. on_converted(chrom, v, rs, re_, d, raw) is called for every kept converted call before its filter is set + (the dev tool annotates through it).""" + conv_vcf = vcflib.VCF(conv_path, 'r') + conv_intervals = load_converted_intervals(conv_vcf, rules['max_converted_loss'], ploidy, + rules['max_converted_gain_haploid']) + support = None + if rules['loss_support_size']: + # every converted loss counts, before the size cap and the haploid rules: did the SV caller see a loss here? + support = load_converted_intervals(conv_vcf) + conv_vcf.close() + conv_vcf = vcflib.VCF(conv_path, 'r') + cnv_vcf = vcflib.VCF(cnv_path, 'r') + raw_vcf = vcflib.VCF(raw_path, 'r') if raw_path is not None else None + seg_vcf = cnv_vcf if raw_vcf is None else raw_vcf # without --raw the records of --cnv are the segments + raw_states = has_hmm_states(seg_vcf) + + # header: CNVscope's, plus the converter's own lines (not its contigs) + conv_extra = [h for h in conv_vcf.headers if h.startswith('##') and not h.startswith(('##contig=', '##fileformat='))] + out_vcf = vcflib.VCF(out_path, 'wb') + out_vcf.copy_header(cnv_vcf, update=tuple(conv_extra) + tuple(headers)) + out_vcf.emit_header() + all_contigs = list(cnv_vcf.contigs.keys()) + [c for c in conv_vcf.contigs if c not in cnv_vcf.contigs] + conv_n = dict.fromkeys(('kept', 'floored', 'uncovered', 'over_cap', 'nodepth', 'haploid_het_loss', + 'haploid_gain_oversize'), 0) + cnv_n = dict.fromkeys(('kept', 'svdup', 'svdup_gain', 'nosv', 'rejected', 'neutral', 'uncovered'), 0) seq = 0 for chrom in all_contigs: records = [] + cnv_records = list(cnv_vcf.range(chrom)) # read once: segments, PASS-call evidence and output + raw = None + if chrom in seg_vcf.contigs: + raw = load_raw_segments(cnv_records if raw_vcf is None else raw_vcf.range(chrom), chrom, ploidy) + if not raw_states: + raw.update(load_pass_calls(cnv_records, chrom, ploidy)) - # Converted SV calls: emit all, add CN from CNDIFF for v in conv_vcf.range(chrom): + if over_loss_cap(v, rules['max_converted_loss']): + conv_n['over_cap'] += 1 + continue + cn0 = ploidy.cn_neutral(chrom, v.pos) + if cn0 < 0: + conv_n['uncovered'] += 1 + continue + why = haploid_drop(v, cn0, rules['max_converted_gain_haploid']) + if why: + conv_n[why] += 1 + continue cndiff = v.info.get('CNDIFF', 0) - v.info['CN'] = args.neutral_cn + cndiff + cn, floored = converted_cn(v, cn0) + if floored: + conv_n['floored'] += 1 + log.info('%s:%d CNDIFF=%d with CN_NEUTRAL=%d: CN floored to 0', chrom, v.pos + 1, cndiff, cn0) + v.info['CN'] = cn + v.info['CN_NEUTRAL'] = cn0 v.info['SOURCE'] = 'SV_CNV' + d = -1 if cndiff < 0 else 1 + rs, re_ = v.info.get('REFSTART', v.pos), v.info.get('REFSTOP', v.end) + if on_converted: + on_converted(chrom, v, rs, re_, d, raw) + if no_depth(cndiff, v.info.get('SVLEN', 0), rs, re_, raw, rules): + v.filter = ['NoDepth'] + conv_n['nodepth'] += 1 + # an unsupported gain must not suppress the CNVscope gain + try: + conv_intervals[(chrom, d)].remove((rs, re_, cndiff, v.info.get('SVLEN', 0))) + except (KeyError, ValueError): + pass + elif not v.filter: + v.filter = ['PASS'] v.line = None - conv_stats['total'] += 1 + conv_n['kept'] += 1 seq += 1 records.append((v.pos, 0, seq, v)) - # CNVscope calls: mark losses overlapping converted REFSTART/REFSTOP - for v in cnv_vcf.range(chrom): - if v.filter and set(v.filter).difference(('PASS',)): - cnv_stats['filtered'] += 1 + for v in cnv_records: + cn0 = v.info.get('CN_NEUTRAL', ploidy.cn_neutral(chrom, v.pos)) + if cn0 < 0: + cnv_n['uncovered'] += 1 continue - cn = v.info.get('CN', args.neutral_cn) - if cn == args.neutral_cn: - cnv_stats['neutral'] += 1 + cn = v.info.get('CN', cn0) + if v.filter and set(v.filter).difference(('PASS',)): + cnv_n['rejected'] += 1 + v.line = with_source(v.line) # a segment the model rejected: CNVModelApply's own line + elif cn == cn0: + cnv_n['neutral'] += 1 continue - v.info['SOURCE'] = 'CNV' - if cn < args.neutral_cn: - key = (chrom, -1) - interval = overlaps_any(conv_intervals.get(key, []), - v.pos, v.end, args.min_overlap) - if interval: + else: + d = -1 if cn < cn0 else 1 + ivs = conv_intervals.get((chrom, d), []) + hit = (loss_duplicate if d < 0 else gain_duplicate)(ivs, v.pos, v.end, rules) + if hit: v.filter = ['SVdup'] - v.info['DEDUP'] = f'{chrom}:{interval[0]}-{interval[1]}' - cnv_stats['marked'] += 1 - v.line = None - seq += 1 - records.append((v.pos, 1, seq, v)) - continue - cnv_stats['kept'] += 1 - v.line = None + v.info['DEDUP'] = f'{chrom}:{hit[0]}-{hit[1]}' + cnv_n['svdup'] += 1 + cnv_n['svdup_gain'] += d > 0 + elif (d < 0 and support is not None and v.end - v.pos < rules['loss_support_size'] + and not loss_supported(support.get((chrom, -1), []), v.pos, v.end)): + v.filter = ['NoSV'] + cnv_n['nosv'] += 1 + else: + cnv_n['kept'] += 1 + v.info['SOURCE'] = 'CNV' + v.line = None seq += 1 records.append((v.pos, 1, seq, v)) records.sort() for _, _, _, v in records: out_vcf.emit(v) + del records, cnv_records # free this chromosome's records before the next is read out_vcf.close() conv_vcf.close() cnv_vcf.close() + if raw_vcf is not None: + raw_vcf.close() + log.info('Converted: %d written (%d NoDepth gains, %d floored to CN 0); dropped %d over the loss cap, %d haploid ' + 'het losses, %d haploid gains over the cap, %d on uncovered sequence', conv_n['kept'], conv_n['nodepth'], + conv_n['floored'], conv_n['over_cap'], conv_n['haploid_het_loss'], conv_n['haploid_gain_oversize'], + conv_n['uncovered']) + log.info('CNVscope: %d PASS calls kept, %d SVdup (%d gains), %d NoSV, %d rejected segments written; left out %d ' + 'neutral PASS records, %d on uncovered sequence', cnv_n['kept'], cnv_n['svdup'], cnv_n['svdup_gain'], + cnv_n['nosv'], cnv_n['rejected'], cnv_n['neutral'], cnv_n['uncovered']) + return conv_n, cnv_n - eprint(f'Converted: {conv_stats["total"]} records') - eprint(f'CNV: {cnv_stats["kept"]} kept, {cnv_stats["marked"]} marked SVdup, ' - f'{cnv_stats["neutral"]} neutral, {cnv_stats["filtered"]} non-PASS') - eprint(f'Written to {args.output}') + +def main(): + parser = build_parser() + args = parser.parse_args() + logging.basicConfig(level=logging.INFO if args.verbose else logging.WARNING, + format='%(asctime)s %(levelname)s %(message)s') + rules, source = preset_rules(parser, args.preset) + log.info('%s, rules from %s: %s', tool_version(), source, describe(rules)) + ploidy, sex = setup(parser, args) + combine(args.cnv, args.converted, args.raw, args.output, rules, ploidy, output_headers(sex)) + log.info('Written to %s', args.output) if __name__ == '__main__': diff --git a/sentieon_cli/scripts/gvcf_clip.py b/sentieon_cli/scripts/gvcf_clip.py new file mode 100644 index 0000000..476798a --- /dev/null +++ b/sentieon_cli/scripts/gvcf_clip.py @@ -0,0 +1,371 @@ +#!/usr/bin/env python +"""Restrict a gVCF to the inside, or the outside, of a set of BED regions + +The DNAscope hybrid pipeline splices its two calling passes together: the +first pass covers the positions outside the realigned regions and the +second pass the positions inside them. Variant records are kept or dropped +by their start position, as with `bcftools view -T`. Reference blocks are +clipped at the region boundaries. A dropped variant record can span +positions on the kept side, such as a deletion that starts inside a region +and ends past it; a reference block with no genotype confidence covers +whichever of those positions no kept record covers. So the spliced gVCF +has neither overlapping blocks nor gaps. + +An output of `-` writes uncompressed VCF to stdout. +""" + +import argparse +import bisect +import collections +import functools +import io +import multiprocessing as mp +import re +import sys + +import vcflib + +REF_BLOCK_ALTS = ([""], ["<*>"]) +# The FORMAT fields of a fill block, in this order, when the header has them +FILL_KEYS = ("GT", "DP", "GQ", "MIN_DP", "PL") +GT_SEP = re.compile(r"[/|]") + +Contig = collections.namedtuple("Contig", "length offset width skip") + + +class Reference: + """Random access to an uncompressed, faidx-indexed FASTA file""" + + def __init__(self, path): + self.path = path + self.index = collections.OrderedDict() + with io.open(self.path + ".fai", "rb") as fp: + for line in fp: + flds = line.rstrip().decode().split("\t") + if len(flds) != 5: + raise RuntimeError("Corrupted index") + self.index[flds[0]] = Contig._make(map(int, flds[1:])) + self.fp = io.open(path, "rb") + + def __getstate__(self): + odict = self.__dict__.copy() + odict["fp"] = self.fp.tell() + return odict + + def __setstate__(self, ndict): + fp = io.open(ndict["path"], "rb") + fp.seek(ndict.pop("fp")) + ndict["fp"] = fp + self.__dict__.update(ndict) + + def get(self, c, s, e): + """The 0-based, half-open reference sequence [s, e) of contig c""" + ci = self.index[c] + e = min(e, ci.length) + if s > e: + raise ValueError + seq = b"" + while s < e: + o = ci.offset + s // ci.width * ci.skip + s % ci.width + n = min(ci.width - s % ci.width, e - s) + self.fp.seek(o) + seq += self.fp.read(n) + s += n + return seq.decode() + + +def load_regions(path): + """Read a BED file into sorted, merged regions per contig + + Each contig maps to a pair of lists, the 0-based region starts and the + exclusive region ends. + """ + raw = collections.defaultdict(list) + with open(path) as fh: + for line in fh: + cols = line.split() + if not cols or cols[0] in ("track", "browser"): + continue + if cols[0].startswith("#"): + continue + raw[cols[0]].append((int(cols[1]), int(cols[2]))) + + regions = {} + for chrom, intervals in raw.items(): + starts, ends = [], [] + for s, e in sorted(intervals): + if starts and s <= ends[-1]: + ends[-1] = max(ends[-1], e) + else: + starts.append(s) + ends.append(e) + regions[chrom] = (starts, ends) + return regions + + +def in_regions(regions, chrom, pos): + """Whether the 0-based position falls inside a region""" + starts, ends = regions.get(chrom, ([], [])) + i = bisect.bisect_right(starts, pos) - 1 + return i >= 0 and pos < ends[i] + + +def pieces(regions, chrom, start, end, include): + """The parts of [start, end) inside, or outside, the regions""" + starts, ends = regions.get(chrom, ([], [])) + out = [] + cur = start + i = bisect.bisect_right(ends, start) # first region ending after start + while i < len(starts) and starts[i] < end: + s, e = max(starts[i], start), min(ends[i], end) + if include: + out.append((s, e)) + elif cur < s: + out.append((cur, s)) + cur = e + i += 1 + if not include and cur < end: + out.append((cur, end)) + return out + + +def set_end(info, end): + """Set the INFO/END value of an INFO column""" + flds = [] if info == "." else info.split(";") + for i, fld in enumerate(flds): + if fld.startswith("END="): + flds[i] = "END=%d" % end + break + else: + flds.append("END=%d" % end) + return ";".join(flds) + + +def clip_block(v, start, end, reference): + """A copy of the reference block `v` covering [start, end)""" + if start == v.pos: + ref = v.ref[0] + else: + ref = reference.get(v.chrom, start, start + 1).upper() + cols = v.line.split("\t") + cols[1] = str(start + 1) + cols[3] = ref + cols[7] = set_end(cols[7], end) + info = dict(v.info, END=end) + return vcflib.Variant( + v.chrom, + start, + v.id, + ref, + v.alt, + v.qual, + v.filter, + info, + v.samples, + end, + "\t".join(cols), + ) + + +def fill_block(v, start, end, reference, alt, keys): + """A reference block over [start, end) for the dropped record `v` + + The block claims no genotype confidence: GQ is 0 and the likelihoods + are flat. DP and MIN_DP are the depth of `v`. + """ + ref = reference.get(v.chrom, start, start + 1).upper() + cols = [v.chrom, str(start + 1), ".", ref, alt, ".", ".", "END=%d" % end] + if v.samples: + cols.append(":".join(keys)) + for sample in v.samples: + ploidy = len(GT_SEP.split(sample.get("GT") or "0/0")) + dp = sample.get("DP") + values = { + "GT": "/".join(["0"] * ploidy), + "DP": "." if dp is None else str(dp), + "GQ": "0", + "MIN_DP": "." if dp is None else str(dp), + "PL": ",".join(["0"] * (ploidy + 1)), + } + cols.append(":".join(values[k] for k in keys)) + return vcflib.Variant( + v.chrom, + start, + ".", + ref, + [alt], + None, + [], + {"END": end}, + [{} for _ in v.samples], + end, + "\t".join(cols), + ) + + +class Writer: + """Writes the records of one shard, or of the whole file, in order + + Kept records and block pieces are written by the shard they start in. + A fill block waits for the records after it, which may cover part of + it, and is cut at the shard boundaries, since the records that could + cover it past a boundary belong to the next shard. + """ + + def __init__(self, outvcf, start, end): + self.outvcf = outvcf + self.start = start + self.end = end + self.chrom = None + # The end of everything written or scheduled on the contig + self.covered = 0 + # Scheduled fill blocks: (start, end, block builder) + self.fills = [] + + def record(self, v): + """Write a kept record or a block piece""" + self.contig(v.chrom) + pending, self.fills = self.fills, [] + for start, end, make in pending: + # The part before `v` is final; `v` covers the part it spans + self.write_fill(start, min(end, v.pos), make) + if max(start, v.end) < end: + self.fills.append((max(start, v.end), end, make)) + self.covered = max(self.covered, v.end) + if v.pos >= self.start: + self.outvcf.emit(v) + + def fill(self, chrom, start, end, make): + """Schedule a fill block over the uncovered part of [start, end)""" + self.contig(chrom) + start = max(start, self.covered) + if start < end: + self.fills.append((start, end, make)) + self.covered = end + + def contig(self, chrom): + if chrom != self.chrom: + self.finish() + self.chrom, self.covered = chrom, 0 + + def finish(self): + """Write the scheduled fill blocks""" + for start, end, make in self.fills: + self.write_fill(start, end, make) + self.fills = [] + + def write_fill(self, start, end, make): + start = max(start, self.start) + if self.end is not None: + end = min(end, self.end) + if start < end: + self.outvcf.emit(make(start, end)) + + +def clip(invcf, outvcf, regions, reference, include, alt, keys): + """Write the records, or parts of records, inside or outside regions + + When sharded, `invcf` yields every record overlapping the shard, so the + records starting before the shard still count toward what is covered. + """ + writer = Writer( + outvcf, getattr(invcf, "start", -1), getattr(invcf, "end", None) + ) + for v in invcf: + if v.alt in REF_BLOCK_ALTS: + for start, end in pieces(regions, v.chrom, v.pos, v.end, include): + if start == v.pos and end == v.end: + writer.record(v) + else: + writer.record(clip_block(v, start, end, reference)) + elif in_regions(regions, v.chrom, v.pos) == include: + writer.record(v) + else: + make = functools.partial( + fill_block, v, reference=reference, alt=alt, keys=keys + ) + for start, end in pieces(regions, v.chrom, v.pos, v.end, include): + writer.fill(v.chrom, start, end, make) + writer.finish() + + +def main(args): + regions = load_regions(args.bed) + reference = Reference(args.reference) + include = not args.exclude + + invcf = vcflib.VCF(args.vcf, "r") + # The symbolic ALT of a fill block: the one the header declares + alt = "" + if "*" in invcf.alts and "NON_REF" not in invcf.alts: + alt = "<*>" + keys = [k for k in FILL_KEYS if k in invcf.formats] + outvcf = vcflib.VCF(args.output, "w") + outvcf.copy_header(invcf) + outvcf.emit_header() + + if args.threads < 2: + clip(invcf, outvcf, regions, reference, include, alt, keys) + else: + sharder = vcflib.Sharder(args.threads) + try: + contig_lengths = [ + (contig, 0, int(d["length"])) + for contig, d in invcf.contigs.items() + ] + except KeyError: + return ( + "This script requires a VCF with contig lengths in" + " the header when using multiple threads" + ) + shards = sharder.cut(contig_lengths, args.step_size) + sharder.run( + shards, + clip, + [], + invcf, + outvcf, + regions, + reference, + include, + alt, + keys, + ) + + outvcf.close() + invcf.close() + return 0 + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + prog="gvcf_clip.py", + usage="%(prog)s [options] -v GVCF -b BED -r REFERENCE output", + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + parser.add_argument("output", help="Output gVCF file name, or -") + parser.add_argument("-v", "--vcf", required=True, help="Input gVCF") + parser.add_argument("-b", "--bed", required=True, help="Region BED file") + parser.add_argument( + "-r", "--reference", required=True, help="Indexed reference FASTA" + ) + parser.add_argument( + "--exclude", + action="store_true", + help="Keep the positions outside the regions rather than inside", + ) + parser.add_argument( + "-t", + "--threads", + type=int, + default=mp.cpu_count(), + help="number of threads", + ) + parser.add_argument( + "--step-size", + type=int, + default=10 * 1000 * 1000, + help=argparse.SUPPRESS, + ) + sys.exit(main(parser.parse_args())) diff --git a/sentieon_cli/scripts/indel2cnv.py b/sentieon_cli/scripts/indel2cnv.py index 68713f9..e3801e1 100755 --- a/sentieon_cli/scripts/indel2cnv.py +++ b/sentieon_cli/scripts/indel2cnv.py @@ -12,13 +12,8 @@ import platform import mappy import os -import pathlib import re - -# When launched as `python indel2cnv.py ...` from an arbitrary CWD, -# `sentieon_cli` is only importable if the package root is on sys.path. -sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[2])) -from sentieon_cli._indel2cnv_utils import find_peaks # noqa: E402 +import datetime # Convert long INDELs (from assembly-based SV calls) into tandem CNV truth events. # Pipeline: find_repeats (detect period) -> match_ref (locate repeat on reference) @@ -27,6 +22,8 @@ # a truncated copy of (SVLEN < PERIOD) an existing tandem array — both are valid CNVs # because reads from the new copy map back to the reference array, increasing depth. +INDEL2CNV_VERSION = '2.3.0' # combine_sv_cnv release this file belongs to + MIN_SEQ_LEN = 400 REF_SEARCH_RNG = 50000 REF_SEARCH_CNT = 6 @@ -34,6 +31,11 @@ MAX_EXTEND_SEQ = 20000 MAX_SEQ_LEN = 200000 MIN_PERIOD_RATIO = 0.1 # skip if matched portion < 10% of period (coincidental similarity) +MAX_N_FRAC = 0.5 # gate 16: a record at least half reference N is a gap artefact, not a CNV + +def n_frac(seq): + """Fraction of N in a sequence (0 for an empty one).""" + return seq.upper().count('N') / len(seq) if seq else 0.0 log = logging.getLogger(__name__) @@ -91,72 +93,342 @@ def copy_variant(v, **overrides): setattr(new_v, k, val) return new_v -# Detect repeat period(s) in an INDEL sequence via autocorrelation-like scanning. -# Compares var_seq[0:n] vs var_seq[n:2n] at increasing n, finds peaks, refines with -# binary search, and explores harmonic multipliers (2x, 3x, ...) as candidate periods. -def find_repeats(var_seq, thresh): - if len(var_seq) < MIN_SEQ_LEN: - return [] - length = len(var_seq) - step = 200 if length > 2000 else 100 - start = step * 2 - scan_limit = min(length//2+step*5, length-step*1) - dist = [Levenshtein.ratio(var_seq[:n], var_seq[n:n*2]) for n in range(start, scan_limit, step)] - if not dist or max(dist) < thresh: - return [len(var_seq)] - peaks = find_peaks(dist, distance=start//step, height=0.95) - if len(peaks) == 0: - mid = length//2 - if Levenshtein.ratio(var_seq[:mid], var_seq[mid:mid*2]) > thresh: - peaks = [mid/step] - else: - peaks = [p+start//step for p in peaks] - if len(peaks) > 0: - est_period = min([peaks[i+1] - peaks[i] for i in range(len(peaks)-1)] + [peaks[0]]) - pidx = [int(round(p/est_period)) for p in peaks] - period = int((sum([pidx[i] * p for i, p in enumerate(peaks)])/sum([pidx[i]**2 for i, p in enumerate(peaks)]) if len(peaks) > 1 else peaks[0]) * step) - cnt = int(round(length/period)) - if cnt < 2: - return [len(var_seq)] - p_min, p_max = int(length/(cnt+0.5)), int(length/(cnt-0.5)) - p_step = max(1, (p_max - p_min)//10) - while True: - ratios = [max([Levenshtein.ratio(var_seq[i*p:(i+1)*p], var_seq[(i+1)*p:(i+2)*p]) for i in range(cnt-1)]) for p in range(p_min, p_max+1, p_step)] - best_ratio = max(ratios) - best_period = p_min + (ratios.index(best_ratio))*p_step - if p_step == 1: +# --------------------------------------------------------------------------- +# Repeat structure of an indel sequence +# --------------------------------------------------------------------------- +# A tandem expansion or contraction is k copies (k need not be an integer) of a +# unit of length P. Earlier versions scored fixed-offset chunks on a coarse grid +# and could neither land on an arbitrary P nor express a fractional k. Here the +# period candidates come from a k-mer gap histogram (single-base resolution) and +# are verified by aligning the sequence against itself shifted by the candidate: +# the alignment absorbs variable-length copies and partial copies. +KMER = 10 # k-mer size for the gap histogram +MAX_KMER_BINS = 24 # strongest gap bins (by support density) that get verified by alignment +MIN_OVERLAP = 50 # a partial tail counts as evidence from this length ... +MIN_OVERLAP_FRAC = 0.3 # ... and this fraction of the period (>= 1.3 copies). A weak fundamental cannot hurt the + # array walk (it tests the unit-step hypothesis too); it still becomes UNIT, the rung base + # when >= MIN_PERIOD, and the fit-ladder base, so it is a hypothesis, not a fact. +MIN_CORE_COVER = 0.5 # k-mer support must span at least half the sequence +SHIFT_WINDOW = 20000 # verify long sequences on a prefix of the core (>= 3 periods) +REFINE_WINDOW = 0.15 # search +-15% around a k-mer gap for the best period (a bin can sit on the shoulder of a broad peak) +LADDER_SLACK = 0.03 # rungs within this of the best self-identity are kept +ARRAY_PERIOD_MARGIN = 0.02 # an array-level period must beat a microsatellite fundamental by this (capped at 0.99) +MAX_RUNG = 16 +MIN_PERIOD = MIN_SEQ_LEN # reporting floor: a smaller unit is reported as one block + +RepeatInfo = collections.namedtuple('RepeatInfo', 'period ratio core_start core_end candidates array_period') + +_BASE_CODE = {ord('A'): 0, ord('C'): 1, ord('G'): 2, ord('T'): 3, + ord('a'): 0, ord('c'): 1, ord('g'): 2, ord('t'): 3} + + +def _kmer_gap_bins(seq, k=KMER): + """Gaps between consecutive occurrences of each k-mer, merged within 2%. + + Returns [(center, support, first, last)] sorted by support, where + [first, last+k) is the span of positions that support the gap (the core). + """ + n = len(seq) + mask = (1 << (2 * k)) - 1 + last_pos = {} + gaps = collections.Counter() + first, last = {}, {} + h = valid = 0 + for i, ch in enumerate(seq.encode()): + c = _BASE_CODE.get(ch) + if c is None: + h = valid = 0 + continue + h = ((h << 2) | c) & mask + valid += 1 + if valid < k: + continue + start = i - k + 1 + j = last_pos.get(h) + if j is not None: + d = start - j + gaps[d] += 1 + if d not in first: + first[d] = j + last[d] = start + last_pos[h] = start + min_support = max(3, n // 500) + bins = [] + for d, c in sorted(gaps.items(), key=lambda t: -t[1]): + if c < min_support: + break + tol = max(2, int(0.02 * d)) + for b in bins: + if abs(b[0] - d) <= tol: + b[1] += c + b[2] = min(b[2], first[d]) + b[3] = max(b[3], last[d]) break - p_min, p_max = best_period - 2 * p_step, best_period + 2 * p_step - p_step = max(1, (p_max - p_min)//10) - ratio_orig = max([Levenshtein.ratio(var_seq[i*period:(i+1)*period], var_seq[(i+1)*period:(i+2)*period]) for i in range(cnt-1)]) - if best_ratio > ratio_orig: - period = best_period - periods = {period: max(best_ratio, ratio_orig)} - i = 2 + else: + if len(bins) < 200: + bins.append([d, c, first[d], last[d]]) + bins.sort(key=lambda b: -b[1] / max(1, n - b[0])) # support per possible position + return [tuple(b) for b in bins] + + +MAX_PAIRS = 40 # adjacent-copy pairs scored per period (evenly subsampled beyond this) + + +def _shift_ratio(core, d, cutoff=0.0): + """Self-similarity of core at period d: length-weighted mean Levenshtein ratio of + adjacent d-length chunks, plus the partial tail against the start of the previous + chunk when the tail is long enough to mean anything. + + Chunks are distinct substrings, so the alignment cannot cheat by sliding one copy + onto the next (aligning core[:L-d] with core[d:] scores 1-d/(L-d) for any small d). + The core is the k-mer-supported span, so junk outside it does not enter the score. + """ + L = len(core) + if d < 1: + return 0.0 + k, t = divmod(L, d) + if k < 1: + return 0.0 + pairs = [(i * d, (i + 1) * d, d) for i in range(k - 1)] + if len(pairs) > MAX_PAIRS: + step = len(pairs) / MAX_PAIRS + pairs = [pairs[int(j * step)] for j in range(MAX_PAIRS)] + if t >= max(MIN_OVERLAP, MIN_OVERLAP_FRAC * d): + pairs.append(((k - 1) * d, k * d, t)) + if not pairs: + return 0.0 + num = den = 0 + for a, b, w in pairs: + num += w * Levenshtein.ratio(core[a:a + w], core[b:b + w], score_cutoff=cutoff) + den += w + return num / den + + +def _refine_period(core, d0, cutoff): + """Local maximum of the shift ratio around d0, coarse-to-fine, re-centred if it + lands on the window edge.""" + best_d, best_r = d0, _shift_ratio(core, d0, cutoff) + w0 = max(2, int(REFINE_WINDOW * d0)) + res = max(1, d0 // 1000) # final resolution: 1 bp up to 1 kb periods + for _ in range(4): + w, step = w0, max(res, w0 // 4) + edge_lo, edge_hi = best_d - w, best_d + w while True: - cur_period = period * i - if cur_period * 2 > length: + for d in range(best_d - w, best_d + w + 1, step): + if d != best_d: + r = _shift_ratio(core, d, cutoff) + if r > best_r: + best_d, best_r = d, r + if step <= res: break - cur_cnt = int(round(length/cur_period)) - if cur_cnt >= 2: - periods[cur_period] = max([Levenshtein.ratio(var_seq[j*cur_period:(j+1)*cur_period], var_seq[(j+1)*cur_period:(j+2)*cur_period]) for j in range(cur_cnt-1)]) - i += 1 - periods = {k: v for k, v in periods.items() if k >= MIN_SEQ_LEN} - if len(periods) > 1: - min_ratio = max(thresh, max(periods.values()) - 0.03) - periods = {k: v for k, v in periods.items() if v > min_ratio} - return list(periods.keys()) - return [len(var_seq)] + w, step = step, max(res, step // 4) + if edge_lo < best_d < edge_hi: + break + return best_d, best_r + + +def _window(seq, lo, hi, d): + """The core span, cut to a prefix of SHIFT_WINDOW (or 3 periods) for long sequences.""" + c = seq[lo:hi] + lim = max(SHIFT_WINDOW, 3 * d) + return c[:lim] if len(c) > lim else c + + +def analyze_repeats(seq, thresh): + """Fundamental period of an indel sequence and the candidate units to try on + the reference. The whole sequence (one block) is always a candidate.""" + n = len(seq) + if n < MIN_SEQ_LEN: + return RepeatInfo(None, 0.0, 0, n, [], None) + fundamental, f_ratio, core = None, 0.0, (0, n) + bins = sorted(_kmer_gap_bins(seq)[:MAX_KMER_BINS]) + for d0, support, lo, hi in bins: + v = _verify_bin(seq, bins, d0, lo, hi + KMER, thresh) + if v: + fundamental, f_ratio, core = v[0], v[1], (v[2], v[3]) + break + # A microsatellite fundamental hides the array-level structure (e.g. 837 bp units that are + # mostly a 9-mer). Look for the verified period >= MIN_PERIOD whose self-identity is clearly + # higher than the fundamental's; proc_variant tries it when the block does not fit the array. + array_period = None + if fundamental is not None and fundamental < MIN_PERIOD: + # must beat the fundamental by the margin; above 0.97 the margin shrinks and at 0.99 a tie is enough, + # since a near-perfect microsatellite leaves no room above it + best = min(f_ratio + ARRAY_PERIOD_MARGIN, 0.99) + for d0, support, lo, hi in bins: + if d0 < MIN_PERIOD: + continue + v = _verify_bin(seq, bins, d0, lo, hi + KMER, thresh) + if v and v[1] >= best: + array_period, best = v[0], v[1] + return RepeatInfo(fundamental, f_ratio, core[0], core[1], + _unit_candidates(seq, fundamental, f_ratio, core, thresh), array_period) + + +def _verify_bin(seq, bins, d0, lo, hi, thresh): + """Refine a k-mer gap bin to the best nearby period and verify it on its core. + Returns (period, ratio, core_lo, core_hi) or None.""" + n = len(seq) + if hi - lo < MIN_CORE_COVER * n: + return None + c = _window(seq, lo, hi, d0) + if _shift_ratio(c, d0, thresh - 0.15) < thresh - 0.1: + return None + d, r = _refine_period(c, d0, thresh - 0.15) + # the refined period may belong to a neighbouring bin (e.g. a bin sitting on the + # shoulder of the true peak): take the union of the supporting spans and rescore + for e0, _, elo, ehi in bins: + if e0 != d0 and abs(e0 - d) <= max(2, 0.02 * d): + lo, hi = min(lo, elo), max(hi, ehi + KMER) + r = _shift_ratio(_window(seq, lo, hi, d), d, thresh - 0.15) + break + return (d, r, lo, hi) if r >= thresh else None + + +def _unit_candidates(seq, period, ratio, core, thresh): + """Units to match against the reference: the smallest rung of the period ladder + within LADDER_SLACK of the best self-identity (a higher-order repeat wins when its + copies are clearly more similar), the best rung, and the whole sequence as one + block. A sub-MIN_PERIOD fundamental is reported as the block only.""" + n = len(seq) + cands = [] + if period is not None and period >= MIN_PERIOD: + c = seq[core[0]:core[1]] + if len(c) > max(SHIFT_WINDOW, 3 * MAX_RUNG * period): + c = c[:max(SHIFT_WINDOW, 3 * MAX_RUNG * period)] + rungs = [(period, ratio)] + for i in range(2, MAX_RUNG + 1): + p = i * period + if len(c) - p < p: + break + r = _shift_ratio(c, p, thresh - 0.05) + if r >= thresh: + rungs.append((p, r)) + best = max(r for p, r in rungs) + keep = [(p, r) for p, r in rungs if r >= best - LADDER_SLACK] + cands.append(keep[0][0]) + best_p = max(keep, key=lambda t: t[1])[0] + if best_p != cands[0]: + cands.append(best_p) + if n not in cands: + cands.append(n) + return cands + + +def find_repeats(var_seq, thresh): + """Candidate repeat units of an indel sequence (see analyze_repeats).""" + return analyze_repeats(var_seq, thresh).candidates + + +# --------------------------------------------------------------------------- +# Tandem array extension on the reference +# --------------------------------------------------------------------------- +EXTEND_SLACK = 0.0 # array walk continues down to thresh - EXTEND_SLACK (set by --extend_slack) +MAX_EXTEND_STEP = 20000 # cap on the chunk compared per step when the fundamental is known (whole periods) +PARTIAL_WINDOW = 100 # local window that must still match at the end of a partial step +MIN_PARTIAL = 50 + + +def _partial_step(ref_chunk, unit_ext, thresh): + """Longest prefix of ref_chunk that still matches the tiled unit, trimmed so the + last PARTIAL_WINDOW bases match on their own.""" + hi = min(len(ref_chunk), len(unit_ext)) + lo = 0 + while hi - lo > 8: + mid = (lo + hi) // 2 + if Levenshtein.ratio(ref_chunk[:mid], unit_ext[:mid], score_cutoff=thresh) >= thresh: + lo = mid + else: + hi = mid + l, w = lo, PARTIAL_WINDOW + while l >= w and Levenshtein.ratio(ref_chunk[l - w:l], unit_ext[l - w:l], + score_cutoff=thresh - 0.05) < thresh - 0.05: + l -= w // 2 + return l if l >= MIN_PARTIAL else 0 + + +def _step_hypotheses(unit, unit_period): + """Ways the reference array may continue past the matched unit: (chunk, E_right, E_left). + The unit itself (step L) is always one hypothesis. When the indel has a fundamental shorter + than the unit, its tiling in whole periods is another: it keeps the phase when the unit is a + fractional number of copies, and it is cheaper for very long units. The walk picks per side + by the first step, so a fundamental that does not describe the reference array is harmless.""" + L = len(unit) + hyps = [(L, unit, unit)] + if unit_period and 0 < unit_period < L: + C = unit_period * max(1, min(L // unit_period, MAX_EXTEND_STEP // unit_period)) + if C != L: + hyps.append((C, unit[L - C:], unit[:C])) + return hyps + + +EXTEND_THRESH = None # set by --match_thresh: the array walk keeps the main threshold + +def _extend_array(ref, chrom, start, stop, unit_period, thresh): + if EXTEND_THRESH is not None: + thresh = EXTEND_THRESH + """Walk outward from the matched reference unit [start, stop) while the reference keeps + tiling it, then a partial step finds the ragged end of the array.""" + unit = ref.get(chrom, start, stop) + L = len(unit) + if L == 0: + return start, stop + thresh = thresh - EXTEND_SLACK + hyps = _step_hypotheses(unit, unit_period) + # right + best = None + for C, E, _ in hyps: + r = Levenshtein.ratio(ref.get(chrom, stop, stop + C), E, score_cutoff=thresh - 0.1) + if best is None or r > best[0]: + best = (r, C, E) + _, C, E = best + pos = stop + for _ in range(MAX_EXTEND_ITERS): + chunk = ref.get(chrom, pos, pos + C) + if len(chunk) < C or Levenshtein.ratio(chunk, E, score_cutoff=thresh) < thresh: + break + pos += C + pos += _partial_step(ref.get(chrom, pos, pos + C), E, thresh) + # left + best = None + for C, _, E in hyps: + r = Levenshtein.ratio(ref.get(chrom, max(0, start - C), start), E, score_cutoff=thresh - 0.1) + if best is None or r > best[0]: + best = (r, C, E) + _, C, E = best + neg = start + for _ in range(MAX_EXTEND_ITERS): + nxt = neg - C + if nxt < 0: + break + chunk = ref.get(chrom, nxt, neg) + if len(chunk) < C or Levenshtein.ratio(chunk, E, score_cutoff=thresh) < thresh: + break + neg = nxt + lo = max(0, neg - C) + neg -= _partial_step(ref.get(chrom, lo, neg)[::-1], E[::-1], thresh) + return neg, pos + # Align one repeat unit (seq) to a reference window around the INDEL position using mappy. # If only part of seq aligned (partial match), extend the alignment boundary via binary # search on Levenshtein ratio to handle split period boundaries. # Then walk outward in steps of matched_ref_length to find the full tandem repeat array. # Returns: (match_start, match_stop), (array_start, array_stop), (seq_start, seq_stop), ratio -def mm2_match(seq, ref, chrom, pos, half_win, thresh): +_COMP = str.maketrans('ACGTacgtNn', 'TGCAtgcaNn') + +def _revcomp(s): + return s.translate(_COMP)[::-1] + + +def mm2_match(seq, ref, chrom, pos, half_win, thresh, unit_period=None): ref_search_start = max(0, pos - half_win) ref_search_stop = pos + half_win - max_nm_ratio = 1 - thresh + inverted = False + # One identity scale for both matching paths. thresh is a Levenshtein ratio, 1 - edits/(2 * length), + # so a ratio of thresh allows an edit rate of 2 * (1 - thresh); mappy's NM/blen is that edit rate. + # (v1 used 1 - thresh here, which made this path twice as strict as the local one.) + max_nm_ratio = 2 * (1 - thresh) min_match_length = MIN_SEQ_LEN ref_seq = ref.get(chrom, ref_search_start, ref_search_stop) if not ref_seq: @@ -168,7 +440,7 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): # filter hits filtered = [] for h in hits: - if h.strand != 1 or h.is_primary == False: + if h.is_primary == False: continue nm = h.NM if hasattr(h, 'NM') else (h.blen - h.mlen) nm_rate = nm / h.blen if h.blen > 0 else 1.0 @@ -177,6 +449,19 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): filtered.append((h, nm_rate)) if not filtered: return None + # an inverted copy (reverse strand) changes read depth like a direct one; use it only when + # no forward-strand hit passes, so existing behaviour is unchanged elsewhere + if any(h.strand == 1 for h, _ in filtered): + filtered = [(h, nr) for h, nr in filtered if h.strand == 1] + else: + seq = _revcomp(seq) + hits = list(aligner.map(seq)) + filtered = [(h, (h.NM if hasattr(h, 'NM') else h.blen - h.mlen) / h.blen) for h in hits + if h.is_primary and h.strand == 1 and h.blen > min_match_length] + filtered = [(h, nr) for h, nr in filtered if nr < max_nm_ratio] + if not filtered: + return None + inverted = True min_nm_rate = min(f[1] for f in filtered) filtered = [(h, nr) for h, nr in filtered if nr <= min(min_nm_rate + 0.01, max_nm_ratio)] best_hit, best_nm_rate = max(filtered, key=lambda x: x[0].blen) @@ -193,7 +478,7 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): # Partial match extension: when mappy aligned only part of seq, the period boundary # may be split (e.g. seq = [tail|head] of repeat). Extend the match into adjacent # reference to recover the full period, using binary search on Levenshtein ratio. - if abs(len(seq) - matched_seq_length) > 50 and (len(seq)/matched_seq_length-1) < thresh and len(seq) <= MAX_EXTEND_SEQ: + if not inverted and abs(len(seq) - matched_seq_length) > 50 and (len(seq)/matched_seq_length-1) < thresh and len(seq) <= MAX_EXTEND_SEQ: if len(seq) - seq_stop < 10: ref_after = ref.get(chrom, ref_stop, ref_stop + len(seq) - matched_seq_length) step = min(100, len(ref_after)//5) @@ -206,7 +491,7 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): last_invalid = -1 best_ext = 0 while True: - for l in range(start_i, end_i, step): + for l in range(max(start_i, seq_len_diff), end_i, step): ratio = Levenshtein.ratio(best_matched_ref + ref_after[:l], seq[seq_start:]+seq[:l-seq_len_diff]) if ratio > thresh: last_valid = l @@ -220,7 +505,7 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): step //= 2 start_i = last_valid end_i = last_invalid - if last_valid_ratio > 0: + if last_valid_ratio > 0 and best_ext > 0: best_ratio = last_valid_ratio ref_stop += best_ext best_matched_ref += ref_after[:best_ext] @@ -253,7 +538,7 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): step //= 2 start_i = last_valid end_i = last_invalid - if last_valid_ratio > 0: + if last_valid_ratio > 0 and best_ext > 0: best_ratio = last_valid_ratio ref_start -= best_ext best_matched_ref = ref_before[-best_ext:] + best_matched_ref @@ -263,37 +548,28 @@ def mm2_match(seq, ref, chrom, pos, half_win, thresh): matched_ref_length += best_ext # Walk outward from the match to find the full tandem repeat array on reference. # REFSTART/REFSTOP define where read-depth change is expected (used by cnv_eval). - ref_start_ext = ref_start - ref_stop_ext = ref_stop - for _ in range(MAX_EXTEND_ITERS): - ref_start_ext -= matched_ref_length - if ref_start_ext < 0 or Levenshtein.ratio(ref.get(chrom, ref_start_ext, ref_start_ext+matched_ref_length), best_matched_ref) < thresh: - break - for _ in range(MAX_EXTEND_ITERS): - ref_stop_ext += matched_ref_length - if Levenshtein.ratio(ref.get(chrom, ref_stop_ext-matched_ref_length, ref_stop_ext), best_matched_ref) < thresh: - break + ref_start_ext, ref_stop_ext = _extend_array(ref, chrom, ref_start, ref_stop, unit_period, thresh) + if inverted: # query coordinates back in the insertion's own orientation + seq_start, seq_stop = len(seq) - seq_stop, len(seq) - seq_start + mm2_match.last_inverted = inverted return (ref_start, ref_stop), (ref_start_ext, ref_stop_ext), (seq_start, seq_stop), best_ratio # Try to match one repeat unit in the reference immediately adjacent to the INDEL. # For insertions: align alt_seq[:period] against ref_after[:p] + ref_before[-(period-p):] # at varying split positions p (binary search). Handles wrap-around period boundaries. # For deletions: the deleted ref sequence itself is the repeat unit. -def local_match(v, ref, period, thresh): +def local_match(v, ref, period, thresh, unit_period=None): if len(v.ref) < len(v.alt[0]): alt_seq = v.alt[0] ref_before = ref.get(v.chrom, v.pos-period, v.pos) ref_after = ref.get(v.chrom, v.pos, v.pos+period) step = max(period//40, 100) - cnt = int(round(len(alt_seq)/period)) - offset = len(alt_seq) - cnt * period - if offset > 0: - ref_after = alt_seq[-offset:] + ref_after - elif offset < 0: - ref_after = ref_after[-offset:] + # v1 prepended the ALT's partial last copy to ref_after when the copy count was fractional. + # That scored ALT bases against themselves (DIST 0.94 where the reference alone gives 0.78) + # and reported a window shifted by the partial copy; only reference bases are compared now. alt_period_seq = alt_seq[:period] - ratios = {p: Levenshtein.ratio(ref_after[:p] + ref_before[-(period-p):], alt_period_seq) for p in range(0, period, step)} - ratios[period] = Levenshtein.ratio(ref_after, alt_period_seq) + ratios = {p: Levenshtein.ratio(ref_after[:p] + ref_before[-(period-p):], alt_period_seq, score_cutoff=0.8) for p in range(0, period, step)} + ratios[period] = Levenshtein.ratio(ref_after, alt_period_seq, score_cutoff=0.8) while True: best_pos = max(ratios, key=ratios.get) max_ratio = ratios[best_pos] @@ -306,41 +582,91 @@ def local_match(v, ref, period, thresh): p_hi = min(best_pos + step, period) for p in (p_lo, p_hi): if p not in ratios: - ratios[p] = Levenshtein.ratio(ref_after[:p] + ref_before[-(period-p):], alt_period_seq) + ratios[p] = Levenshtein.ratio(ref_after[:p] + ref_before[-(period-p):], alt_period_seq, score_cutoff=0.8) if max_ratio < thresh: return None ref_start = v.pos - (period - best_pos) ref_end = v.pos + best_pos - ref_seq = ref_before[-(period-best_pos):] + alt_seq[:best_pos] if ref_start != v.pos else alt_period_seq else: - ref_start = v.pos - ref_end = v.pos + period - ref_seq = v.ref[:period] + ref_start = v.pos + 1 # first deleted base (v.pos is the anchor) + ref_end = v.pos + 1 + period max_ratio = 1. - ref_start_ext = ref_start - ref_end_ext = ref_end - for _ in range(MAX_EXTEND_ITERS): - ref_start_ext -= period - if ref_start_ext < 0 or Levenshtein.ratio(ref.get(v.chrom, ref_start_ext, ref_start_ext+period), ref_seq) < thresh: - break - for _ in range(MAX_EXTEND_ITERS): - ref_end_ext += period - if Levenshtein.ratio(ref.get(v.chrom, ref_end_ext-period, ref_end_ext), ref_seq) < thresh: - break + # the array walk compares the reference itself over [ref_start, ref_end) (v1 used a hybrid of + # reference and ALT bases for insertions; they agree wherever the inserted copy matches) + ref_start_ext, ref_end_ext = _extend_array(ref, v.chrom, ref_start, ref_end, unit_period, thresh) return (ref_start, ref_end), (ref_start_ext, ref_end_ext), (0, period), max_ratio # Try local_match first (fast, exact adjacent match). Fall back to mm2_match # (mappy alignment in a wider window) for segdups where the repeat may be distant. -def match_ref(v, ref, period, thresh): - local_result = local_match(v, ref, period, thresh) +def match_ref(v, ref, period, thresh, unit_period=None): + mm2_match.last_inverted = False + local_result = local_match(v, ref, period, thresh, unit_period) if local_result is not None: return local_result - return mm2_match(v.alt[0][:period], ref, v.chrom, v.pos, min(max(period, len(v.alt[0]), len(v.ref))*REF_SEARCH_CNT, REF_SEARCH_RNG), thresh) + return mm2_match(v.alt[0][:period], ref, v.chrom, v.pos, min(max(period, len(v.alt[0]), len(v.ref))*REF_SEARCH_CNT, REF_SEARCH_RNG), thresh, unit_period) + +RATIO_SLACK = 0.02 # match ratios within this of the best are treated as equal +NET_RECIPROCAL = 0.8 # a DUP and a DEL whose arrays overlap reciprocally by this much are one array and net +EXACT_FRAC = 0.05 # |round(n/p)*p - n| within this fraction of n counts as a whole copy count + +def select_candidate(results, periods, indel_len): + """Pick one matched unit. Among matches within RATIO_SLACK of the best ratio prefer the + longest reference array, then a whole copy count (the block is exact by construction), + then the smallest unit. results: (idx, match, array, seq_range, ratio).""" + best_ratio = max(r[4] for r in results) + pool = [r for r in results if best_ratio - r[4] <= RATIO_SLACK] + max_ext = max(r[2][1] - r[2][0] for r in pool) + ext_tol = 0.5 * min(periods[r[0]] for r in pool) + def key(r): + p = periods[r[0]] + ext = r[2][1] - r[2][0] + err = abs(int(indel_len / p + 0.5) * p - indel_len) + return (max_ext - ext > ext_tol, err > EXACT_FRAC * indel_len, p, -r[4]) + return min(pool, key=key) + +def _indel_len(v, is_insert): + return len(v.alt[0]) if is_insert else len(v.ref) - len(v.alt[0]) + + +def _fit_unit(sv, ref, unit, n, thresh, max_len=None): + """An insertion longer than the reference array cannot match as one block. Find the + largest multiple of the fundamental unit (below n, at least MIN_SEQ_LEN) that matches + the reference next to the insertion. The match is roughly monotone in the unit length + (it fails once the unit outgrows the array), but local_match's split search is a + heuristic, so probe a coarse ladder first and binary search from the largest success. + Returns (period, match) or None.""" + top = n - 1 if not max_len else min(n - 1, max_len) # a unit cannot be longer than the array + lo, hi = -(-MIN_SEQ_LEN // unit), top // unit + if lo > hi: + return None + probes = sorted({lo + (hi - lo) * i // 6 for i in range(7)}, reverse=True) + best = None + fail = hi + 1 + for k in probes: + r = local_match(sv, ref, k * unit, thresh, unit) + if r is not None: + best, lo = (k * unit, r), k + break + fail = k + if best is None: + return None + while fail - lo > 1: + mid = (lo + fail) // 2 + r = local_match(sv, ref, mid * unit, thresh, unit) + if r is not None: + lo, best = mid, (mid * unit, r) + else: + fail = mid + return best + # Process one input INDEL variant into CNV call(s). # Multi-allelic variants are split; nearby het INS with opposite phase are merged. # For each allele: find_repeats -> match_ref -> pick best period -> emit DUP/DEL. -def proc_variant(v, ref, thresh): +def proc_variant(v, ref, thresh, match_thresh=None): + # thresh: period verification (and, unless --match_thresh is given, matching and the array walk); + # match_thresh: the identity a local or mappy match must reach to be accepted + mt = thresh if match_thresh is None else match_thresh stime = time.time() result_vs = [] # Split multi-allelic into per-allele variants, trim shared prefix/suffix @@ -348,7 +674,11 @@ def proc_variant(v, ref, thresh): vs = [] is_insert = [] if split_var: + # Only process alleles referenced by GT + gt_indices = set(int(g) for g in re.split(r'[|/]', v.samples[0].get('GT', '0/0')) if g not in ('.', '0')) for i, a in enumerate(v.alt): + if (i + 1) not in gt_indices: + continue if a == '*': continue if len(a) < len(v.ref): @@ -357,7 +687,7 @@ def proc_variant(v, ref, thresh): if len(a) > 1: if a == v1.ref[:len(a)] or len(a) < 10 or Levenshtein.ratio(a, v1.ref[:len(a)]) > thresh: v1.pos += len(a) - 1 - v1.ref = v.ref[len(a):] + v1.ref = v.ref[len(a) - 1:] # the new anchor is the allele's last base v1.alt = [a[-1]] elif Levenshtein.ratio(a, v.ref[-len(a):]) > thresh: v1.ref = v.ref[:-len(a)] @@ -383,6 +713,8 @@ def proc_variant(v, ref, thresh): vs = [v] is_insert = [len(v.alt[0]) > len(v.ref)] periods = [] + units = [] + aperiods = [] for si, vv in enumerate(vs): if not is_insert[si]: del_len = len(vv.ref) - len(vv.alt[0]) @@ -393,35 +725,51 @@ def proc_variant(v, ref, thresh): end_pos = vv.pos + len(vv.ref) vv.info = {'SVTYPE': 'DEL', 'CNDIFF': -gt_sum, 'POSSHIFT': 0, 'POSSHIFT2': 0, 'DIST': 1.0, 'SVLEN': del_len, 'PERIOD': del_len, - 'REFSTART': vv.pos + 1, 'REFSTOP': end_pos, 'END': end_pos, + 'REFSTART': vv.pos + 2, 'REFSTOP': end_pos, 'END': end_pos, 'INDEL': '%s:%d' % (vv.chrom, vv.pos + 1)} - vv.end = end_pos - 1 + vv.end = end_pos vv.ref = ref.get(vv.chrom, vv.pos, vv.pos + 1) vv.alt = [''] vv.line = None result_vs.append(vv) periods.append([]) + units.append(None) + aperiods.append(None) else: - periods.append(find_repeats(vv.ref[1:], thresh)) + info = analyze_repeats(vv.ref[1:], thresh) + periods.append(info.candidates) + units.append(info.period) + aperiods.append(info.array_period) else: periods.append([]) + units.append(None) + aperiods.append(None) else: if len(vv.alt[0]) - len(vv.ref) >= MIN_SEQ_LEN and len(vv.alt[0]) <= MAX_SEQ_LEN: - periods.append(find_repeats(vv.alt[0], thresh)) + info = analyze_repeats(vv.alt[0], thresh) + periods.append(info.candidates) + units.append(info.period) + aperiods.append(info.array_period) else: periods.append([]) + units.append(None) + aperiods.append(None) if len(sum(periods, [])) == 0: - return [], time.time() - stime - # If both alleles are same direction with matching periods, merge as homozygous - if (all(is_insert) or not any(is_insert)) and len(periods) > 1 and min(len(p) for p in periods) > 0: + return result_vs, time.time() - stime + # If exactly 2 alleles, same direction with matching periods, merge as homozygous + if (all(is_insert) or not any(is_insert)) and len(periods) == 2 and min(len(p) for p in periods) > 0 \ + and abs(_indel_len(vs[0], is_insert[0]) / _indel_len(vs[1], is_insert[1]) - 1) <= 0.05: updated_periods = [] for p0 in periods[0]: p1s = [p1 for p1 in periods[1] if abs(p0/p1-1) <= 0.05] - if p1s and Levenshtein.ratio(*(a[:min(p0, p1s[0])] for a in v.alt)) > 1 - (1-thresh)/2: + if p1s and Levenshtein.ratio(vs[0].alt[0][:min(p0, p1s[0])], vs[1].alt[0][:min(p0, p1s[0])]) > 1 - (1-thresh)/2: updated_periods.append(min(p0, p1s[0])) if updated_periods: periods = [updated_periods] - vs = [vs[1]] if len(vs[1].alt[0]) > len(vs[0].alt[0]) else [vs[0]] + keep = 1 if len(vs[1].alt[0]) > len(vs[0].alt[0]) else 0 + vs = [vs[keep]] + units = [units[keep] or units[1 - keep]] + aperiods = [aperiods[keep] or aperiods[1 - keep]] vs[0].samples[0]['GT'] = '1/1' for si, sv in enumerate(vs): @@ -429,45 +777,52 @@ def proc_variant(v, ref, thresh): if not period: continue results = [] + inverted = {} for j, p in enumerate(period): try: - result = match_ref(sv, ref, p, thresh) + result = match_ref(sv, ref, p, mt, units[si]) except Exception as e: log.warning("match_ref failed for %s:%d period=%d: %s", sv.chrom, sv.pos, p, e) result = None if result: results.append((j, *result)) + inverted[j] = mm2_match.last_inverted + indel_len = len(sv.alt[0]) if is_insert[si] else len(sv.ref) - len(sv.alt[0]) + # The whole insertion did not match locally as one block (longer than the array, + # or only a mappy partial match): try the largest unit that fits the array. + if is_insert[si] and units[si] and not any(period[r[0]] == indel_len and r[3] == (0, indel_len) for r in results): + fit = None + if aperiods[si]: # the array-level period of a microsatellite-rich indel + for p in (aperiods[si], 2 * aperiods[si]): + if p < indel_len: + m = local_match(sv, ref, p, mt, p) + if m: + fit = (p, m) + break + if fit is None: + # bound the ladder by twice the array already measured (a larger unit can + # extend further than the smaller candidates did, but not without limit) + max_ext = max((r[2][1] - r[2][0] for r in results), default=None) + fit = _fit_unit(sv, ref, units[si], indel_len, mt, 2 * max_ext if max_ext else None) + if fit: + period.append(fit[0]) + results.append((len(period) - 1, *fit[1])) + # Filter: matched portion must be large enough and a meaningful fraction of period + results = [r for r in results if r[3][1] - r[3][0] >= MIN_SEQ_LEN + and (r[3][1] - r[3][0]) / period[r[0]] >= MIN_PERIOD_RATIO] if not results: continue - # Pick the best period: highest Levenshtein ratio, break ties by matched length - if len(results) == 1: - result = results[0] - else: - best_ratios = [r[4] for r in results] - best_ratio = max(best_ratios) - best_i = best_ratios.index(best_ratio) - result = results[best_i] - best_seq_length = result[3][1] - result[3][0] - for ri, r in enumerate(results): - if r and best_i != ri and best_ratio - r[4] < 0.01: - seq_length = r[3][1] - r[3][0] - if seq_length > best_seq_length: - result = r - best_i = ri - best_seq_length = seq_length + result = select_candidate(results, period, indel_len) ri, best_match, ref_range, seq_range, best_ratio = result best_period = period[ri] if is_insert[si]: period_cnt = len(sv.alt[0]) / best_period else: - overlap = max(0, min(sv.pos + len(sv.ref), ref_range[1]) - max(sv.pos, ref_range[0])) + overlap = max(0, min(sv.pos + len(sv.ref), ref_range[1]) - max(sv.pos + 1, ref_range[0])) period_cnt = overlap / best_period - period_cnt = int(round(period_cnt)) + period_cnt = max(1, int(period_cnt + 0.5)) seq_length = seq_range[1] - seq_range[0] - # Filter: matched portion must be large enough and a meaningful fraction of period - if seq_length < MIN_SEQ_LEN or seq_length / best_period < MIN_PERIOD_RATIO: - continue gt = re.split(r"\||/", sv.samples[0]['GT'].replace('.', '0')) gt_sum = sum(int(g) for g in gt) is_dup = is_insert[si] @@ -482,11 +837,14 @@ def proc_variant(v, ref, thresh): 'DIST': best_ratio, 'SVLEN': seq_length, 'PERIOD': best_period, + 'UNIT': units[si] or best_period, 'REFSTART': ref_range[0] + 1, 'REFSTOP': ref_range[1] + 1, - 'END': best_match[1] + 1, + 'END': best_match[1], 'INDEL': indel_pos} - sv.pos = best_match[0] + if inverted.get(ri): + sv.info['INV'] = 1 + sv.pos = best_match[0] if is_dup else best_match[0] - 1 # a DEL's POS is the base before the deletion sv.end = best_match[1] sv.ref = ref.get(sv.chrom, sv.pos, sv.pos+1) sv.alt = [''] if is_dup else [''] @@ -494,11 +852,15 @@ def proc_variant(v, ref, thresh): result_vs.append(sv) return result_vs, time.time() - stime -def proc_batch(ref, variants, thresh): +def proc_batch(ref, variants, thresh, extend_slack=0.0, match_thresh=None): + global EXTEND_THRESH + EXTEND_THRESH = thresh if match_thresh is not None else None + global EXTEND_SLACK + EXTEND_SLACK = extend_slack results = [] for v in variants: try: - results.append(proc_variant(v, ref, thresh)) + results.append(proc_variant(v, ref, thresh, match_thresh)) except Exception as e: log.error("Failed processing %s:%d: %s", v.chrom, v.pos, e) results.append(([], 0.0)) @@ -510,15 +872,18 @@ def merge(ref, v1, v2): return None if set([v1.samples[0]['GT'], v2.samples[0]['GT']]) != set(['1|0', '0|1']): return None - if len(v1.alt) != len(v2.alt) or len(v2.alt) != 1: - return None - if len(v1.ref) > len(v1.alt[0]) or len(v1.ref) > len(v2.alt[0]): + if len(v1.alt) != 1 or len(v2.alt) != 1: return None + if len(v1.ref) > len(v1.alt[0]) or len(v2.ref) > len(v2.alt[0]): + return None # both must be insertions ref_gap = ref.get(v1.chrom, v1.end, v2.pos+1) v1.ref += ref_gap v1.alt[0] += ref_gap v1.alt.append(v1.ref[:-1] + v2.alt[0]) - v1.samples[0]['GT'] = '1|1' + # allele 1 sits on the haplotype v1 was on, allele 2 on the other (v1 used to write 1|1, + # which made proc_variant drop the second allele) + v1.samples[0]['GT'] = '1|2' if v1.samples[0]['GT'] == '1|0' else '2|1' + v1.end = v2.pos + 1 return v1 def gt_set(gt): @@ -531,8 +896,10 @@ def extend_v(v): period = v.info['PERIOD'] if not cndiff or not period: return v - num_periods = round(v.info['SVLEN'] / period) + num_periods = max(1, round(v.info['SVLEN'] / period)) gt_sum = sum(gt_set(v.samples[0]['GT'])) or 1 + if cndiff % gt_sum: + gt_sum = 1 # a merged count that is not per-haplotype: re-express it as a whole diff = cndiff // gt_sum max_num_period = (v.info['REFSTOP'] - v.info['REFSTART']) / period if max_num_period - round(max_num_period) < 0.1: @@ -548,10 +915,12 @@ def extend_v(v): if abs(diff) == best_d or best_d == -1: return v svlen_periods = round(num_periods * abs(diff) / best_d) - v.info['CNDIFF'] = best_d if cndiff > 0 else -best_d + # keep the total change in base pairs: CNDIFF counts copies over all haplotypes + v.info['CNDIFF'] = best_d * gt_sum if cndiff > 0 else -best_d * gt_sum v.info['SVLEN'] = svlen_periods * period - v.pos = v.info['REFSTART'] + v.pos = v.info['REFSTART'] - 1 # 0-based start of the array v.end = v.pos + v.info['SVLEN'] + v.info['END'] = v.end v.line = None return v @@ -559,25 +928,37 @@ def extend_v(v): # Handles: (1) opposite-phase same-type pairs -> hom, (2) opposite-type pairs -> net CNDIFF, # (3) same-type DEL/DUP pairs with compatible periods -> combined CNDIFF. def merge_output_vars(variants): - if not variants: - return [] + """Records arrive sorted by POS. Each one is tried against every kept record whose array it + overlaps, not only the previous one (a DUP formed by merging two DUPs must still net against a DEL + that came before them), and a merged record is tried again until nothing more merges.""" result = [] - last = variants[0] - for v in variants[1:]: - if last is None: - last = v - continue - merged = _try_merge_output(last, v) - if merged == 'drop': - last = None - elif merged: - last = merged - else: - result.append(last) - last = v - if last is not None: - result.append(last) - return result + for v in variants: + cur = v + while cur is not None: + for i in range(len(result) - 1, -1, -1): + merged = _try_merge_output(result[i], cur) + if merged is None: + continue + del result[i] + cur = None if merged == 'drop' else merged + break + else: + result.append(cur) + cur = None + # An inverted copy is used only to cancel the deletion side of an inversion (assemblies write an + # inversion as a deletion plus an inverted insertion, sometimes in several pieces). On its own it is + # not reported: read depth supports 19 of 45 such records on NA21102, i.e. none beyond chance. + for v in result: + v.info.pop('_BP', None) + return [v for v in sorted(result, key=lambda v: v.pos) if not v.info.get('INV')] + +def _event_bp(v): + """Signed change in reference-matching base pairs over all haplotypes. CNDIFF counts copies + of the matched unit: SVLEN, which equals PERIOD for a full local match, the matched part of + the unit for a partial mappy match (the rest of the insertion is novel sequence and changes no + depth), and the re-expressed unit after extend_v.""" + return v.info.get('_BP', v.info['CNDIFF'] * v.info['SVLEN']) + def _try_merge_output(v1, v2): if v1.chrom != v2.chrom: @@ -591,47 +972,73 @@ def _try_merge_output(v1, v2): if len1 <= 0 or len2 <= 0: return None r1, r2 = overlap / len1, overlap / len2 - if max(r1, r2) < 0.9: - return None - # Opposite types (DUP+DEL or DEL+DUP): net the CNDIFF - if v1.alt[0] != v2.alt[0] and min(r1, r2) > 0.9: - net = v1.info['CNDIFF'] + v2.info['CNDIFF'] - if net == 0: + # Opposite types (DUP+DEL or DEL+DUP) on one array net their change in base pairs. CNDIFF counts + # copies of each record's own unit, so counts are only comparable through the unit. "One array" + # is a reciprocal overlap of NET_RECIPROCAL: the two records' arrays were extended independently + # and differ by more than 10% at ~30 pairs genome-wide that are plainly the same locus. + if v1.alt[0] != v2.alt[0]: + if min(r1, r2) < NET_RECIPROCAL: + return None + bp1, bp2 = _event_bp(v1), _event_bp(v2) + net = bp1 + bp2 + if bool(v1.info.get('INV')) != bool(v2.info.get('INV')) and abs(net) > 0.1 * max(abs(bp1), abs(bp2)): + return None # an inverted copy cancels only a deletion of its own size + # cancel when the residual is below the event floor or 5% of the larger event (an inversion + # represented as a deletion plus an inverted insertion leaves a few unmatched bases per copy) + if abs(net) < max(MIN_SEQ_LEN, EXACT_FRAC * max(abs(bp1), abs(bp2))): return 'drop' # cancel out - keep = v1 if abs(v1.info['CNDIFF']) > abs(v2.info['CNDIFF']) else v2 - keep.info['CNDIFF'] = net + keep = v1 if abs(bp1) >= abs(bp2) else v2 + # count the residual in the finer unit of the two, bounded by the array: a one-copy deletion must + # not vanish by rounding against a duplication re-expressed as one copy of seven units + unit = min(v1.info['PERIOD'], v1.info['SVLEN'], v2.info['PERIOD'], v2.info['SVLEN'], + keep.info['REFSTOP'] - keep.info['REFSTART']) + keep.info['PERIOD'] = keep.info['SVLEN'] = unit + keep.info['CNDIFF'] = (1 if net > 0 else -1) * max(1, int(abs(net) / unit + 0.5)) + keep.info['_BP'] = net # exact total, so a chain of merges rounds once + keep = extend_v(keep) keep.line = None return keep - # Same type DEL + # Same-type merges also need one array. A record nested inside a much larger array is a distinct event + # (a 700 bp deletion of one sub-array inside a 720 kb satellite deletion); adding it would drive the + # merged unit down to the smaller record's and count the large deletion in hundreds of tiny units. + if min(r1, r2) < NET_RECIPROCAL: + return None + # Same type DEL: two deletions on one array add up in base pairs (the two alleles of a 1|2 + # site are two events, one per haplotype; the old rule kept one of them when the GTs matched) if v1.alt[0] == '' and v2.alt[0] == '': - keep = v1 if len1 >= len2 else v2 - if keep.info['CNDIFF'] == -2 or v1.samples[0]['GT'] == v2.samples[0]['GT']: - return keep - if min(r1, r2) > 0.9: - gt1, gt2 = gt_set(v1.samples[0]['GT']), gt_set(v2.samples[0]['GT']) - dot = '.' in v1.samples[0]['GT'] or '.' in v2.samples[0]['GT'] - if not dot: - gt = [gt1[i] | gt2[i] for i in range(min(len(gt1), len(gt2)))] - keep.samples[0]['GT'] = '|'.join(str(g) for g in gt) - else: - keep.samples[0]['GT'] = str(gt1[0] | gt2[0]) + '/.' - keep.info['CNDIFF'] = min(-2, v1.info['CNDIFF'] + v2.info['CNDIFF']) - keep.line = None - return keep - # Same type DUP with compatible periods + keep, other = (v1, v2) if len1 >= len2 else (v2, v1) + # count in the finer unit so a nested unequal deletion is not rounded away + unit = min(v1.info['PERIOD'], v1.info['SVLEN'], v2.info['PERIOD'], v2.info['SVLEN'], + keep.info['REFSTOP'] - keep.info['REFSTART']) + total = _event_bp(keep) + _event_bp(other) + gt1, gt2 = gt_set(v1.samples[0]['GT']), gt_set(v2.samples[0]['GT']) + dot = '.' in v1.samples[0]['GT'] or '.' in v2.samples[0]['GT'] + if not dot and len(gt1) == len(gt2): + keep.samples[0]['GT'] = '|'.join(str(g) for g in (gt1[i] | gt2[i] for i in range(len(gt1)))) + elif gt1 and gt2: + keep.samples[0]['GT'] = str(gt1[0] | gt2[0]) + '/.' + keep.info['PERIOD'] = keep.info['SVLEN'] = unit + keep.info['CNDIFF'] = -max(1, int(abs(total) / unit + 0.5)) + keep.info['_BP'] = total + keep.line = None + return keep + # Same type DUP with compatible periods (multiples of each other, or the same fundamental unit) if v1.alt[0] == '' and v2.alt[0] == '': + if bool(v1.info.get('INV')) != bool(v2.info.get('INV')): + return None # an inverted copy is not reported; keep it out of reported records p1, p2 = v1.info['PERIOD'], v2.info['PERIOD'] r = max(p1, p2) / min(p1, p2) if min(p1, p2) > 0 else 999 - if abs(r - round(r)) > 0.1: - log.warning("Incompatible periods %d vs %d at %s:%d", p1, p2, v1.chrom, v1.pos) + u1, u2 = v1.info.get('UNIT', 0), v2.info.get('UNIT', 0) + if abs(r - round(r)) > 0.1 and not (u1 and u2 and abs(u1 / u2 - 1) <= 0.05): + log.debug("Incompatible periods %d vs %d at %s:%d", p1, p2, v1.chrom, v1.pos) return None - period = min(p1, p2) keep, other = (v1, v2) if len1 >= len2 else (v2, v1) - # Normalize CNDIFF to common period; use round to avoid truncating partial matches - n1 = max(1, round(keep.info['SVLEN'] / period)) - n2 = max(1, round(other.info['SVLEN'] / period)) - cd1 = keep.info['CNDIFF'] * n1 - cd2 = other.info['CNDIFF'] * n2 + # the common unit: the finer of the two (SVLEN where a partial match carries PERIOD > SVLEN), or + # the array when a partially matched block exceeds it (depth over an array smaller than the unit + # changes by base pairs / array) + period = min(p1, v1.info['SVLEN'], p2, v2.info['SVLEN'], keep.info['REFSTOP'] - keep.info['REFSTART']) + cd1 = int(_event_bp(keep) / period + 0.5) + cd2 = int(_event_bp(other) / period + 0.5) gt1, gt2 = gt_set(v1.samples[0]['GT']), gt_set(v2.samples[0]['GT']) if gt1 == gt2 or (cd1 == cd2 and [1, 1] not in (gt1, gt2)): dot = '.' in v1.samples[0]['GT'] or '.' in v2.samples[0]['GT'] @@ -640,15 +1047,43 @@ def _try_merge_output(v1, v2): keep.samples[0]['GT'] = '|'.join(str(g) for g in gt) else: keep.samples[0]['GT'] = str(gt1[0] | gt2[0]) + '/.' - keep.info['PERIOD'] = period - keep.info['SVLEN'] = period - keep.info['CNDIFF'] = cd1 + cd2 + total = _event_bp(keep) + _event_bp(other) # before the unit changes + keep.info['PERIOD'] = keep.info['SVLEN'] = period + keep.info['CNDIFF'] = max(1, int(total / period + 0.5)) + keep.info['_BP'] = total keep = extend_v(keep) keep.line = None return keep return None +def tool_version(): + """COMBINE_SV_CNV_VERSION when a driver set it, else combine_sv_cnv-v, followed by + ' (sentieon-cli-)' when this file sits in the sentieon_cli package (sentieon_cli/scripts/).""" + if os.environ.get('COMBINE_SV_CNV_VERSION'): + return os.environ['COMBINE_SV_CNV_VERSION'] + ver = 'combine_sv_cnv-v' + INDEL2CNV_VERSION + if os.path.basename(os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) == 'sentieon_cli': + try: + from sentieon_cli import __version__ as scli_version + except Exception: + try: + from importlib.metadata import version + scli_version = version('sentieon_cli') + except Exception: + scli_version = 'unknown' + ver += ' (sentieon-cli-%s)' % scli_version + return ver + +def provenance_line(): + """The ##CommandLine header line, in the form of Sentieon's ##SentieonCommandLine.""" + date = datetime.datetime.now(datetime.timezone.utc).strftime('%Y-%m-%dT%H:%M:%SZ') + cmd = ' '.join(sys.argv).replace('\\', '\\\\').replace('"', '\\"') + return '##CommandLine.indel2cnv=' % (tool_version(), date, cmd) + def main(args): + global EXTEND_SLACK, EXTEND_THRESH + EXTEND_SLACK = args.extend_slack + EXTEND_THRESH = args.thresh if args.match_thresh is not None else None log_level = logging.DEBUG if args.verbose else logging.INFO logging.basicConfig(level=log_level, format='%(asctime)s %(levelname)s %(message)s', datefmt='%H:%M:%S') @@ -658,8 +1093,9 @@ def main(args): update = ('##INFO=', '##INFO=', '##INFO=', - '##INFO=', + '##INFO=', '##INFO=', + '##INFO=', '##INFO=', '##INFO=', '##INFO=', @@ -668,6 +1104,7 @@ def main(args): '##ALT=', '##INFO=', '##INFO=', + provenance_line(), ) fout.copy_header(input_vcf, update=update) fout.emit_header() @@ -682,10 +1119,20 @@ def main(args): contigs = all_contigs t0 = time.time() variants = [] + n_skip_in = n_skip_out = 0 + skipped_spans = collections.defaultdict(list) # gate 16: what the skipped records covered for c in contigs: vcf_c = input_vcf.range(c) for v in vcf_c: if (len(v.filter) == 0 or 'PASS' in v.filter) and max(len(a) for a in [v.ref] + v.alt) >= MIN_SEQ_LEN: + if len(v.ref) >= max((len(a) for a in v.alt), default=0): + nf = n_frac(v.ref) + if nf >= MAX_N_FRAC: + log.info('skipped N-gap record %s:%d len=%d nfrac=%.2f', v.chrom, v.pos + 1, len(v.ref), nf) + n_skip_in += 1 + # the deleted span: the REF allele without its anchor base, 0-based half-open + skipped_spans[v.chrom].append((v.pos + 1, v.pos + len(v.ref), 'in:%.2f' % nf)) + continue if variants: v_new = merge(ref, variants[-1], v) if v_new: @@ -696,7 +1143,11 @@ def main(args): if args.threads == 1: all_vars = [] for v in variants: - vs, t = proc_variant(v, ref, args.thresh) + try: + vs, t = proc_variant(v, ref, args.thresh, args.match_thresh) + except Exception as e: + log.error("Failed processing %s:%d: %s", v.chrom, v.pos, e) + continue all_vars += vs if vs: for sv in vs: @@ -709,15 +1160,16 @@ def main(args): n_threads = args.threads or os.cpu_count() # Sort by expected difficulty (sequence length) descending so that # expensive variants are dispatched first and don't pile up at the end. - # Use small batches (5) for better load balancing across workers. + # Hand out one variant at a time: with larger batches, or starmap's + # default chunksize, the most expensive variants all land on one worker. variants.sort(key=lambda v: max(len(a) for a in [v.ref] + v.alt), reverse=True) - batch_size = 5 - args_in = [] - for i in range(0, len(variants), batch_size): - args_in.append((ref, variants[i:i+batch_size], args.thresh)) + args_in = [ + (ref, [v], args.thresh, args.extend_slack, args.match_thresh) + for v in variants + ] all_vars = [] with Pool(n_threads) as p: - for result in p.starmap(proc_batch, args_in): + for result in p.starmap(proc_batch, args_in, chunksize=1): for vs, t in result: all_vars += vs # Merge overlapping output calls (same array, different haplotypes) and write @@ -725,6 +1177,13 @@ def main(args): out_c = sorted([v for v in all_vars if v.chrom == c], key=lambda v: v.pos) out_c = merge_output_vars(out_c) for v in sorted(out_c, key=lambda v: (v.pos, v.end)): + nf = n_frac(ref.get(v.chrom, v.info['REFSTART'] - 1, v.info['REFSTOP'] - 1)) + if nf >= MAX_N_FRAC: + log.info('skipped N-gap record %s:%d len=%d nfrac=%.2f', v.chrom, v.pos + 1, + v.info['REFSTOP'] - v.info['REFSTART'], nf) + n_skip_out += 1 + skipped_spans[v.chrom].append((v.info['REFSTART'] - 1, v.info['REFSTOP'] - 1, 'out:%.2f' % nf)) + continue fout.emit(v) t1 = time.time() mm, ut, st = 0, 0, 0 @@ -734,6 +1193,13 @@ def main(args): mm += ru.ru_maxrss * mem_scale ut += ru.ru_utime st += ru.ru_stime + log.info('N-gap records skipped: %d on input, %d on output', n_skip_in, n_skip_out) + if args.skipped_bed: + with open(args.skipped_bed, 'w') as fb: + for c in contigs: + for bs, be, name in sorted(skipped_spans.get(c, ())): + fb.write('%s\t%d\t%d\t%s\n' % (c, bs, be, name)) + log.info('wrote the skipped spans to %s', args.skipped_bed) log.info('overall: %d mem %.3f user %.3f sys %.3f real', mm, ut, st, t1-t0) fout.close() input_vcf.close() @@ -745,7 +1211,13 @@ def main(args): parser.add_argument('out_vcf', help='Output VCF file name') parser.add_argument('--contig', help="Contigs to process (comma-separated)") parser.add_argument('-t', '--threads', help="Concurrent processes (default: nproc)", default=None, type=int) - parser.add_argument('--thresh', help="Minimum Levenshtein similarity ratio", default=0.9, type=float) + parser.add_argument('--thresh', help="Minimum Levenshtein similarity ratio of a match, 1 - edits/(2 * length), applied identically to the local and mappy paths (0.9 is about 80%% identity)", default=0.9, type=float) + parser.add_argument('--extend_slack', help="Array extension continues down to thresh minus this", default=0.0, type=float) + parser.add_argument('--match_thresh', help="Identity a match must reach (default: thresh); period verification and the array walk keep thresh", default=None, type=float) + parser.add_argument('--skipped_bed', help="Write the spans of the records gate 16 skipped to this bed, " + "0-based half-open, name in|out:nfrac: the deleted span for an input skip and the " + "array for an output skip. Subtracting it from the high-confidence bed makes a " + "mostly-N record's real flank count as not assessed") parser.add_argument('-v', '--verbose', action='store_true', help="Verbose logging") args = parser.parse_args() main(args) diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index ed445e8..b41c547 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -26,7 +26,9 @@ from .logging import get_logger from .shell_pipeline import Command, Pipeline from .stages.alignment import ( + lane_path, BwaExtractStage, + PANGENOME_BWA_K, find_unzip, ) from .stages.base import StageContext @@ -50,6 +52,7 @@ parse_rg_line, path_arg, require_versions, + split_alignment, total_memory, ) from .shard import ( @@ -57,6 +60,13 @@ parse_fai, ) +# Peak memory, in GB, of vg haplotypes on the HPRC v2.0 and v2.1 graphs +# (about 60 with vg 1.78.0, up to 92 with a build that samples in parallel) +# and of one bwa-extract lane (about 20), measured on HG002 and HG005. +# vg-haplotypes runs alongside bwa only if the two fit together. +VG_HAPLOTYPES_MEM_GB = 100 +BWA_LANE_MEM_GB = 32 + SENT_PANGENOME_MIN_VERSIONS = { "kmc": None, "sentieon driver": packaging.version.Version("202503.02"), @@ -131,6 +141,10 @@ class SentieonPangenome(BasePangenome): "help": argparse.SUPPRESS, "action": "store_true", }, + "no_split_alignment": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, } ) @@ -149,6 +163,7 @@ def __init__(self) -> None: self.skip_multiqc = False self.skip_pop_vcf_id_check: bool = False self.skip_model_apply = False + self.no_split_alignment = False self.skip_small_variants = False def _cnv_in_second_dag(self) -> bool: @@ -334,8 +349,12 @@ def validate_bundle(self) -> None: else: self.extract_model_name = extract_candidate + # `dnascope.model` is only used for small-variant calling if ( - "dnascope.model" not in bundle_members + ( + "dnascope.model" not in bundle_members + and not self.skip_small_variants + ) or self.extract_model_name not in bundle_members or "minimap2.model" not in bundle_members ): @@ -345,6 +364,7 @@ def validate_bundle(self) -> None: sys.exit(2) self.has_cnv_model = "cnv.model" in bundle_members + self.has_combine_preset = "combine_cnv.json" in bundle_members if self.call_svs and not self.has_cnv_model: self.logger.warning( "The model bundle does not contain a 'cnv.model' file. " @@ -408,7 +428,10 @@ def collect_readgroups(self) -> None: def configure(self) -> None: """Configure pipeline parameters""" - pass + # bwa and minimap2 run as one lane per NUMA node (pieces of at most + # 64 CPUs, as in DNAscope) on large machines with fastq input + if self.r1_fastq and not self.no_split_alignment: + self.numa_nodes = split_alignment(self.cores) def build_dag(self) -> DAG: """Build the first DAG for the Sentieon pangenome pipeline""" @@ -463,17 +486,24 @@ def build_dag(self) -> DAG: dag.add_job(kmc_job) haplotype_dependencies.add(kmc_job) - # BWA alignment and extraction - bwa_result = self.bwa_extract_stage( - ctx, bwa_bam, ext_fastq - ).add_to(dag) - bwa_job = bwa_result.jobs[0] - mm2_dependencies.add(bwa_job) - bwa_lc_dependencies.add(bwa_job) - # Do not run vg-haplotypes with bwa in low-mem environments - if total_mem_gb < 70: - haplotype_dependencies.add(bwa_job) + # BWA alignment and extraction, one lane per NUMA node + bwa_stage = self.bwa_extract_stage(ctx, bwa_bam, ext_fastq) + lanes = bwa_stage.lane_paths() + bwa_result = bwa_stage.add_to(dag) + bwa_jobs = list(bwa_result.jobs) + mm2_dependencies.update(bwa_jobs) + bwa_lc_dependencies.update(bwa_jobs) + # Run vg-haplotypes alongside bwa when both fit in memory + vg_with_bwa = ( + total_mem_gb + >= VG_HAPLOTYPES_MEM_GB + BWA_LANE_MEM_GB * len(bwa_jobs) + ) + if not vg_with_bwa: + haplotype_dependencies.update(bwa_jobs) else: + lanes = [(bwa_bam, ext_fastq)] + bwa_jobs = [] + vg_with_bwa = False dnascope_bams = copy.deepcopy(self.sample_input) # ReadWriter cannot write to /dev/stdout directly; pre-create a # symlink so the driver writes to a real path that resolves to @@ -506,7 +536,14 @@ def build_dag(self) -> DAG: haplotype_dependencies.add(extract_kmc_job) # vg haplotypes - create a sample-specific pangenome - haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) + # Alongside bwa, vg-haplotypes reserves no threads, as kmc does. + # Otherwise the scheduler holds it until bwa releases its threads, + # whatever its dependencies say. + haplotypes_job = self.build_haplotypes_job( + sample_pangenome, + kmer_file, + job_threads=0 if vg_with_bwa else self.cores, + ) dag.add_job(haplotypes_job, haplotype_dependencies) # Confirm the sampled pangenome kept its reference paths before @@ -527,14 +564,24 @@ def build_dag(self) -> DAG: # minimap2 alignment of the extracted fastq dnascope_dependencies = set() - mm2_job = self.build_minimap2_lift_job( - mm2_bam, - ext_fastq, - sample_fasta, - sample_gfa, - ) - dag.add_job(mm2_job, mm2_dependencies | {check_gfa_job, fasta_job}) - dnascope_dependencies.add(mm2_job) + n_lanes = len(lanes) + mm2_jobs: List[Job] = [] + mm2_bams: List[pathlib.Path] = [] + for j, (_, lane_fastq) in enumerate(lanes): + lane_bam = mm2_bam if n_lanes == 1 else lane_path(mm2_bam, j) + lane_job = self.build_minimap2_lift_job( + lane_bam, + lane_fastq, + sample_fasta, + sample_gfa, + lane=j if n_lanes > 1 else None, + n_lanes=n_lanes, + ) + lane_deps = {bwa_jobs[j]} if n_lanes > 1 else mm2_dependencies + dag.add_job(lane_job, lane_deps | {check_gfa_job, fasta_job}) + dnascope_dependencies.add(lane_job) + mm2_jobs.append(lane_job) + mm2_bams.append(lane_bam) # With fastq input, perform dedup and metrics sr_alignments: List[pathlib.Path] = [] @@ -554,7 +601,7 @@ def build_dag(self) -> DAG: bwa_dedup = DedupStage( ctx=ctx, tag="bwa", - inputs=[bwa_bam], + inputs=[bam for bam, _ in lanes], output=out_bwa_aln, score_file=self.tmp_dir.joinpath("sample-bwa-score.txt.gz"), dedup_metrics=dedup_metrics, @@ -564,14 +611,14 @@ def build_dag(self) -> DAG: mm2_dedup = DedupStage( ctx=ctx, tag="mm2", - inputs=[mm2_bam], + inputs=mm2_bams, output=out_mm2_aln, score_file=self.tmp_dir.joinpath("sample-mm2-score.txt.gz"), read_filters=[ "IndelLeftAlignReadTransform,rgid=" f"{self.fastq_readgroup['ID']}-mm2" ], - ).add_to(dag, {mm2_job}) + ).add_to(dag, set(mm2_jobs)) dnascope_dependencies.add(mm2_dedup.dedup_job) sr_alignments = [out_bwa_aln] @@ -704,6 +751,7 @@ def build_dag(self) -> DAG: gvcf=out_gvcf, output=ctx.output_vcf, interval=self.bed, + dbsnp=self.dbsnp, ).add_to(dag, {ad_update_job}) # CNV calling is sex-aware and runs in the second DAG @@ -737,10 +785,15 @@ def bwa_extract_stage( extract_model=bundle.joinpath(self.extract_model_name), bwa_model=bundle.joinpath("bwa.model"), unzip=find_unzip(self.logger), + numa_nodes=self.numa_nodes, + bwa_k=PANGENOME_BWA_K, ) def build_haplotypes_job( - self, output_gbz: pathlib.Path, kmer_file: pathlib.Path + self, + output_gbz: pathlib.Path, + kmer_file: pathlib.Path, + job_threads: Optional[int] = None, ) -> Job: """Build vg haplotypes job""" hapl_file = self.required(self.hapl, "hapl") @@ -761,7 +814,7 @@ def build_haplotypes_job( ], ), "vg-haplotypes", - self.cores, + self.cores if job_threads is None else job_threads, task_name="pangenome", ) return haplotypes_job @@ -804,8 +857,12 @@ def build_minimap2_lift_job( ext_fastq: pathlib.Path, sample_fasta: pathlib.Path, sample_gfa: pathlib.Path, + lane: Optional[int] = None, + n_lanes: int = 1, ) -> Job: - """Build minimap2 alignment with pgutil lift job""" + """Build minimap2 alignment with pgutil lift job; with ``lane`` it is + that NUMA lane's share, pinned to its node""" + threads = max(1, int(self.cores / n_lanes)) bundle = self.required(self.model_bundle, "model_bundle") reference = self.required(self.reference, "reference") @@ -828,11 +885,13 @@ def build_minimap2_lift_job( reference, "@RG\\t" + "\\t".join([f"{x[0]}:{x[1]}" for x in rg2.items()]), mm2_model, - threads=self.cores, + threads=threads, lift_prefix=self.contig_prefix(), + numa=self.numa_nodes[lane] if lane is not None else None, ), - "mm2-lift", - self.cores, + "mm2-lift" if lane is None else f"mm2-lift-{lane}", + threads, + resources={f"node{lane}": 1} if lane is not None else None, task_name="pangenome-alignment", ) return mm2_job diff --git a/sentieon_cli/stages/alignment.py b/sentieon_cli/stages/alignment.py index 32736df..4656bbd 100644 --- a/sentieon_cli/stages/alignment.py +++ b/sentieon_cli/stages/alignment.py @@ -13,7 +13,7 @@ import logging import pathlib import shutil -from typing import Iterable, List, Optional, Union +from typing import Iterable, List, Optional, Tuple, Union import packaging.version @@ -377,6 +377,19 @@ def build(self) -> AlignResult: ) +# `sentieon fqidx extract` rejects -K below 10,000,000, so a split run needs +# at least this. dnascope-pangenome uses it whether or not it splits: each +# lane then forms the batches an unsplit run would, and results match +# across machines +PANGENOME_BWA_K = "10000000" + + +def lane_path(path: pathlib.Path, lane: int) -> pathlib.Path: + """`sample-bwa.fq.gz` -> `sample-bwa-.fq.gz`""" + stem, _, ext = path.name.partition(".") + return path.with_name(f"{stem}-{lane}.{ext}" if ext else f"{stem}-{lane}") + + @dataclass(kw_only=True) class BwaExtractStage(AlignmentStage): """Align short-read fastq with bwa and extract the reads a pangenome @@ -396,29 +409,54 @@ class BwaExtractStage(AlignmentStage): unzip: str = "gzip" name: str = "bwa-extract" task_name: str = "alignment" + numa_nodes: List[str] = field(default_factory=list) + # None keeps cmd_bwa_extract's default, or PANGENOME_BWA_K when split + bwa_k: Optional[str] = None + + def lane_paths(self) -> List[Tuple[pathlib.Path, pathlib.Path]]: + """(bam, fastq) per lane: one lane per NUMA node when there are + several, else the single output pair""" + if len(self.numa_nodes) <= 1: + return [(self.output_bam, self.output_fastq)] + return [ + (lane_path(self.output_bam, j), lane_path(self.output_fastq, j)) + for j in range(len(self.numa_nodes)) + ] def build(self) -> AlignResult: - job = Job( - cmds.cmd_bwa_extract( - self.output_bam, - self.output_fastq, - self.ctx.reference, - self.r1_fastq, - self.r2_fastq, - self.readgroup, - self.extract_model, - self.bwa_model, - self.ctx.cores, - unzip=self.unzip, - ), - self.name, - self.ctx.cores, - task_name=self.task_name, - ) + lanes = self.lane_paths() + n = len(lanes) + cores = max(1, int(self.ctx.cores / n)) + bwa_k = self.bwa_k or (PANGENOME_BWA_K if n > 1 else None) + jobs: List[Job] = [] + for j, (bam, fastq) in enumerate(lanes): + jobs.append( + Job( + cmds.cmd_bwa_extract( + bam, + fastq, + self.ctx.reference, + self.r1_fastq, + self.r2_fastq, + self.readgroup, + self.extract_model, + self.bwa_model, + cores, + unzip=self.unzip, + **({"bwa_k": bwa_k} if bwa_k else {}), + numa=self.numa_nodes[j] if n > 1 else None, + split=f"{j}/{n}" if n > 1 else None, + ), + f"{self.name}-{j}" if n > 1 else self.name, + cores, + resources={f"node{j}": 1} if n > 1 else None, + task_name=self.task_name, + ) + ) return AlignResult( - jobs=[job], - terminal={job}, - outputs=[self.output_bam, self.output_fastq], + jobs=jobs, + terminal=set(jobs), + outputs=[x for lane in lanes for x in lane], cleanup_paths=[], ) diff --git a/sentieon_cli/stages/cnv.py b/sentieon_cli/stages/cnv.py index d02e48e..f68e00e 100644 --- a/sentieon_cli/stages/cnv.py +++ b/sentieon_cli/stages/cnv.py @@ -13,11 +13,12 @@ from ..dag import DAG from ..driver import CNVModelApply, CNVscope from ..job import Job -from ..util import SampleSex, cnvscope_sex_args +from ..shard import parse_fai +from ..util import SampleSex, cnvscope_sex_args, unknown_sex_blocks from .base import Stage, StageResult, driver_job CNV_MIN_VERSIONS = { - # 202503.04 adds the CNVscope `--sex` and `--par` arguments + # 202503.04 adds the CNVscope `--sex`, `--par` and `--blocks` arguments "sentieon driver": packaging.version.Version("202503.04"), } @@ -29,6 +30,11 @@ class CNVResult(StageResult): cnvscope_job: Job apply_job: Job cnv_vcf: pathlib.Path + # The CNVscope `--sex` and `--par` arguments, None when omitted + sex: Optional[str] = None + par: Optional[pathlib.Path] = None + # The CNVscope `--blocks` BED file, used when the sex is not known + blocks: Optional[pathlib.Path] = None @dataclass(kw_only=True) @@ -38,7 +44,11 @@ class CNVscopeStage(Stage): The sample sex has to be known before the jobs are built, so callers run this stage after ploidy estimation. Version checks stay with the pipelines, which run them during validation. ``dump_probes`` is the - path CNVscope writes its per-probe debugging output to. + path CNVscope writes its per-probe debugging output to. When the sex + is not known, CNVscope is given a `--blocks` BED file of assumed sex + chromosome ploidies (`unknown_sex_blocks`) in place of `--sex/--par`. + CNVModelApply runs with `--emit_conf 0`, so `cnv_vcf` keeps the + segments the model rejects too, with FILTER `MLrejected`. """ inputs: List[pathlib.Path] @@ -57,6 +67,15 @@ class CNVscopeStage(Stage): def add_to(self, dag: DAG, upstream: Iterable[Job] = ()) -> CNVResult: deps = set(upstream) sex, par = cnvscope_sex_args(self.sample_sex, self.par_bed) + blocks: Optional[pathlib.Path] = None + if sex is None: + blocks = self.ctx.tmp_dir.joinpath(f"{self.name}-blocks.bed") + if not self.ctx.dry_run: + fai = parse_fai(pathlib.Path(f"{self.ctx.reference}.fai")) + lengths = {ctg: d["length"] for ctg, d in fai.items()} + with open(blocks, "w") as fh: + for block in unknown_sex_blocks(lengths, self.par_bed): + print(*block, sep="\t", file=fh) cnvscope = driver_job( self.ctx, @@ -66,6 +85,7 @@ def add_to(self, dag: DAG, upstream: Iterable[Job] = ()) -> CNVResult: self.model, sex=sex, par=par, + blocks=blocks, dump_probes=self.dump_probes, ) ], @@ -82,6 +102,8 @@ def add_to(self, dag: DAG, upstream: Iterable[Job] = ()) -> CNVResult: self.cnv_vcf, self.model, vcf=self.cnvscope_vcf, + # Downstream tools need the rejected segments too + emit_conf=0, ) ], name=self.apply_name, @@ -97,6 +119,9 @@ def add_to(self, dag: DAG, upstream: Iterable[Job] = ()) -> CNVResult: cnvscope_job=cnvscope, apply_job=apply_job, cnv_vcf=self.cnv_vcf, + sex=sex, + par=par, + blocks=blocks, ) @@ -120,6 +145,10 @@ class PangenomeCNVStage(Stage): sets are combined by `combine_cnv.py` into `_cnv.vcf.gz`. Both inputs -- the alignments and the SV VCF -- come from an earlier DAG, and the sample sex has to be known before the jobs are built. + `combine_cnv.py` gets the same sex and PAR BED file (or `--blocks` + BED file) as CNVscope. `combine_preset` is the JSON file of its rule + overrides, the model bundle's `combine_cnv.json`; without one, + `combine_cnv.py` uses its defaults, the paired-end values. """ inputs: List[pathlib.Path] @@ -127,6 +156,7 @@ class PangenomeCNVStage(Stage): sv_vcf: pathlib.Path sample_sex: Optional[SampleSex] = None par_bed: Optional[pathlib.Path] = None + combine_preset: Optional[pathlib.Path] = None interval: Optional[pathlib.Path] = None replace_rg: Optional[List[List[str]]] = None task_name: str = "cnv" @@ -178,7 +208,8 @@ def add_to( ) dag.add_job(indel2cnv_job, deps) - # Combine the CNVModelApply output with the converted SVs + # Combine the CNVModelApply output, rejected segments included, + # with the converted SVs combine_script = pathlib.Path( str(files("sentieon_cli.scripts").joinpath("combine_cnv.py")) ) @@ -188,6 +219,10 @@ def add_to( apply_vcf, converted_vcf, combine_script, + preset=self.combine_preset, + sex=cnv_result.sex, + par=cnv_result.par, + blocks=cnv_result.blocks, ), "combine-cnv", 0, diff --git a/sentieon_cli/stages/small_variants.py b/sentieon_cli/stages/small_variants.py index 4880b55..ac20fd3 100644 --- a/sentieon_cli/stages/small_variants.py +++ b/sentieon_cli/stages/small_variants.py @@ -120,6 +120,7 @@ class GVCFtyperStage(Stage): gvcf: pathlib.Path output: pathlib.Path interval: Optional[pathlib.Path] = None + dbsnp: Optional[pathlib.Path] = None name: str = "gvcftyper" task_name: str = "gvcftyper" @@ -128,7 +129,7 @@ def add_to( ) -> GVCFtyperResult: job = driver_job( self.ctx, - [GVCFtyper(output=self.output, vcf=self.gvcf)], + [GVCFtyper(output=self.output, vcf=self.gvcf, dbsnp=self.dbsnp)], interval=self.interval, name=self.name, task_name=self.task_name, diff --git a/sentieon_cli/stages/transfer.py b/sentieon_cli/stages/transfer.py index b6b9fff..480df46 100644 --- a/sentieon_cli/stages/transfer.py +++ b/sentieon_cli/stages/transfer.py @@ -94,10 +94,12 @@ def build_transfer_jobs( print(f"{shard.contig}\t0\t{ctg_len}", file=fh) view_job = Job( - cmds.cmd_bcftools_view_regions( + cmds.cmd_bcftools_subset_remove( subset_vcf, raw_vcf, regions_file=subset_bed, + # The merged shards drop FILTER too + remove="FILTER", ), f"{name_base}-extra", 1, @@ -124,14 +126,9 @@ def build_transfer_jobs( pop_vcf, trim_script, subset_bed, + shard.contig, merge_rules=merge_rules, - merge_xargs=[ - "--no-version", - "--regions-overlap", - "pos", - "-m", - "all", - ], + merge_xargs=["--no-version", "-m", "all"], view_xargs=["--no-version"], ), f"{name_base}-{i}", diff --git a/sentieon_cli/util.py b/sentieon_cli/util.py index 7d79905..d55eb38 100644 --- a/sentieon_cli/util.py +++ b/sentieon_cli/util.py @@ -73,21 +73,61 @@ def cnvscope_sex_args( Male samples are called with the pseudo-autosomal regions (PAR) BED file. Female samples do not need a PAR BED file. When the sample sex - is not known, both arguments are omitted and CNVscope assumes a - diploid genome, matching the behavior of previous releases. + is not known, both arguments are omitted and the caller supplies + assumed sex chromosome ploidies instead (`unknown_sex_blocks`). """ if sample_sex is SampleSex.MALE: return ("M", par_bed) if sample_sex is SampleSex.FEMALE: return ("F", None) logger.warning( - "The sample sex is not known. CNVscope will assume a diploid " - "genome, matching the behavior of previous releases. Supply " + "Sex estimation failed: the sample sex could not be determined " + "from read coverage. CNV calling is proceeding with assumed sex " + "chromosome ploidies: two copies of chrX, and one copy of chrY " + "outside the pseudo-autosomal regions (PAR). Supply " "`--sample_sex` for sex-aware CNV calling." ) return (None, None) +def unknown_sex_blocks( + contig_lengths: Mapping[str, int], + par_bed: Optional[pathlib.Path], +) -> List[Tuple[str, int, int, int]]: + """CNVscope `--blocks` intervals for a sample of unknown sex. + + Each block is (contig, start, end, neutral copy number), 0-based and + half-open. chrX is diploid and chrY is haploid outside the PAR + intervals of `par_bed`. Regions without a block, including the chrY + PAR, keep CNVscope's diploid default. Contigs are named with or + without the `chr` prefix, following the reference. + """ + blocks: List[Tuple[str, int, int, int]] = [] + x_contig = next((c for c in ("chrX", "X") if c in contig_lengths), None) + if x_contig: + blocks.append((x_contig, 0, contig_lengths[x_contig], 2)) + y_contig = next((c for c in ("chrY", "Y") if c in contig_lengths), None) + if not y_contig: + return blocks + + par: List[Tuple[int, int]] = [] + if par_bed: + with open(par_bed) as fh: + for line in fh: + flds = line.split() + if len(flds) >= 3 and flds[0] == y_contig: + par.append((int(flds[1]), int(flds[2]))) + pos = 0 + for start, end in sorted(par): + if start > pos: + blocks.append((y_contig, pos, start, 1)) + pos = max(pos, end) + y_len = contig_lengths[y_contig] + if pos < y_len: + blocks.append((y_contig, pos, y_len, 1)) + return blocks + + def caller_sex_arg(sample_sex: Optional[SampleSex]) -> str: """The `--sex male|female` value for ExpansionHunter and segdup-caller. diff --git a/tests/unit/test_cnv_sex.py b/tests/unit/test_cnv_sex.py index 3fc4f1d..ab8126f 100644 --- a/tests/unit/test_cnv_sex.py +++ b/tests/unit/test_cnv_sex.py @@ -28,6 +28,7 @@ caller_sex_arg, cnvscope_sex_args, sample_sex_arg, + unknown_sex_blocks, ) PACKAGE_LOGGER = "sentieon_cli" @@ -208,7 +209,61 @@ def test_cnvscope_sex_args_for_a_female_sample(): @pytest.mark.parametrize("sex", [None, SampleSex.UNKNOWN]) def test_cnvscope_sex_args_without_a_known_sex(sex, messages): assert cnvscope_sex_args(sex, None) == (None, None) - assert any("diploid" in msg for msg in messages) + assert any( + "Sex estimation failed" in msg + and "assumed sex chromosome ploidies" in msg + for msg in messages + ) + + +# The CNVscope `--blocks` intervals for a sample of unknown sex + +HG38_LENGTHS = {"chr20": 64444167, "chrX": 156040895, "chrY": 57227415} + + +def test_unknown_sex_blocks_haploid_chry_outside_the_par(tmp_path): + par_bed = tmp_path / "par.bed" + par_bed.write_text( + "chrX\t10000\t2781479\n" + "chrX\t155701382\t156030895\n" + "chrY\t56887902\t57217415\n" + "chrY\t10000\t2781479\n" + ) + + assert unknown_sex_blocks(HG38_LENGTHS, par_bed) == [ + ("chrX", 0, 156040895, 2), + ("chrY", 0, 10000, 1), + ("chrY", 2781479, 56887902, 1), + ("chrY", 57217415, 57227415, 1), + ] + + +def test_unknown_sex_blocks_without_a_par_bed(): + assert unknown_sex_blocks(HG38_LENGTHS, None) == [ + ("chrX", 0, 156040895, 2), + ("chrY", 0, 57227415, 1), + ] + + +def test_unknown_sex_blocks_with_b37_names_and_overlapping_par(tmp_path): + par_bed = tmp_path / "par.bed" + par_bed.write_text( + "X\t60000\t2699520\tPAR1\n" + "Y\t10000\t2649520\tPAR1\n" + "Y\t2000000\t2700000\tPAR1b\n" + "Y\t59034049\t59373566\tPAR2\n" + ) + lengths = {"20": 63025520, "X": 155270560, "Y": 59373566} + + assert unknown_sex_blocks(lengths, par_bed) == [ + ("X", 0, 155270560, 2), + ("Y", 0, 10000, 1), + ("Y", 2700000, 59034049, 1), + ] + + +def test_unknown_sex_blocks_without_sex_chromosomes(): + assert unknown_sex_blocks({"chr20": 64444167}, None) == [] # The ExpansionHunter and segdup-caller `--sex` argument @@ -253,6 +308,7 @@ def test_cnvscope_algo_omits_the_sex_arguments_by_default(): assert "--sex" not in cmd assert "--par" not in cmd + assert "--blocks" not in cmd assert "--dump_probes" not in cmd assert cmd[-1] == "out.vcf.gz" @@ -412,3 +468,4 @@ def test_call_cnvs_without_a_sample_sex(): cmd = _cnvscope_cmd() assert "--sex" not in cmd assert "--par" not in cmd + assert "--blocks /tmp/cnv/cnvscope-blocks.bed" in cmd diff --git a/tests/unit/test_dag_construction.py b/tests/unit/test_dag_construction.py index b7faa77..6320a4f 100644 --- a/tests/unit/test_dag_construction.py +++ b/tests/unit/test_dag_construction.py @@ -47,9 +47,12 @@ def setup_method(self): self.mock_bundle = pathlib.Path(self.temp_dir) / "model.bundle" # Create empty files - for file_path in [self.mock_ref, self.mock_bam, self.mock_bundle]: + for file_path in [self.mock_ref, self.mock_bam]: file_path.touch() + # An empty ar archive: a bundle without a bundle_info.json + self.mock_bundle.write_bytes(b"!\n") + # Create reference index (pathlib.Path(str(self.mock_ref) + ".fai")).touch() diff --git a/tests/unit/test_dnascope_bundle.py b/tests/unit/test_dnascope_bundle.py new file mode 100644 index 0000000..bd5eb16 --- /dev/null +++ b/tests/unit/test_dnascope_bundle.py @@ -0,0 +1,329 @@ +""" +Unit tests for the DNAscope pipeline's model-bundle handling: the +`--pop_vcf` annotation transfer and the platform's read filter +""" + +import json +import pathlib +from typing import Dict, List, Optional +from unittest.mock import MagicMock + +import pytest + +from sentieon_cli.dnascope import DNAscopePipeline + +BUNDLE_VCF_ID = "pop-vcf-v1" +CHR1_LENGTH = 25 * 1000 * 1000 + + +def write_bundle(path: pathlib.Path, bundle_info: Optional[Dict]) -> None: + """Write an ar archive holding at most a bundle_info.json member""" + data = b"!\n" + if bundle_info is not None: + body = json.dumps(bundle_info).encode() + name = b"bundle_info.json" # exactly 16 bytes, so no "/" suffix + header = ( + name.ljust(16) + + b"0".ljust(12) + + b"0".ljust(6) + + b"0".ljust(6) + + b"644".ljust(8) + + str(len(body)).encode().ljust(10) + + b"`\n" + ) + data += header + body + (b"\n" if len(body) % 2 else b"") + path.write_bytes(data) + + +def write_pop_vcf( + path: pathlib.Path, + vcf_id: Optional[str], + contigs: Dict[str, int], +) -> None: + """Write a header-only population VCF""" + lines = ["##fileformat=VCFv4.2"] + if vcf_id is not None: + lines.append(f"##SentieonVcfID={vcf_id}") + for ctg, length in contigs.items(): + lines.append(f"##contig=") + lines.append("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO") + path.write_text("\n".join(lines) + "\n") + + +@pytest.fixture +def files(tmp_path): + """The reference, alignment, bundle and pop VCF of a run""" + ref = tmp_path / "reference.fa" + ref.touch() + pathlib.Path(str(ref) + ".fai").write_text( + f"chr1\t{CHR1_LENGTH}\t6\t60\t61\nchrM\t16569\t25416684\t60\t61\n" + ) + bam = tmp_path / "sample.bam" + bam.touch() + bundle = tmp_path / "model.bundle" + write_bundle(bundle, {"SentieonVcfID": BUNDLE_VCF_ID}) + pop_vcf = tmp_path / "pop.vcf.gz" + write_pop_vcf(pop_vcf, BUNDLE_VCF_ID, {"chr1": CHR1_LENGTH}) + tmp_dir = tmp_path / "tmp" + tmp_dir.mkdir() + return { + "ref": ref, + "bam": bam, + "bundle": bundle, + "pop_vcf": pop_vcf, + "tmp_dir": tmp_dir, + "output_vcf": tmp_path / "output.vcf.gz", + } + + +def make_pipeline(files, **overrides) -> DNAscopePipeline: + """A dry-run DNAscope pipeline over a single BAM input""" + pipeline = DNAscopePipeline() + pipeline.logger = MagicMock() + pipeline.output_vcf = files["output_vcf"] + pipeline.reference = files["ref"] + pipeline.model_bundle = files["bundle"] + pipeline.sample_input = [files["bam"]] + pipeline.pop_vcf = files["pop_vcf"] + pipeline.cores = 2 + pipeline.dry_run = True + pipeline.skip_version_check = True + pipeline.tmp_dir = files["tmp_dir"] + for key, value in overrides.items(): + setattr(pipeline, key, value) + return pipeline + + +def build(pipeline: DNAscopePipeline): + """Validate, configure and build the DAG""" + pipeline.validate() + pipeline.configure() + dag = pipeline.build_dag() + return dag + + +def all_jobs(dag) -> List: + return list(dag.waiting_jobs) + list(dag.ready_jobs) + + +def job_named(dag, name: str): + matches = [job for job in all_jobs(dag) if job.name == name] + assert len(matches) == 1, (name, [job.name for job in all_jobs(dag)]) + return matches[0] + + +def dep_names(dag, job) -> List[str]: + return sorted(dep.name for dep in dag.waiting_jobs.get(job, set())) + + +def exit_with_error(pipeline: DNAscopePipeline, message: str) -> None: + """`validate()` exits 2 and logs an error containing `message`""" + with pytest.raises(SystemExit) as excinfo: + pipeline.validate() + assert excinfo.value.code == 2 + logged = [str(call.args[0]) for call in pipeline.logger.error.call_args_list] + assert any(message in line for line in logged), logged + + +class TestPopVcfValidation: + """The `--pop_vcf` is checked against the model bundle""" + + def test_matching_pop_vcf_is_accepted(self, files): + pipeline = make_pipeline(files) + pipeline.validate() + pipeline.logger.error.assert_not_called() + + def test_bundle_requiring_a_pop_vcf_rejects_a_run_without_one( + self, files + ): + pipeline = make_pipeline(files, pop_vcf=None) + exit_with_error(pipeline, "requires a population VCF") + + def test_mismatched_pop_vcf_id_is_rejected(self, files): + write_pop_vcf(files["pop_vcf"], "other-id", {"chr1": CHR1_LENGTH}) + pipeline = make_pipeline(files) + exit_with_error(pipeline, "does not match the population VCF") + + def test_pop_vcf_without_an_id_is_rejected(self, files): + write_pop_vcf(files["pop_vcf"], None, {"chr1": CHR1_LENGTH}) + pipeline = make_pipeline(files) + exit_with_error(pipeline, "does not match the population VCF") + + def test_skip_pop_vcf_id_check_accepts_a_mismatch(self, files): + write_pop_vcf(files["pop_vcf"], "other-id", {"chr1": CHR1_LENGTH}) + pipeline = make_pipeline(files, skip_pop_vcf_id_check=True) + pipeline.validate() + pipeline.logger.error.assert_not_called() + + @pytest.mark.parametrize("bundle_info", [None, {"minScriptVersion": "1"}]) + def test_bundle_without_an_id_rejects_a_pop_vcf(self, files, bundle_info): + write_bundle(files["bundle"], bundle_info) + pipeline = make_pipeline(files) + exit_with_error(pipeline, "does not require a population VCF") + + @pytest.mark.parametrize("bundle_info", [None, {"minScriptVersion": "1"}]) + def test_bundle_without_an_id_runs_without_a_pop_vcf( + self, files, bundle_info + ): + write_bundle(files["bundle"], bundle_info) + pipeline = make_pipeline(files, pop_vcf=None) + pipeline.validate() + pipeline.logger.error.assert_not_called() + + def test_transfer_inputs_are_collected(self, files): + pipeline = make_pipeline(files) + pipeline.validate() + assert pipeline.pop_vcf_contigs == {"chr1": CHR1_LENGTH} + assert set(pipeline.fai_data) == {"chr1", "chrM"} + # 10 Mb shards: three for chr1, one for chrM + assert [s.contig for s in pipeline.shards] == ["chr1"] * 3 + ["chrM"] + + def test_no_reference_index_is_read_without_a_pop_vcf(self, files): + write_bundle(files["bundle"], None) + pathlib.Path(str(files["ref"]) + ".fai").unlink() + pipeline = make_pipeline(files, pop_vcf=None) + pipeline.validate() + assert pipeline.fai_data == {} + assert pipeline.shards == [] + + +class TestPopVcfDag: + """Annotations are transferred before DNAModelApply""" + + def test_transfer_runs_between_dnascope_and_model_apply(self, files): + dag = build(make_pipeline(files)) + + names = [job.name for job in all_jobs(dag)] + assert sorted(n for n in names if n.startswith("merge-trim")) == [ + "merge-trim-0", + "merge-trim-1", + "merge-trim-2", + "merge-trim-concat", + # chrM is not in the pop VCF, so it is copied, not merged + "merge-trim-extra", + ] + for name in ("merge-trim-0", "merge-trim-extra"): + assert dep_names(dag, job_named(dag, name)) == ["dnascope"] + + apply_job = job_named(dag, "model-apply") + assert dep_names(dag, apply_job) == ["merge-trim-concat"] + transfer_vcf = files["tmp_dir"] / "sample-dnascope_transfer.vcf.gz" + assert str(transfer_vcf) in str(job_named(dag, "merge-trim-concat").shell) + assert str(transfer_vcf) in str(apply_job.shell) + assert str(files["output_vcf"]) in str(apply_job.shell) + + assert dep_names(dag, job_named(dag, "rm-tmp-vcf")) == ["model-apply"] + + def test_no_transfer_without_a_pop_vcf(self, files): + write_bundle(files["bundle"], None) + dag = build(make_pipeline(files, pop_vcf=None)) + + names = [job.name for job in all_jobs(dag)] + assert not any(n.startswith("merge-trim") for n in names) + apply_job = job_named(dag, "model-apply") + assert dep_names(dag, apply_job) == ["dnascope"] + + def test_gvcf_transfer(self, files): + dag = build(make_pipeline(files, gvcf=True)) + + transfer_gvcf = files["tmp_dir"] / "sample-dnascope_transfer.g.vcf.gz" + apply_job = job_named(dag, "model-apply") + assert str(transfer_gvcf) in str(apply_job.shell) + out_gvcf = str(files["output_vcf"]).replace(".vcf.gz", ".g.vcf.gz") + assert out_gvcf in str(apply_job.shell) + assert dep_names(dag, job_named(dag, "gvcftyper")) == ["model-apply"] + + def test_gvcftyper_annotates_dbsnp(self, files): + dbsnp = files["tmp_dir"] / "dbsnp.vcf.gz" + dag = build(make_pipeline(files, gvcf=True, dbsnp=dbsnp)) + assert f"--dbsnp {dbsnp}" in str(job_named(dag, "gvcftyper").shell) + + dag = build(make_pipeline(files, gvcf=True, dbsnp=None)) + assert "--dbsnp" not in str(job_named(dag, "gvcftyper").shell) + + def test_skip_small_variants_has_no_transfer(self, files): + dag = build(make_pipeline(files, skip_small_variants=True)) + names = [job.name for job in all_jobs(dag)] + assert not any(n.startswith("merge-trim") for n in names) + + +class TestUltimaReadFilter: + """An Ultima bundle adds `UltimaReadFilter` to the DNAscope call""" + + @staticmethod + def build_with_bundle(files, bundle_info, **overrides): + write_bundle(files["bundle"], bundle_info) + return build(make_pipeline(files, pop_vcf=None, **overrides)) + + @pytest.mark.parametrize("platform", ["Ultima", "ULTIMA", "ultima"]) + def test_ultima_bundle_filters_the_dnascope_call(self, files, platform): + dag = self.build_with_bundle(files, {"platform": platform}) + + cmd = str(job_named(dag, "dnascope").shell) + assert cmd.count("--read_filter UltimaReadFilter") == 1 + # Duplicate marking and the rest of the pipeline are unfiltered + others = [j for j in all_jobs(dag) if j.name != "dnascope"] + assert others + for job in others: + assert "UltimaReadFilter" not in str(job.shell), job.name + + def test_ultima_filter_covers_sv_calling(self, files): + dag = self.build_with_bundle(files, {"platform": "Ultima"}) + + # Small variants and SVs come from the same driver call + cmd = str(job_named(dag, "dnascope").shell) + assert "--var_type BND" in cmd + assert "--read_filter UltimaReadFilter" in cmd + + def test_ultima_bundle_with_a_pop_vcf(self, files): + write_bundle( + files["bundle"], + {"platform": "Ultima", "SentieonVcfID": BUNDLE_VCF_ID}, + ) + dag = build(make_pipeline(files)) + + cmd = str(job_named(dag, "dnascope").shell) + assert "--read_filter UltimaReadFilter" in cmd + assert dep_names(dag, job_named(dag, "model-apply")) == [ + "merge-trim-concat" + ] + + @pytest.mark.parametrize("platform", ["Ultima", "ultima"]) + def test_short_read_platform_key(self, files, platform): + # SentieonUltima1.1.bundle names the platform `shortReadPlatform` + dag = self.build_with_bundle(files, {"shortReadPlatform": platform}) + + cmd = str(job_named(dag, "dnascope").shell) + assert cmd.count("--read_filter UltimaReadFilter") == 1 + + def test_platform_key_takes_precedence(self, files): + dag = self.build_with_bundle( + files, {"platform": "Illumina", "shortReadPlatform": "Ultima"} + ) + assert "UltimaReadFilter" not in str(job_named(dag, "dnascope").shell) + + @pytest.mark.parametrize( + "bundle_info", + [ + None, + {}, + {"platform": "Illumina"}, + {"shortReadPlatform": "Illumina"}, + ], + ) + def test_other_bundles_are_unfiltered(self, files, bundle_info): + dag = self.build_with_bundle(files, bundle_info) + + for job in all_jobs(dag): + assert "--read_filter" not in str(job.shell), job.name + + def test_platform_is_read_from_the_bundle(self, files): + write_bundle(files["bundle"], {"platform": "Ultima"}) + pipeline = make_pipeline(files, pop_vcf=None) + pipeline.validate() + assert pipeline.tech == "ULTIMA" + + write_bundle(files["bundle"], None) + pipeline = make_pipeline(files, pop_vcf=None) + pipeline.validate() + assert pipeline.tech == "" diff --git a/tests/unit/test_dnascope_hybrid_cnv.py b/tests/unit/test_dnascope_hybrid_cnv.py index 99287e4..34f6ef3 100644 --- a/tests/unit/test_dnascope_hybrid_cnv.py +++ b/tests/unit/test_dnascope_hybrid_cnv.py @@ -183,6 +183,22 @@ def test_cnvscope_dumps_the_probes_next_to_the_cnv_vcf(self): assert f"--dump_probes {probes}" in str(jobs["cnvscope"].shell) assert "--dump_probes" not in str(jobs["cnv-model-apply"].shell) + def test_cnv_model_apply_writes_every_segment(self): + # The CNV VCF keeps the rejected segments too, as MLrejected + pipeline = self.create_pipeline() + pipeline.sample_sex = SampleSex.FEMALE + self.build_dag(pipeline) + _dag, jobs = self.build_second_dag(pipeline) + + cnv_vcf = self.mock_dir / "output.cnv.vcf.gz" + assert str(jobs["cnv-model-apply"].shell) == ( + f"sentieon driver --reference {self.mock_ref} " + "--thread_count 2 --algo CNVModelApply " + f"--model {self.mock_bundle}/cnv.model " + f"--vcf {self.mock_dir}/cnvscope.vcf.gz " + f"--emit_conf 0 {cnv_vcf}" + ) + def test_a_dry_run_assumes_a_male_sample(self): # Validation has already confirmed a PAR BED file is available pipeline = self.create_pipeline() diff --git a/tests/unit/test_dnascope_hybrid_gvcf.py b/tests/unit/test_dnascope_hybrid_gvcf.py new file mode 100644 index 0000000..867685b --- /dev/null +++ b/tests/unit/test_dnascope_hybrid_gvcf.py @@ -0,0 +1,318 @@ +""" +Unit tests for the DNAscope hybrid pipeline's `--gvcf` output +""" + +from unittest.mock import MagicMock, patch + +import pytest + +from sentieon_cli.dnascope_hybrid import DNAscopeHybridPipeline + + +@pytest.fixture +def pipeline(tmp_path): + """A DNAscopeHybridPipeline ready for `build_dag`""" + for name in [ + "reference.fa", + "short.bam", + "long.bam", + "population.vcf.gz", + "interval.bed", + ]: + (tmp_path / name).touch() + (tmp_path / "model.bundle").touch() + + with patch("sys.exit"): + pipeline = DNAscopeHybridPipeline() + pipeline.logger = MagicMock() + + pipeline.output_vcf = tmp_path / "output.vcf.gz" + pipeline.reference = tmp_path / "reference.fa" + pipeline.model_bundle = tmp_path / "model.bundle" + pipeline.bed = tmp_path / "interval.bed" + pipeline.cores = 2 + pipeline.dry_run = True + pipeline.skip_version_check = True + pipeline.tmp_dir = tmp_path + pipeline.pop_vcf = None + pipeline.gvcf = True + + pipeline.fai_data = {"chr1": {"length": 1000}} + shard = MagicMock() + shard.contig = "chr1" + shard.start = 1 + shard.stop = 1000 + pipeline.shards = [shard] + pipeline.pop_vcf_contigs = {"chr1": 1000} + + pipeline.lr_aln_readgroups = [[{"ID": "lr_rg1", "SM": "sample1"}]] + pipeline.sr_aln_readgroups = [[{"ID": "sr_rg1", "SM": "sample1"}]] + pipeline.hybrid_rg_sm = "sample1" + pipeline.hybrid_set_rg = False + pipeline.shortread_tech = "Illumina" + pipeline.longread_tech = "ONT" + pipeline.sr_aln = [tmp_path / "short.bam"] + pipeline.lr_aln = [tmp_path / "long.bam"] + return pipeline + + +def _jobs(pipeline): + """Build the DAG, returning {job name: (shell, dependency names)}""" + dag = pipeline.build_dag() + jobs = {} + for job, deps in dag.waiting_jobs.items(): + assert job.name not in jobs, f"duplicate job name {job.name}" + jobs[job.name] = (str(job.shell), {d.name for d in deps}) + return jobs + + +def _out(pipeline, suffix): + """A path next to the output VCF""" + return str(pipeline.output_vcf).replace(".vcf.gz", suffix) + + +class TestHybridGvcf: + def test_calling_passes_emit_gvcfs(self, pipeline): + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + call2, _ = jobs["dnascope-2"] + assert "--emit_mode gvcf" in call2 + assert call2.endswith(str(tmp / "hybrid_pass2.g.vcf.gz")) + + def test_pass1_also_writes_a_vcf_for_region_selection(self, pipeline): + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + call1, _ = jobs["dnascope-1"] + # One driver run, two DNAscope algos: the gVCF, then the VCF + algos = call1.split(" --algo ")[1:] + assert len(algos) == 2 + assert algos[0].startswith("DNAscope ") + assert "--emit_mode gvcf" in algos[0] + assert algos[0].endswith(str(tmp / "initial.g.vcf.gz")) + assert algos[1].startswith("DNAscope ") + assert "--emit_mode" not in algos[1] + assert algos[1].endswith(str(tmp / "initial.vcf.gz")) + # Both algos otherwise take the same arguments + assert algos[0].replace(" --emit_mode gvcf", "").rsplit(" ", 1)[0] == ( + algos[1].rsplit(" ", 1)[0] + ) + + select, select_deps = jobs["hybrid-select"] + assert f"-v {tmp / 'initial.vcf.gz'}" in select + assert "initial.g.vcf.gz" not in select + assert select_deps == {"dnascope-1"} + + cleanup, _ = jobs["rm-tmp5"] + assert str(tmp / "initial.vcf.gz") in cleanup + assert str(tmp / "initial.g.vcf.gz") in cleanup + + def test_both_passes_are_clipped_before_the_concat(self, pipeline): + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + recall_bed = str(tmp / "hybrid_recall.bed") + + pass1, pass1_deps = jobs["subset-calls"] + assert "gvcf_clip.py" in pass1 + assert f"-v {tmp / 'initial.g.vcf.gz'}" in pass1 + assert f"-b {recall_bed}" in pass1 + assert f"-r {pipeline.reference}" in pass1 + assert "--exclude - | sentieon util vcfconvert - " in pass1 + assert pass1.endswith(str(tmp / "mix_subset.g.vcf.gz")) + assert pass1_deps == {"intersect-stage2-bed"} + + pass2, pass2_deps = jobs["subset-pass2-calls"] + assert "gvcf_clip.py" in pass2 + assert f"-v {tmp / 'hybrid_pass2.g.vcf.gz'}" in pass2 + assert f"-b {recall_bed}" in pass2 + assert "--exclude" not in pass2 + assert " - | sentieon util vcfconvert - " in pass2 + assert pass2.endswith(str(tmp / "hybrid_pass2_subset.g.vcf.gz")) + assert pass2_deps == {"dnascope-2"} + + concat, concat_deps = jobs["concat-calls"] + assert "-aD --threads 2 " in concat + assert concat.endswith( + f"{tmp / 'mix_subset.g.vcf.gz'} " + f"{tmp / 'hybrid_pass2_subset.g.vcf.gz'}" + ) + assert concat_deps == {"subset-calls", "subset-pass2-calls"} + + cleanup, _ = jobs["rm-tmp5"] + assert str(tmp / "hybrid_pass2_subset.g.vcf.gz") in cleanup + + def test_annotation_streams_through_vcfconvert(self, pipeline): + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + anno, anno_deps = jobs["anno-calls"] + assert f"-v {tmp / 'combined_tmp.g.vcf.gz'}" in anno + assert anno.endswith( + f" - | sentieon util vcfconvert - " + f"{tmp / 'combined_tmp_anno.g.vcf.gz'}" + ) + assert anno_deps == {"concat-calls"} + + def test_model_apply_writes_the_output_gvcf(self, pipeline): + jobs = _jobs(pipeline) + apply, apply_deps = jobs["model-apply"] + assert f"--vcf {pipeline.tmp_dir / 'combined_tmp_anno.g.vcf.gz'}" in ( + apply + ) + assert apply.endswith(_out(pipeline, ".g.vcf.gz")) + assert apply_deps == {"anno-calls"} + + def test_gvcftyper_then_field_removal_then_final_norm(self, pipeline): + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + gvcftyper_vcf = str(tmp / "combined_gvcftyper.vcf.gz") + typed_vcf = str(tmp / "combined_typed.vcf.gz") + + typer, typer_deps = jobs["gvcftyper"] + assert f"--interval {pipeline.bed}" in typer + assert f"--vcf {_out(pipeline, '.g.vcf.gz')}" in typer + assert "--dbsnp" not in typer + assert typer.endswith(gvcftyper_vcf) + assert typer_deps == {"model-apply"} + + drop, drop_deps = jobs["drop-allele-fields"] + assert drop == ( + "bcftools annotate -x FORMAT/LAD,FORMAT/LPL,FORMAT/SAD,FORMAT/SPL " + f"-Oz -W=tbi -o {typed_vcf} {gvcftyper_vcf}" + ) + assert drop_deps == {"gvcftyper"} + + norm, norm_deps = jobs["final-norm"] + assert f"""-e 'GT="0/0"' {typed_vcf}""" in norm + assert norm.endswith(str(pipeline.output_vcf)) + assert norm_deps == {"drop-allele-fields"} + + def test_gvcftyper_annotates_dbsnp(self, pipeline): + pipeline.dbsnp = pipeline.tmp_dir / "dbsnp.vcf.gz" + jobs = _jobs(pipeline) + typer, _ = jobs["gvcftyper"] + assert f"--dbsnp {pipeline.dbsnp}" in typer + + def test_pop_vcf_transfer_runs_on_the_gvcf(self, pipeline): + pipeline.pop_vcf = pipeline.tmp_dir / "population.vcf.gz" + jobs = _jobs(pipeline) + transfer_vcf = str(pipeline.tmp_dir / "combined_tmp_transfer.g.vcf.gz") + + concat, _ = jobs["merge-trim-concat"] + assert f"--output {transfer_vcf}" in concat + apply, apply_deps = jobs["model-apply"] + assert f"--vcf {transfer_vcf}" in apply + assert apply.endswith(_out(pipeline, ".g.vcf.gz")) + assert apply_deps == {"merge-trim-concat"} + assert jobs["gvcftyper"][1] == {"model-apply"} + + def test_skip_model_apply(self, pipeline): + pipeline.skip_model_apply = True + jobs = _jobs(pipeline) + gvcftyper_vcf = str(pipeline.tmp_dir / "combined_gvcftyper.vcf.gz") + + anno, _ = jobs["anno-calls"] + assert anno.endswith(_out(pipeline, ".g.vcf.gz")) + typer, typer_deps = jobs["gvcftyper"] + assert f"--vcf {_out(pipeline, '.g.vcf.gz')}" in typer + assert typer.endswith(gvcftyper_vcf) + assert typer_deps == {"anno-calls"} + drop, drop_deps = jobs["drop-allele-fields"] + assert drop.endswith(f"-o {pipeline.output_vcf} {gvcftyper_vcf}") + assert drop_deps == {"gvcftyper"} + assert "model-apply" not in jobs + assert "final-norm" not in jobs + + def test_skip_model_apply_with_pop_vcf(self, pipeline): + pipeline.skip_model_apply = True + pipeline.pop_vcf = pipeline.tmp_dir / "population.vcf.gz" + jobs = _jobs(pipeline) + + concat, _ = jobs["merge-trim-concat"] + assert f"--output {_out(pipeline, '.g.vcf.gz')}" in concat + assert jobs["gvcftyper"][1] == {"merge-trim-concat"} + drop, _ = jobs["drop-allele-fields"] + assert f"-o {pipeline.output_vcf} " in drop + assert "final-norm" not in jobs + + +@pytest.mark.parametrize("gvcf", [True, False]) +class TestRecallBed: + """The second pass and the merge stay inside --bed""" + + def test_stage2_bed_is_intersected_with_the_bed(self, pipeline, gvcf): + pipeline.gvcf = gvcf + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + stage2_bed = str(tmp / "hybrid_stage2.bed") + recall_bed = str(tmp / "hybrid_recall.bed") + + intersect, intersect_deps = jobs["intersect-stage2-bed"] + assert ( + f"bedtools intersect -a {stage2_bed} -b {pipeline.bed}" + in intersect + ) + assert "bedtools merge" in intersect + assert recall_bed in intersect.rsplit(">", 1)[1] + assert intersect_deps == {"second-stage"} + + call2, call2_deps = jobs["dnascope-2"] + assert f"--interval {recall_bed}" in call2 + assert call2_deps == {"third-stage", "intersect-stage2-bed"} + assert recall_bed in jobs["subset-calls"][0] + assert jobs["subset-calls"][1] == {"intersect-stage2-bed"} + if gvcf: + assert recall_bed in jobs["subset-pass2-calls"][0] + + # Stage 3 keeps the whole stage-2 BED + stage3, _ = jobs["third-stage"] + assert f"--interval {stage2_bed}" in stage3 + assert recall_bed not in stage3 + + def test_without_a_bed(self, pipeline, gvcf): + pipeline.gvcf = gvcf + pipeline.bed = None + jobs = _jobs(pipeline) + stage2_bed = str(pipeline.tmp_dir / "hybrid_stage2.bed") + + assert "intersect-stage2-bed" not in jobs + call2, call2_deps = jobs["dnascope-2"] + assert f"--interval {stage2_bed}" in call2 + assert call2_deps == {"third-stage"} + assert stage2_bed in jobs["subset-calls"][0] + assert jobs["subset-calls"][1] == {"second-stage"} + if gvcf: + assert stage2_bed in jobs["subset-pass2-calls"][0] + assert "--interval" not in jobs["gvcftyper"][0] + + +class TestHybridWithoutGvcf: + def test_vcf_path(self, pipeline): + pipeline.gvcf = False + jobs = _jobs(pipeline) + tmp = pipeline.tmp_dir + + call1, _ = jobs["dnascope-1"] + assert call1.count("--algo DNAscope") == 1 + assert "--emit_mode" not in call1 + assert call1.endswith(str(tmp / "initial.vcf.gz")) + assert "--emit_mode" not in jobs["dnascope-2"][0] + assert f"-v {tmp / 'initial.vcf.gz'}" in jobs["hybrid-select"][0] + subset, _ = jobs["subset-calls"] + assert "bcftools view -T" in subset + assert "gvcf_clip.py" not in subset + assert "subset-pass2-calls" not in jobs + assert jobs["concat-calls"][1] == {"subset-calls", "dnascope-2"} + assert "--threads" not in jobs["concat-calls"][0] + assert "vcfconvert" not in jobs["anno-calls"][0] + assert jobs["anno-calls"][0].endswith( + str(tmp / "combined_tmp_anno.vcf.gz") + ) + assert "gvcftyper" not in jobs + assert "drop-allele-fields" not in jobs + + apply, _ = jobs["model-apply"] + assert apply.endswith(str(tmp / "combined_apply.vcf.gz")) + norm, norm_deps = jobs["final-norm"] + assert f"""-e 'GT="0/0"' {tmp / 'combined_apply.vcf.gz'}""" in norm + assert norm.endswith(str(pipeline.output_vcf)) + assert norm_deps == {"model-apply"} + assert not any(".g.vcf.gz" in shell for shell, _ in jobs.values()) diff --git a/tests/unit/test_gvcf_clip.py b/tests/unit/test_gvcf_clip.py new file mode 100644 index 0000000..75df0e5 --- /dev/null +++ b/tests/unit/test_gvcf_clip.py @@ -0,0 +1,410 @@ +""" +Unit tests for the `gvcf_clip.py` helper script. + +The script is not importable as part of the package (the `scripts` +directory has no `__init__.py`), so its helpers are loaded from the file +path, and the end-to-end tests run it the way the pipeline does, with +`sys.executable`. +""" + +import gzip +import importlib.util +import pathlib +import subprocess +import sys + +import pytest +import vcflib + +SCRIPT = ( + pathlib.Path(__file__).resolve().parents[2] + / "sentieon_cli" + / "scripts" + / "gvcf_clip.py" +) + + +def _load_script(): + spec = importlib.util.spec_from_file_location("gvcf_clip", SCRIPT) + assert spec is not None and spec.loader is not None + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +gvcf_clip = _load_script() + +SEQ = "ACGT" * 25 +HEADER = [ + "##fileformat=VCFv4.2", + "##contig=", + '##ALT=', + '##INFO=', + '##INFO=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS1", +] +BLOCK_SAMPLE = "GT:PL\t0/0:0,30,300" +VARIANT_SAMPLE = "GT:PL\t0/1:500,0,600,500,600,900" +# Records in 1-based coordinates, around the BED regions 21-30 and 61-70 +RECORDS = [ + f"chr1\t1\t.\tA\t\t.\t.\tEND=25;DP=10\t{BLOCK_SAMPLE}", + f"chr1\t26\t.\tC\t\t.\t.\t.\t{BLOCK_SAMPLE}", + f"chr1\t27\t.\tG\tC,\t50\t.\t.\t{VARIANT_SAMPLE}", + f"chr1\t28\t.\tT\t\t.\t.\tEND=40\t{BLOCK_SAMPLE}", + f"chr1\t41\t.\tA\tT,\t50\t.\t.\t{VARIANT_SAMPLE}", + f"chr1\t42\t.\tC\t\t.\t.\tEND=100\t{BLOCK_SAMPLE}", +] +REGIONS = "chr1\t20\t30\nchr1\t60\t70\n" + + +def _write_inputs(tmp_path, records=RECORDS, regions=REGIONS): + """Write the reference, the region BED, and a bgzipped, indexed gVCF""" + ref = tmp_path / "ref.fa" + # Two FASTA lines, so lookups cross a line break + ref.write_text(f">chr1\n{SEQ[:60]}\n{SEQ[60:]}\n") + (tmp_path / "ref.fa.fai").write_text("chr1\t100\t6\t60\t61\n") + + bed = tmp_path / "regions.bed" + bed.write_text(regions) + + gvcf = tmp_path / "in.g.vcf.gz" + out = vcflib.VCF(str(gvcf), "w") + out.headers = list(HEADER) + out.parse_header() + out.emit_header() + for line in records: + out.emit(out.parse(line)) + out.close() + return gvcf, bed, ref + + +def _run(tmp_path, exclude, threads=1, step_size=None, **inputs): + gvcf, bed, ref = _write_inputs(tmp_path, **inputs) + name = "exclude" if exclude else "include" + output = tmp_path / f"{name}-{threads}.g.vcf.gz" + cmd = [ + sys.executable, + str(SCRIPT), + "-v", + str(gvcf), + "-b", + str(bed), + "-r", + str(ref), + "-t", + str(threads), + str(output), + ] + if exclude: + cmd.append("--exclude") + if step_size: + cmd.extend(["--step-size", str(step_size)]) + subprocess.run(cmd, check=True) + assert pathlib.Path(str(output) + ".tbi").exists() + with gzip.open(output, "rt") as fh: + return [line.rstrip("\n") for line in fh if not line.startswith("#")] + + +def _columns(lines): + """The CHROM, POS, REF, ALT and INFO columns of each record""" + return [tuple(line.split("\t")[i] for i in (0, 1, 3, 4, 7)) for line in lines] + + +def _covered(lines): + """Every 1-based position the records span, with repeats""" + positions = [] + for line in lines: + cols = line.split("\t") + start = int(cols[1]) + end = start + len(cols[3]) - 1 + for fld in cols[7].split(";"): + if fld.startswith("END="): + end = int(fld[4:]) + positions.extend(range(start, end + 1)) + return positions + + +class TestRegions: + def test_load_regions_sorts_and_merges(self, tmp_path): + bed = tmp_path / "regions.bed" + bed.write_text( + "track name=test\n" + "# comment\n" + "chr1\t50\t60\textra\n" + "chr1\t10\t20\n" + "chr1\t15\t30\n" + "chr1\t30\t35\n" + "chr2\t0\t5\n" + ) + regions = gvcf_clip.load_regions(str(bed)) + assert regions == { + "chr1": ([10, 50], [35, 60]), + "chr2": ([0], [5]), + } + + def test_in_regions(self): + regions = {"chr1": ([10, 50], [20, 60])} + assert not gvcf_clip.in_regions(regions, "chr1", 9) + assert gvcf_clip.in_regions(regions, "chr1", 10) + assert gvcf_clip.in_regions(regions, "chr1", 19) + assert not gvcf_clip.in_regions(regions, "chr1", 20) + assert gvcf_clip.in_regions(regions, "chr1", 55) + assert not gvcf_clip.in_regions(regions, "chr2", 15) + + @pytest.mark.parametrize( + "start, end, include, expected", + [ + # Outside every region + (0, 5, True, []), + (0, 5, False, [(0, 5)]), + # Across the start of a region + (5, 15, True, [(10, 15)]), + (5, 15, False, [(5, 10)]), + # Across the end of a region + (15, 25, True, [(15, 20)]), + (15, 25, False, [(20, 25)]), + # Within a region + (12, 18, True, [(12, 18)]), + (12, 18, False, []), + # Spanning two regions + (0, 70, True, [(10, 20), (50, 60)]), + (0, 70, False, [(0, 10), (20, 50), (60, 70)]), + # Ending at a region start, starting at a region end + (0, 10, True, []), + (20, 50, False, [(20, 50)]), + ], + ) + def test_pieces(self, start, end, include, expected): + regions = {"chr1": ([10, 50], [20, 60])} + assert ( + gvcf_clip.pieces(regions, "chr1", start, end, include) == expected + ) + + def test_pieces_without_regions_on_contig(self): + assert gvcf_clip.pieces({}, "chr1", 0, 10, False) == [(0, 10)] + assert gvcf_clip.pieces({}, "chr1", 0, 10, True) == [] + + @pytest.mark.parametrize( + "info, expected", + [ + ("END=40", "END=50"), + ("DP=10;END=40;MQ=60", "DP=10;END=50;MQ=60"), + ("BLOCKEND=3", "BLOCKEND=3;END=50"), + (".", "END=50"), + ], + ) + def test_set_end(self, info, expected): + assert gvcf_clip.set_end(info, 50) == expected + + +class TestClip: + def test_exclude_clips_blocks_at_region_boundaries(self, tmp_path): + records = _run(tmp_path, exclude=True) + assert _columns(records) == [ + ("chr1", "1", "A", "", "END=20;DP=10"), + ("chr1", "31", "G", "", "END=40"), + ("chr1", "41", "A", "T,", "."), + ("chr1", "42", "C", "", "END=60"), + ("chr1", "71", "G", "", "END=100"), + ] + + def test_include_keeps_the_region_parts(self, tmp_path): + records = _run(tmp_path, exclude=False) + assert _columns(records) == [ + ("chr1", "21", "A", "", "END=25;DP=10"), + ("chr1", "26", "C", "", "."), + ("chr1", "27", "G", "C,", "."), + ("chr1", "28", "T", "", "END=30"), + ("chr1", "61", "A", "", "END=70"), + ] + + def test_clipped_blocks_keep_the_sample_columns(self, tmp_path): + records = _run(tmp_path, exclude=True) + for line in records: + cols = line.split("\t") + if cols[4] == "": + assert "\t".join(cols[8:]) == BLOCK_SAMPLE + + def test_the_two_sides_tile_the_input(self, tmp_path): + outside = _run(tmp_path, exclude=True) + inside = _run(tmp_path, exclude=False) + positions = _covered(outside) + _covered(inside) + assert sorted(positions) == list(range(1, 101)) + + @pytest.mark.parametrize("exclude", [True, False]) + def test_sharded_matches_single_thread(self, tmp_path, exclude): + single = _run(tmp_path, exclude) + # Shard boundaries at 30 and 60 fall inside blocks + sharded = _run(tmp_path, exclude, threads=2, step_size=30) + assert sharded == single + + def test_empty_bed(self, tmp_path): + assert _run(tmp_path, exclude=True, regions="") == RECORDS + assert _run(tmp_path, exclude=False, regions="") == [] + + +def _block(pos, end): + return ( + f"chr1\t{pos}\t.\t{SEQ[pos - 1]}\t\t.\t.\tEND={end}\t" + "GT:DP:GQ:MIN_DP:PL\t0/0:10:30:8:0,30,300" + ) + + +def _deletion(pos, length, gt="0/1", dp=12): + ref = SEQ[pos - 1 : pos - 1 + length] + return ( + f"chr1\t{pos}\t.\t{ref}\t{ref[0]},\t50\t.\t.\t" + f"GT:DP:PL\t{gt}:{dp}:500,0,600,500,600,900" + ) + + +def _snv(pos): + ref = SEQ[pos - 1] + alt = "C" if ref != "C" else "T" + return f"chr1\t{pos}\t.\t{ref}\t{alt},\t50\t.\t.\t{VARIANT_SAMPLE}" + + +def _fill(pos, end, dp=12): + return ( + "chr1", + str(pos), + SEQ[pos - 1], + "", + f"END={end}", + "GT:DP:GQ:MIN_DP:PL", + f"0/0:{dp}:0:{dp}:0,0,0", + ) + + +def _summary(lines): + """POS, END or REF span, and whether each record is a fill block""" + out = [] + for line in lines: + cols = line.split("\t") + out.append(tuple(cols[i] for i in (0, 1, 3, 4, 7)) + tuple(cols[8:])) + return out + + +OUTSIDE = list(range(1, 21)) + list(range(31, 61)) + list(range(71, 101)) +INSIDE = list(range(21, 31)) + list(range(61, 71)) + + +class TestFill: + """Positions spanned only by a dropped variant record are filled""" + + def test_tail_of_a_dropped_deletion(self, tmp_path): + # The deletion starts at 29, inside 21-30, and spans 29-33 + records = [_block(1, 28), _deletion(29, 5), _block(34, 60)] + records.append(_block(61, 100)) + out = _run(tmp_path, exclude=True, records=records) + + summary = _summary(out) + assert _fill(31, 33) in summary + assert [x[1] for x in summary] == ["1", "31", "34", "71"] + assert sorted(_covered(out)) == OUTSIDE + + def test_records_inside_the_tail_split_the_fill(self, tmp_path): + # The deletion spans 29-36; a SNV at 33 is inside its span + records = [_block(1, 28), _deletion(29, 8), _snv(33)] + records += [_block(37, 60), _block(61, 100)] + out = _run(tmp_path, exclude=True, records=records) + + summary = _summary(out) + assert [x[1] for x in summary] == ["1", "31", "33", "34", "37", "71"] + assert summary[1] == _fill(31, 32) + assert summary[3] == _fill(34, 36) + assert sorted(_covered(out)) == OUTSIDE + + def test_a_block_inside_the_tail_trims_the_fill(self, tmp_path): + records = [_block(1, 28), _deletion(29, 8), _block(33, 60)] + records.append(_block(61, 100)) + out = _run(tmp_path, exclude=True, records=records) + + summary = _summary(out) + assert [x[1] for x in summary] == ["1", "31", "33", "71"] + assert summary[1] == _fill(31, 32) + assert sorted(_covered(out)) == OUTSIDE + + def test_no_fill_where_a_kept_record_covers_the_tail(self, tmp_path): + # The deletion at 19 starts outside and spans 19-34 + records = [_block(1, 18), _deletion(19, 16), _deletion(29, 5)] + records += [_block(35, 60), _block(61, 100)] + out = _run(tmp_path, exclude=True, records=records) + + assert [x[1] for x in _summary(out)] == ["1", "19", "35", "71"] + + def test_include_side(self, tmp_path): + # The deletion at 59 starts outside 61-70 and spans 59-63 + records = [_block(1, 58), _deletion(59, 5, gt="0/0", dp=7)] + records.append(_block(64, 100)) + out = _run(tmp_path, exclude=False, records=records) + + summary = _summary(out) + assert [x[1] for x in summary] == ["21", "61", "64"] + assert summary[1] == _fill(61, 63, dp=7) + assert sorted(_covered(out)) == INSIDE + + def test_the_two_sides_tile_the_contig(self, tmp_path): + records = [_block(1, 28), _deletion(29, 8), _snv(33)] + records += [_block(37, 58), _deletion(59, 5), _block(64, 100)] + outside = _run(tmp_path, exclude=True, records=records) + inside = _run(tmp_path, exclude=False, records=records) + positions = _covered(outside) + _covered(inside) + assert set(positions) == set(range(1, 101)) + # The deletion kept at 59 spans 61-63, inside 61-70, as in the input + assert sorted(_covered(outside)) == sorted(OUTSIDE + [61, 62, 63]) + + def test_sharded_matches_single_thread(self, tmp_path): + records = [_block(1, 28), _deletion(29, 8), _snv(33)] + records += [_block(37, 60), _block(61, 100)] + single = _run(tmp_path, exclude=True, records=records) + # A shard boundary between positions 32 and 33 + sharded = _run( + tmp_path, exclude=True, threads=2, step_size=32, records=records + ) + assert sharded == single + + def test_fill_is_cut_at_a_shard_boundary(self, tmp_path): + records = [_block(1, 28), _deletion(29, 5), _block(34, 60)] + records.append(_block(61, 100)) + out = _run( + tmp_path, exclude=True, threads=2, step_size=32, records=records + ) + + summary = _summary(out) + assert summary[1] == _fill(31, 32) + assert summary[2] == _fill(33, 33) + assert sorted(_covered(out)) == OUTSIDE + + def test_stdout_output(self, tmp_path): + records = [_block(1, 28), _deletion(29, 5), _block(34, 60)] + records.append(_block(61, 100)) + expected = _run(tmp_path, exclude=True, records=records) + gvcf, bed, ref = _write_inputs(tmp_path, records=records) + result = subprocess.run( + [ + sys.executable, + str(SCRIPT), + "-v", + str(gvcf), + "-b", + str(bed), + "-r", + str(ref), + "-t", + "2", + "--exclude", + "-", + ], + check=True, + stdout=subprocess.PIPE, + text=True, + ) + lines = result.stdout.splitlines() + assert lines[0] == "##fileformat=VCFv4.2" + assert [x for x in lines if not x.startswith("#")] == expected diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index fbda6eb..21c4d70 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -3,6 +3,7 @@ """ from importlib.resources import files +import json import logging import os import pathlib @@ -25,6 +26,7 @@ from sentieon_cli.command_strings import LONGREAD_SV_BED_AWK from sentieon_cli.dag import DAG from sentieon_cli.util import SampleSex +from tests.utils.test_helpers import write_bundle PACKAGE_LOGGER = "sentieon_cli" @@ -214,6 +216,25 @@ def test_kmc_command(self): assert str(self.mock_r1) in cmd_str assert str(self.mock_r2) in cmd_str + @pytest.mark.parametrize("aligned", [False, True]) + def test_kmc_threads_capped(self, aligned): + """KMC is capped at 128 threads; the rest of the job is not""" + if aligned: + pipeline = self.create_aligned_pipeline() + name = "extract-kmc" + else: + pipeline = self.create_pipeline() + name = "kmc" + pipeline.cores = 256 + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, name).shell) + assert " -t128 " in cmd_str + assert "-t256" not in cmd_str + assert "samtools fasta" in cmd_str + assert "-@ 256" in cmd_str + def test_bwa_readgroup_lr0(self): """The bwa alignment carries the LR:0 readgroup attribute""" pipeline = self.create_pipeline() @@ -1455,7 +1476,8 @@ def test_cnvscope_without_a_bed(self): assert "--interval" not in str(job.shell) def test_cnv_model_apply_command(self): - """CNVModelApply filters the CNVscope output""" + """CNVModelApply filters the CNVscope output, keeping the rejected + segments""" pipeline = self.enable_cnv(self.create_pipeline()) job = self._second_dag_job(pipeline, "cnv-model-apply") @@ -1464,6 +1486,7 @@ def test_cnv_model_apply_command(self): f"--thread_count {pipeline.cores} " f"--algo CNVModelApply --model {self.mock_bundle}/cnv.model " f"--vcf {self.mock_dir}/sample-cnvscope.vcf.gz " + "--emit_conf 0 " f"{self.mock_dir}/sample-cnv_model_apply.vcf.gz" ) # The BED restricts CNVscope, not the model apply @@ -1490,7 +1513,8 @@ def test_cnvscope_female_sample(self): assert "--par" not in str(job.shell) def test_cnvscope_unknown_sex(self, messages): - """An unknown sex calls a diploid genome, with a warning""" + """An unknown sex calls assumed sex chromosome ploidies, with a + warning""" pipeline = self.enable_cnv(self.create_pipeline()) job = self._second_dag_job( pipeline, "cnvscope", sample_sex=SampleSex.UNKNOWN @@ -1498,7 +1522,8 @@ def test_cnvscope_unknown_sex(self, messages): assert "--sex" not in str(job.shell) assert "--par" not in str(job.shell) - assert any("diploid" in msg for msg in messages) + assert "--blocks" in str(job.shell) + assert any("Sex estimation failed" in msg for msg in messages) def test_indel2cnv_command(self): """The PangenomeSV output is converted to CNVs""" @@ -1522,9 +1547,42 @@ def test_combine_cnv_command(self): f"{sys.executable} {self.script_path('combine_cnv.py')} " f"--cnv {self.mock_dir}/sample-cnv_model_apply.vcf.gz " f"--converted {self.mock_dir}/sample-sv_cnv.vcf.gz " + "--sex F " f"-o {cnv_vcf}" ) + def test_combine_cnv_male_sample(self): + """combine_cnv.py gets the same sex and PAR BED file as CNVscope""" + pipeline = self.enable_cnv(self.create_pipeline()) + job = self._second_dag_job( + pipeline, "combine-cnv", sample_sex=SampleSex.MALE + ) + + assert f"--sex M --par {pipeline.cnv_par_bed}" in str(job.shell) + + def bundle_combine_cnv_shell(self, members): + """The combine-cnv command for a bundle of Ultima short reads""" + self.write_bundle( + members + ["cnv.model"], {"shortReadPlatform": "Ultima"} + ) + pipeline = self.enable_cnv(self.create_pipeline()) + pipeline.validate_bundle() + return str(self._second_dag_job(pipeline, "combine-cnv").shell) + + def test_combine_cnv_bundle_preset(self): + """The bundle's `combine_cnv.json` is the preset""" + shell = self.bundle_combine_cnv_shell( + self.BUNDLE_MEMBERS + ["combine_cnv.json"] + ) + + assert f"--preset {self.mock_bundle}/combine_cnv.json " in shell + + def test_combine_cnv_without_a_bundle_preset(self): + """Without one the script uses its defaults, whatever the platform""" + shell = self.bundle_combine_cnv_shell(self.BUNDLE_MEMBERS) + + assert "--preset" not in shell + def test_expansion_hunter_command(self): """ExpansionHunter genotypes the short reads""" catalog = self.mock_dir / "catalog.json" @@ -1794,19 +1852,31 @@ def test_validate_cnv_does_not_warn_with_a_bed(self): # Model bundle validation - def bundle_pipeline(self, monkeypatch, members): - """A pipeline whose bundle holds `members`""" + def bundle_pipeline(self, monkeypatch, members, bundle_info=None): + """A pipeline whose bundle holds `members` and `bundle_info`""" pipeline = self.create_pipeline() pipeline.skip_pop_vcf_id_check = True + info = {"pipeline": "Hybrid pangenome", **(bundle_info or {})} def fake_ar_load(path): if str(path).endswith("bundle_info.json"): - return b'{"pipeline": "Hybrid pangenome"}' + return json.dumps(info).encode() return list(members) monkeypatch.setattr(hybrid_pangenome, "ar_load", fake_ar_load) return pipeline + def write_bundle(self, members, bundle_info=None): + """Write the model bundle with `members` and a bundle_info.json""" + info = {"pipeline": "Hybrid pangenome", **(bundle_info or {})} + write_bundle( + self.mock_bundle, + { + "bundle_info.json": json.dumps(info).encode(), + **{name: b"" for name in members}, + }, + ) + BUNDLE_MEMBERS = [ "dnascope.model", "longreadsv.model", @@ -1846,3 +1916,59 @@ def test_validate_bundle_diploid_model_not_required_by_default( ): pipeline = self.bundle_pipeline(monkeypatch, self.BUNDLE_MEMBERS) pipeline.validate_bundle() # no SystemExit + + def test_validate_bundle_dnascope_model_with_skip_small_variants( + self, monkeypatch + ): + """`dnascope.model` is only needed for small-variant calling""" + members = [m for m in self.BUNDLE_MEMBERS if m != "dnascope.model"] + + pipeline = self.bundle_pipeline(monkeypatch, members) + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_bundle() + assert excinfo.value.code == 2 + + pipeline = self.bundle_pipeline(monkeypatch, members) + pipeline.skip_small_variants = True + pipeline.validate_bundle() # no SystemExit + + def test_validate_bundle_bwa_model_with_aligned_input(self, monkeypatch): + """`bwa.model` is only needed with FASTQ input""" + members = [m for m in self.BUNDLE_MEMBERS if m != "bwa.model"] + + pipeline = self.bundle_pipeline(monkeypatch, members) + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_bundle() + assert excinfo.value.code == 2 + + pipeline = self.bundle_pipeline(monkeypatch, members) + pipeline.r1_fastq = [] + pipeline.validate_bundle() # no SystemExit + + @pytest.mark.parametrize( + "bundle_info", + [ + {}, + {"shortReadPlatform": "Illumina"}, + {"shortReadPlatform": "Ultima"}, + ], + ) + def test_validate_bundle_combine_preset(self, bundle_info): + """The CNV combine preset is the bundle's `combine_cnv.json`""" + self.write_bundle( + self.BUNDLE_MEMBERS + ["combine_cnv.json"], bundle_info + ) + pipeline = self.create_pipeline() + pipeline.validate_bundle() + assert pipeline.has_combine_preset is True + assert ( + pipeline.cnv_combine_preset() + == self.mock_bundle / "combine_cnv.json" + ) + + # Without one, whatever the platform, the script's defaults apply + self.write_bundle(self.BUNDLE_MEMBERS, bundle_info) + pipeline = self.create_pipeline() + pipeline.validate_bundle() + assert pipeline.has_combine_preset is False + assert pipeline.cnv_combine_preset() is None diff --git a/tests/unit/test_indel2cnv_utils.py b/tests/unit/test_indel2cnv_utils.py deleted file mode 100644 index 9722b96..0000000 --- a/tests/unit/test_indel2cnv_utils.py +++ /dev/null @@ -1,167 +0,0 @@ -"""Golden-value tests for ``sentieon_cli._indel2cnv_utils``. - -Expected values were captured from ``Levenshtein.ratio`` (v0.27.3) and -``scipy.signal.find_peaks`` (v1.17.1). Keeping the deps as test-time -imports would defeat the point of this module, so the reference -outputs are frozen as literals here. -""" - -import os -import sys - -sys.path.insert( - 0, - os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")), -) - -import pytest - -from sentieon_cli._indel2cnv_utils import find_peaks, ratio - - -LONG_A = ( - "AAGCCCAATAAACCACTCTGACTGGCCGAATAGGGATATAGGCAACGACATGTGCGGCGACC" - "CTTGCGACAGTGACGCTTTCGCCGTTGCCTAAACCTATTTGAAGGAGTCTAGCAGCCGCAGT" - "AAGGCACAATACCTCGTCCGTGTTACCAGACCAAACAAGACGTCCTCTTCAATGTTTAAATG" - "ACCCTCTCGTCATAAAACCTTTCTACTATGTGTTCCGCAAGAATCAACAACTACAATGGCGC" - "GTCGTGAATAACGCGACGGCTGAGACGAACGGCGCGTGAATGAAGCGCTTAAACAGCTCAGG" - "AGCCAGTCCCCTACGTCGCATATCCTGGCCACTGGAGGTGAAGCGAATGGTATCGATACGTA" - "GGAGGTGTGCCTTCGTAGGCTGTTTCTCAGGACGCCCAACTATTCTTTCCAATCCTACATCT" - "GTTTCTTGCGTCGTAGCGGGACCCTCCATTGTTACTTATTAGGTTCTCGTTATGTCTCATAA" - "TCTC" -) -LONG_B = ( - "AAGCCCAATAAACCACTCTGCCTGGGGGAATAGGGTTATAGGCAACGACATGTGCGGGCTAC" - "CCTTGCGACAGTGACGCTTTCGCCAGTTGCCTAAACCTATTTAAAGGACTCTAGTAGCCGCA" - "ATAAGGCCCAATACCTCGTCGGTGTTACCAGACCAAACAAGACGTGCTCTTCAATGTTTAAA" - "TGACCCTCTCCTCATAAAACCTTTCTACTATGTGTTCCGCAAGAATCAACAACTGACCAGGG" - "CGCGTCGTGAATAACGCGACGGCTGAGACGAACGGCGCGTGAATGATAGCGCTTAAACAGCT" - "CAGGAGCCAGTCACCTACGTCGCATATCCTGGCCACTGGAGGTGAAGCGAATGGTATCGATA" - "CGTAGGAGGTGTGCCGTCGTACGCTGTTGTCTCAGGACGCCCAACTATTCTTTCCAATCCTA" - "CATCTGTTTCTTGCGTCGTAGCGGGACCCTCCATTGTTACTTATTAGGTTCTCGTTCAGTCT" - "CATAATCTC" -) -LONG_C = ( - "TCACAGAAGTGAGATTATGTCTCGTTTGGCAGTCTTGATGCTCGGGGGACACTTCTTTAAGC" - "TCGGTGTGGTGGGCACGACCCTGGACGCGCGACGAAGCTAAGTTTGCAGTAATTAACCGACA" - "TCTTTGTGAACCGACCCACATTTGACGGTACGCTACCGCAACGGTATGTGTTAATGGAACAG" - "ACTTGCTTATGTGGACGTTGTATAGGGATATTACGTTACGCGTTAACCGATACATACTGGTT" - "TCTCTCCAGTGGAGGTCTTGGTTGCCTCTAGTTTCTACGATATACTCATGGTAGTGTAACGC" - "ATAATCGAAGAGGGTCCTCCCATCTCCTGTGATGCATGGTGTGCTTACTGGGATGAATGCGC" - "CGCAAGTAGCAGGTCCCGGCGTGGATACCTGATAGATGGTGACTAGCATGTACAAGTAACCT" - "TGTCTATTGAGCTTCGAGGATGCATACAAGCCCACCCGCAGCCGCAACAGCGACGACTAATT" - "GATC" -) -RATIO_AB = 0.9572139303482587 -RATIO_AC = 0.636 - -TANDEM_UNIT = "ACGTACGAT" -TANDEM = TANDEM_UNIT * 6 -TANDEM_VAR = ( - TANDEM[: len(TANDEM) // 2] + "T" + TANDEM[len(TANDEM) // 2 + 1:] -) -TANDEM_RATIO = 0.9814814814814815 - - -class TestRatio: - """``ratio`` must match ``Levenshtein.ratio`` on the captured cases.""" - - @pytest.mark.parametrize( - "a,b,expected", - [ - ("ACGTACGT", "ACGTACGT", 1.0), - ("", "", 1.0), - ("ACGT", "", 0.0), - ("", "ACGT", 0.0), - ("AAAA", "CCCC", 0.0), - ("ACGT", "ACGGT", 0.8888888888888888), - ("ACGGT", "ACGT", 0.8888888888888888), - ("ACGT", "AGGT", 0.75), - ("ACGTACGT" * 4, "ACGTACGT" * 4, 1.0), - ("ACGTACGT" * 4, "CGTACGTA" * 4, 0.96875), - ("AACCGGTT", "AACCGGTG", 0.875), - (TANDEM, TANDEM_VAR, TANDEM_RATIO), - ], - ) - def test_golden_values(self, a, b, expected): - assert ratio(a, b) == pytest.approx(expected, abs=1e-12) - - def test_long_similar_dna(self): - assert ratio(LONG_A, LONG_B) == pytest.approx(RATIO_AB, abs=1e-12) - - def test_long_unrelated_dna(self): - assert ratio(LONG_A, LONG_C) == pytest.approx(RATIO_AC, abs=1e-12) - - def test_symmetry(self): - assert ratio(LONG_A, LONG_B) == ratio(LONG_B, LONG_A) - - def test_threshold_pass_returns_exact_value(self): - # Real ratio is above threshold — must return the exact value. - assert ratio(LONG_A, LONG_B, threshold=0.9) == pytest.approx( - RATIO_AB, abs=1e-12 - ) - - def test_threshold_fail_returns_below_threshold(self): - # Real ratio is below threshold — returned value must be < threshold. - result = ratio("AAAA", "CCCC", threshold=0.5) - assert result < 0.5 - - def test_threshold_zero_disables_abort(self): - # Whatever value threshold=0 returns must equal the no-threshold call. - assert ratio(LONG_A, LONG_C, threshold=0.0) == ratio(LONG_A, LONG_C) - - def test_threshold_very_asymmetric_lengths(self): - # Length-based upper bound rejects cheaply. - result = ratio("A" * 10, "A" * 10000, threshold=0.9) - assert result < 0.9 - - -class TestFindPeaks: - """``find_peaks`` must match ``scipy.signal.find_peaks(x, distance=d, - height=h)`` for the call pattern used by indel2cnv.py.""" - - @pytest.mark.parametrize( - "values,distance,height,expected", - [ - ([0.0, 0.5, 1.0, 0.5, 0.0], 1, 0.9, [2]), - ( - [0.0, 0.95, 0.0, 0.0, 0.0, 0.97, 0.0], - 2, 0.9, [1, 5], - ), - ([0.0, 0.95, 0.0, 0.97, 0.0], 3, 0.9, [3]), - ([0.0, 0.5, 0.0, 0.6, 0.0], 1, 0.9, []), - ([0.0, 0.95, 0.95, 0.0], 1, 0.9, [1]), - ([1.0, 0.9, 0.8], 1, 0.9, []), - ([0.8, 0.9, 1.0], 1, 0.9, []), - ( - [0.3, 0.5, 0.6, 0.95, 0.7, 0.6, 0.5, 0.96, 0.4, 0.3], - 2, 0.9, [3, 7], - ), - ], - ) - def test_golden_values(self, values, distance, height, expected): - assert find_peaks( - values, distance=distance, height=height - ) == expected - - def test_empty_input(self): - assert find_peaks([], distance=1, height=0.9) == [] - - def test_two_element_input(self): - # Neither element has both neighbors; never a peak. - assert find_peaks([0.0, 1.0], distance=1, height=0.5) == [] - - def test_left_edge_plateau_not_a_peak(self): - # scipy ignores plateaus without a strictly lower left neighbor. - assert find_peaks([0.95, 0.95, 0.0], distance=1, height=0.9) == [] - - def test_right_edge_plateau_not_a_peak(self): - # scipy ignores plateaus without a strictly lower right neighbor. - assert find_peaks([0.0, 0.95, 0.95], distance=1, height=0.9) == [] - - def test_wider_plateau_midpoint(self): - # Plateau of width 3 at indices 2..4; midpoint = 3. - assert find_peaks( - [0.0, 0.3, 0.95, 0.95, 0.95, 0.3, 0.0], - distance=1, height=0.9, - ) == [3] diff --git a/tests/unit/test_pangenome_numa_lanes.py b/tests/unit/test_pangenome_numa_lanes.py new file mode 100644 index 0000000..4a3fb72 --- /dev/null +++ b/tests/unit/test_pangenome_numa_lanes.py @@ -0,0 +1,137 @@ +"""bwa-extract and mm2-lift run as one lane per NUMA node""" + +import os +import sys +from unittest.mock import patch + +sys.path.insert( + 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +) + +import test_sentieon_pangenome # noqa: E402 + +TWO_NODES = ["0-95", "96-191"] +LANES = ["0-47", "48-95", "96-143", "144-191"] + + +class TestNumaLanes: + # Use the pangenome tests' fixtures without inheriting, and so + # re-running, their tests + def setup_method(self): + self.base = test_sentieon_pangenome.TestSentieonPangenome() + self.base.setup_method() + + def lane_pipeline(self, cores=192, nodes=TWO_NODES, no_split=False): + pipeline = self.base.create_fastq_pipeline() + pipeline.cores = cores + pipeline.no_split_alignment = no_split + with patch("sentieon_cli.util.find_numa_nodes", return_value=nodes): + pipeline.configure() + return pipeline + + def jobs(self, dag, name): + """The step's jobs: `name`, or `name-` when split""" + everything = list(dag.waiting_jobs) + list(dag.ready_jobs) + return [ + j + for j in everything + if j.name == name or j.name.rpartition("-")[0] == name + ] + + def test_configure_splits_like_dnascope(self): + assert self.lane_pipeline().numa_nodes == LANES + + def test_bwa_extract_lanes(self): + dag = self.lane_pipeline().build_dag() + bwa = self.jobs(dag, "bwa-extract") + assert len(bwa) == 4 + for j, cpus in enumerate(LANES): + job = next(b for b in bwa if f"taskset -c {cpus} " in str(b.shell)) + cmd = str(job.shell) + assert job.name == f"bwa-extract-{j}" + assert job.threads == 48 + assert job.resources == {f"node{j}": 1} + assert f"fqidx extract -F {j}/4 -K 10000000" in cmd + assert "bwa mem" in cmd and cmd.count("-K 10000000") == 2 + assert " -p " in cmd and "-t 48" in cmd + assert ( + f"sample-bwa-{j}.bam" in cmd and f"sample-bwa-{j}.fq.gz" in cmd + ) + # pgutil extract is pinned too + assert cmd.count(f"taskset -c {cpus}") == 2 + + def test_mm2_lanes_follow_their_bwa_lane(self): + pipeline = self.lane_pipeline() + dag = pipeline.build_dag() + mm2 = self.jobs(dag, "mm2-lift") + assert len(mm2) == 4 + for j, cpus in enumerate(LANES): + job = next( + m for m in mm2 if f"sample-bwa-{j}.fq.gz" in str(m.shell) + ) + cmd = str(job.shell) + assert job.name == f"mm2-lift-{j}" + assert job.threads == 48 and job.resources == {f"node{j}": 1} + assert cmd.count(f"taskset -c {cpus}") == 3 # minimap2, lift, sort + assert f"sample-mm2-{j}.bam" in cmd + deps = self.base._get_deps(dag, job) + lane_bwa = [d for d in deps if d.name.startswith("bwa-extract")] + assert [d.name for d in lane_bwa] == [f"bwa-extract-{j}"] + assert f"sample-bwa-{j}.fq.gz" in str(lane_bwa[0].shell) + assert {"check-sample-gfa", "vg-paths-fasta"} <= { + d.name for d in deps + } + + def test_dedups_read_every_lane(self): + dag = self.lane_pipeline().build_dag() + everything = list(dag.waiting_jobs) + list(dag.ready_jobs) + lc = [ + str(j.shell) + for j in everything + if "LocusCollector" in str(j.shell) + ] + bwa_lc = next(c for c in lc if "sample-bwa-0.bam" in c) + mm2_lc = next(c for c in lc if "sample-mm2-0.bam" in c) + for j in range(4): + assert f"sample-bwa-{j}.bam" in bwa_lc + assert f"sample-mm2-{j}.bam" in mm2_lc + + def test_no_split_alignment_keeps_one_job(self): + dag = self.lane_pipeline(no_split=True).build_dag() + bwa = self.jobs(dag, "bwa-extract") + mm2 = self.jobs(dag, "mm2-lift") + assert len(bwa) == 1 and len(mm2) == 1 + assert "taskset" not in str(bwa[0].shell) + str(mm2[0].shell) + assert "fqidx" not in str(bwa[0].shell) + assert "sample-bwa.bam" in str(bwa[0].shell) + assert bwa[0].name == "bwa-extract" and mm2[0].name == "mm2-lift" + # the same -K as a split run, so the alignments match + assert "-K 10000000" in str(bwa[0].shell) + assert "-K 2000000" not in str(bwa[0].shell) + + def test_small_machine_keeps_one_job(self): + dag = self.lane_pipeline(cores=48, nodes=["0-47"]).build_dag() + assert len(self.jobs(dag, "bwa-extract")) == 1 + assert len(self.jobs(dag, "mm2-lift")) == 1 + + def test_job_names_are_unique(self): + dag = self.lane_pipeline().build_dag() + everything = list(dag.waiting_jobs) + list(dag.ready_jobs) + names = [ + j.name + for j in everything + if j.name.startswith(("bwa-extract", "mm2-lift")) + ] + assert len(names) == 8 and len(set(names)) == 8 + + def test_bam_input_is_unchanged(self): + pipeline = self.base.create_pipeline() + pipeline.cores = 192 + with patch( + "sentieon_cli.util.find_numa_nodes", return_value=TWO_NODES + ): + pipeline.configure() + assert pipeline.numa_nodes == [] + dag = pipeline.build_dag() + mm2 = self.jobs(dag, "mm2-lift") + assert len(mm2) == 1 and "taskset" not in str(mm2[0].shell) diff --git a/tests/unit/test_pangenome_vg_scheduling.py b/tests/unit/test_pangenome_vg_scheduling.py new file mode 100644 index 0000000..d09c2ce --- /dev/null +++ b/tests/unit/test_pangenome_vg_scheduling.py @@ -0,0 +1,68 @@ +"""vg-haplotypes runs alongside bwa-extract when both fit in memory""" + +import os +import sys +from unittest.mock import patch + +sys.path.insert( + 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +) + +import test_sentieon_pangenome # noqa: E402 + +GB = 1024**3 +TWO_NODES = ["0-95", "96-191"] + + +class TestVgHaplotypesScheduling: + # Use the pangenome tests' fixtures without inheriting, and so + # re-running, their tests + def setup_method(self): + self.base = test_sentieon_pangenome.TestSentieonPangenome() + self.base.setup_method() + + def build(self, memory_gb, cores=64, nodes=None): + pipeline = self.base.create_fastq_pipeline() + pipeline.cores = cores + with patch( + "sentieon_cli.util.find_numa_nodes", return_value=nodes or [] + ): + pipeline.configure() + with patch( + "sentieon_cli.sentieon_pangenome.total_memory", + return_value=memory_gb * GB, + ): + dag = pipeline.build_dag() + everything = list(dag.waiting_jobs) + list(dag.ready_jobs) + vg = next(j for j in everything if j.name == "vg-haplotypes") + return pipeline, dag, vg + + def test_runs_alongside_bwa_with_enough_memory(self): + """No bwa dependency, and no threads reserved from the scheduler""" + _, dag, vg = self.build(384, cores=192, nodes=TWO_NODES) + deps = {d.name for d in self.base._get_deps(dag, vg)} + assert not any(d.startswith("bwa-extract") for d in deps) + assert "kmc" in deps + assert vg.threads == 0 + # vg itself still uses every core + assert "-t 192" in str(vg.shell) + + def test_waits_for_bwa_when_memory_is_short(self): + """A 128 GB host has no room for vg next to a bwa lane""" + pipeline, dag, vg = self.build(128) + deps = {d.name for d in self.base._get_deps(dag, vg)} + assert "bwa-extract" in deps + assert vg.threads == pipeline.cores + + def test_waits_for_every_bwa_lane(self): + """The memory check counts each lane, and vg waits for all of them""" + _, dag, vg = self.build(200, cores=192, nodes=TWO_NODES) + bwa = [ + d + for d in self.base._get_deps(dag, vg) + if d.name.startswith("bwa-extract") + ] + assert sorted(d.name for d in bwa) == [ + f"bwa-extract-{j}" for j in range(4) + ] + assert vg.threads == 192 diff --git a/tests/unit/test_pipeline_validation.py b/tests/unit/test_pipeline_validation.py index ba5a546..c4b408c 100644 --- a/tests/unit/test_pipeline_validation.py +++ b/tests/unit/test_pipeline_validation.py @@ -48,9 +48,12 @@ def setup_method(self): self.mock_bundle = pathlib.Path(self.temp_dir) / "model.bundle" # Create empty files - for file_path in [self.mock_ref, self.mock_bam, self.mock_fastq, self.mock_bundle]: + for file_path in [self.mock_ref, self.mock_bam, self.mock_fastq]: file_path.touch() + # An empty ar archive: a bundle without a bundle_info.json + self.mock_bundle.write_bytes(b"!\n") + # Create BWA index files next to the reference for suf in (".amb", ".ann", ".bwt", ".pac", ".sa"): pathlib.Path(str(self.mock_ref) + suf).touch() diff --git a/tests/unit/test_sentieon_pangenome.py b/tests/unit/test_sentieon_pangenome.py index 9d96fe2..18ef369 100644 --- a/tests/unit/test_sentieon_pangenome.py +++ b/tests/unit/test_sentieon_pangenome.py @@ -18,6 +18,7 @@ from sentieon_cli.sentieon_pangenome import SentieonPangenome from sentieon_cli.util import SampleSex from sentieon_cli.dag import DAG +from tests.utils.test_helpers import write_bundle class TestSentieonPangenome: @@ -139,6 +140,17 @@ def test_dedup_metrics_output(self): mm2_dedup = next(j for j in all_jobs if j.name == "dedup-mm2") assert "--metrics" not in str(mm2_dedup.shell) + @pytest.mark.parametrize("cores,kmc_threads", [(2, 2), (256, 128)]) + def test_kmc_threads(self, cores, kmc_threads): + """KMC runs with the pipeline's cores, capped at 128 threads""" + pipeline = self.create_fastq_pipeline() + pipeline.cores = cores + dag = pipeline.build_dag() + + _, all_jobs = self._get_all_job_names(dag) + kmc_job = next(j for j in all_jobs if j.name == "kmc") + assert f" -t{kmc_threads} " in str(kmc_job.shell) + def test_metrics_input_bwa_only(self): """The metrics job reads only the bwa alignment; including the mm2 alignment would count the extracted reads it re-aligns twice.""" @@ -267,6 +279,18 @@ def test_gvcf_adds_gvcftyper_and_routes_outputs(self): assert apply_gvcf not in cmd_str assert str(pipeline.output_vcf) in cmd_str assert dag.waiting_jobs[gvcftyper_job] == {ad_update_job} + assert "--dbsnp" not in cmd_str + + def test_gvcftyper_annotates_dbsnp(self): + pipeline = self.create_pipeline() + pipeline.gvcf = True + pipeline.dbsnp = pathlib.Path(self.temp_dir) / "dbsnp.vcf.gz" + dag = pipeline.build_dag() + + gvcftyper_job = next( + j for j in dag.waiting_jobs if j.name == "gvcftyper" + ) + assert f"--dbsnp {pipeline.dbsnp}" in str(gvcftyper_job.shell) def snv_apply_vcf(self, pipeline): """The small-variant model-apply intermediate for a pipeline""" @@ -933,6 +957,8 @@ def test_cnv_model_apply_command(self): cmd_str = str(job.shell) assert "--algo CNVModelApply" in cmd_str assert "cnv.model" in cmd_str + # The rejected segments are written too + assert "--emit_conf 0 " in cmd_str def test_indel2cnv_command(self): """Test indel2cnv script command""" @@ -968,10 +994,80 @@ def test_combine_cnv_command(self): assert "combine_cnv.py" in cmd_str assert "--cnv" in cmd_str assert "--converted" in cmd_str + # CNVModelApply keeps the rejected segments, so there is no raw VCF + assert "--raw" not in cmd_str + # A bundle without `combine_cnv.json` leaves the script its defaults + assert "--preset" not in cmd_str + # The same ploidy as CNVscope, a female sample + assert "--sex F" in cmd_str + assert "--par" not in cmd_str # Output should use _cnv.vcf.gz suffix cnv_vcf = str(pipeline.output_vcf).replace(".vcf.gz", "_cnv.vcf.gz") assert cnv_vcf in cmd_str + BUNDLE_MEMBERS = [ + "dnascope.model", + "extract.model", + "minimap2.model", + "cnv.model", + ] + + def write_bundle(self, members, platform): + """Write the model bundle with `members` and a bundle_info.json""" + bundle_info = { + "minScriptVersion": "1.0.0", + "pipeline": "Sentieon pangenome", + "platform": platform, + "SentieonVcfID": "population", + } + write_bundle( + self.mock_bundle, + { + "bundle_info.json": json.dumps(bundle_info).encode(), + **{name: b"" for name in members}, + }, + ) + + def bundle_combine_cnv_shell(self): + """The combine-cnv command for the bundle on disk""" + pipeline = self.create_pipeline() + pipeline.call_svs = True + pipeline.validate_bundle() + dag = self._build_cnv_dag(pipeline) + + _, all_jobs = self._get_all_job_names(dag) + job = next(j for j in all_jobs if j.name == "combine-cnv") + return str(job.shell) + + @pytest.mark.parametrize("platform", ["Illumina", "Ultima"]) + def test_combine_cnv_bundle_preset(self, platform): + """The bundle's `combine_cnv.json` is the preset""" + self.write_bundle(self.BUNDLE_MEMBERS + ["combine_cnv.json"], platform) + + shell = self.bundle_combine_cnv_shell() + assert f"--preset {self.mock_bundle}/combine_cnv.json " in shell + + @pytest.mark.parametrize("platform", ["Illumina", "Ultima"]) + def test_combine_cnv_without_a_bundle_preset(self, platform): + """Without one the script uses its defaults, whatever the platform""" + self.write_bundle(self.BUNDLE_MEMBERS, platform) + + assert "--preset" not in self.bundle_combine_cnv_shell() + + def test_combine_cnv_male_sample(self): + """A male sample is combined with the PAR BED file""" + par_bed = self.mock_dir / "par.bed" + par_bed.touch() + + pipeline = self.create_pipeline() + pipeline.call_svs = True + pipeline.cnv_par_bed = par_bed + dag = self._build_cnv_dag(pipeline, sample_sex=SampleSex.MALE) + + _, all_jobs = self._get_all_job_names(dag) + job = next(j for j in all_jobs if j.name == "combine-cnv") + assert f"--sex M --par {par_bed}" in str(job.shell) + def test_cnv_with_skip_model_apply(self): """Test CNV jobs are added even with skip_model_apply""" pipeline = self.create_pipeline() @@ -988,3 +1084,29 @@ def test_cnv_with_skip_model_apply(self): assert "cnv-model-apply" in job_names assert "indel2cnv" in job_names assert "combine-cnv" in job_names + + def test_validate_bundle_dnascope_model_with_skip_small_variants(self): + """`dnascope.model` is only needed for small-variant calling""" + bundle_info = { + "minScriptVersion": "1.0.0", + "pipeline": "Sentieon pangenome", + "platform": "Illumina", + "SentieonVcfID": "population", + } + + def fake_ar_load(path): + if str(path).endswith("bundle_info.json"): + return json.dumps(bundle_info).encode() + return ["extract.model", "minimap2.model", "cnv.model"] + + with patch( + "sentieon_cli.sentieon_pangenome.ar_load", fake_ar_load + ): + pipeline = self.create_pipeline() + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_bundle() + assert excinfo.value.code == 2 + + pipeline = self.create_pipeline() + pipeline.skip_small_variants = True + pipeline.validate_bundle() # no SystemExit diff --git a/tests/unit/test_stage_alignment.py b/tests/unit/test_stage_alignment.py index 6e4829e..787c118 100644 --- a/tests/unit/test_stage_alignment.py +++ b/tests/unit/test_stage_alignment.py @@ -396,6 +396,26 @@ def test_unzip_reaches_the_command(self, tmp_path): assert "igzip -dc" in str(result.jobs[0].shell) + def test_bwa_k_default_and_override(self, tmp_path): + default = str(self.make(tmp_path).build().jobs[0].shell) + override = str( + self.make(tmp_path, bwa_k="10000000").build().jobs[0].shell + ) + + assert "-K 2000000" in default + assert "-K 10000000" in override and "-K 2000000" not in override + + def test_split_lanes_have_unique_names(self, tmp_path): + result = self.make(tmp_path, numa_nodes=["0-47", "48-95"]).build() + + assert [j.name for j in result.jobs] == [ + "bwa-extract-0", + "bwa-extract-1", + ] + # fqidx extract needs -K of at least 10,000,000 + for job in result.jobs: + assert str(job.shell).count("-K 10000000") == 2 + class TestReadWriterStage: """`sentieon driver --algo ReadWriter` merges and extracts""" diff --git a/tests/unit/test_stage_cnv.py b/tests/unit/test_stage_cnv.py index e0ee30c..ae29ffd 100644 --- a/tests/unit/test_stage_cnv.py +++ b/tests/unit/test_stage_cnv.py @@ -2,7 +2,10 @@ Unit tests for the CNVscope stage """ +import dataclasses +import importlib.util from importlib.resources import files +import json import logging import pathlib import sys @@ -19,6 +22,7 @@ PangenomeCNVStage, ) from sentieon_cli.util import SampleSex +from tests.utils.test_helpers import write_bundle PACKAGE_LOGGER = "sentieon_cli" @@ -90,18 +94,21 @@ def test_cnvscope_command(self, tmp_path): f"sentieon driver --input {tmp_path}/sample.cram " f"--reference {tmp_path}/ref.fa --thread_count 4 " f"--algo CNVscope --model {tmp_path}/bundle/cnv.model " + f"--blocks {tmp_path}/cnvscope-blocks.bed " f"{tmp_path}/cnvscope.vcf.gz" ) def test_cnv_model_apply_command(self, tmp_path): result = make_stage(tmp_path).add_to(DAG()) - # CNVModelApply reads the CNVscope VCF; it takes no alignment + # CNVModelApply reads the CNVscope VCF; it takes no alignment. + # `--emit_conf 0` writes the rejected segments too assert str(result.apply_job.shell) == ( f"sentieon driver --reference {tmp_path}/ref.fa " f"--thread_count 4 --algo CNVModelApply " f"--model {tmp_path}/bundle/cnv.model " f"--vcf {tmp_path}/cnvscope.vcf.gz " + "--emit_conf 0 " f"{tmp_path}/output.cnv.vcf.gz" ) @@ -153,6 +160,8 @@ def test_male_sample(self, tmp_path): shell = str(result.cnvscope_job.shell) assert "--sex M" in shell assert f"--par {par_bed}" in shell + # The result records the ploidy CNVscope was called with + assert (result.sex, result.par) == ("M", par_bed) def test_female_sample(self, tmp_path): # A female sample does not need the PAR regions @@ -171,9 +180,42 @@ def test_without_a_known_sex(self, tmp_path, sample_sex, messages): result = make_stage(tmp_path, sample_sex=sample_sex).add_to(DAG()) shell = str(result.cnvscope_job.shell) + blocks = tmp_path / "cnvscope-blocks.bed" assert "--sex" not in shell assert "--par" not in shell - assert any("diploid" in msg for msg in messages) + assert f"--blocks {blocks}" in shell + assert result.blocks == blocks + # A dry run does not write the file + assert not blocks.exists() + assert any("Sex estimation failed" in msg for msg in messages) + + def test_without_a_known_sex_writes_the_blocks(self, tmp_path): + (tmp_path / "ref.fa.fai").write_text( + "chr20\t64444167\t7\t60\t61\n" + "chrX\t156040895\t65518245\t60\t61\n" + "chrY\t57227415\t224159863\t60\t61\n" + ) + par_bed = tmp_path / "par.bed" + par_bed.write_text( + "chrX\t10000\t2781479\nchrY\t10000\t2781479\n" + "chrX\t155701382\t156030895\nchrY\t56887902\t57217415\n" + ) + ctx = dataclasses.replace(make_ctx(tmp_path), dry_run=False) + + result = make_stage( + tmp_path, + ctx=ctx, + sample_sex=SampleSex.UNKNOWN, + par_bed=par_bed, + ).add_to(DAG()) + + assert result.blocks is not None + assert result.blocks.read_text() == ( + "chrX\t0\t156040895\t2\n" + "chrY\t0\t10000\t1\n" + "chrY\t2781479\t56887902\t1\n" + "chrY\t57217415\t57227415\t1\n" + ) class TestJobMetadata: @@ -280,6 +322,16 @@ def test_derived_paths(self, tmp_path): assert result.cnv_vcf == tmp_path / "output_cnv.vcf.gz" + def test_cnv_model_apply_writes_every_segment(self, tmp_path): + """combine_cnv.py reads the rejected segments from CNVModelApply""" + result = make_pangenome_stage(tmp_path).add_to(DAG()) + + assert "--emit_conf 0 " in str(result.apply_job.shell) + # So it needs no raw CNVscope VCF + combine_shell = str(result.combine_job.shell) + assert "--raw" not in combine_shell + assert "sample-cnvscope.vcf.gz" not in combine_shell + def test_indel2cnv_command(self, tmp_path): result = make_pangenome_stage(tmp_path).add_to(DAG()) @@ -290,15 +342,73 @@ def test_indel2cnv_command(self, tmp_path): ) def test_combine_cnv_command(self, tmp_path): + """Without a preset combine_cnv.py uses its defaults""" result = make_pangenome_stage(tmp_path).add_to(DAG()) assert str(result.combine_job.shell) == ( f"{sys.executable} {script_path('combine_cnv.py')} " f"--cnv {tmp_path}/sample-cnv_model_apply.vcf.gz " f"--converted {tmp_path}/sample-sv_cnv.vcf.gz " + f"--blocks {tmp_path}/cnvscope-blocks.bed " f"-o {tmp_path}/output_cnv.vcf.gz" ) + def test_combine_cnv_preset(self, tmp_path): + """The preset, a model bundle member, is passed as it is""" + preset = tmp_path / "model.bundle" / "combine_cnv.json" + stage = make_pangenome_stage(tmp_path, combine_preset=preset) + result = stage.add_to(DAG()) + + assert str(result.combine_job.shell) == ( + f"{sys.executable} {script_path('combine_cnv.py')} " + f"--cnv {tmp_path}/sample-cnv_model_apply.vcf.gz " + f"--converted {tmp_path}/sample-sv_cnv.vcf.gz " + f"--preset {preset} " + f"--blocks {tmp_path}/cnvscope-blocks.bed " + f"-o {tmp_path}/output_cnv.vcf.gz" + ) + + def test_combine_cnv_male_sample(self, tmp_path): + """combine_cnv.py gets the same sex and PAR BED file as CNVscope""" + par_bed = tmp_path / "par.bed" + result = make_pangenome_stage( + tmp_path, sample_sex=SampleSex.MALE, par_bed=par_bed + ).add_to(DAG()) + + assert str(result.combine_job.shell) == ( + f"{sys.executable} {script_path('combine_cnv.py')} " + f"--cnv {tmp_path}/sample-cnv_model_apply.vcf.gz " + f"--converted {tmp_path}/sample-sv_cnv.vcf.gz " + f"--sex M --par {par_bed} " + f"-o {tmp_path}/output_cnv.vcf.gz" + ) + + def test_combine_cnv_female_sample(self, tmp_path): + result = make_pangenome_stage( + tmp_path, + sample_sex=SampleSex.FEMALE, + par_bed=tmp_path / "par.bed", + ).add_to(DAG()) + + shell = str(result.combine_job.shell) + assert "--sex F" in shell + assert "--par" not in shell + + def test_combine_cnv_unknown_sex(self, tmp_path): + """Without a known sex both tools get the same assumed ploidies""" + result = make_pangenome_stage( + tmp_path, + sample_sex=SampleSex.UNKNOWN, + par_bed=tmp_path / "par.bed", + ).add_to(DAG()) + + blocks = tmp_path / "cnvscope-blocks.bed" + for job in (result.cnvscope_job, result.combine_job): + shell = str(job.shell) + assert "--sex" not in shell + assert "--par" not in shell + assert f"--blocks {blocks}" in shell + def test_upstream_lands_on_the_entry_jobs(self, tmp_path): dag = DAG() upstream = rm_job([tmp_path / "upstream"], "upstream") @@ -372,3 +482,41 @@ def test_result_fields(self, tmp_path): ] assert result.terminal == {result.combine_job} assert result.cnv_vcf == tmp_path / "output_cnv.vcf.gz" + + +def load_combine_cnv(): + """The packaged combine_cnv.py, imported as a module""" + spec = importlib.util.spec_from_file_location( + "combine_cnv", script_path("combine_cnv.py") + ) + assert spec is not None and spec.loader is not None + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +class TestCombinePreset: + """combine_cnv.py reads the preset the pipelines pass""" + + def test_bundle_preset_loads_with_combine_cnv(self, tmp_path): + """`--preset /combine_cnv.json` names a bundle member""" + preset = { + "description": "Single-end CNVscope models", + "max_converted_loss": 3000, + "loss_support_size": 10000, + } + bundle = write_bundle( + tmp_path / "model.bundle", + { + "cnv.model": b"", + "combine_cnv.json": json.dumps(preset).encode(), + }, + ) + combine_cnv = load_combine_cnv() + + rules = combine_cnv.load_preset(bundle / "combine_cnv.json") + assert rules == dict( + combine_cnv.DEFAULTS, + max_converted_loss=3000, + loss_support_size=10000, + ) diff --git a/tests/unit/test_stage_small_variants.py b/tests/unit/test_stage_small_variants.py index 56ae89f..e7af083 100644 --- a/tests/unit/test_stage_small_variants.py +++ b/tests/unit/test_stage_small_variants.py @@ -268,6 +268,23 @@ def test_interval(self, tmp_path): assert f"--interval {bed}" in str(result.job.shell) assert result.output == tmp_path / "sample.vcf.gz" + def test_dbsnp(self, tmp_path): + dbsnp = tmp_path / "dbsnp.vcf.gz" + result = GVCFtyperStage( + ctx=make_ctx(tmp_path), + gvcf=tmp_path / "sample.g.vcf.gz", + output=tmp_path / "sample.vcf.gz", + dbsnp=dbsnp, + ).add_to(DAG()) + + assert str(result.job.shell) == ( + f"sentieon driver --reference {tmp_path}/ref.fa " + f"--thread_count 4 --algo GVCFtyper " + f"--vcf {tmp_path}/sample.g.vcf.gz " + f"--dbsnp {dbsnp} " + f"{tmp_path}/sample.vcf.gz" + ) + class TestTransferApplyStage: """Transfer, apply, or both""" diff --git a/tests/unit/test_stage_transfer.py b/tests/unit/test_stage_transfer.py index d9c03f7..7f4de84 100644 --- a/tests/unit/test_stage_transfer.py +++ b/tests/unit/test_stage_transfer.py @@ -2,13 +2,24 @@ Unit tests for the annotation-transfer stage """ +import asyncio import pathlib +import shutil +import subprocess +from importlib.resources import files from types import SimpleNamespace import pytest +from sentieon_cli.command_strings import ( + KEEP_CONTIG_AWK, + bcftools_concat, + cmd_bcftools_merge_trim, + cmd_bcftools_subset_remove, +) from sentieon_cli.dag import DAG from sentieon_cli.shard import Shard +from sentieon_cli.shell_pipeline import Context, InputProcSub from sentieon_cli.stages.base import StageContext, rm_job from sentieon_cli.stages.transfer import TransferConfig, TransferStage @@ -92,6 +103,345 @@ def test_concat_writes_the_output(self, tmp_path): assert result.out_vcf == out_vcf +class TestMergeCommand: + """The shard merge keeps the input VCF's QUAL and drops FILTER""" + + @staticmethod + def shard_job(tmp_path, name): + result = make_stage(tmp_path).add_to(DAG()) + (job,) = [job for job in result.shard_jobs if job.name == name] + return job + + def merge_cmd(self, tmp_path): + return self.shard_job(tmp_path, "merge-trim-0").shell.nodes[0] + + def test_merge_reads_both_inputs_unindexed(self, tmp_path): + merge = self.merge_cmd(tmp_path) + + assert merge.args[0] == "merge" + assert "--no-index" in merge.args + # `-R` cannot be combined with `--no-index`; the inputs are + # subset as they are read instead + assert "--regions-file" not in merge.args + assert "--regions-overlap" not in merge.args + raw_sub, pop_sub = merge.args[-2:] + assert isinstance(raw_sub, InputProcSub) + assert isinstance(pop_sub, InputProcSub) + + def test_both_inputs_are_subset_to_the_shard(self, tmp_path): + merge = self.merge_cmd(tmp_path) + + beds = set() + for sub, vcf in zip(merge.args[-2:], ("raw.vcf.gz", "pop.vcf.gz")): + annotate, _ = sub.node.nodes + assert annotate.args[0] == "annotate" + assert annotate.args[-1] == str(tmp_path / vcf) + i = annotate.args.index("--regions-file") + beds.add(annotate.args[i + 1]) + i = annotate.args.index("--regions-overlap") + assert annotate.args[i + 1] == "pos" + # The same shard BED for both inputs + (bed,) = beds + assert bed.endswith("sample-dnascope_transfer-subset0.bed") + + def test_removed_annotations(self, tmp_path): + merge = self.merge_cmd(tmp_path) + removed = [] + for sub in merge.args[-2:]: + annotate, _ = sub.node.nodes + removed.append(annotate.args[annotate.args.index("-x") + 1]) + + # QUAL goes from the pop VCF only; FILTER from both + assert removed == ["FILTER", "QUAL,FILTER"] + + def test_both_inputs_keep_only_the_shard_contig(self, tmp_path): + merge = self.merge_cmd(tmp_path) + for sub in merge.args[-2:]: + _, keep = sub.node.nodes + assert keep.executable == "awk" + assert keep.args == ["-v", "contig=chr1", KEEP_CONTIG_AWK] + + def test_extra_contig_drops_filter_too(self, tmp_path): + (annotate,) = self.shard_job(tmp_path, "merge-trim-extra").shell.nodes + + assert annotate.args[0] == "annotate" + assert annotate.args[annotate.args.index("-x") + 1] == "FILTER" + assert annotate.args[-1] == str(tmp_path / "raw.vcf.gz") + i = annotate.args.index("--regions-file") + assert annotate.args[i + 1].endswith( + "sample-dnascope_transfer-subset2.bed" + ) + + +HEADER = ( + "##fileformat=VCFv4.2\n" + "##contig=\n" + '##FILTER=\n' + '##FILTER=\n' +) +RAW_VCF = ( + HEADER + + '##INFO=\n' + + '##FORMAT=\n' + + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tS\n" + + "chr1\t100\t.\tC\tG\t30\tPASS\tAF=0.5\tGT\t0/1\n" + + "chr1\t200\t.\tC\tG\t150\tLowQual\tAF=0.5\tGT\t0/1\n" + + "chr1\t300\t.\tA\tC\t.\tPASS\tAF=0.5\tGT\t0/1\n" + + "chr1\t500\t.\tA\tC\t25\tPASS\tAF=0.5\tGT\t0/1\n" + + "chr1\t600\t.\tA\tC\t35\tPASS\tAF=0.5\tGT\t0/1\n" +) +POP_VCF = ( + HEADER + + '##INFO=\n' + + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + + "chr1\t100\trs100\tC\tG\t100\tAC0\tAF_v20=0.1\n" + + "chr1\t200\trs200\tC\tG\t50\tPASS\tAF_v20=0.2\n" + + "chr1\t300\trs300\tA\tC\t60\tPASS\tAF_v20=0.3\n" + + "chr1\t550\trs550\tG\tC\t999\tPASS\tAF_v20=0.05\n" + + "chr1\t600\trs600\tA\tG\t500\tAC0\tAF_v20=0.6\n" +) + + +@pytest.mark.skipif( + shutil.which("bcftools") is None, reason="bcftools is not installed" +) +class TestMergeTrimRun: + """The transfer's shard commands run through the shell pipeline""" + + @staticmethod + def bgzip(tmp_path, name, text): + plain = tmp_path / f"{name}.vcf" + plain.write_text(text) + out = tmp_path / f"{name}.vcf.gz" + subprocess.run( + ["bcftools", "view", "-Oz", "-W=tbi", "-o", str(out), str(plain)], + check=True, + ) + return out + + @staticmethod + async def run(pipeline): + context = Context() + try: + proc = await pipeline.run(context) + return await proc.async_wait() + finally: + await context.cleanup() + + def test_keeps_the_input_qual_and_drops_filter(self, tmp_path): + raw_vcf = self.bgzip(tmp_path, "raw", RAW_VCF) + pop_vcf = self.bgzip(tmp_path, "pop", POP_VCF) + subset_bed = tmp_path / "shard.bed" + subset_bed.write_text("chr1\t0\t1000\n") + shard_vcf = tmp_path / "shard.vcf.gz" + trim_script = pathlib.Path( + str(files("sentieon_cli.scripts").joinpath("trimalt.py")) + ) + + pipeline = cmd_bcftools_merge_trim( + shard_vcf, + raw_vcf, + pop_vcf, + trim_script, + subset_bed, + "chr1", + merge_rules="AF_v20:sum", + merge_xargs=["--no-version", "-m", "all"], + view_xargs=["--no-version"], + ) + assert asyncio.run(self.run(pipeline)) == 0 + + query = subprocess.run( + [ + "bcftools", + "query", + "-f", + "%POS\t%ID\t%ALT\t%QUAL\t%FILTER\t%INFO/AF_v20\n", + str(shard_vcf), + ], + check=True, + capture_output=True, + text=True, + ) + rows = [line.split("\t") for line in query.stdout.splitlines()] + assert rows == [ + # QUAL is the input's, whichever is higher, and FILTER is + # dropped from both inputs + ["100", "rs100", "G", "30", ".", "0.1"], + ["200", "rs200", "G", "150", ".", "0.2"], + # A missing input QUAL is not filled from the pop VCF + ["300", "rs300", "C", ".", ".", "0.3"], + ["500", ".", "C", "25", ".", "."], + # The pop allele is trimmed and its QUAL does not survive + ["600", "rs600", "C", "35", ".", "0"], + ] + + def merge_shard(self, tmp_path, raw_vcf, pop_vcf, contig, i): + """Run the merge of one whole-contig shard""" + subset_bed = tmp_path / f"shard{i}.bed" + subset_bed.write_text(f"{contig}\t0\t1000\n") + shard_vcf = tmp_path / f"shard{i}.vcf.gz" + trim_script = pathlib.Path( + str(files("sentieon_cli.scripts").joinpath("trimalt.py")) + ) + pipeline = cmd_bcftools_merge_trim( + shard_vcf, + raw_vcf, + pop_vcf, + trim_script, + subset_bed, + contig, + merge_rules="AF_v20:sum", + merge_xargs=["--no-version", "-m", "all"], + view_xargs=["--no-version"], + ) + assert asyncio.run(self.run(pipeline)) == 0 + return shard_vcf + + @staticmethod + def contig_lines(vcf): + header = subprocess.run( + ["bcftools", "view", "-h", str(vcf)], + check=True, + capture_output=True, + text=True, + ) + return [ + line + for line in header.stdout.splitlines() + if line.startswith("##contig") + ] + + def test_contig_lists_that_differ(self, tmp_path): + # The input VCF lists the reference's contigs; the pop VCF lists + # other contigs, in another order + raw_vcf = self.bgzip( + tmp_path, + "raw", + RAW_VCF.replace( + "##contig=\n", + "##contig=\n" + "##contig=\n" + "##contig=\n", + ), + ) + pop_vcf = self.bgzip( + tmp_path, + "pop", + POP_VCF.replace( + "##contig=\n", + "##contig=\n" + "##contig=\n" + "##contig=\n", + ), + ) + shard_vcf = self.merge_shard(tmp_path, raw_vcf, pop_vcf, "chr1", 0) + + query = subprocess.run( + ["bcftools", "query", "-f", "%POS\t%ID\t%QUAL\n", str(shard_vcf)], + check=True, + capture_output=True, + text=True, + ) + assert query.stdout.splitlines() == [ + "100\trs100\t30", + "200\trs200\t150", + "300\trs300\t.", + "500\t.\t25", + "600\trs600\t35", + ] + # The shard lists its own contig, with the input VCF's line + assert self.contig_lines(shard_vcf) == [ + "##contig=" + ] + + def test_concat_lists_the_contigs_in_reference_order(self, tmp_path): + # chrM, which the pop VCF lacks, sits between the merged contigs + contigs = "".join( + f"##contig=\n" for c in ("chr1", "chrM", "chr2") + ) + raw_vcf = self.bgzip( + tmp_path, + "raw", + RAW_VCF.replace("##contig=\n", contigs) + + "chrM\t10\t.\tA\tC\t40\tPASS\tAF=0.5\tGT\t0/1\n" + + "chr2\t10\t.\tA\tC\t40\tPASS\tAF=0.5\tGT\t0/1\n", + ) + pop_vcf = self.bgzip( + tmp_path, + "pop", + POP_VCF.replace( + "##contig=\n", + "##contig=\n" + "##contig=\n", + ), + ) + shard0 = self.merge_shard(tmp_path, raw_vcf, pop_vcf, "chr1", 0) + extra_bed = tmp_path / "extra.bed" + extra_bed.write_text("chrM\t0\t1000\n") + shard1 = tmp_path / "extra.vcf.gz" + assert ( + asyncio.run( + self.run( + cmd_bcftools_subset_remove( + shard1, raw_vcf, regions_file=extra_bed, remove="FILTER" + ) + ) + ) + == 0 + ) + shard2 = self.merge_shard(tmp_path, raw_vcf, pop_vcf, "chr2", 2) + + out_vcf = tmp_path / "transfer.vcf.gz" + concat = bcftools_concat( + out_vcf, [shard0, shard1, shard2], xargs=["--no-version"] + ) + assert asyncio.run(self.run(concat)) == 0 + assert self.contig_lines(out_vcf) == [ + "##contig=", + "##contig=", + "##contig=", + ] + query = subprocess.run( + ["bcftools", "query", "-f", "%CHROM:%POS\n", str(out_vcf)], + check=True, + capture_output=True, + text=True, + ) + assert query.stdout.splitlines()[-2:] == ["chrM:10", "chr2:10"] + + def test_extra_contig_subset_drops_filter(self, tmp_path): + raw_vcf = self.bgzip(tmp_path, "raw", RAW_VCF) + subset_bed = tmp_path / "shard.bed" + subset_bed.write_text("chr1\t150\t1000\n") + out_vcf = tmp_path / "subset.vcf.gz" + + pipeline = cmd_bcftools_subset_remove( + out_vcf, raw_vcf, regions_file=subset_bed, remove="FILTER" + ) + assert asyncio.run(self.run(pipeline)) == 0 + + query = subprocess.run( + [ + "bcftools", + "query", + "-f", + "%POS\t%QUAL\t%FILTER\n", + str(out_vcf), + ], + check=True, + capture_output=True, + text=True, + ) + assert query.stdout.splitlines() == [ + "200\t150\t.", + "300\t.\t.", + "500\t25\t.", + "600\t35\t.", + ] + assert (tmp_path / "subset.vcf.gz.tbi").exists() + + class TestTag: """`tag` disambiguates two fan-outs in one run""" diff --git a/tests/utils/test_helpers.py b/tests/utils/test_helpers.py index af59da6..8aba479 100644 --- a/tests/utils/test_helpers.py +++ b/tests/utils/test_helpers.py @@ -26,6 +26,42 @@ def create_mock_args(loglevel="WARNING"): return args +def write_bundle( + path: pathlib.Path, members: Dict[str, bytes] +) -> pathlib.Path: + """Write a model bundle: a GNU ar archive of `members`. + + As in GNU ar, names of 16 or more characters, `combine_cnv.json` + among them, go to the archive's `//` name table. + """ + + def entry(name: bytes, body: bytes) -> bytes: + header = ( + name.ljust(16) + + b"0".ljust(12) + + b"0".ljust(6) + + b"0".ljust(6) + + b"644".ljust(8) + + str(len(body)).encode().ljust(10) + + b"`\n" + ) + return header + body + (b"\n" if len(body) % 2 else b"") + + names = b"" + entries = b"" + for name, body in members.items(): + if len(name) < 16: + entries += entry(f"{name}/".encode(), body) + else: + entries += entry(f"/{len(names)}".encode(), body) + names += f"{name}/\n".encode() + data = b"!\n" + if names: + data += entry(b"//", names) + path.write_bytes(data + entries) + return path + + class MockFileSystem: """Helper class for creating mock file systems for testing""" @@ -71,13 +107,14 @@ def create_bam_header(self, path: str, readgroups: Optional[List[Dict[str, str]] def create_model_bundle(self, path: str, pipeline_type: str = "DNAscope", longread_tech: str = "HiFi", shortread_tech: str = "Illumina") -> pathlib.Path: - """Create a mock model bundle file (ar archive will be mocked)""" - # Create an empty file - the ar_load function will be mocked + """Create a mock model bundle file (an empty ar archive)""" bundle_path = pathlib.Path(self.temp_dir) / path bundle_path.parent.mkdir(parents=True, exist_ok=True) - # Just create an empty file - the actual content will be provided by mocked ar_load - bundle_path.write_bytes(b"mock_ar_archive") + # An empty ar archive: an unmocked `ar_load` reads it as a bundle + # without a bundle_info.json. Tests that need bundle members or a + # bundle_info.json mock `ar_load`. + bundle_path.write_bytes(b"!\n") self.files[path] = bundle_path return bundle_path