Refine DP filter to require intact genotype and valid DP value

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.
 

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

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.
 

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.
 

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",

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

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)])
 

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,

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

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:

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