Clarification on recommended filtering thresholds for confident 5mC calls from modkit pileup

[Alexkortsi]

Hi developers,

I’m using modkit pileup (v0.6.0) to extract confident 5mC calls in BED format from Nanopore data. I’ve successfully generated the .bed.gz output, and I understand the meaning of the columns based on your documentation. However, I’m still unsure about how to properly filter this output to retain high-confidence methylated sites.

My questions is what are the recommended filtering criteria for downstream use of 5mC calls from modkit pileup?
Are there any published or internal best practices for balancing sensitivity and specificity in calling modified bases (especially 5mC in CpG context)?
Should thresholds be chosen dynamically per-sample, e.g., based on modkit summary output?

For example in my .bed.gz i have:

contig_23 28713 28714 m 23 + 28713 28714 255,0,0 23 4.35 1 22 0 1 9 1 1

From what I understand:

Nmod = 1
Nvalid_cov = 23
fraction_modified = 4.35%
Nfail = 1

So this doesn’t pass strict thresholds, but I’m trying to define the optimal cutoff.

Any clarification on which values are most robust and biologically meaningful would be much appreciated. I’d also love to hear how others in the community are doing this in practice.

Thanks again!

[ArtRand]

Hello @Alexkortsi,

The recommended best practice is to let Modkit estimate the 10th percentile of base modification call probabilities and discard (/fail) calls with probability below this number. We also recommend estimating this threshold on a per-sample basis. One note, since the 10th percentile value is an estimate, if you change --interval-size, --threads, --region, or --num-reads the estimate may change slightly (I'm working on removing this). One way around this is to run modkit sample-probs once per sample and use the value it outputs for that sample.

To your questions specifically:

My questions is what are the recommended filtering criteria for downstream use of 5mC calls from modkit pileup?

It depends what you're trying to do. Using the defaults:

$ modkit pileup $bam $bed --modified-bases 5mC 5hmC --ref $genome --bgzf

Is the recommended approach.

Should thresholds be chosen dynamically per-sample, e.g., based on modkit summary output?

Yes, dynamically per-sample.

Are there any published or internal best practices for balancing sensitivity and specificity in calling modified bases (especially 5mC in CpG context)?

We've published some information regrading the accuracy of the modified base classifiers here.

contig_23 28713 28714 m 23 + 28713 28714 255,0,0 23 4.35 1 22 0 1 9 1 1
From what I understand:
Nmod = 1
Nvalid_cov = 23
fraction_modified = 4.35%
Nfail = 1
So this doesn’t pass strict thresholds, but I’m trying to define the optimal cutoff.

What do you mean "doesn't pass strict thresholds". If you tell me more about what you're trying to do maybe I can advise you.