estimated pass threshold with default settings very high (v0.5.0 and v.0.6.1)

[AmberVerhasselt]

Hi!

I ran modkit pileup (v.0.5.0) on one of my ONT samples for 5mC and noticed that the pass threshold with the default settings was put very high. A lot of methylated calls were filtered out giving very low percentages. This is what the logfile shows:
estimated threshold too high 1, using 0.9941406

I also tried runnning it with the newer version (v0.6.1), but got the same issue. This is what the log file shows:
estimated threshold too high 1, using 0.9746094 (maximum explicit canonical probability)

I also ran the sample probs to see what the thresholds and distribution looks like. The threshold was the same as for the pileup for both versions, put at 1:
base percentile threshold
C 10 1
C 50 1
C 90 1

I understood from the documentation that with the default settings, the lowest 10% of probabilities are filtered out. But when i look at the probability distributions of this sample it does not seem to put the pass threshold at the 10th lowest percentile, but a lot higher (see attached picture)? This leads to the majority of data not passing the threshold and being filtered out.

For reference, all my ONT samples get the exact same threshold put. I've played around with the filter settings myself, but there is no best practice for what threshold to put it at (yet) so this is not the best option for me at this moment.

Could you tell me why this happens and how I could fix it?

[ArtRand]

Hello @AmberVerhasselt,

Sorry about the delay getting back to you.

When Modkit estimates the pass threshold it accumulates all of the call-probabilities* for each primary base (A, C, G, T) together. With the all-context models, there is a space-saving optimization where base modification calls with probability of modification $\lt 5\%$ (by default) are "implicit unmodified" and Modkit interprets these as unmodified with 100% probability. If you have many more unmodified calls than modified ones it's possible that all of the modified calls are in the lowest 10%, leaving only the 100%-confident unmodified calls. There is a check in pileup v0.6.z that fall back to the highest explicit probability, often times that will still fail out a lot of the modified calls.

The way to handle this, as you've started to see, is to use the --mod-threshold option with a modification-specific threshold. I would look at the probabilities.tsv table and find the m probability that is roughly the 10th percentile rank and use that value.

* A call probability is $\arg\max_{x \in \text{mods and unmodified}} P_x$, the most likely base modification state.

[AmberVerhasselt]

Hi @ArtRand,

Thanks for the explanation, that really helps! If I look at all Cs, that seems to explain the threshold setting indeed. I'll look at probabilities for C:m specifically and see if I can create my own mod-threshold that way.

Thanks a lot!