Does --single-code h perform 5hmC-vs-all or 5hmC-only DMR analysis?

[danammak]

Hello Modkit team,

Thank you for developing and maintaining Modkit.

I have a question regarding the behaviour of the "--single-code" option in differential methylation analysis.
Suppose a modBAM contains both 5mC ("m") and 5hmC ("h") calls from a Dorado model, and I run:

"modkit dmr pair --single-code h"

Could you please clarify how non-target modification calls are handled internally?

For example, at a CpG site with the following read-level composition:

• 10 reads called as 5hmC ("h")
• 80 reads called as 5mC ("m")
• 10 reads called as canonical C

When using "--single-code h":

1. Are 5mC calls treated as canonical/unmodified bases?
2. Are 5mC calls excluded from the coverage denominator?
3. Are 5mC calls retained as an "other modification" category?
4. How exactly is the modification proportion calculated for DMR testing?

Similarly, when using:
modkit dmr pair --single-code m

how are 5hmC calls handled?

My main biological question is whether "--single-code h" can be interpreted as a differential analysis of only 5hmC sites, independent of 5mC. In other words, if a CpG is called as 5mC in some reads and 5hmC in other reads, does the presence of 5mC contribute to the denominator or influence the estimated 5hmC frequency?

I am interested in generating separate 5mC and 5hmC landscapes from Nanopore data and would like to understand whether:

1. "--single-code h" effectively measures 5hmC versus everything else, or
2. it measures 5hmC while excluding 5mC calls from the analysis, or
3. it uses some other counting model internally.

Any clarification on the underlying counting model, treatment of other modification classes, and denominator used for DMR testing would be greatly appreciated.

Thank you very much for your help.

[ArtRand]

Hello @danammak,

When you use --single-code h you make two categories: 5hmC calls and ¬5hmC calls. For example:

Condition	5hmC	5mC	Canonical
A	2	5	3
B	2	3	5

This will become:

Condition	5hmC	Not 5hmC
A	2	8
B	2	8

So effect size is (20% - 20% = 0%), and similarly the other statistics will show small/no change.

On the other hand, if you used --single-code m:

Condition	5mC	Not 5mC
A	5	5
B	3	7

Now the effect size is 50% - 30% = 20%.

So to answer your questions concretely:

When using "--single-code h":

Are 5mC calls treated as canonical/unmodified bases?

Yes, although I'd say they are treated as NOT 5hmC, and that canonical and 5mC are both NOT 5hmC.

Are 5mC calls excluded from the coverage denominator?

No

Are 5mC calls retained as an "other modification" category?

Yes.

How exactly is the modification proportion calculated for DMR testing?

Shown above.

My main biological question is whether "--single-code h" can be interpreted as a differential analysis of only 5hmC sites, independent of 5mC.

Yes, that's the whole idea!

In other words, if a CpG is called as 5mC in some reads and 5hmC in other reads, does the presence of 5mC contribute to the denominator or influence the estimated 5hmC frequency?

The 5mC reads contribute to the denominator, but they don't influence the estimated 5hmC stoichiometry.

I'll add this example to the documentation if you think it makes sense.