From f083b34b3c71f0b24c47406706909dcb054a5730 Mon Sep 17 00:00:00 2001 From: Eric Talevich <52723+etal@users.noreply.github.com> Date: Wed, 14 Oct 2026 16:55:29 -0700 Subject: [PATCH] Replace snpfilter.sh with a pysam-based snpfilter.py scripts/snpfilter.sh called samtools mpileup --VCF, an output mode that samtools removed in 1.10, and it never called genotypes, so its output gave CNVkit no heterozygous SNPs. Its pre-filter also discarded every record of a sites-only VCF, the input that doc/baf.rst recommends. The replacement, installed as the console script snpfilter.py, genotypes a sample at known SNP sites with the bcftools bundled in pysam: mpileup with per-sample AD and DP and a depth cap of 8000 reads rather than bcftools' default of 250, then call -m. It needs no external binary, so it also runs in the Docker image, which never had samtools or bcftools. A matched normal BAM, if given, is genotyped jointly, and a PEDIGREE header pairs the two samples for CNVkit. Sites on chromosomes missing from the BAM or the FASTA are skipped, an input error stops the run before the output file is opened, and an output file that is also an input is refused. doc/baf.rst documents this route beside the Mutect2 germline options, including why tumor-only genotyping loses LOH and strongly imbalanced sites, and how --zygosity-freq helps. Fixes #656. Refs #364. --- Dockerfile | 3 - cnvlib/cli/snpfilter.py | 249 ++++++++++++++++++++++++++++ doc/baf.rst | 45 ++++++ doc/docker.rst | 2 +- doc/scripts.rst | 8 + pyproject.toml | 4 +- scripts/snpfilter.sh | 40 ----- test/test_snpfilter.py | 351 ++++++++++++++++++++++++++++++++++++++++ 8 files changed, 655 insertions(+), 47 deletions(-) create mode 100755 cnvlib/cli/snpfilter.py delete mode 100755 scripts/snpfilter.sh create mode 100644 test/test_snpfilter.py diff --git a/Dockerfile b/Dockerfile index e8744aa8..a0ce0beb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -43,9 +43,6 @@ RUN if [ -n "$CNVKIT_VERSION" ]; then \ # Verify installation and build matplotlib font cache RUN cnvkit.py version -# Copy helper scripts to PATH -COPY scripts/snpfilter.sh /usr/local/bin/ - # Create non-root user for security RUN useradd -m -s /bin/bash cnvkit USER cnvkit diff --git a/cnvlib/cli/snpfilter.py b/cnvlib/cli/snpfilter.py new file mode 100755 index 00000000..df2518c7 --- /dev/null +++ b/cnvlib/cli/snpfilter.py @@ -0,0 +1,249 @@ +#!/usr/bin/env python +"""Genotype a sample at known SNP sites, for allele-frequency (BAF) analysis. + +Pile up the reads at each site listed in a VCF -- for example, common SNPs from +dbSNP or the 1000 Genomes Project, or a normal sample's germline calls -- and +call genotypes there with bcftools. The output VCF carries each sample's +genotype (GT) and allele depths (AD, DP), which is what CNVkit reads to find +heterozygous germline SNPs and their B-allele frequencies. + +This is useful when the available VCF holds only somatic variants, as from +GATK Mutect2 run without germline output, and so has no heterozygous SNPs. + +Give the matched normal BAM, if there is one, with -n. Both samples are then +called jointly, and the output's PEDIGREE header pairs them so that CNVkit +judges zygosity from the normal sample, as it does for a tumor-normal VCF from +a variant caller. + +Sample names in the output are the BAM filenames without directory or +extension. Chromosome names must match between the BAM, the FASTA and the +sites VCF; sites on chromosomes missing from the BAM or the FASTA are skipped. +""" + +import argparse +import logging +import os +import sys +import tempfile +from collections import Counter +from collections.abc import Iterable, Iterator + +import pysam +import pysam.bcftools + +from skgenome import tabio + +from ..core import fbase +from ..samutil import ensure_bam_index, get_bam_chroms + +# Sites per bcftools run. pysam returns each run's output as one string, so +# a whole-genome site list is processed in blocks to bound memory use. +# Splitting the list changes no allele depth (AD, DP). At a site deeper than +# 255 reads, though, bcftools computes the genotype likelihoods from a random +# subsample of the bases, drawn from a generator that runs on from site to +# site. PL and QUAL there -- and so the genotype, at a borderline site -- +# therefore depend on which sites preceded it in the same run, and so on where +# the blocks begin. +CHUNK_SIZE = 50_000 +# Per-BAM read depth cap. bcftools' default of 250 would subsample the reads +# at a deeply sequenced site and so coarsen its allele frequency; this matches +# the default of samtools mpileup. +MAX_DEPTH = 8000 + + +def argument_parsing() -> argparse.Namespace: + AP = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + AP.add_argument("sample_bam", help="Sample (e.g. tumor) BAM or CRAM file.") + AP.add_argument( + "sites", + help="""VCF of the SNP sites to genotype. Only positions are used, so + a sites-only VCF without sample columns is fine.""", + ) + AP.add_argument( + "-f", + "--fasta", + metavar="FILENAME", + required=True, + help="""Reference genome, FASTA format: the one the reads were aligned + to. An index (.fai) is created beside it if absent.""", + ) + AP.add_argument( + "-n", + "--normal", + metavar="BAM", + help="Matched normal BAM or CRAM file, to be genotyped jointly.", + ) + AP.add_argument( + "-o", "--output", metavar="FILENAME", help="Output VCF. [Default: stdout]" + ) + return AP.parse_args() + + +def snpfilter( + sample_bam: str, + sites_fname: str, + fasta: str, + normal_bam: str | None = None, + chunk_size: int = CHUNK_SIZE, +) -> Iterator[str]: + """Genotype a sample, and optionally its matched normal, at the given sites. + + Yield the lines of the output VCF. + """ + bams = [sample_bam] if normal_bam is None else [sample_bam, normal_bam] + names = [fbase(bam) for bam in bams] + if len(set(names)) < len(names): + raise ValueError( + f"The sample and normal BAMs would both get the same sample name " + f"'{names[0]}'; rename one of the files" + ) + for bam in bams: + ensure_bam_index(bam) + bam_contigs = set.intersection(*(get_bam_chroms(bam, fasta) for bam in bams)) + with pysam.FastaFile(fasta) as fa: + fasta_contigs = set(fa.references) + contigs = bam_contigs & fasta_contigs + if not contigs: + raise ValueError( + "The BAM and the FASTA have no chromosome name in common " + f"(BAM: {', '.join(_abbreviate(sorted(bam_contigs)))}; " + f"FASTA: {', '.join(_abbreviate(sorted(fasta_contigs)))})" + ) + n_sites = n_records = 0 + with tempfile.TemporaryDirectory() as tmpdir: + regions_fname = os.path.join(tmpdir, "sites.tsv") + pileup_fname = os.path.join(tmpdir, "pileup.vcf") + positions = _read_positions(sites_fname, contigs) + for i, chunk in enumerate(_chunks(positions, chunk_size)): + n_sites += len(chunk) + with open(regions_fname, "w") as handle: + handle.writelines(f"{chrom}\t{pos}\n" for chrom, pos in chunk) + # Without split_lines, pysam's dispatcher returns standard output + # as one string + pileup: str = pysam.bcftools.mpileup( # type: ignore[assignment] + "--fasta-ref", fasta, + "--regions-file", regions_fname, + "--annotate", "FORMAT/AD,FORMAT/DP", + "--max-depth", str(MAX_DEPTH), + "--skip-indels", + "--count-orphans", + "--ignore-RG", + "--output-type", "v", + *bams, + ) # fmt: skip + with open(pileup_fname, "w") as handle: + handle.write(pileup) + # Keep sites where no sample carries an ALT allele (no --variants- + # only), so that the output covers every site that had reads. + calls: str = pysam.bcftools.call( # type: ignore[assignment] + "--multiallelic-caller", "--output-type", "v", pileup_fname + ) + for line in calls.splitlines(keepends=True): + if not line.startswith("#"): + n_records += 1 + yield line + elif i == 0: + if line.startswith("#CHROM"): + if normal_bam is not None: + yield f"##PEDIGREE=\n" + fields = line.rstrip("\n").split("\t") + line = "\t".join(fields[:9] + names) + "\n" + yield line + if n_records: + logging.info("%d of %d sites had reads to genotype", n_records, n_sites) + else: + logging.warning( + "WARNING: None of the %d sites had any reads; the VCF has no records", + n_sites, + ) + + +def _read_positions(sites_fname: str, contigs: set[str]) -> Iterator[tuple[str, int]]: + """Yield the distinct (chromosome, 1-based position) pairs of the sites. + + Sites on chromosomes outside ``contigs`` are skipped: bcftools fails on a + region it cannot find in the BAM, and where the FASTA lacks the chromosome + it piles up against a reference base of N, which the caller then drops + without notice. + Consecutive records at one position, such as the split alleles of a + multiallelic site, give a single position, since bcftools would otherwise + genotype the site once for each. + """ + skipped: Counter[str] = Counter() + n_kept = 0 + last = None + with pysam.VariantFile(sites_fname, drop_samples=True) as sites: + for record in sites: + if record.chrom not in contigs: + skipped[record.chrom] += 1 + continue + position = (record.chrom, record.pos) + if position != last: + n_kept += 1 + yield position + last = position + if skipped: + logging.info( + "Skipped %d sites on %d chromosomes not in both the BAM and the FASTA: %s", + skipped.total(), + len(skipped), + ", ".join(_abbreviate(skipped)), + ) + if not n_kept: + if not skipped: + raise ValueError(f"No sites found in {sites_fname}") + raise ValueError( + f"None of the sites in {sites_fname} is on a chromosome in both the " + "BAM and the FASTA. Check that the chromosome names match " + f"(sites VCF: {', '.join(_abbreviate(skipped))}; " + f"BAM and FASTA: {', '.join(_abbreviate(sorted(contigs)))})" + ) + + +def _abbreviate(names: Iterable[str], limit: int = 5) -> list[str]: + names = list(names) + return names if len(names) <= limit else [*names[:limit], "..."] + + +def _chunks( + positions: Iterable[tuple[str, int]], size: int +) -> Iterator[list[tuple[str, int]]]: + chunk: list[tuple[str, int]] = [] + for position in positions: + chunk.append(position) + if len(chunk) == size: + yield chunk + chunk = [] + if chunk: + yield chunk + + +def main() -> None: + logging.basicConfig(level=logging.INFO, format="%(message)s") + args = argument_parsing() + if args.output and os.path.exists(args.output): + # The inputs are still being read while the output is written. An input + # that does not exist cannot be compared, and pysam reports it instead. + for input_fname in (args.sites, args.sample_bam, args.normal, args.fasta): + if ( + input_fname + and os.path.exists(input_fname) + and os.path.samefile(args.output, input_fname) + ): + raise ValueError( + f"Output file {args.output} is also an input file ({input_fname})" + ) + lines = snpfilter(args.sample_bam, args.sites, args.fasta, args.normal) + # snpfilter() is a generator: its input checks and its first block of sites + # run at this next(), before the output file is opened. A failure in a + # later block truncates the output, as for other CNVkit commands. + first_line = next(lines) + with tabio.safe_write(args.output or sys.stdout) as handle: # type: ignore[arg-type] + handle.write(first_line) + handle.writelines(lines) + + +if __name__ == "__main__": + main() diff --git a/doc/baf.rst b/doc/baf.rst index d7e75fd8..0696d298 100644 --- a/doc/baf.rst +++ b/doc/baf.rst @@ -76,6 +76,9 @@ genotypes -- a somatic-only VCF gives CNVkit no heterozygous SNPs to use: somatic VCF contains only somatic variants and omits the ``GT`` field, so it has no usable germline SNPs. +If the caller cannot be re-run, genotype the known SNP sites directly from the +BAM files instead; see :ref:`baf-snpfilter` below. + An `example VCF `_ constructed from the 1000 Genomes samples NA12878 and NA12882 is included in @@ -83,6 +86,48 @@ CNVkit's test suite. See also: :ref:`vcfformat`. +.. _baf-snpfilter: + +Genotyping known SNP sites with ``snpfilter.py`` +```````````````````````````````````````````````` + +The script ``snpfilter.py``, installed with CNVkit, builds a VCF that CNVkit +can use for BAF analysis from the sample's BAM file and a list of known SNP +sites. It piles up the reads at each site and calls genotypes there with +bcftools -- the copy bundled with pysam, so no separate installation is +needed -- and records each sample's genotype (``GT``) and allele depths +(``AD``, ``DP``):: + + snpfilter.py Sample.bam common_snps.vcf.gz -f hg38.fa -n Normal.bam -o Sample.snps.vcf + cnvkit.py call Sample.cns -v Sample.snps.vcf -o Sample.call.cns + +- **Sites.** Any VCF of known SNP positions will do: common SNPs from dbSNP or + the 1000 Genomes Project, or a normal sample's germline calls. Only the + positions are read, so a sites-only VCF without sample columns is fine. +- **Reference.** ``-f`` is required, and must be the FASTA the reads were + aligned to. If it has no index (``.fai``), one is created beside it. +- **Chromosome names** must match between the BAM, the FASTA, and the sites + VCF (e.g. ``chr1`` versus ``1``). Sites on chromosomes absent from the BAM + or the FASTA are skipped and counted in the log; if no site remains, the + script stops with an error and writes no VCF. +- **Sample names** in the output are the BAM filenames without directory or + extension, as in CNVkit's own output files. Given a matched normal BAM + (``-n``), the script genotypes both BAMs together and declares the pair in + the output's PEDIGREE header, so :ref:`call`, :ref:`scatter` and + :ref:`segment` read it as a tumor-normal VCF without ``--sample-id`` and + ``--normal-id``. +- **Without a matched normal**, the tumor's own genotype decides which sites + are heterozygous. A site where the tumor has lost one allele entirely is + called homozygous and excluded. At high read depth, a site with strong but + incomplete allelic imbalance -- the signal BAF analysis looks for -- can be + called homozygous too (for example, an alternate-allele fraction of 0.8 at + 300 reads), because bcftools' diploid genotype model finds a 50% allele + fraction implausible there. To classify sites by allele frequency instead, + pass ``--zygosity-freq`` (``-z``) with the lowest minor-allele fraction to + keep: ``-z 0.1`` treats sites with an alternate-allele fraction from 0.1 up + to 0.9 as heterozygous. With a matched normal, the normal's genotype selects + the heterozygous sites, so neither problem arises. + Sample identification in VCFs ----------------------------- diff --git a/doc/docker.rst b/doc/docker.rst index c824ec3b..916aafd7 100644 --- a/doc/docker.rst +++ b/doc/docker.rst @@ -124,7 +124,7 @@ The production image includes: - Python 3.11+ with CNVkit and all dependencies - R with the DNAcopy package for segmentation - A Conda environment with all required packages -- Helper scripts (e.g. ``snpfilter.sh``) +- The :doc:`additional scripts ` (e.g. ``snpfilter.py``) - A non-root user (``cnvkit``) for security - A working directory at ``/data`` for mounted volumes diff --git a/doc/scripts.rst b/doc/scripts.rst index 1885568b..b49520b1 100644 --- a/doc/scripts.rst +++ b/doc/scripts.rst @@ -85,3 +85,11 @@ Additional scripts supported by skgenome.tabio, including BED and UCSC RefFlat (e.g. refFlat.txt from `UCSC Genome Bioinformatics `_). + +``snpfilter.py`` + Genotype a sample, and optionally its matched normal, at a list of known + SNP sites, producing a VCF of germline SNPs for :doc:`BAF analysis `. + Useful when the available VCF has only somatic variants, as from GATK + Mutect2 run without germline output. See :ref:`baf-snpfilter`. :: + + snpfilter.py Sample.bam common_snps.vcf.gz -f hg38.fa -n Normal.bam -o Sample.snps.vcf diff --git a/pyproject.toml b/pyproject.toml index 1de230cb..a6fe4c12 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,9 +49,7 @@ changelog = "https://github.com/etal/cnvkit/releases" "guess_baits.py" = "cnvlib.cli.guess_baits:main" "reference2targets.py" = "cnvlib.cli.reference2targets:main" "skg_convert.py" = "skgenome.cli.skg_convert:main" -# XXX TODO -- figure out whether this script should and could be included -# in the installation -#"snpfilter.sh" = "scripts/snpfilter.sh" +"snpfilter.py" = "cnvlib.cli.snpfilter:main" [build-system] # Minimum requirements for the build system to execute (PEP 508) diff --git a/scripts/snpfilter.sh b/scripts/snpfilter.sh deleted file mode 100755 index 7df365e7..00000000 --- a/scripts/snpfilter.sh +++ /dev/null @@ -1,40 +0,0 @@ -#!/bin/bash -# Call the given germline-heterozygous SNPs in a tumor sample -# Output filename is "calls.vcf" -# NB: Chromosome names must match between the BAM, VCF, and (optional) FASTA. -# Without the FASTA genome sequence, the VCF is a little weird but still usable. -# If given, FASTA should be bgzipped and indexed with 'samtools faidx'. - -if [ -z "$1" ] -then - echo "usage: $0 bam vcf [fasta]" - exit 1 -fi - -bam=$1 -vcf=$2 -fasta=$3 -# TODO make sure chroms match - -if [ -z "$fasta" ] -then - fa_opts="" -else - fa_opts="--fasta-ref $fasta" - bgzip -r $fasta - if [ ! -e "${fasta}.fai" ] - then - samtools faidx $fasta - fi -fi - -# TODO/ENH: Can bcftools emit just the sites, no samples? Should it? -bcftools view --exclude-uncalled --trim-alt-alleles --genotype ^miss \ - --output-type v --output-file _snps.vcf $vcf -grep -v '^#' _snps.vcf | cut -f1,2 > sites.txt -samtools mpileup \ - --ignore-RG --skip-indels --count-orphans --output-tags DP,AD \ - $fa_opts \ - --positions sites.txt \ - --VCF --uncompressed --output calls.vcf \ - $bam diff --git a/test/test_snpfilter.py b/test/test_snpfilter.py new file mode 100644 index 00000000..5a66fe6c --- /dev/null +++ b/test/test_snpfilter.py @@ -0,0 +1,351 @@ +#!/usr/bin/env python +"""Tests for snpfilter.py, which genotypes a sample at known SNP sites. + +Each test builds a small synthetic genome in ``tmp_path``: a random reference +sequence, a sites VCF, and tumor and normal BAMs whose allele counts at every +site are set exactly. The pipeline runs on the bcftools bundled with pysam, so +no external binary is needed. The output is then read back through +``cmdutil.load_het_snps``, the loader that the BAF-consuming commands use. +""" + +import logging +import random +import sys + +import numpy as np +import pysam +import pytest + +from cnvlib import cmdutil +from cnvlib.cli import snpfilter + +READ_LEN = 100 +# Above the 250 reads per file to which bcftools mpileup subsamples by default +DEPTH = 300 +CONTIGS = {"chr1": 8000, "chr2": 6000} +# (normal ALT fraction, tumor ALT fraction, number of sites) +DESIGN = { + "het": (0.5, 0.5, 50), + "imbalanced": (0.5, 0.8, 6), + "loh": (0.5, 1.0, 6), + "hom_ref": (0.0, 0.0, 3), + "hom_alt": (1.0, 1.0, 3), +} + + +@pytest.fixture(scope="module") +def genome(tmp_path_factory): + """Write the synthetic reference, sites VCF and BAMs; return their paths.""" + out = tmp_path_factory.mktemp("snpfilter") + rng = random.Random(0) + seqs = { + c: "".join(rng.choice("ACGT") for _ in range(n)) for c, n in CONTIGS.items() + } + fasta = out / "ref.fa" + with open(fasta, "w") as handle: + for chrom, seq in seqs.items(): + handle.write(f">{chrom}\n{seq}\n") + pysam.faidx(str(fasta)) + + # Sites 200 bp apart, so that no read reaches a neighboring site + slots = [ + (chrom, pos) + for chrom, size in CONTIGS.items() + for pos in range(150, size - 150, 200) + ] + kinds = [kind for kind, (_, _, n) in DESIGN.items() for _ in range(n)] + rng.shuffle(kinds) + sites = [] # (chrom, 0-based position, ref, alt, kind) + for (chrom, pos), kind in zip(slots, kinds, strict=True): + ref = seqs[chrom][pos] + alt = rng.choice([b for b in "ACGT" if b != ref]) + sites.append((chrom, pos, ref, alt, kind)) + + sites_vcf = out / "sites.vcf" + write_sites(sites_vcf, [(c, p + 1, r, a) for c, p, r, a, _ in sites]) + bams = { + name: write_bam(out / f"{name}.bam", seqs, sites, column) + for column, name in enumerate(("normal", "tumor")) + } + return { + "fasta": str(fasta), + "sites": str(sites_vcf), + "tumor": bams["tumor"], + "normal": bams["normal"], + "dir": out, + "kinds": {(c, p + 1): kind for c, p, _, _, kind in sites}, + } + + +def write_sites(path, records): + """Write a sites-only VCF of (chrom, 1-based position, ref, alt) records.""" + with open(path, "w") as handle: + handle.write("##fileformat=VCFv4.2\n") + for chrom, size in CONTIGS.items(): + handle.write(f"##contig=\n") + handle.write("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n") + for chrom, pos, ref, alt in records: + handle.write(f"{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\t.\t.\n") + + +def write_bam(path, seqs, sites, column): + """Write a sorted, indexed BAM with exactly the designed allele counts.""" + # Both BAMs name the same sample in their read groups, as a pipeline might + # for one patient. Without --ignore-RG, bcftools would merge the two files + # into a single sample column. + header = { + "HD": {"VN": "1.6", "SO": "coordinate"}, + "SQ": [{"SN": c, "LN": n} for c, n in CONTIGS.items()], + "RG": [{"ID": "rg1", "SM": "patient"}], + } + tids = {c: i for i, c in enumerate(CONTIGS)} + reads = [] + for chrom, pos, _ref, alt, kind in sites: + n_alt = round(DESIGN[kind][column] * DEPTH) + for i in range(DEPTH): + start = pos - (i * 37 % READ_LEN) + bases = list(seqs[chrom][start : start + READ_LEN]) + if i < n_alt: + bases[pos - start] = alt + reads.append((tids[chrom], start, "".join(bases))) + reads.sort() + with pysam.AlignmentFile(str(path), "wb", header=header) as bam: + for i, (tid, start, seq) in enumerate(reads): + read = pysam.AlignedSegment(bam.header) + read.query_name = f"r{i}" + read.reference_id = tid + read.reference_start = start + read.mapping_quality = 60 + read.cigarstring = f"{READ_LEN}M" + read.query_sequence = seq + read.query_qualities = pysam.qualitystring_to_array("I" * READ_LEN) + read.set_tag("RG", "rg1") + bam.write(read) + pysam.index(str(path)) + return str(path) + + +def run(path, *args, **kwargs): + """Run the helper and write its output VCF to ``path``.""" + with open(path, "w") as handle: + handle.writelines(snpfilter.snpfilter(*args, **kwargs)) + return str(path) + + +def records(vcf_fname): + with open(vcf_fname) as handle: + return [line for line in handle if not line.startswith("#")] + + +def sample_fields(record, *keys): + """The named FORMAT fields of each sample in a VCF record line.""" + fields = record.rstrip("\n").split("\t") + names = fields[8].split(":") + return [ + tuple(dict(zip(names, sample.split(":"), strict=True))[key] for key in keys) + for sample in fields[9:] + ] + + +def het_kinds(varr, kinds): + return sorted( + kinds[(c, s + 1)] for c, s in zip(varr.chromosome, varr.start, strict=True) + ) + + +def test_tumor_only_round_trip(genome, tmp_path, caplog): + """Germline hets come back as hets with BAF near 0.5.""" + vcf = run(tmp_path / "calls.vcf", genome["tumor"], genome["sites"], genome["fasta"]) + with open(vcf) as handle: + sample_line = next(line for line in handle if line.startswith("#CHROM")) + assert sample_line.rstrip("\n").split("\t")[9:] == ["tumor"] + + # Only load_het_snps's own warnings count, not snpfilter's log records + caplog.clear() + with caplog.at_level(logging.WARNING): + varr = cmdutil.load_het_snps(vcf) + # Enough hets for load_het_snps to check the median, and it did not object + assert len(varr) >= cmdutil._BAF_INPUT_MIN_HET + assert not caplog.records + assert np.isfinite(varr["alt_freq"]).all() + assert abs(np.median(varr["alt_freq"]) - 0.5) < 0.05 + # Genotyped alone, the tumor's imbalanced and LOH sites are called + # homozygous at this depth. Classifying by allele frequency instead + # recovers the imbalanced ones, as doc/baf.rst advises. + assert het_kinds(varr, genome["kinds"]) == ["het"] * 50 + varr = cmdutil.load_het_snps(vcf, zygosity_freq=0.1) + assert het_kinds(varr, genome["kinds"]) == ["het"] * 50 + ["imbalanced"] * 6 + + +def test_paired_round_trip(genome, tmp_path): + """A matched normal is paired automatically and recovers the LOH sites.""" + vcf = run( + tmp_path / "calls.vcf", + genome["tumor"], + genome["sites"], + genome["fasta"], + normal_bam=genome["normal"], + ) + # Each sample's own depth, not the joint INFO/DP that the reader would + # otherwise fall back to, and not subsampled + for record in records(vcf): + assert sample_fields(record, "DP") == [(str(DEPTH),)] * 2 + assert all(gt != "./." for (gt,) in sample_fields(record, "GT")) + varr = cmdutil.load_het_snps(vcf) + assert "n_zygosity" in varr + assert (varr["n_zygosity"] == 0.5).all() + assert het_kinds(varr, genome["kinds"]) == ( + ["het"] * 50 + ["imbalanced"] * 6 + ["loh"] * 6 + ) + loh = varr[varr["alt_freq"] > 0.9] + assert len(loh) == 6 + + +def test_chunking_does_not_change_allele_depths(genome, tmp_path): + """Splitting the site list into blocks changes no site, genotype or depth. + + Above 255 reads, bcftools subsamples the bases for its genotype likelihoods + with a generator whose state carries over from site to site, so PL and QUAL + can differ between the runs; at these unambiguous sites, GT cannot. + """ + args = (genome["tumor"], genome["sites"], genome["fasta"]) + whole = run(tmp_path / "whole.vcf", *args, normal_bam=genome["normal"]) + split = run( + tmp_path / "split.vcf", *args, normal_bam=genome["normal"], chunk_size=7 + ) + + def contents(vcf_fname): + with open(vcf_fname) as handle: + lines = handle.readlines() + # bcftools' command lines name temporary files + header = [line for line in lines if line.startswith("#")] + header = [line for line in header if not line.startswith("##bcftools")] + calls = [ + (line.split("\t")[:5], sample_fields(line, "GT", "AD", "DP")) + for line in lines + if not line.startswith("#") + ] + return header, calls + + assert contents(split) == contents(whole) + assert len(records(whole)) == sum(n for _, _, n in DESIGN.values()) + + +def test_sites_off_the_bam_contigs(genome, tmp_path): + """Sites on contigs the BAM lacks are skipped; repeated positions collapse.""" + lines = (genome["dir"] / "sites.vcf").read_text().splitlines(keepends=True) + meta = [line for line in lines if line.startswith("##")] + header = [line for line in lines if line.startswith("#CHROM")] + body = [line for line in lines if not line.startswith("#")] + # A second allele at the first site, as in a split multiallelic record + first = body[0].split("\t") + other_alt = next(b for b in "ACGT" if b not in first[3:5]) + second_allele = "\t".join([*first[:4], other_alt, *first[5:]]) + extra = tmp_path / "extra.vcf" + extra.write_text( + "".join( + [ + *meta, + "##contig=\n", + *header, + body[0], + second_allele, + *body[1:], + "chrUn\t500\t.\tA\tG\t.\t.\t.\n", + ] + ) + ) + plain = run( + tmp_path / "plain.vcf", genome["tumor"], genome["sites"], genome["fasta"] + ) + padded = run(tmp_path / "padded.vcf", genome["tumor"], str(extra), genome["fasta"]) + assert records(padded) == records(plain) + + +def test_no_sites_on_the_bam_contigs(genome, tmp_path): + """Mismatched chromosome naming is reported, not written as an empty VCF.""" + with open(genome["sites"]) as handle: + renamed = handle.read().replace("chr", "") + sites = tmp_path / "nochr.vcf" + sites.write_text(renamed) + with pytest.raises(ValueError, match="chromosome names"): + list(snpfilter.snpfilter(genome["tumor"], str(sites), genome["fasta"])) + + empty = tmp_path / "empty.vcf" + write_sites(empty, []) + with pytest.raises(ValueError, match="No sites"): + list(snpfilter.snpfilter(genome["tumor"], str(empty), genome["fasta"])) + + +def test_fasta_contigs(genome, tmp_path, caplog): + """Chromosomes the FASTA lacks are skipped and reported.""" + with pysam.FastaFile(genome["fasta"]) as fa: + chr1 = fa.fetch("chr1") + partial = tmp_path / "chr1.fa" + partial.write_text(f">chr1\n{chr1}\n") + with caplog.at_level(logging.INFO): + vcf = run(tmp_path / "chr1.vcf", genome["tumor"], genome["sites"], str(partial)) + n_chr1 = sum(1 for chrom, _pos in genome["kinds"] if chrom == "chr1") + assert [rec.split("\t")[0] for rec in records(vcf)] == ["chr1"] * n_chr1 + n_chr2 = len(genome["kinds"]) - n_chr1 + assert f"Skipped {n_chr2} sites on 1 chromosomes" in caplog.text + + renamed = tmp_path / "nochr.fa" + renamed.write_text(f">1\n{chr1}\n") + with pytest.raises(ValueError, match="no chromosome name in common"): + list(snpfilter.snpfilter(genome["tumor"], genome["sites"], str(renamed))) + + +def test_refused_run_leaves_files_alone(genome, tmp_path, monkeypatch): + """With -o, an input error or an input named as the output changes nothing.""" + + # A copy, since a failed guard would overwrite it, and the genome fixture is + # shared with the other tests in this module + reference = (genome["dir"] / "ref.fa").read_bytes() + fasta = tmp_path / "ref.fa" + fasta.write_bytes(reference) + + def run_main(sites, output): + argv = ["snpfilter.py", genome["tumor"], str(sites), "-f", str(fasta)] + monkeypatch.setattr(sys, "argv", [*argv, "-o", str(output)]) + snpfilter.main() + + output = tmp_path / "calls.vcf" + run_main(genome["sites"], output) + complete = output.read_text() + assert len(records(output)) == len(genome["kinds"]) + + nochr = tmp_path / "nochr.vcf" + nochr.write_text((genome["dir"] / "sites.vcf").read_text().replace("chr", "")) + with pytest.raises(ValueError, match="chromosome names"): + run_main(nochr, output) + # The earlier output, spelled differently, is a valid sites VCF that would be + # truncated while still being read. The ".." keeps the two paths textually + # distinct, so that os.path.samefile must decide; pathlib removes a ".". + with pytest.raises(ValueError, match="also an input"): + run_main(tmp_path / ".." / tmp_path.name / "calls.vcf", output) + assert output.read_text() == complete + + alias = tmp_path / "alias.fa" + alias.symlink_to(fasta) + with pytest.raises(ValueError, match="also an input"): + run_main(genome["sites"], alias) + assert fasta.read_bytes() == reference + + +def test_same_sample_name_refused(genome, tmp_path): + """Two BAMs that would get the same sample name are refused.""" + normal = tmp_path / "tumor.bam" + for suffix in ("", ".bai"): + (tmp_path / f"tumor.bam{suffix}").write_bytes( + (genome["dir"] / f"normal.bam{suffix}").read_bytes() + ) + with pytest.raises(ValueError, match="same sample name"): + list( + snpfilter.snpfilter( + genome["tumor"], + genome["sites"], + genome["fasta"], + normal_bam=str(normal), + ) + )