Explanation for modkit stats output

[dietvin]

Hi,

I'm currently trying the modkit stats command and wanted to clarify a couple of things.

From what I understand, count_valid_x corresponds to the number of canonical bases that correspond to a given modification. So when the bedmethyl file contains only CpG sites, the count_valid_m entry corresponds to the number of CpG sites in a given region, right?

Then count_x corresponds to the number of bases that are considered modified. How is the modification status of a single site determined here? Is there a set threshold for the percent modified value from the bedmethyl file?

Thank you!

Best
Vincent

[ArtRand]

Hello @dietvin Sorry about the delay in getting back to you.

From what I understand, count_valid_x corresponds to the number of canonical bases that correspond to a given modification. So when the bedmethyl file contains only CpG sites, the count_valid_m entry corresponds to the number of CpG sites in a given region, right?

Correct. The output will have groups of 3 columns count_x, count_valid_x, percent_x. count_x is the sum of the reads reporting modification code x in the region, count_valid_x is the sum of valid coverage for this modification code in the region. Finally, percent_x is count_x / count_valid_x. So if you have run pileup with --cpg then count_valid_m will be the number of CpGs in the region with at least one valid base modification call (not necessarily the number of CpG sites, since some could have zero valid coverage).

How is the modification status of a single site determined here?

There isn't any logic around determining the modification status of single sites. The calls from each read, at each site, are summed together. If you wanted something that, say, calculated the number of sites with >x% 5mC, that would be a different calculation. Don't quote me on this, but I think Bedtools can do this.

Is there a set threshold for the percent modified value from the bedmethyl file?

Based on the above, it's probably clear that there isn't, but if this is a feature you have a use case for I can look into adding it (or you can submit a PR!).