Native methylation-aware phasing: wire the ported algorithm into make_examples#6
Open
nh13 wants to merge 4 commits into
Open
Native methylation-aware phasing: wire the ported algorithm into make_examples#6nh13 wants to merge 4 commits into
nh13 wants to merge 4 commits into
Conversation
ac32b5f to
c4b614a
Compare
0cd41f8 to
2724a87
Compare
|
Heads-up @nh13: I had to squash Side effect: this rewrote the base of your 4 open PRs (#3, #4, #6, #7), so they now show an inflated diff against Sorry for the churn, and thanks again for all the fixes. Your authorship is preserved via a |
This was referenced Jun 24, 2026
2724a87 to
a21e4d4
Compare
--parse_sam_aux_fields was set on MakeExamplesOptions but never on the SamReaderOptions used to open the BAM, and sam_reader.cc's ParseAuxFields is a no-op unless aux_field_handling == PARSE_ALL_AUX_FIELDS. So the flag did nothing: no aux fields (MM/ML base modifications, HP, ...) were ever parsed into the Read proto. Set aux_field_handling = PARSE_ALL_AUX_FIELDS on sam_opts when the flag is on; default off keeps the byte-identical baseline (WGS chr20 unchanged at 182 examples; 14/14 ctest green). Found while validating methylation-aware phasing on real PacBio HiFi 5mC data: without this, read.base_modifications() is always empty so the 5mC path can never engage. (It also gates BAM HP-tag parsing used by the phasing fallback.)
The single-sample path used vcf_candidate_importer::VariantCaller
(variant_calling.cc), which does not emit methylated reference sites
(alt='.'); only multi_sample::VariantCaller (variant_calling_multisample.cc)
does. Upstream always routes through the multi-sample VerySensitiveCaller for
every sample count — vcf_candidate_importer is reserved for the special
--variant_caller=vcf_candidate_importer import mode. So methylation-aware
phasing could never get methylated_ref_sites in the single-sample path.
Switch the single-sample path to multi_sample::VariantCaller (matching the
trio path and upstream), calling CallsFromAlleleCounts with a single-element
{sample_name -> AlleleCounter*} map and role 'sample'. This is the change that
makes methylation-aware phasing actually engage.
Validated on the chr20 HG002 WGS fixture (full chr20, b37, vs GIAB v4.2.1):
F1 0.9974 -> 0.9975, recall 0.9955 -> 0.9957, precision 0.9993 unchanged
(57 FP) — the broader candidate set recovers a few true positives the
narrower importer caller missed, with no precision cost, i.e. it is also more
upstream-faithful. On real PacBio HiFi 5mC data the methylated-ref-site path
now engages (6 methylated ref sites/region; PerformMethylationAwarePhasing
runs). 14/14 ctest green. Apple-only; nothing ports upstream.
The methylation_aware_phasing algorithm was compiled + unit-tested but never invoked by the native pipeline; --enable_methylation_aware_phasing did not exist, so the feature was unreachable. Port upstream's integration (make_examples_core.py): when the flag is set, AlleleCounter extracts 5mC levels and the VariantCaller emits methylated reference sites (alt='.') as candidates; in the single-sample phasing block those sites are split out of the SNP phasing graph, DirectPhasing runs on the SNP candidates, and PerformMethylationAwarePhasing then assigns a haplotype (Wilcoxon rank-sum vote on 5mC) to reads DirectPhasing left unphased. The refined per-read phases feed read_hp_tags / the haplotype small_model exactly as the DirectPhasing output does. The option is set on AlleleCounterOptions, VariantCallerOptions, and MakeExamplesOptions; dv_methylation_aware_phasing is now linked into dv_make_examples_lib. Scope: single-sample path only (the primary PacBio/ONT methylation use case), matching how the other phasing features are wired; the per-site p-values are dropped (upstream uses them only for --phasing_error_stats_output, not ported). Like --enable_methylation_calling, the flag is exposed at the make_examples stage and not forwarded through the `run` orchestrator. Validation: build + 14/14 ctest green (incl. the methylation_aware_phasing unit test). Flag off (default) is byte-identical — the option is false on all three option protos and the phasing block reduces to phasing the full candidate set as before. Flag on, chr20 HG002 (Illumina, no 5mC): runs clean, 182 examples (same as off) — the methylation path is correctly inert without modification data. NOTE: the active 5mC phasing path is NOT validated on aligned methylation data (no PacBio/ONT methylation BAM fixture is available); its correctness rests on the algorithm's unit tests + structural parity with upstream. Opt-in and gated. All changes are Apple-only native code.
When --enable_methylation_aware_phasing engages the 5mC vote, log the number of methylated reference sites and reads rephased per region, so the active path is observable in run logs. Fires only when the flag is on and methylated ref sites are present (no per-region noise otherwise).
c4b614a to
ec06b4b
Compare
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
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 #7 (which is stacked on #5) — review after they merge. Opt-in (
--enable_methylation_aware_phasing, default off); default behavior is byte-identical. Apple-only native code; nothing ports upstream.Depends on #7 for its prerequisites (aux-field parsing + the single-sample multi-sample VariantCaller). This PR is now just the methylation-aware phasing feature itself.
What this does
Wires the already-ported, unit-tested
methylation_aware_phasingalgorithm into the make_examples single-sample phasing block (mirroringmake_examples_core.py): methylated reference sites (alt='.') are split out of the SNP phasing graph, DirectPhasing runs on the SNP candidates, andPerformMethylationAwarePhasingthen assigns a haplotype (Wilcoxon rank-sum vote on 5mC) to reads DirectPhasing left unphased. The refined per-read phases feedread_hp_tags/ the haplotype small_model exactly as the DirectPhasing output does.AlleleCounterOptions/VariantCallerOptions/MakeExamplesOptions.dv_methylation_aware_phasingis linked intodv_make_examples_lib.Scope: single-sample path (the primary PacBio/ONT methylation use case); per-site p-values are dropped (upstream uses them only for
--phasing_error_stats_output). Like--enable_methylation_calling, the flag is exposed at the make_examples stage, not forwarded throughrun.Validation
ctestgreen (incl. the upstreammethylation_aware_phasingunit test run against the native build). Flag off (default): byte-identical.HG002-CpG-methylation-202202BAM, chr20 ~40×, streamed via a.bai-guided byte-range fetch): with Native single-sample correctness: multi-sample VariantCaller + aux-field parsing fix #7's prerequisites in place, the active path engages — 6 methylated reference sites per region,PerformMethylationAwarePhasingruns on real 5mC. (At 40× HiFi DirectPhasing already phases all reads, so the vote correctly finds 0 additional reads to rephase; the algorithm executing on real methylation data is the end-to-end proof. A sparser-SNP locus would exercise the rephase-count path.)Mm/Mltags (pre-2021 SAM convention); DeepVariant expects uppercaseMM/ML. Renaming them (semantically identical) is required for that specific file — modern PacBio/ONT outputs use the uppercase tags.