Native port feature parity: haploid calling, select_variant_types, configurable extract geometry, CoreML strides, methylation build#5
Merged
BenjaminDEMAILLE merged 22 commits intoJun 24, 2026
Conversation
…delity Several native make_examples paths diverged from the upstream Python reference, changing emitted values (not just internals): - haplotypes.cc NormalizeLog10Probs subtracted the max instead of log10sumexp, so resolved-variant GL/PL diverged from genomics_math.normalize_log10_probs whenever probability mass was spread across genotypes. Use the existing Log10SumExp helper. - make_examples_main.cc skipped zero-coverage regions before gVCF emission; upstream candidates_in_region only aborts early when gVCF is disabled, so reference blocks went missing under --gvcf. - small_model_features.cc computed is_multiallelic from the raw alt count; upstream is_multiallelic excludes the symbolic <*>, <NON_REF> and . alleles by default. - gvcf_emit.cc MED_DP used the upper-middle element; statistics.median averages the two middle values on even-length blocks. - regions.cc rejected comma-grouped coordinates that ranges.parse_literal accepts, silently dropping the region.
… EOF Both TFRecord readers (the standalone dv_tfrecord and the macOS nucleus patch) skipped the length and payload masked-CRC32C checks and treated a short payload read as a clean end-of-stream. A truncated or corrupt examples/CVO shard therefore lost records silently (missing variant calls). Read and verify both CRCs against the writer's MaskedCrc32c, and distinguish a partial read (truncation -> error) from gcount()==0 at a record boundary (clean EOF). Add F_FULLFSYNC + close-status on the writer's Close() and a sector-alignment static_assert. The native lib stays dependency-light by reporting via stderr rather than absl::log.
…ites postprocess_main.cc and call_variants_main.cc discarded the bool status of WriteRecord/Close on the output and gVCF-temp paths, so a failed final flush was reported as success (silently short VCF/gVCF). Check them and fail with a non-zero exit. Also make the site comparator a strict total order (chrom, pos, ref, alts) so equal-position records sort deterministically before MaybeResolveConflictingVariants instead of relying on std::sort's unspecified order for equal keys.
…tput MetalInception::Predict committed and waited on the deterministic-path command buffers without inspecting status/error, so a GPU fault (OOM, timeout) silently fed stale/garbage features into the classifier. Guard both det-path buffers on MTLCommandBufferStatusError and return false.
- cli.cc only matched the model_type string "ONT", so upstream's canonical ONT_R104 was misrouted to WGS defaults; normalize aliases via CanonicalModelType at every flag read site. - Germline RunAll never created --intermediate_results_dir (the trio/ somatic/pangenome dispatchers do), so a clean machine failed opaquely inside make_examples; create it with std::filesystem. - The somatic/pangenome Stage-2.5 CVO merge shelled out to an unquoted cat > file via std::system (breaks on spaces, injects); replace with a shared MergeCvoFiles byte-copy helper matching the germline path. - Validate the numeric position arg in extract_pileup_at_pos and the region ints in dump_allele_counts instead of defaulting bad input to 0. - Fix an operator-precedence out-of-bounds read in the cigar microtest.
tensor_bundle_reader mapped DT_BFLOAT16 to (np.float32, 2) and then read the 2-byte data with np.frombuffer(dtype=float32), reinterpreting pairs of bf16 values as one float32. Widen bf16 to fp32 via uint16<<16, and add a size guard asserting the decoded element count matches the tensor shape for every dtype so a mismatch fails loudly.
…rotests - dv_metal_conv_serial and dv_metal_conv_kahan PUBLIC-linked each other, a static-library cycle ld reports as "ignoring duplicate libraries". Only serial->kahan is a real edge (MetalConvKahan::Create); kahan only needs ConvDesc from the header, so drop the spurious back-edge. - Resolve libdeflate via 'brew --prefix libdeflate' + find_library instead of a hardcoded /opt/homebrew/lib/libdeflate.a. - Replace the deprecated single-arg FetchContent_Populate(libssw) with FetchContent_MakeAvailable. - Add CONFIGURE_DEPENDS to the channel/proto GLOBs. - Register the five self-contained determinism/fidelity microtests with ctest; each exits non-zero on mismatch.
…ve fs 'cmake -B build' collides with the repo's Bazel 'BUILD' file on macOS's default case-insensitive filesystem and clobbers it. README, docs, and release/build_release.sh already use 'build-macos'; align build-prereq-macos.sh and CLAUDE.md.
- third_party/tf_example was a committed symlink to /Users/benjamin/... that dangles on every other machine; the build uses the unrelated proto_tf_example cmake alias. - Untrack .claude/settings.json (it disabled local gateguard hooks) and add .claude/ to .gitignore.
…ifacts - Add curl --fail (+retry/timeout, atomic .partial->mv) to the GIAB and GLnexus downloads so a 404/HTML error page isn't consumed by the next step. - Guard the unconditional examples_chr20.tfrecord read in capture_linux_x86.sh and the shape-detection command-substitution in convert_via_docker.sh. - trap-clean the mktemp dir in diff_filter_classes.sh on the early-exit path, and add a liveness check for the detached Docker baseline in run_hg002_wg_compare.sh.
…idation - tfrecord_writer_macos Close() returned on flush failure before close(), leaking the open stream in the impls map; always close and aggregate the flush+close status. - postprocess gVCF temp path now removes the partial .variants.tmp.tfrecord on a WriteRecord/Close failure instead of leaving it on disk. - dump_allele_counts rejects non-positive 1-based coordinates and inverted ranges (s >= 1, e >= s).
The Core ML backend copied the input batch into the MLMultiArray and the prediction result out of it with flat memcpy, assuming both buffers were densely packed row-major. MLMultiArray makes no such layout guarantee: initWithShape: allocations and especially prediction-result arrays can carry padded strides, which would silently scramble the pileup tensor on copy-in and the genotype probabilities on copy-out. Check the array's strides against the packed row-major layout and keep the fast memcpy when they match; otherwise gather/scatter element-by-element honoring the actual strides. Also bounds-check the result array's element count before reading. Affects the non-default coreml backend only.
debug_metal links dv_metal_inference, which already exposes dv_weights as a PUBLIC dependency, so listing dv_weights again directly made the linker emit 'ld: warning: ignoring duplicate libraries: lib/libdv_weights.a'. Remove the redundant direct entry; dv_weights still reaches the link line transitively and debug_metal links cleanly.
The conversion README still described the abandoned tensorflow-metal voie and a TF-based flow that no longer exists: it listed a convert_metal.py and venv-metal (neither is created by setup_venvs.sh, which builds only venv-coreml/venv-mlx and hard-fails if tensorflow is importable), claimed TF 2.16.2 / tensorflow-metal / tf2onnx dependencies, a coremltools 7.2 pin, a --via-onnx fallback, and a coreml '--saved-model' invocation. Rewrite Layout, Pipeline, the venv table, and stop conditions to match the scripts as they actually exist: a two-way Core ML vs MLX bench, TF-free weight ingest via the in-repo pure-protobuf readers, coremltools 9.0+ with the direct MIL path, and convert_coreml.py's real '--bundle' CLI. Every file named in the layout now exists on disk.
extract_pileup_npy and extract_pileup_at_pos hardcoded the WGS pileup geometry (100x221x7) and rejected every other size as a 'bad image size', so they could not be used on WES/PacBio/ONT examples whose pileup tensors have different height/width/channel counts. Decode the example's image/shape int64 feature (packed or unpacked) and derive H/W/C per record, falling back to WGS 100x221x7 only when the feature is absent. The npy tool fixes geometry from the first record and rejects any later record whose shape differs (it packs the whole batch into one array); the at-pos tool sizes each matched image independently. Verified on a real WGS shard: both still emit correct 100x221x7 FP32 .npy.
The native make_examples accepted no --select_variant_types flag, so the DeepVariantOptions.select_variant_types field (proto field 30) was silently ignored: a user asking for snps-only examples got every candidate type. The make_examples_core.filter_candidates step had no native equivalent. Add the --select_variant_types flag (whitespace-separated subset of snps, indels, insertions, deletions, multi-allelics, all) to make_examples_main, validated up front so an unknown selector fails fast. Port the nucleus variant_utils type predicates faithfully — excluding the '<*>', '<NON_REF>', and '.' alleles from classification, and bi-allelic-gating snps/indels/insertions/deletions exactly as VARIANT_TYPE_SELECTORS does — and drop candidates matching no selector at both candidate-generation sites (single-sample and multi-sample). Forward the flag from the run/trio/somatic/ pangenome orchestrators in cli.cc. Validated on chr20 HG002 (20:1.0-1.2Mb): of 326 candidates (253 SNP, 46 indel, 27 multi-allelic), 'snps' keeps exactly the 253 SNPs, 'indels' the 46 indels, 'snps indels' their 299-record union, 'multi-allelics' the 27, and 'all' all 326 — matching upstream selector semantics.
deepvariant/methylation_aware_phasing.cc — a faithful port of upstream's Wilcoxon-rank-sum 5mC haplotype-voting phasing — was built by no native CMake target, so it was dead source that could silently bitrot against the native proto definitions, and its upstream unit test never ran here. Add a dv_methylation_aware_phasing static library and register the upstream methylation_aware_phasing_test.cc verbatim in the native ctest harness (it needs only protos + gtest, no heavy TF I/O). 13/13 ctest green. Verification result, documented in the build file: no native code path calls PerformMethylationAwarePhasing — the native make_examples phasing step uses DirectPhasing only and there is no native --enable_methylation_aware_phasing flag, so methylation-aware phasing is not reachable by a native user. Full pipeline integration (methylated-ref-site partitioning + a post-DirectPhasing call, mirroring make_examples_core.py) is deferred follow-up work; the library is intentionally not linked into dv_make_examples_lib yet. This commit ensures the ported algorithm stays compiling and correct in the meantime. Touches only Apple-only build files; the shared .cc/.h are unchanged.
…ions) The native port had no haploid calling: chrX/chrY in an XY sample were called diploid, so males got spurious 0/1 calls and diploid gVCF reference confidence on the sex chromosomes. Upstream's --haploid_contigs / --par_regions_bed were unimplemented. Port both halves of the upstream haploid path: - postprocess: in CombineLikelihoods, for a variant on a --haploid_contig outside the PAR, zero every heterozygous genotype probability and renormalize before deriving GT/QUAL/GQ/GL (mirror of correct_nonautosome_probabilities applied in merge_predictions). - make_examples gVCF: ReferenceConfidence drives the het likelihood to the impossible value (IMPOSSIBLE_PROBABILITY_LOG10 = 999) on haploid contigs outside the PAR, matching variant_caller.py's is_haploid branch. Shared --haploid_contigs (comma/space list) + --par_regions_bed (BED) flags are defined once in postprocess_main.cc and declared in make_examples_main.cc (same multi-call binary); cli.cc forwards both to the make_examples and postprocess workers in all four run modes. The PAR/contig gating and the het-zeroing math live in a shared header (haploid_regions.h) with a unit test (14/14 ctest green). A malformed/missing PAR BED fails fast. Validated on chr20 HG002 (treating contig 20 as haploid): 62,197 het calls -> 0 (reclassified to hom/no-call, every record preserved); a PAR BED over 20:0-30Mb keeps the 29,444 hets inside it diploid and zeroes only those outside; the gVCF reference-block het PL is driven from 147 to 9989 (ref/alt unchanged); and the default (no-flag) path is byte-consistent with the established 97,527-PASS baseline (no regression). All changes are Apple-only native code.
An adversarial review pass (each fix checked against the upstream Python) found two parity bugs in variant-affecting code plus robustness gaps: - haploid (MAJOR): the postprocess PAR exemption tested whole-variant-span overlap, but upstream ranges.RangeSet.variant_overlaps tests only the single 0-based point variant.start. A multi-base variant whose start is outside every PAR interval but whose tail reaches into one was wrongly kept diploid. IsHaploidVariant now probes [start, start+1), matching upstream. (The gVCF side already used the single-point form.) - select_variant_types (MAJOR): the single-sample path pruned candidates BEFORE DirectPhasing, but upstream phases the full candidate set in candidates_in_region and only filters afterward in process(). Pruning early drops the SNP backbone DirectPhasing relies on, changing per-read HP tags (hence the haplotype small_model's 106-feature vector and genotypes) under PacBio/ONT or --use_direct_phasing. The filter now runs after phasing, before small_model dispatch. The multi-sample/trio path keeps the filter before small_model (required so pruned types aren't emitted via the small_model CVO; the trio small_model uses no HP features) with a comment documenting the residual is_phased/PS annotation gap. - CoreML (MINOR): the copy-out bounds check rejected only a too-small output array; an unexpectedly larger shape would silently scramble probabilities. Now requires exact logical-element equality. - extract tools (MINOR): a present-but-non-3-D image/shape silently fell back to WGS geometry; now it's an error (a genuinely absent feature still falls back). The npy tool no longer writes a bogus (0,0,0,0) array on empty input. - PAR BED (MINOR): LoadParRegions now detects a gzipped BED and fails with a clear message (upstream's htslib reader accepts .bed.gz; this plain parser does not) instead of reporting confusing malformed lines. Re-validated on chr20 HG002: haploid het->0 with PAR exemption intact, select_variant_types still type-pure, extract happy-path unchanged; 14/14 ctest green. All changes remain Apple-only native code.
Applied the still-valid findings from a CodeRabbit --agent pass and a local review; skipped the ones that would diverge from the upstream reference (the campaign's goal is parity, so 'more correct than upstream' is a regression). Applied: - CorrectNonautosomeProbabilities: guard the caller's size contract ((n_alts+1)(n_alts+2)/2 entries) before indexing, so a mismatched n_alts is a no-op instead of an out-of-bounds write. Added a unit test. - extract tools: track image/shape *presence* separately from its decoded arity, so a present-but-empty/wrong-type shape is an error rather than a silent WGS fallback (a genuinely absent feature still falls back). - extract tools: bound the packed-varint inner loop by the packed-blob end (pend) rather than the whole int64_list length, and range-check each dimension (1..100000) before narrowing int64->int and guard the batch allocation against int64 overflow. - cli.cc: corrected the AppendHaploidFlags comment, which wrongly claimed the flags were not forwarded to make_examples (they are, for haploid gVCF reference confidence). Skipped, with reasons documented in code: - is_haploid NOT added to the gVCF block-grouping key: upstream variant_caller.make_gvcfs groups only on (quantized_gq, has_valid_gl); keying on ploidy would emit more blocks than upstream and break parity (and haploid sites usually fall in a different GQ bin anyway). Comment added. - CorrectNonautosomeProbabilities kept void (not bool) and the zero-mass case kept as all-zeros: upstream correct_nonautosome_probabilities never errors (returns p/(sum or 1.0)); the OOB concern is fully handled by the size guard. - VarHasInsertion/VarHasDeletion keep iterating ALL alternate_bases: upstream has_insertion/has_deletion do exactly that (not the non-excluded subset). Rebuilt, 14/14 ctest green, extract happy-path and haploid/select e2e re-validated. All changes remain Apple-only native code.
Upstream loads --par_regions_bed through htslib's BedReader, which transparently handles .bed.gz; the native loader was plaintext-only and failed fast on a gzipped file. Move LoadParRegions out of the header-only haploid_regions.h into haploid_regions.cc and read via zlib's gzgets, which reads both plaintext and gzipped/bgzipped input through one path. The header stays dependency-free (so the unit test keeps linking only absl + gtest); a small dv_haploid_regions static lib carries the zlib dependency and is linked by both make_examples and postprocess. Validated on chr20 HG002: a gzipped PAR BED and its plaintext source produce byte-identical PAR exemptions (29444 hets kept diploid inside the region, 0 outside). 14/14 ctest green.
The multi-sample/trio path ran DirectPhasing over big_candidates (the post-small_model, post-select_variant_types subset), so its is_phased/PS annotations were built from a degraded SNP phasing graph. Upstream phases the full candidate set in candidates_in_region *before* filter_candidates and before any small_model dispatch. Capture the pre-prune candidate set (only when --use_direct_phasing is on) and phase that; select_variant_types pruning and small_model dispatch still run on the filtered set, and is_phased/PS is applied to the surviving big_candidates by position. This removes the documented select_variant_types+trio phasing gap and also aligns normal trio phasing with upstream's phase-the-full-set behavior. No effect when --use_direct_phasing is off (default), and the single-sample path is untouched. Validation: build + 14/14 ctest green; single-sample chr20 make_examples unchanged (separate code path); a trio make_examples smoke (chr20 HG002 BAM as all three samples) with --use_direct_phasing + --select_variant_types=snps runs clean and emits SNP-pure child examples. is_phased/PS *correctness* on real trio data is not validated here — no DeepTrio fixture is available; the change is opt-in (--use_direct_phasing) and moves toward the upstream reference.
This was referenced Jun 23, 2026
Author
|
@BenjaminDEMAILLE — this is the base of a 3-PR stack of follow-on work on the native port:
All validated locally on macOS arm64 (clean build, 14/14 ctest, chr20 HG002 GIAB concordance; #6 also exercised on real PacBio HiFi 5mC data). If it's useful, you might reference these on google#1085 so the upstream reviewers can see the port is getting active follow-on work. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stacked on #2 — review after it merges, when this diff reduces to its unique commits. Closes feature-parity gaps found while reviewing the Apple Silicon native port. Every change is Apple-only native/conversion code (zero edits to shared
deepvariant/*.cc|*.py); each realigns the port to upstream's already-correct behavior. Build green; 14/14ctestpass.Variant-affecting parity (new features)
--haploid_contigs/--par_regions_bed, previously unimplemented). Both upstream halves are ported: postprocess zeroes heterozygous genotype probabilities and renormalizes on a haploid contig outside the PAR (mirror ofcorrect_nonautosome_probabilities), and make_examples drives the gVCF reference-confidence het likelihood to the impossible value (mirror ofvariant_caller.py'sis_haploidbranch). Shared flags are defined once in the multi-call binary and forwarded to every run mode. Validated on chr20 HG002 (treating contig 20 as haploid): 62,197 het calls → 0, a PAR BED keeps the 29,444 hets inside it diploid, the gVCF reference-block het PL goes 147 → 9989, and the default (no-flag) path is byte-consistent with the established 97,527-PASS baseline.--select_variant_typescandidate filtering (previously accepted-but-ignored). Faithful port ofmake_examples_core.filter_candidates+ the nucleusvariant_utilspredicates (<*>/<NON_REF>/.exclusion, bi-allelic gating), applied at both candidate-generation sites after DirectPhasing so the phasing graph sees the full candidate set (matching upstream's filter-in-process()ordering). chr20:snpskeeps exactly the 253 SNPs,indelsthe 46 indels,multi-allelicsthe 27,allall 326.Tooling / build / docs
extract_pileup_*now read H/W/C from the example'simage/shape(WGS fallback only when the feature is genuinely absent; a present-but-malformed shape is an error), so the profiling tools work on WES/PacBio/ONT examples instead of rejecting them.MLMultiArray.strideson copy-in/out (keeping the fastmemcpywhen packed) and requires an exact output-element count — closes a latent probability-scrambling bug in the non-default backend.methylation_aware_phasing.ccis now compiled as a library and its upstream unit test runs in ctest (it was built by no target). Pipeline integration is deferred and documented; the algorithm can no longer bitrot.debug_metalno longer emits the duplicate-libdv_weightslinker warning.tools/conversion/README.mdrewritten to match the shipped TF-free two-venv (Core ML / MLX) pipeline; the prior text referenced a nonexistentconvert_metal.py/venv-metaland dead TF/tf2onnx deps.Review
Each fix was checked against the upstream Python in an adversarial pass, then a second CodeRabbit + local review pass. Together these caught two parity bugs (haploid PAR point-vs-span overlap; select-before-phasing) and several robustness gaps, all fixed here. Three review suggestions were intentionally declined because they would diverge from the upstream reference (gVCF block-grouping key, the
has_insertion/has_deletionallele set, and erroring on zero-mass genotype probabilities) — each is documented in a code comment. New shared logic (haploid_regions.h) has unit tests.Suggested reading order: haploid (
haploid_regions.h→postprocess_main.cc→gvcf_emit.cc), thenselect_variant_typesinmake_examples_main.cc, then the tooling commits.