Add --variants-vcf support for alignments

README.md

@@ -25,7 +25,7 @@ python -m pip install "git+https://github.com/samuk-lab/sprite.git"
 ```bash
 sprite from-alignments \
   --samples samples.tsv \
-  --min-dp 10 \
+  --variants-vcf variants.vcf.gz \
   --out results \
   --threads 4 \
   --jobs 2
@@ -41,6 +41,14 @@ sample_2   popB        /path/sample_2.cram
 
 If BAM/CRAM read groups include sample names, they must match the corresponding `sample_id`.
 
+When `--variants-vcf` is provided, `sprite` uses the variants-only VCF to fill in omitted
+alignment thresholds: `--min-dp`/`--max-dp` from per-sample `FORMAT/DP` and `--min-mapq`
+from `INFO/MQ` when those fields are available. Manually supplied threshold flags take
+precedence. The same VCF is also scanned for indels, symbolic structural variants, breakends,
+and multi-nucleotide polymorphisms; those reference spans are removed from every sample pass BED,
+so they are omitted from the sparse final BED and interpreted as zero passing samples in every
+population.
+
 ### From an all-sites VCF
 
 ```bash

docs/about.rst

@@ -35,6 +35,13 @@ the passing intervals, optionally clips them to a mask BED, intersects all
 sample pass BEDs with ``bedtools multiinter``, and assembles them into a
 population count mask.
 
+An optional variants-only VCF can modify this alignment workflow. ``sprite``
+can estimate omitted depth and mapping-quality thresholds from the VCF, and
+it subtracts indel, structural-variant, breakend, and multi-nucleotide
+polymorphism spans from every sample pass BED. Because the final BED is
+sparse, these excluded spans are represented as absent intervals, meaning
+zero passing samples in every population.
+
 All-sites VCF mode
 ------------------
 

docs/arguments.rst

@@ -15,13 +15,15 @@ Commands
     Build a population count mask from a prefiltered all-sites VCF with
     per-sample ``FORMAT/DP`` values.
 
-Required for every run
-======================
+Core arguments
+==============
 
 **--min-dp INTEGER**
     Minimum depth for a sample to pass a site. Must be non-negative.
-    ``--min-dp 0`` requires ``--mask`` because no genome-wide coordinate
-    file is supplied.
+    Required unless ``from-alignments`` is run with ``--variants-vcf`` and
+    the VCF contains per-sample ``FORMAT/DP`` values for threshold estimation.
+    ``--min-dp 0`` requires ``--mask`` because no genome-wide coordinate file
+    is supplied.
 
 **--out PATH**
     Output directory for the final ``sprite.bed.gz`` and tabix index.
@@ -33,6 +35,17 @@ Input-specific arguments
     Sample metadata TSV for BAM/CRAM mode. Must provide sample ID,
     population, and alignment path columns. See :doc:`inputs`.
 
+**--variants-vcf PATH**
+    Optional variants-only VCF for BAM/CRAM mode. When present, omitted
+    ``--min-dp`` and ``--max-dp`` values are estimated from positive
+    per-sample ``FORMAT/DP`` values among selected samples and omitted
+    ``--min-mapq`` is estimated from ``INFO/MQ``. Manual threshold flags
+    override these estimates.
+    Indels, symbolic structural variants, breakends, and multi-nucleotide
+    polymorphisms are converted to exclusion intervals and removed from all
+    sample pass BEDs, so those spans are omitted from the sparse final BED
+    and interpreted as zero passing samples in every population.
+
 **--all-sites-vcf PATH**
     All-sites VCF for VCF mode. Must include per-sample ``FORMAT/DP``
     values for every sample in ``--popfile``.
@@ -82,10 +95,12 @@ BAM/CRAM mode options
     is ``--jobs × --threads``.
 
 **--min-mapq INTEGER**
-    Minimum read mapping quality.
+    Minimum read mapping quality. In BAM/CRAM mode, defaults from
+    ``--variants-vcf`` when available.
 
 **--max-dp INTEGER**
-    Maximum depth to pass a site.
+    Maximum depth to pass a site. In BAM/CRAM mode, defaults from
+    ``--variants-vcf`` when available.
 
 **--exclude-flag INTEGER**
     SAM FLAG bits to exclude reads.
@@ -108,6 +123,7 @@ BAM/CRAM mode:
    sprite from-alignments \
      --samples tests/test_data/1000g_5sample_chr20_smoke/samples.tsv \
      --min-dp 10 \
+     --variants-vcf validation/cohort.variants.vcf.gz \
      --mask tests/test_data/1000g_5sample_chr20_smoke/targets.bed \
      --out results \
      --work work \

docs/changelog.rst

@@ -9,6 +9,9 @@ Initial documented release.
 Highlights:
 
 * Build sparse population-count BEDs from BAM/CRAM cohorts.
+* Use an optional variants-only VCF in BAM/CRAM mode to estimate omitted
+  depth/MAPQ thresholds and exclude indel, structural-variant, breakend, and
+  multi-nucleotide polymorphism spans.
 * Build the same output from prefiltered all-sites VCF ``FORMAT/DP`` values.
 * Write bgzipped and tabix-indexed ``sprite.bed.gz`` output.
 * Include JSON metadata and population column headers in the output BED.

docs/examples.rst

@@ -40,6 +40,36 @@ Add ``--keep-work`` to inspect the mosdepth outputs, sample pass BEDs, and
      --work work/chr20_debug \
      --keep-work
 
+BAM/CRAM with a variants-only VCF
+=================================
+
+Use ``--variants-vcf`` when you have a variants-only VCF from the same callset
+and want ``sprite`` to estimate omitted alignment thresholds and mask
+non-SNP variant spans:
+
+.. code-block:: console
+
+   sprite from-alignments \
+     --samples tests/test_data/1000g_5sample_chr20_smoke/samples.tsv \
+     --variants-vcf validation/cohort.variants.vcf.gz \
+     --mask tests/test_data/1000g_5sample_chr20_smoke/targets.bed \
+     --out results/chr20_variants_vcf \
+     --work work/chr20_variants_vcf \
+     --threads 2 \
+     --jobs 2
+
+Manual threshold flags override VCF-derived estimates:
+
+.. code-block:: console
+
+   sprite from-alignments \
+     --samples samples.tsv \
+     --variants-vcf cohort.variants.vcf.gz \
+     --min-dp 8 \
+     --max-dp 80 \
+     --min-mapq 30 \
+     --out results/manual_thresholds
+
 All-sites VCF run
 =================
 

docs/index.rst

@@ -55,7 +55,9 @@ but works equally well on its own. Source code is on
 The tool produces the same ``sprite.bed.gz`` output from two input modes:
 
 * BAM/CRAM alignments, using ``mosdepth`` to quantize each sample and
-  ``bedtools multiinter`` to combine samples.
+  ``bedtools multiinter`` to combine samples. In this mode, an optional
+  variants-only VCF can estimate omitted thresholds and exclude non-SNP
+  variant spans from the final sparse mask.
 * A prefiltered all-sites VCF, using per-sample ``FORMAT/DP`` values directly.
 
 The output is bgzip-compressed and tabix-indexed. Large cohorts and large

docs/inputs.rst

@@ -66,6 +66,43 @@ been filtered as desired. Duplicate ``CHROM:POS`` records are merged with OR
 semantics per sample, but duplicates must be contiguous, as in a
 coordinate-sorted VCF.
 
+Variants-only VCF for BAM/CRAM mode
+===================================
+
+``--variants-vcf`` is optional in BAM/CRAM mode. It should point to a
+coordinate-sorted variants-only VCF for the same samples and reference
+coordinate system as the alignments. When sample columns are present, every
+sample ID in ``--samples`` must appear in the VCF; extra VCF samples are
+ignored for threshold estimation.
+
+``sprite`` uses this VCF in two ways.
+
+Threshold estimation
+--------------------
+
+If ``--min-dp`` or ``--max-dp`` is omitted, ``sprite`` estimates it from
+positive per-sample ``FORMAT/DP`` values at variant records. ``--min-dp``
+uses the smallest observed positive DP and ``--max-dp`` uses the largest
+observed positive DP. If ``--min-mapq`` is omitted, ``sprite`` estimates it
+from the smallest ``INFO/MQ`` value, rounded down to an integer.
+
+Any threshold supplied manually on the command line takes precedence over
+the VCF-derived estimate. If ``--min-dp`` is omitted and the VCF has no usable
+``FORMAT/DP`` values, the run is rejected.
+
+Variant exclusions
+------------------
+
+The same VCF is scanned for variant classes that should not contribute
+callable single-base denominators: indels, symbolic structural variants,
+breakends, and multi-nucleotide polymorphisms. SNP-only records are retained.
+
+Exclusion spans are emitted in BED coordinates. For symbolic structural
+variants, ``INFO/END`` is preferred when present; otherwise ``SVLEN`` is used
+when available. For ordinary sequence alleles, the reference allele length is
+used. The resulting intervals are sorted, merged, and subtracted from every
+sample pass BED before population counts are built.
+
 Mask BED
 ========
 

docs/output.rst

@@ -60,6 +60,9 @@ The ``#sprite_mask_metadata`` line is JSON. It includes:
 * ``population_sample_counts``
 * input paths such as ``samples_path``, ``popfile``, ``all_sites_vcf``, and
   ``mask_bed`` when applicable
+* ``variants_vcf``, ``threshold_sources``, and
+  ``variant_vcf_threshold_estimates`` when ``--variants-vcf`` is used in
+  BAM/CRAM mode
 
 Sparse interpretation
 =====================
@@ -68,6 +71,11 @@ The mask is sparse by design. A missing interval means zero passing samples
 in every population — not that the interval was skipped or that counts are
 unknown.
 
+When ``--variants-vcf`` excludes an indel, structural variant, breakend, or
+multi-nucleotide polymorphism span in BAM/CRAM mode, that span is removed from
+all sample pass BEDs. In the final sparse BED, this is represented the same
+way as any other all-zero region: no data row is written for the span.
+
 Intermediate files
 ==================
 
@@ -84,7 +92,14 @@ When ``--keep-work`` is set, BAM/CRAM mode retains files like:
    <sample>.d<min-dp>.mosdepth.stderr.log
    <sample>.d<min-dp>.pass.bed
    <sample>.d<min-dp>.pass.targets.bed
+   variants_vcf.excluded.raw.bed
+   variants_vcf.excluded.sorted.merged.bed
+   <sample>.d<min-dp>.pass.variants.bed
+   <sample>.d<min-dp>.pass.targets.variants.bed
    cohort.d<min-dp>.multiinter.tsv
    cohort.d<min-dp>.population_count_quantized.bed
 
+The ``variants_vcf.*`` and ``*.variants.bed`` files are present only when
+``--variants-vcf`` finds non-SNP exclusion intervals.
+
 VCF mode retains the uncompressed population count mask in the work directory.

src/sprite_mask/bedtools.py

@@ -31,6 +31,18 @@ def intersect_sort_merge(a_bed: Path, b_bed: Path, out_bed: Path) -> Path:
     return out_bed
 
 
+def subtract_sort_merge(a_bed: Path, b_bed: Path, out_bed: Path) -> Path:
+    run_pipeline(
+        [
+            ["bedtools", "subtract", "-a", str(a_bed), "-b", str(b_bed)],
+            ["bedtools", "sort", "-i", "-"],
+            ["bedtools", "merge", "-i", "-"],
+        ],
+        out_bed,
+    )
+    return out_bed
+
+
 def run_multiinter(pass_beds: Sequence[Path], names: Sequence[str], out_tsv: Path) -> Path:
     command = build_multiinter_command(pass_beds, names)
     out_tsv.parent.mkdir(parents=True, exist_ok=True)

src/sprite_mask/cli.py

@@ -67,6 +67,7 @@ def _cmd_from_alignments(args: argparse.Namespace) -> int:
         threads=args.threads,
         jobs=args.jobs,
         mask_bed=Path(args.mask) if args.mask else None,
+        variants_vcf=Path(args.variants_vcf) if args.variants_vcf else None,
         min_mapq=args.min_mapq,
         max_dp=args.max_dp,
         exclude_flag=args.exclude_flag,
@@ -117,8 +118,13 @@ def build_parser() -> argparse.ArgumentParser:
     return parser
 
 
-def _add_common_run_args(p: argparse.ArgumentParser) -> None:
-    p.add_argument("--min-dp", required=True, type=int, help="minimum depth to pass a site")
+def _add_common_run_args(p: argparse.ArgumentParser, *, min_dp_required: bool = True) -> None:
+    p.add_argument(
+        "--min-dp",
+        required=min_dp_required,
+        type=int,
+        help="minimum depth to pass a site",
+    )
     p.add_argument("--out", required=True, help="output directory")
     p.add_argument(
         "--output-prefix",
@@ -149,7 +155,7 @@ def _build_from_alignments_parser(subparsers: argparse._SubParsersAction) -> Non
         required=True,
         help="sample metadata TSV (sample_id, population, alignment)",
     )
-    _add_common_run_args(p)
+    _add_common_run_args(p, min_dp_required=False)
     p.add_argument(
         "--threads",
         type=int,
@@ -162,8 +168,23 @@ def _build_from_alignments_parser(subparsers: argparse._SubParsersAction) -> Non
         default=1,
         help="samples to process concurrently; total parallelism = --jobs × --threads",
     )
-    p.add_argument("--min-mapq", type=int, help="minimum read mapping quality")
-    p.add_argument("--max-dp", type=int, help="maximum depth to pass a site")
+    p.add_argument(
+        "--variants-vcf",
+        help=(
+            "variants-only VCF used to estimate omitted depth/MAPQ thresholds "
+            "and mask non-SNP variant spans"
+        ),
+    )
+    p.add_argument(
+        "--min-mapq",
+        type=int,
+        help="minimum read mapping quality; defaults from --variants-vcf when available",
+    )
+    p.add_argument(
+        "--max-dp",
+        type=int,
+        help="maximum depth to pass a site; defaults from --variants-vcf when available",
+    )
     p.add_argument("--exclude-flag", type=int, help="SAM FLAG bits to exclude reads")
     p.add_argument("--reference", help="FASTA reference for CRAM inputs")
     p.add_argument(

src/sprite_mask/config.py

@@ -7,12 +7,13 @@
 @dataclass(frozen=True)
 class AlignmentRunConfig:
     samples_path: Path
-    min_dp: int
+    min_dp: int | None
     out_dir: Path
     work_dir: Path | None = None
     threads: int = 1
     jobs: int = 1
     mask_bed: Path | None = None
+    variants_vcf: Path | None = None
     min_mapq: int | None = None
     max_dp: int | None = None
     exclude_flag: int | None = None

src/sprite_mask/mosdepth.py

@@ -9,6 +9,8 @@
 
 
 def run_mosdepth(sample: Sample, config: AlignmentRunConfig) -> MosdepthOutputs:
+    if config.min_dp is None:
+        raise ValueError("--min-dp is required before running mosdepth")
     prefix = config.resolved_work_dir / f"{sample.sample_id}.d{config.min_dp}"
     outputs = mosdepth_outputs_for_prefix(prefix)
     outputs.stderr_log.parent.mkdir(parents=True, exist_ok=True)
@@ -54,6 +56,8 @@ def run_mosdepth(sample: Sample, config: AlignmentRunConfig) -> MosdepthOutputs:
 def build_mosdepth_command(sample: Sample, config: AlignmentRunConfig, prefix: Path) -> list[str]:
     if sample.alignment is None:
         raise ValueError(f"alignment for sample {sample.sample_id!r} is required")
+    if config.min_dp is None:
+        raise ValueError("--min-dp is required before running mosdepth")
 
     quantize = (
         f"0:{config.min_dp}:{config.max_dp}:"

src/sprite_mask/validation.py

@@ -41,6 +41,13 @@ def validate_vcf_inputs(all_sites_vcf: Path, popfile_path: Path) -> None:
         raise ValueError(f"--popfile is a directory: {popfile_path}")
 
 
+def validate_variants_vcf_input(variants_vcf: Path) -> None:
+    if not variants_vcf.exists():
+        raise ValueError(f"--variants-vcf does not exist: {variants_vcf}")
+    if variants_vcf.is_dir():
+        raise ValueError(f"--variants-vcf is a directory: {variants_vcf}")
+
+
 def validate_alignment_sample_headers(samples: list[Sample]) -> None:
     for sample in samples:
         if sample.alignment is None: