From c8ef1ffd10330ee35981b10bdf5ead8bd56c1d7b Mon Sep 17 00:00:00 2001 From: Kieran Samuk Date: Tue, 2 Jun 2026 21:07:28 -0700 Subject: [PATCH] allow max-dp flag in from-vcf mode --- README.md | 4 +++- docs/about.rst | 8 +++++--- docs/arguments.rst | 6 ++++-- docs/examples.rst | 1 + docs/inputs.rst | 9 +++++---- src/sprite_mask/cli.py | 6 ++++++ src/sprite_mask/config.py | 1 + src/sprite_mask/validation.py | 8 +++++++- src/sprite_mask/vcf.py | 8 +++++++- src/sprite_mask/workflow.py | 11 +++++++++-- tests/test_cli.py | 3 +++ tests/test_validation.py | 7 +++++++ tests/test_vcf.py | 33 +++++++++++++++++++++++++++++++++ tests/test_workflow_commands.py | 3 +++ tests/test_workflow_unit.py | 10 +++++++--- 15 files changed, 101 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index 45580e5..b318f8e 100644 --- a/README.md +++ b/README.md @@ -56,6 +56,7 @@ sprite from-vcf \ --all-sites-vcf all_sites.vcf.gz \ --popfile populations.tsv \ --min-dp 10 \ + --max-dp 80 \ --out results ``` @@ -71,7 +72,8 @@ A few things to know about VCF mode: - Every sample in the population file must appear in the VCF. VCF samples absent from the population file are ignored (with a warning). - Records may carry any `FILTER` value — the input is assumed to have been filtered as desired before running `sprite`. -- At duplicate `CHROM:POS` records, a sample passes a site if any duplicate has `FORMAT/DP ≥ --min-dp`. Duplicates must be contiguous, as in a coordinate-sorted VCF. +- A sample passes a site when `FORMAT/DP >= --min-dp` and, if supplied, `FORMAT/DP <= --max-dp`. +- At duplicate `CHROM:POS` records, a sample passes a site if any duplicate passes the depth thresholds. Duplicates must be contiguous, as in a coordinate-sorted VCF. - By default, all record types are used (SNPs, indels, symbolic alleles, invariant sites). Pass `--snps-only` to exclude indel sites while retaining invariant sites. ## Output diff --git a/docs/about.rst b/docs/about.rst index 836fec3..3a12cec 100644 --- a/docs/about.rst +++ b/docs/about.rst @@ -47,9 +47,11 @@ All-sites VCF mode In VCF mode, ``sprite`` reads ``FORMAT/DP`` values directly from an all-sites VCF. A sample passes a base when its DP value is greater than or equal to -``--min-dp``. Duplicate records at the same ``CHROM:POS`` are merged with OR -semantics per sample: if any duplicate passes, the sample passes that base. -Duplicate records must be contiguous, as in a coordinate-sorted VCF. +``--min-dp`` and, if ``--max-dp`` is supplied, less than or equal to +``--max-dp``. Duplicate records at the same ``CHROM:POS`` are merged with OR +semantics per sample: if any duplicate passes the depth thresholds, the sample +passes that base. Duplicate records must be contiguous, as in a +coordinate-sorted VCF. Sparse output ============= diff --git a/docs/arguments.rst b/docs/arguments.rst index 52dc386..df24a43 100644 --- a/docs/arguments.rst +++ b/docs/arguments.rst @@ -99,8 +99,9 @@ BAM/CRAM mode options ``--variants-vcf`` when available. **--max-dp INTEGER** - Maximum depth to pass a site. In BAM/CRAM mode, defaults from - ``--variants-vcf`` when available. + Maximum depth to pass a site. In VCF mode, omitted means no upper + depth bound. In BAM/CRAM mode, defaults from ``--variants-vcf`` when + available. **--exclude-flag INTEGER** SAM FLAG bits to exclude reads. @@ -139,6 +140,7 @@ All-sites VCF mode: --all-sites-vcf validation/cohort.all_sites.vcf.gz \ --popfile validation/sample_populations.tsv \ --min-dp 10 \ + --max-dp 80 \ --mask validation/targets.bed \ --out results \ --work work diff --git a/docs/examples.rst b/docs/examples.rst index cbcef36..f372271 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -81,6 +81,7 @@ When you have a prefiltered all-sites VCF with per-sample DP values: --all-sites-vcf validation/1000g_20sample_highcov_4chrom_subset/1000g_20samples_highcov_4chroms.all_sites.bam_call.trim_alt.vcf.gz \ --popfile validation/1000g_20sample_highcov_4chrom_subset/sample_populations.tsv \ --min-dp 10 \ + --max-dp 80 \ --mask validation/1000g_20sample_highcov_4chrom_subset/targets.bed \ --out results/vcf_mode \ --work work/vcf_mode diff --git a/docs/inputs.rst b/docs/inputs.rst index 47f6f8f..771cc08 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -61,10 +61,11 @@ absent from ``--popfile`` produce a warning and are ignored. For each record, ``sprite`` reads the ``DP`` field from the sample's ``FORMAT`` value. Missing DP values do not pass. Non-integer DP values are -rejected. Records may carry any ``FILTER`` value; the file is assumed to have -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. +rejected. A sample passes when DP is at least ``--min-dp`` and, if +``--max-dp`` is supplied, no greater than ``--max-dp``. Records may carry any +``FILTER`` value; the file is assumed to have 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 =================================== diff --git a/src/sprite_mask/cli.py b/src/sprite_mask/cli.py index b221d80..9db1193 100644 --- a/src/sprite_mask/cli.py +++ b/src/sprite_mask/cli.py @@ -87,6 +87,7 @@ def _cmd_from_vcf(args: argparse.Namespace) -> int: popfile_path=Path(args.popfile), min_dp=args.min_dp, out_dir=Path(args.out), + max_dp=args.max_dp, output_prefix=args.output_prefix, work_dir=Path(args.work) if args.work else None, mask_bed=Path(args.mask) if args.mask else None, @@ -210,6 +211,11 @@ def _build_from_vcf_parser(subparsers: argparse._SubParsersAction) -> None: # t ) p.add_argument("--popfile", required=True, help="sample/population TSV (sample_id, population)") _add_common_run_args(p) + p.add_argument( + "--max-dp", + type=int, + help="maximum depth to pass a site", + ) p.add_argument( "--snps-only", action="store_true", diff --git a/src/sprite_mask/config.py b/src/sprite_mask/config.py index c2c96dd..b658e6e 100644 --- a/src/sprite_mask/config.py +++ b/src/sprite_mask/config.py @@ -35,6 +35,7 @@ class VcfRunConfig: popfile_path: Path min_dp: int out_dir: Path + max_dp: int | None = None work_dir: Path | None = None mask_bed: Path | None = None keep_work: bool = False diff --git a/src/sprite_mask/validation.py b/src/sprite_mask/validation.py index 51e22b6..f428d12 100644 --- a/src/sprite_mask/validation.py +++ b/src/sprite_mask/validation.py @@ -11,11 +11,17 @@ def validate_threshold( threshold: int, *, + max_threshold: int | None = None, targets_bed: Path | None = None, all_sites_vcf: Path | None = None, ) -> None: if threshold < 0: - raise ValueError("--threshold must be a non-negative integer") + raise ValueError("--min-dp must be a non-negative integer") + if max_threshold is not None: + if max_threshold < 0: + raise ValueError("--max-dp must be a non-negative integer") + if max_threshold < threshold: + raise ValueError("--max-dp must be greater than or equal to --min-dp") if threshold == 0 and targets_bed is None and all_sites_vcf is None: raise ValueError("--min-dp 0 requires --mask because no genome.txt is required") diff --git a/src/sprite_mask/vcf.py b/src/sprite_mask/vcf.py index 3799057..46c273b 100644 --- a/src/sprite_mask/vcf.py +++ b/src/sprite_mask/vcf.py @@ -32,6 +32,7 @@ def build_population_counts_from_all_sites_vcf( output_bed: Path, *, threshold: int, + max_depth: int | None = None, targets_bed: Path | None = None, depth_field: str = "DP", snps_only: bool = False, @@ -120,6 +121,7 @@ def build_population_counts_from_all_sites_vcf( selected_sample_columns, depth_index, threshold, + max_depth, all_sites_vcf, line_number, ) @@ -571,6 +573,7 @@ def _update_sample_passes( selected_sample_columns: list[int], depth_index: int, threshold: int, + max_depth: int | None, path: Path, line_number: int, ) -> None: @@ -585,6 +588,7 @@ def _update_sample_passes( fields[field_index], depth_index, threshold, + max_depth, path, line_number, ) @@ -594,6 +598,7 @@ def _sample_depth_passes( sample_field: str, depth_index: int, threshold: int, + max_depth: int | None, path: Path, line_number: int, ) -> bool: @@ -608,9 +613,10 @@ def _sample_depth_passes( return False try: - return int(depth_text) >= threshold + depth = int(depth_text) except ValueError as error: raise ValueError(f"{path}:{line_number} has non-integer sample DP") from error + return depth >= threshold and (max_depth is None or depth <= max_depth) def _sample_integer_format_value( diff --git a/src/sprite_mask/workflow.py b/src/sprite_mask/workflow.py index 6931615..eba9ca8 100644 --- a/src/sprite_mask/workflow.py +++ b/src/sprite_mask/workflow.py @@ -51,6 +51,7 @@ def run_workflow(config: RunConfig) -> WorkflowOutputs: logger.info("Analysis validation: checking threshold and VCF input paths") validate_threshold( config.min_dp, + max_threshold=config.max_dp, targets_bed=config.mask_bed, all_sites_vcf=config.all_sites_vcf, ) @@ -66,7 +67,11 @@ def run_workflow(config: RunConfig) -> WorkflowOutputs: config, alignment_metadata = _resolve_alignment_thresholds(config, samples) logger.info("Analysis validation: checking threshold") assert config.min_dp is not None - validate_threshold(config.min_dp, targets_bed=config.mask_bed) + validate_threshold( + config.min_dp, + max_threshold=config.max_dp, + targets_bed=config.mask_bed, + ) _log_sample_summary(samples) required_tools = _required_tools(config) @@ -208,8 +213,9 @@ def _build_from_all_sites_vcf( config: VcfRunConfig, generated_work_files: list[Path], ) -> None: - metadata = { + metadata: dict[str, object] = { "min_dp": config.min_dp, + "max_dp": config.max_dp, "sample_count": len(samples), "popfile": str(config.popfile_path), "all_sites_vcf": str(config.all_sites_vcf), @@ -231,6 +237,7 @@ def _build_from_all_sites_vcf( config.all_sites_vcf, population_count_bed, threshold=config.min_dp, + max_depth=config.max_dp, targets_bed=config.mask_bed, snps_only=config.snps_only, metadata=metadata, diff --git a/tests/test_cli.py b/tests/test_cli.py index 5326dda..70ddb6c 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -283,6 +283,8 @@ def fake_run_workflow(config: object) -> WorkflowOutputs: str(tmp_path / "popfile.tsv"), "--min-dp", "30", + "--max-dp", + "80", "--out", str(tmp_path / "out"), "--snps-only", @@ -292,6 +294,7 @@ def fake_run_workflow(config: object) -> WorkflowOutputs: assert status == 0 assert seen_config is not None assert seen_config.snps_only is True + assert seen_config.max_dp == 80 def test_main_reports_subprocess_errors( diff --git a/tests/test_validation.py b/tests/test_validation.py index 01ba3e9..c68877e 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -21,6 +21,13 @@ def test_validate_threshold_rejects_negative() -> None: with pytest.raises(ValueError, match="non-negative"): validate_threshold(-1) + with pytest.raises(ValueError, match="non-negative"): + validate_threshold(1, max_threshold=-1) + + +def test_validate_threshold_rejects_max_below_min() -> None: + with pytest.raises(ValueError, match="greater than or equal"): + validate_threshold(10, max_threshold=9) def test_validate_threshold_zero_requires_targets_or_vcf(tmp_path: Path) -> None: diff --git a/tests/test_vcf.py b/tests/test_vcf.py index 5453cb6..4914c82 100644 --- a/tests/test_vcf.py +++ b/tests/test_vcf.py @@ -63,6 +63,39 @@ def test_build_population_counts_from_all_sites_vcf_merges_site_counts( ] +def test_build_population_counts_from_all_sites_vcf_respects_max_depth( + tmp_path: Path, +) -> None: + samples = [ + Sample("s1", "popA"), + Sample("s2", "popA"), + ] + vcf = tmp_path / "all_sites.vcf" + vcf.write_text( + "##fileformat=VCFv4.2\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\n" + "chr1\t1\t.\tA\t.\t.\t.\t.\tGT:DP\t0/0:5\t0/0:6\n" + "chr1\t2\t.\tC\t.\t.\t.\t.\tGT:DP\t0/0:8\t0/0:6\n" + "chr1\t3\t.\tG\t.\t.\t.\t.\tGT:DP\t0/0:9\t0/0:9\n" + "chr1\t3\t.\tG\tT\t.\t.\t.\tGT:DP\t0/1:7\t0/1:4\n" + ) + out = tmp_path / "population_counts.bed" + + build_population_counts_from_all_sites_vcf( + samples, + vcf, + out, + threshold=5, + max_depth=7, + ) + + assert out.read_text().splitlines()[1:] == [ + "#chrom\tstart\tend\tpopA", + "chr1\t0\t1\t2", + "chr1\t1\t3\t1", + ] + + def test_build_population_counts_from_all_sites_vcf_requires_popfile_samples( tmp_path: Path, ) -> None: diff --git a/tests/test_workflow_commands.py b/tests/test_workflow_commands.py index 6797bc4..6369904 100644 --- a/tests/test_workflow_commands.py +++ b/tests/test_workflow_commands.py @@ -188,6 +188,8 @@ def test_cli_parser_from_vcf_subcommand(tmp_path: Path) -> None: str(tmp_path / "popfile.tsv"), "--min-dp", "30", + "--max-dp", + "80", "--out", str(tmp_path / "out"), "--snps-only", @@ -197,4 +199,5 @@ def test_cli_parser_from_vcf_subcommand(tmp_path: Path) -> None: assert args.all_sites_vcf == str(tmp_path / "all_sites.vcf.gz") assert args.popfile == str(tmp_path / "popfile.tsv") assert args.min_dp == 30 + assert args.max_dp == 80 assert args.snps_only is True diff --git a/tests/test_workflow_unit.py b/tests/test_workflow_unit.py index 4f9a472..b69e635 100644 --- a/tests/test_workflow_unit.py +++ b/tests/test_workflow_unit.py @@ -661,6 +661,7 @@ def test_build_from_all_sites_vcf_writes_population_bed_and_metadata( popfile_path=popfile, min_dp=5, out_dir=tmp_path / "out", + max_dp=12, work_dir=tmp_path / "work", mask_bed=targets, ) @@ -673,6 +674,7 @@ def fake_build( output_bed: Path, *, threshold: int, + max_depth: int | None, targets_bed: Path | None, snps_only: bool, metadata: dict[str, object], @@ -682,6 +684,7 @@ def fake_build( all_sites_vcf, output_bed, threshold, + max_depth, targets_bed, snps_only, metadata, @@ -704,10 +707,11 @@ def fake_sort(in_bed: Path, out_bed_gz: Path) -> None: output_bed = tmp_path / "work" / "cohort.d5.population_count_quantized.bed" assert generated == [output_bed] build_call = calls["build"] - assert build_call[:5] == (samples, vcf, output_bed, 5, targets) - assert build_call[5] is False - assert build_call[6] == { + assert build_call[:6] == (samples, vcf, output_bed, 5, 12, targets) + assert build_call[6] is False + assert build_call[7] == { "min_dp": 5, + "max_dp": 12, "sample_count": 1, "popfile": str(popfile), "all_sites_vcf": str(vcf),