From d85aee50ad7db935669028b5c0471e0cc073c632 Mon Sep 17 00:00:00 2001 From: Kieran Samuk Date: Wed, 3 Jun 2026 08:36:55 -0700 Subject: [PATCH 1/2] Condition DP filters on a present genotype - Its possible to have a missing genotype with a DP value, which are interpreted as valid sites - Now the DP filter requires an intact genotype AND a passing DP format value. --- README.md | 2 +- docs/about.rst | 7 ++++--- docs/inputs.rst | 7 ++++--- src/sprite_mask/vcf.py | 28 +++++++++++++++++++++++++ tests/test_vcf.py | 47 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 84 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index b318f8e..48e5c15 100644 --- a/README.md +++ b/README.md @@ -72,7 +72,7 @@ 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`. -- A sample passes a site when `FORMAT/DP >= --min-dp` and, if supplied, `FORMAT/DP <= --max-dp`. +- A sample passes a site when `FORMAT/DP >= --min-dp`, if supplied `FORMAT/DP <= --max-dp`, and, when `FORMAT/GT` is present, the genotype is not missing. - 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. diff --git a/docs/about.rst b/docs/about.rst index 3a12cec..a977552 100644 --- a/docs/about.rst +++ b/docs/about.rst @@ -48,9 +48,10 @@ 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`` 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 +``--max-dp``. When ``FORMAT/GT`` is present, the genotype must also be +non-missing. Duplicate records at the same ``CHROM:POS`` are merged with OR +semantics per sample: if any duplicate passes the depth and genotype checks, +the sample passes that base. Duplicate records must be contiguous, as in a coordinate-sorted VCF. Sparse output diff --git a/docs/inputs.rst b/docs/inputs.rst index 771cc08..6c767a9 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -61,9 +61,10 @@ 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. 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. +rejected. A sample passes when DP is at least ``--min-dp``, if ``--max-dp`` is +supplied no greater than ``--max-dp``, and, when ``GT`` is present in +``FORMAT``, the genotype is not missing. 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. diff --git a/src/sprite_mask/vcf.py b/src/sprite_mask/vcf.py index 46c273b..61394f4 100644 --- a/src/sprite_mask/vcf.py +++ b/src/sprite_mask/vcf.py @@ -115,11 +115,13 @@ def build_population_counts_from_all_sites_vcf( current_passes = [False] * len(samples) depth_index = _depth_format_index(fields, depth_field, all_sites_vcf, line_number) + genotype_index = _optional_genotype_format_index(fields) _update_sample_passes( current_passes, fields, selected_sample_columns, depth_index, + genotype_index, threshold, max_depth, all_sites_vcf, @@ -507,6 +509,17 @@ def _optional_depth_format_index(fields: list[str], depth_field: str = "DP") -> return None +def _optional_genotype_format_index(fields: list[str]) -> int | None: + format_text = fields[8] + if format_text in {"", "."}: + return None + format_fields = format_text.split(":") + try: + return format_fields.index("GT") + except ValueError: + return None + + def _parse_info_values(info_text: str) -> dict[str, list[str]]: if info_text in {"", "."}: return {} @@ -572,6 +585,7 @@ def _update_sample_passes( fields: list[str], selected_sample_columns: list[int], depth_index: int, + genotype_index: int | None, threshold: int, max_depth: int | None, path: Path, @@ -587,6 +601,7 @@ def _update_sample_passes( passes[selected_index] = _sample_depth_passes( fields[field_index], depth_index, + genotype_index, threshold, max_depth, path, @@ -597,6 +612,7 @@ def _update_sample_passes( def _sample_depth_passes( sample_field: str, depth_index: int, + genotype_index: int | None, threshold: int, max_depth: int | None, path: Path, @@ -616,9 +632,21 @@ def _sample_depth_passes( depth = int(depth_text) except ValueError as error: raise ValueError(f"{path}:{line_number} has non-integer sample DP") from error + if genotype_index is not None and _sample_genotype_is_missing(parts, genotype_index): + return False return depth >= threshold and (max_depth is None or depth <= max_depth) +def _sample_genotype_is_missing(parts: list[str], genotype_index: int) -> bool: + if genotype_index >= len(parts): + return True + genotype_text = parts[genotype_index] + if genotype_text in {"", "."}: + return True + alleles = genotype_text.replace("|", "/").split("/") + return any(allele in {"", "."} for allele in alleles) + + def _sample_integer_format_value( sample_field: str, value_index: int, diff --git a/tests/test_vcf.py b/tests/test_vcf.py index 4914c82..24257bc 100644 --- a/tests/test_vcf.py +++ b/tests/test_vcf.py @@ -423,6 +423,53 @@ def test_build_population_counts_from_all_sites_vcf_omits_missing_and_zero_depth assert out.read_text().splitlines()[1:] == ["#chrom\tstart\tend\tpopA"] +def test_build_population_counts_from_all_sites_vcf_omits_missing_genotypes_with_depth( + tmp_path: Path, +) -> None: + samples = [ + Sample("s1", "popA"), + Sample("s2", "popA"), + Sample("s3", "popB"), + ] + vcf = tmp_path / "all_sites.vcf" + vcf.write_text( + "##fileformat=VCFv4.2\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\ts3\n" + "chr1\t1\t.\tA\t.\t.\t.\t.\tGT:DP\t./.:10\t.|.:11\t0/0:12\n" + "chr1\t2\t.\tC\t.\t.\t.\t.\tGT:DP\t0/.:10\t0|0:11\t.:12\n" + "chr1\t3\t.\tG\t.\t.\t.\t.\tGT:DP\t./.:10\t.:11\t.|.:12\n" + ) + out = tmp_path / "population_counts.bed" + + build_population_counts_from_all_sites_vcf(samples, vcf, out, threshold=5) + + assert out.read_text().splitlines()[1:] == [ + "#chrom\tstart\tend\tpopA\tpopB", + "chr1\t0\t1\t0\t1", + "chr1\t1\t2\t1\t0", + ] + + +def test_build_population_counts_from_all_sites_vcf_supports_depth_only_format( + tmp_path: Path, +) -> None: + samples = [Sample("s1", "popA")] + vcf = tmp_path / "all_sites.vcf" + vcf.write_text( + "##fileformat=VCFv4.2\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\n" + "chr1\t1\t.\tA\t.\t.\t.\t.\tDP\t7\n" + ) + out = tmp_path / "population_counts.bed" + + build_population_counts_from_all_sites_vcf(samples, vcf, out, threshold=5) + + assert out.read_text().splitlines()[1:] == [ + "#chrom\tstart\tend\tpopA", + "chr1\t0\t1\t1", + ] + + def test_build_population_counts_from_all_sites_vcf_with_no_records_writes_only_header( tmp_path: Path, ) -> None: From dd74d8a1f48de29879c7fa0aaad2f922629de4fe Mon Sep 17 00:00:00 2001 From: Kieran Samuk Date: Wed, 3 Jun 2026 11:36:18 -0700 Subject: [PATCH 2/2] expose mosdepth's include flag --- docs/arguments.rst | 4 ++++ src/sprite_mask/cli.py | 2 ++ src/sprite_mask/config.py | 1 + src/sprite_mask/mosdepth.py | 2 ++ tests/test_cli.py | 3 +++ tests/test_workflow_commands.py | 3 +++ 6 files changed, 15 insertions(+) diff --git a/docs/arguments.rst b/docs/arguments.rst index df24a43..81b35de 100644 --- a/docs/arguments.rst +++ b/docs/arguments.rst @@ -106,6 +106,10 @@ BAM/CRAM mode options **--exclude-flag INTEGER** SAM FLAG bits to exclude reads. +**--include-flag INTEGER** + SAM FLAG bits required to include reads. Passed through to + ``mosdepth --include-flag``. + **--reference PATH** FASTA reference for CRAM inputs. diff --git a/src/sprite_mask/cli.py b/src/sprite_mask/cli.py index 9db1193..bcb3a33 100644 --- a/src/sprite_mask/cli.py +++ b/src/sprite_mask/cli.py @@ -71,6 +71,7 @@ def _cmd_from_alignments(args: argparse.Namespace) -> int: min_mapq=args.min_mapq, max_dp=args.max_dp, exclude_flag=args.exclude_flag, + include_flag=args.include_flag, reference=Path(args.reference) if args.reference else None, fast_mode=args.fast_mode, keep_work=args.keep_work, @@ -187,6 +188,7 @@ def _build_from_alignments_parser(subparsers: argparse._SubParsersAction) -> Non 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("--include-flag", type=int, help="SAM FLAG bits required to include reads") p.add_argument("--reference", help="FASTA reference for CRAM inputs") p.add_argument( "--fast-mode", diff --git a/src/sprite_mask/config.py b/src/sprite_mask/config.py index b658e6e..194f389 100644 --- a/src/sprite_mask/config.py +++ b/src/sprite_mask/config.py @@ -17,6 +17,7 @@ class AlignmentRunConfig: min_mapq: int | None = None max_dp: int | None = None exclude_flag: int | None = None + include_flag: int | None = None reference: Path | None = None fast_mode: bool = False keep_work: bool = False diff --git a/src/sprite_mask/mosdepth.py b/src/sprite_mask/mosdepth.py index f322570..28eb316 100644 --- a/src/sprite_mask/mosdepth.py +++ b/src/sprite_mask/mosdepth.py @@ -78,6 +78,8 @@ def build_mosdepth_command(sample: Sample, config: AlignmentRunConfig, prefix: P command.extend(["--mapq", str(config.min_mapq)]) if config.exclude_flag is not None: command.extend(["--flag", str(config.exclude_flag)]) + if config.include_flag is not None: + command.extend(["--include-flag", str(config.include_flag)]) if config.reference is not None: command.extend(["--fasta", str(config.reference)]) diff --git a/tests/test_cli.py b/tests/test_cli.py index 70ddb6c..1a3c37e 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -158,6 +158,8 @@ def fake_run_workflow(config: object) -> WorkflowOutputs: "20", "--exclude-flag", "1796", + "--include-flag", + "2", "--reference", str(tmp_path / "ref.fa"), "--fast-mode", @@ -177,6 +179,7 @@ def fake_run_workflow(config: object) -> WorkflowOutputs: assert seen_config.variants_vcf == tmp_path / "variants.vcf.gz" assert seen_config.min_mapq == 20 assert seen_config.exclude_flag == 1796 + assert seen_config.include_flag == 2 assert seen_config.reference == tmp_path / "ref.fa" assert seen_config.fast_mode is True assert seen_config.keep_work is True diff --git a/tests/test_workflow_commands.py b/tests/test_workflow_commands.py index 6369904..932a5d7 100644 --- a/tests/test_workflow_commands.py +++ b/tests/test_workflow_commands.py @@ -20,6 +20,7 @@ def test_build_mosdepth_command_default_omits_fast_mode(tmp_path: Path) -> None: threads=4, min_mapq=20, exclude_flag=1796, + include_flag=2, reference=tmp_path / "ref.fa", ) @@ -36,6 +37,8 @@ def test_build_mosdepth_command_default_omits_fast_mode(tmp_path: Path) -> None: "20", "--flag", "1796", + "--include-flag", + "2", "--fasta", str(tmp_path / "ref.fa"), str(tmp_path / "work" / "s1.d30"),