Large no. of reads go into "Nfail" column in pileup

[gsukrit]

Dear authors,

I ran the modkit pileup command for modified bases a, h, and c. It appears that when I specify the 10th percentile threshold, I get an extremely low overall modification frequency (0.09%), since all mod probabilities are ~1. When I add a slightly conservative modification threshold, a large chunk of reads are going to the Nfail column.
I have added the log files for both runs of cytosine modifications.
Can you please advise as to what is the best way to move ahead in such a case?

P54log.txt

P54log.txt

Thank you!

[ArtRand]

Hello @gsukrit

Thanks for sending the logs over and sorry about the slow reply. One gotcha I noticed is that you have --mod-threshold C:0.78 isn't going to do exactly what you want. It's going to match to "C" as the mod code, which although I appreciate that means "any C-mod" the code currently will not match it to "h" or "m". I can add this logic as it's a bit of an oversight. If you want to use 0.78 as your "h" and "m" threshold add this to you command line: --mod-threshold m:0.78 --mod-threshold h:0.78.

What depth do you have? If you have depth higher than 30-50X I would look at the results from --no-filtering, in my experience as the depth gets high the estimates of stoichiometry at each site depend less on the threshold values. If you don't have very high depth, try modifying the command line as suggested above.