Using modkit pileup and --modified-bases

[danammak]

Hi Team,

modkit pileup
--cpg
--combine-strands
--modified-bases C:h
--bgzf
--ref ${ref}
--threads ${threads}
--log-filepath "$logfile"
"$bam_file"
"$output_pileup"

I used the above options to get pileup data purely for only 5hmC using modkit latest version 0.6.0. But, I see that the output file has both 'm' and 'h'. The earlier "--ignore" option was good to avoid this.
chr1 10468 10469 h 6 . 10468 10469 255,0,0 6 33.33 2 4 0 0 0 3 11
chr1 10468 10469 m 6 . 10468 10469 255,0,0 6 0.00 0 4 2 0 0 3 11
chr1 10470 10471 h 4 . 10470 10471 255,0,0 4 25.00 1 3 0 0 3 3 10
chr1 10470 10471 m 4 . 10470 10471 255,0,0 4 0.00 0 3 1 0 3 3 10
chr1 10488 10489 h 1 . 10488 10489 255,0,0 1 0.00 0 1 0 1 0 0 18
chr1 10488 10489 m 1 . 10488 10489 255,0,0 1 0.00 0 1 0 1 0 0 18
chr1 10496 10497 h 5 . 10496 10497 255,0,0 5 20.00 1 4 0 0 1 0 14
chr1 10496 10497 m 5 . 10496 10497 255,0,0 5 0.00 0 4 1 0 1 0 14
chr1 10524 10525 h 1 . 10524 10525 255,0,0 1 0.00 0 1 0 0 1 0 18
chr1 10524 10525 m 1 . 10524 10525 255,0,0 1 0.00 0 1 0 0 1 0 18

Below are the modbam summaries for both 5mC and 5hmC:

1. For 5mC modbam:

bases C

total_reads_used 9085

count_reads_C 9085

pass_threshold_C 0.80859375

base code pass_count pass_frac all_count all_frac
C - 249191 0.52520305 281094 0.53459847
C m 225275 0.47479692 244710 0.46540156
C h 0 0 0 0

2. For 5hmC modbam:

bases C

total_reads_used 9085

count_reads_C 9085

pass_threshold_C 0.78125

base code pass_count pass_frac all_count all_frac
C - 259450 0.9509585 282582 0.9333439
C h 13380 0.04904153 20181 0.0666561
C m 0 0 0 0

Below are the options used for generating 5mC and 5hmC only modbams.

1. For 5mC:
   modkit adjust-mods -t 64 --filter-probs --mod-threshold h:1.1 "$bam_file" "$output_bam"

2. For 5hmC:
   modkit adjust-mods -t 64 --filter-probs --mod-threshold m:1.1 "$bam_file" "$output_bam"

I did not use --ignore in modkit adjust-mods because I saw a discrepancy in my results, and so came up with the idea to use m:1.1 or h:1.1, which gave logical results when I ran modkit summary as shown above.

Could you please suggest with modkit pileup for 5hmC using version 0.6.0?

[ArtRand]

Hello @danammak,

This is a bug in Modkit v0.6.0 that has been raised here before. I've attached a build that has a fix, you can also pass the --use-dynamic flag with the released v0.6.0 to get the correct results. I would not pass a threshold greater than 1 in order to shift the 5hmC or 5mC calls into the failed category. I recommend doing one of the following:

1. Use pileup with --modified-bases [5hmC|5mC] to generate records with only 5mC or 5hmC counts. Use the build here or pass --use-dynamic.
2. If you want to combine modification counts pass the --combine-mods flag.
3. The --ignore flag isn't supported in Pileup as of 0.6.0, if you want to keep using this functionality (although --combine-mods is more robust), use modkit modbam adjust-mods --ignore docs here on the modBAM before running pileup.

modkit_v0.6.1_4d52624.zip

[danammak]

Thank you so much for your reply. I will give it a try with this.

[ArteenTorabiMarashi]

Hello, hope all is well. I experienced the same issue and tried both --use-dynamic and the new build you attached, but neither are working. They both still produce bedmethyl files with both h and m flags, and when trying to split and use them in dmr pair I am running into the encountered invalid bedMethyl record(s) error where there is a mismatch between N valid coverage and the sum of N canonical and N modified counts (because I assume the N other needs to be added to N canonical).