Bio Long Read Sequencing Nanopore Methylation

Version Compatibility

Reference examples tested with: modkit 0.3+, dorado 1.0+, minimap2 2.28+, samtools 1.19+, pb-CpG-tools 2.3+.

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags

Inputs that determine what is even possible - record them:

• The basecaller MODIFICATION model (e.g. 5mCG_5hmCG) fixes which mods can ever be piled up; it must be requested at basecall time and cannot be added later.
• The MM/ML tags must survive every fastq/alignment step or methylation is silently lost.
• modkit auto-estimates the pass threshold from the data (per run); fix it for cross-sample comparisons.

If code throws an error, introspect the installed tool (modkit pileup --help, modkit --help) and adapt the example to the actual API rather than retrying.

Nanopore Methylation

"Call methylation from my long reads" -> First confirm the MM/ML tags exist and survived alignment, then pile them into per-site bedMethyl - because methylation is a basecalling decision, not something that can be added now.

• CLI: modkit pileup aligned.bam out.bed --ref ref.fa --cpg --combine-strands

The Single Most Important Modern Insight -- Methylation Is a Basecalling Decision, and the Tags Silently Die in Alignment

Two facts gate the entire skill:

1. If the reads were not basecalled with a modification model, the signal is already gone. Mods are inferred from raw signal at basecall time (ONT Remora model in Dorado; PacBio kinetics model in jasmine) and written into the unaligned BAM as MM/ML tags. A plain BAM or a FASTQ cannot yield methylation - there is no post-hoc tool. The only fix is re-basecalling from POD5 (dorado basecaller sup,5mCG_5hmCG pod5/). The agent's FIRST move is to check the tags exist: samtools view in.bam | head | grep -o 'MM:Z:[^\t]*'.
2. The MM/ML tags silently die in a normal alignment workflow. samtools fastq drops auxiliary tags unless given -T MM,ML; minimap2 ignores them unless given -y; hard-clipping breaks MM's per-base skip counting unless -Y is set. Miss any one and the aligned BAM still sorts, indexes, and looks fine, but modkit pileup returns an empty/all-canonical bedMethyl with no error. Use dorado aligner (carries tags natively) or samtools fastq -T MM,ML | minimap2 -y -Y, and re-grep for MM:Z AFTER alignment.

End-to-End modBAM Pipeline

# 1. MODS-BASECALL (from POD5; the only step that can ever produce methylation)
dorado basecaller sup,5mCG_5hmCG pod5/ > calls.bam      # unaligned BAM, has MM/ML

# 2. TAG-PRESERVING ALIGN (route a is simplest)
dorado aligner ref.mmi calls.bam > aligned.bam                                  # a) native
samtools fastq -T MM,ML calls.bam | minimap2 -y -Y -ax lr:hq ref.fa - \
  | samtools sort -o aligned.bam && samtools index aligned.bam                  # b) manual
samtools view aligned.bam | head | grep -q 'MM:Z' && echo 'tags survived'       # verify!

# 3. PILEUP -> bedMethyl (auto-thresholds at the 10th percentile of ML; NOT 0.5)
modkit pileup aligned.bam out.bed --ref ref.fa --cpg --combine-strands
bgzip out.bed && tabix -p bed out.bed.gz

# 4. (optional) DIFFERENTIAL methylation, long-read native
modkit dmr pair -a A.bed.gz -b B.bed.gz --ref ref.fa --regions cpgislands.bed -o dmr.tsv

MM / ML Tag Spec (what the numbers mean)

• MM:Z encodes modification positions: <canonical base><strand><mod code><. or ?>,<skip counts>;. Mod codes: m=5mC, h=5hmC, a=6mA, c=4mC. The ./? modifier is load-bearing: . = skipped bases are implicitly canonical (count toward the unmodified denominator); ? = skipped bases are no-call/unknown (land in Nnocall, outside the denominator). Misreading ? as . inflates the canonical denominator and deflates methylation.
• ML:B:C is a uint8 per call: value N means probability in [N/256, (N+1)/256), so 255 is ~0.998, never exactly 1.0. Do not threshold == 1.0.

bedMethyl Columns (modkit, 18 columns)

Cols 1-9 are tab-delimited BED9; cols 10-18 are space-delimited (a parsing gotcha). The ones that matter:

For count-based DMR (DSS/methylKit) hand over Nmod and Nvalid_cov, never percent_modified.

Decision Tree by Scenario

Phased Allele-Specific Methylation

# Order is strict: align (tags preserved) -> phase+haplotag -> pileup partitioned by HP
# 1-2. Clair3/DeepVariant -> whatshap/longphase phase + haplotag (adds HP:i:1/2) -> haplotype-phasing
modkit pileup aligned.haplotagged.bam asm_out/ --ref ref.fa --cpg --combine-strands --partition-tag HP
# --partition-tag writes one UNCOMPRESSED bedMethyl per HP value into asm_out/, named by the tag
# value (e.g. 1.bed, 2.bed). bgzip + tabix each before dmr, which requires indexed inputs:
bgzip asm_out/1.bed && tabix -p bed asm_out/1.bed.gz
bgzip asm_out/2.bed && tabix -p bed asm_out/2.bed.gz
modkit dmr pair -a asm_out/1.bed.gz -b asm_out/2.bed.gz --ref ref.fa -o asm.tsv

Each haplotype gets ~half the coverage, so the per-site 10x floor effectively wants ~20x total. Imprinted loci (one haplotype ~fully methylated) are the canonical positive control.

Per-Method Failure Modes

No methylation in a plain BAM

Trigger: BAM basecalled without a mods model. Mechanism: mods are a basecall-time decision. Symptom: no MM:Z tags; modkit returns nothing. Fix: re-basecall from POD5 with a mods model; there is no post-hoc tool.