modkit questions about filtering

[weishwu]

Hi @ArtRand

I have a couple of questions about modkit (v0.6.1):

1. I used to run modkit pileup with --ignore h since I need 5mCG only. I noticed that this argument is gone in the latest modkit version (v0.6.1). Instead, it has --modified-bases. So I can use --modified-bases 5mC to do the same thing as the previous --ignore h, right?

2. When I run modkit pileup with the default setting, I get this in the log: "Threshold of 0.6621094 for base C is low. Consider increasing the filter-percentile or specifying a higher threshold." My ONT sequencing was targeted to a 34Mb panel of ROI enriched for imprinted regions, so I wonder if this could make the default estimate unstable (which can potentially affect sample comparison?). Is it OK to set an explicit threshold for the filter like this:

modkit pileup \
      -t {threads} \
      --cpg \
      --modified-bases 5mC \
      --filter-threshold C:0.8 \
      --ref {params.genome_fasta} \
      --combine-strands \
      --include-bed {params.roi_full} \
      --bgzf \
      --phased \
      --prefix {params.out_prefix} \
      {input} \
      {params.out_dir}

3. When I run modkit dmr pair with the outputs from above, I get this error:

batch failed: invalid data, valid coverage (51) is not equal to the sum of canonical and modified counts (49), [BedMethylLine { chrom: "chr1", interval: Interval { start: 1008712, stop: 1008713, val: () }, raw_mod_code: Code('m'), strand: Both, count_methylated: 31, valid_coverage: 51, count_canonical: 18, count_other: 2, count_delete: 0, count_fail: 1, count_diff: 0, count_nocall: 0 }] chrom: chr1 starting at 1008712, stopping > 6 batches processed > Error! invalid-bedmethyl-data

It is because in some lines, there are reads categorized as Nother_mod, making Nmod + Ncanonical != Nvalid_cov. Is this because of the --filter-threshold C:0.8 parameter I used in modkit pileup?

[ArtRand]

Hello @weishwu,

1. I used to run modkit pileup with --ignore h since I need 5mCG only. I noticed that this argument is gone in the latest modkit version (v0.6.1). Instead, it has --modified-bases. So I can use --modified-bases 5mC to do the same thing as the previous --ignore h, right?

Not quite. The --ignore option has been removed as you noticed, but --modified-bases 5mC does not perform the collapse algorithm that --ignore used to. If you have a modBAM with 5hmC and 5mC calls and you use --modified-bases 5mC this is equivalent to filtering out the 5hmC records not removing/ignoring the 5hmC probabilities. The recommended way to process 5mC, 5hmC data is to either

1. Use --modified-bases 5mC 5hmC --bgzf with pileup then use dmr pair on the results or
2. Use --combine-mods --modified-bases C (or --modified-bases 5mC 5hmC) and use the result in modkit dmr pair.

2. When I run modkit pileup with the default setting, I get this in the log: "Threshold of 0.6621094 for base C is low. Consider increasing the filter-percentile or specifying a higher threshold." My ONT sequencing was targeted to a 34Mb panel of ROI enriched for imprinted regions, so I wonder if this could make the default estimate unstable (which can potentially affect sample comparison?). Is it OK to set an explicit threshold for the filter like this..

Yes using a manual threshold is OK, we've generally found the probabilities to be well calibrated. I may suggest running modkit modbam sample-probs --hist on your targeted sample and inspecting the output distributions as a sanity check.

When I run modkit dmr pair with the outputs from above, I get this error [...] It is because in some lines, there are reads categorized as Nother_mod, making Nmod + Ncanonical != Nvalid_cov. Is this because of the --filter-threshold C:0.8 parameter I used in modkit pileup?

This is because for DMR to work you need to specify all of the modifications that were originally present in the modBAM (--modified-bases 5hmC 5mC). Lots of people have been encountering this error. There is a check in the DMR routine that makes sure all of the modified counts (N_mod and N_other) have codes assigned to them. When you have a bedMethyl that is missing the h records - it will complain. I need to come up with a robust solution for this. But in the mean time I'd use one of the two suggestions at the top of this message.

[weishwu]

Thanks @ArtRand for your help!
Attached are the output files from running modkit modbam sample-probs --hist on a sample. Do you think --filter-threshold C:0.8 would be a good cutoff based on the histogram?
Another option is to increase --filter-percentile (e.g. to 0.2). Which way is better?

thresholds.tsv
counts.html
proportion.html
probabilities.tsv