Add --Methylation mode for bisulfite/EM-seq contamination estimation - #90
Open
tfenne wants to merge 6 commits into
Open
Add --Methylation mode for bisulfite/EM-seq contamination estimation#90tfenne wants to merge 6 commits into
tfenne wants to merge 6 commits into
Conversation
SimplePileup() filled baseInfo and qualInfo in two loops kept in sync only by an identical min-baseQ test. pileup_seq() emits no base for a read with a deletion at the marker (is_del), but the quality loop pushed unconditionally, so every base after such a read at that marker paired with the wrong quality. Merge into one loop; pileup_seq() now returns whether it pushed a base so the quality follows only when it did. Prerequisite for methylation mode, which adds a second per-observation drop over these same vectors.
In WGBS/EM-seq, unmethylated C reads as T (and G as A on the opposite strand), so raw pileups disagree with the reference in ways unrelated to genotype and would inflate contamination estimates. New MethylationModel provides, as small pure/testable functions: the per-read conversion strand (aligner tag YD/ZS/XG, else a proper-pair FLAG inference reading read 1's orientation), and a usability predicate that keeps an observation only when it maps unambiguously to an allele on that strand -- the no-collapse rule (A/T on both strands, other transversions and transitions on one, C/G on neither). --Methylation wires this in: SimplePileup drops bisulfite-ambiguous observations and reads with no reliable strand; supplementary alignments are excluded; and MAPQ capping (sam_cap_mapq) and BAQ are disabled, since both score against the unconverted reference and were measured to hard-reject ~99% of EM-seq reads. --PileupFile warns (strand cannot be recovered from pileup text). Non-methyl behavior is unchanged. Exhaustive unit test (TestMethylationModel, 33 checks); regression 17/17. Validated by an in-silico titration on EM-seq NA12878+NA24385 recovering known contamination 0-10% (e.g. 5%->4.9%, 10%->9.5%).
MethylateTestBam synthesizes a directional methyl-seq BAM from resource/test/test.bam at test time (no committed test data; strand rule written independently of MethylationModel), and a Run+Diff CTest pair checks VerifyBamID --Methylation against a committed golden Ancestry file. Also add an explicit <cstdint> include to the unit test rather than relying on transitive inclusion.
Adds the --Methylation option entry and a Methylation-sequencing section covering the per-class strand rule, the conversion-efficiency independence, MAPQ/BAQ and strand-reliability handling, and the directional-only / --PileupFile limitations.
… test converter - SimplePileupViewer: only copy the panel allele table into bedTable when methylation mode is on (it is unused otherwise), avoiding a panel-sized copy on ordinary runs. - MethylateTestBam: after rewriting SEQ, drop the now-inconsistent MD tag; NM is left as-is since a methyl-aware aligner would not count the introduced C>T/G>A conversions as mismatches (VerifyBamID uses neither).
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds an opt-in --Methylation mode to support contamination estimation on directional bisulfite/WGBS and EM-seq data by making pileup observation selection strand-aware (dropping bisulfite-ambiguous observations) and disabling read-scoring features (MAPQ adjustment/BAQ) that are incompatible with converted reads.
Changes:
- Introduces
MethylationModel.{h,cpp}to infer conversion strand (tags or FLAG) and decide whether a marker is usable on that strand. - Updates
SimplePileupViewerto (a) keep base/quality vectors synchronized (fixing a deletion-related desync) and (b) apply methylation-mode read/observation filtering using the marker allele table. - Adds unit + end-to-end tests and documents the new mode in
README.md.
Reviewed changes
Copilot reviewed 11 out of 12 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| TestMethylationModel.cpp | New unit tests covering strand inference, tag parsing, and usability truth-table. |
| SimplePileupViewer.h | Adds methylation flag and updates constructor to pass marker allele table for methylation filtering. |
| SimplePileupViewer.cpp | Implements strand gating in read callback and per-marker observation filtering; fixes base/qual desync on deletions. |
| resource/test/expected/result.methyl.Ancestry | Golden output for the new end-to-end methylation test. |
| README.md | Documents --Methylation, its behavior, and limitations. |
| MethylationModel.h | Declares conversion-strand inference and strand/allele usability functions. |
| MethylationModel.cpp | Implements tag/FLAG parsing for conversion strand and strand-based usability rule. |
| MethylateTestBam.cpp | Test-only helper to synthesize a converted BAM for end-to-end testing. |
| main.cpp | Adds --Methylation CLI flag and disables MAPQ adjustment/BAQ + filters supplementary alignments in this mode. |
| ContaminationEstimator.cpp | Passes marker allele table into SimplePileupViewer for methylation mode. |
| CMakeLists.txt | Builds new helper/test executables and wires unit + end-to-end methylation tests into CTest. |
| .gitignore | Ignores new test executables and synthesized methylation BAM/index. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+189
to
+193
| if (ma->conf->methylation && | ||
| conversionStrandOf(b) == ConversionStrand::Unknown) { | ||
| skip = 1; | ||
| continue; | ||
| } |
Comment on lines
+478
to
+482
| // Methylation mode: skip observations that cannot be | ||
| // unambiguously assigned to an allele on this read's | ||
| // conversion strand (e.g. a T at a C/T marker on a | ||
| // C->T-converted read). mRef==0 means the marker's | ||
| // alleles were not found, which fails closed. |
Comment on lines
+37
to
+41
| // Whether an observation from a read in channel `cs` unambiguously identifies | ||
| // one of the marker's two alleles. An observation is unusable exactly when the | ||
| // channel's conversion could turn one allele into the other: | ||
| // - Ct is ambiguous whenever C is an allele (a T could be a converted C); | ||
| // - Ga is ambiguous whenever G is an allele (an A could be a converted G). |
Comment on lines
+173
to
+175
| add_test(NAME testMethylationDiff | ||
| WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} | ||
| COMMAND sh -c "diff resource/test/expected/result.methyl.Ancestry result.Methylation.Ancestry && rm result.Methylation.Ancestry result.Methylation.selfSM") |
- mplp_func: check BED overlap before the methylation strand gate so aux tags are only parsed for reads that cover a marker (no behavior change; both filters just drop the read). - MethylationModel.h: correct observationUsable's doc -- a marker is unusable on a channel whenever one allele IS the convertible base (C on Ct, G on Ga), independent of the other allele; e.g. A/C is unusable on Ct. - SimplePileupViewer.cpp: clarify that the per-observation filter decides per (marker alleles, strand), not per observed base. - CMakeLists: remove the generated test.methyl.bam(.bai) after the end-to-end test, matching the other golden tests' cleanup.
Contributor
Author
|
Thanks for the review! I've pushed a commit that addresses all four comments. |
Owner
|
@tfenne if I understand correctly, the "C/G" SNP dropped all the reads, meaning the marker is completely dropped, is that correct?
Otherwise, the code looks good to me, thank you very much! |
Griffan
approved these changes
Aug 3, 2026
Contributor
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 11 out of 12 changed files in this pull request and generated 1 comment.
Suppressed comments (1)
MethylationModel.h:14
- This header comment references docs/methylation-design.md, but there is no docs/ directory in the repository. Either add that document or update the reference to an existing location (e.g., the README methylation section) so readers aren’t sent to a dead path.
// See docs/methylation-design.md for the full model and derivation.
| // quality for them would shift every later quality at this | ||
| // marker onto the wrong base. | ||
| for (j = 0; j < n_plp[i]; ++j) {//each covered read in ith bam file | ||
| const bam_pileup1_t *p = plp[i] + j; |
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.
Add
--Methylationmode for bisulfite / EM-seq contamination estimationMotivation
VerifyBamID2 cannot currently be used on whole-genome bisulfite (WGBS) or enzymatic methyl-seq (EM-seq) data. In these libraries unmethylated cytosines are read as
T(and, on reads from the opposite strand, unmethylatedGpositions read asA), so a regular pileup disagrees with the reference at a large fraction of positions for reasons unrelated to genotype. Left uncorrected these conversions are counted as alternate-allele or error observations and inflate the contamination estimate — on our EM-seq test data an unfiltered run reported ~0.20 FREEMIX on an in-silico-clean sample.No published tool estimates contamination from WGBS/EM-seq; the closest prior art (a GRAIL patent, and the Xu et al. 2023 sample-tagging work) keeps only A/T SNPs, which on the shipped panels is 1,398 markers for 1000g and zero for HGDP (HGDP excludes palindromic SNPs). This PR adds a strand-aware mode that keeps essentially all markers.
What it does
--Methylationuses, at each marker, only the observations that still identify an allele unambiguously on the read's bisulfite-conversion strand, and drops the rest (the "no-collapse" rule):Because ambiguous observations are dropped rather than disambiguated, no methylation-rate or conversion-rate parameter enters the model — the estimate is independent of methylation level and of bisulfite/enzymatic conversion efficiency, and CpG context is irrelevant. Selection depends only on the read's conversion strand and the marker's alleles, never on the sample genotype, so the selection step itself adds no bias; its cost is reduced effective depth (~half the observations at transition and most transversion markers). (A small reference-mapping underestimate does appear at high contamination — see Validation.)
The mode also adapts read handling to converted reads:
--adjust-MQ) and BAQ are disabled. Both score reads against the unconverted reference, so an EM-seq read's C→T conversions look like dense mismatches. Measured on EM-seq NA12878, the default MAPQ cap (sam_cap_mapq) hard-rejects ~99% of usable reads (70 vs 12,118 base observations with the cap off). Both are htslib functions we call rather than vendor, so the mode disables them rather than forking htslib to make them conversion-aware.YD— accepting both bwa-meth's:Z:and BISCUIT's:A:—ZS, orXG) when present, else a proper-pair FLAG inference (paired reads must be properly paired; single-end reads use their own strand). Reads with no reliable strand are skipped, since a misassigned strand turns converted bases into apparent alternate alleles and inflates the estimate.Design
The conversion logic is isolated in a new, dependency-light
MethylationModel.h/.cppas small pure functions over tiny domains —conversionStrandOf(bam1_t*)/strandFromFlag(uint16_t)andobservationUsable(ref, alt, strand)— so it is exhaustively unit-testable.SimplePileupViewerconsumes them:mplp_funcgates reads on a reliable strand, and the (single, post-desync-fix) pileup loop drops bisulfite-ambiguous observations using the marker allele table passed in fromContaminationEstimator. The likelihood, optimizer, and SVD/ancestry model are untouched.Scope and limitations
--help.--PileupFilecannot carry read-pair/tag information, so--Methylation --PileupFileis allowed but warns: the input must already contain only methylation-safe observations (e.g. from a prior--Methylation --BamFile … --OutputPileuprun).Validation
Unit test (
TestMethylationModel, 33 checks): the full 6-class × 2-strand usability truth table, exhaustive FLAG combinations for strand inference, and tag parsing/precedence (YD:Z,YD:A,YD:A:u,ZS,XG, tag-beats-FLAG).Regression: with
--Methylationoff, output is unchanged; the existing suite passes 17/17.In-silico titration on real EM-seq (SEQC2 EpiQC, ~10×, GRCh38, 1000g 100k panel): host HG001/NA12878 mixed with an unrelated contaminant at known α, recovered FREEMIX vs spiked —
Clean at 0% and accurate/monotonic across two contaminant populations. A mild underestimate grows with α and with contaminant divergence (up to −3.6 pp at 25% cross-population), consistent with reference-mapping bias against the contaminant's non-reference alleles amplified by three-letter bisulfite alignment; the QC-relevant low-α regime is well calibrated.
Compatibility
Opt-in via
--Methylation; default behavior and output format are unchanged (the existing test suite passes untouched). The first commit fixes a pre-existing base/quality desync in the pileup loop — a read with a deletion at a marker pushed a quality with no base — which the per-observation filter then builds on.