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), + ) + )