--motif % 0 and --ignore % not working as expected?

[CarolinaA09]

Hello, I basecalled all 4 modifications with Dorado and then used this modkit command: modkit pileup --ref %s/dmel-all-chromosome-current.fasta --motif %s 0 --filter-threshold T:%s --filter-threshold C:%s --filter-threshold A:%s --mod-threshold a:%s --ignore 17596 --with-header %s.final.sorted.bam %s.%s.%s.pileup.final.sorted.bam.bed --log-filepath %s%s.%s.pileup.header.lo_g. It was of my understanding that I will be getting only the m6A modifications (a) in the output table but I am also getting m5C and pseudouridine. Example:

mitochondrion_genome 736 737 mitochondrion_genome 736 737 a 113 + 736 737 255,0,0 113 28.32 32 81 0 132 7282 7 0
mitochondrion_genome 1416 1417 mitochondrion_genome 1416 1417 a 124 + 1416 1417 255,0,0 124 12.90 16 108 0 66 1716 3 0
mitochondrion_genome 2585 2586 mitochondrion_genome 2585 2586 17802 5 + 2585 2586 255,0,0 5 80.00 4 1 0 57 5222 3188 0
mitochondrion_genome 2963 2964 mitochondrion_genome 2963 2964 a 49 + 2963 2964 255,0,0 49 12.24 6 43 0 63 8564 1 0
mitochondrion_genome 3090 3091 mitochondrion_genome 3090 3091 a 116 + 3090 3091 255,0,0 116 99.14 115 1 0 4 2981 5 0
mitochondrion_genome 13299 13300 mitochondrion_genome 13299 13300 17802 8 - 13299 13300 255,0,0 8 12.50 1 7 0 38 7460 38 0
mitochondrion_genome 13299 13300 mitochondrion_genome 13299 13300 a 31 - 13299 13300 255,0,0 31 19.35 6 25 0 38 7460 15 0
mitochondrion_genome 14146 14148 mitochondrion_genome 14146 14147 a 173 - 14146 14147 255,0,0 173 23.12 40 133 0 57 2076 6 0
mitochondrion_genome 14146 14148 mitochondrion_genome 14147 14148 a 290 - 14147 14148 255,0,0 290 68.28 198 92 0 103 1881 39 0
mitochondrion_genome 14146 14148 mitochondrion_genome 14147 14148 m 34 - 14147 14148 255,0,0 34 100.00 34 0 0 103 1881 295 0
rDNA 71960 71962 rDNA 71960 71961 m 8 + 71960 71961 255,0,0 8 50.00 4 4 0 13 1557 55 0
Why? is there a way I can only get the "a" modifications? I am not filtering those myself because then modkit dmr pair will fail. Am I doing something wrong here? Thank you for your help and attention!

[ArtRand]

Hello @CarolinaA09

Sorry for the delay getting back to you. It looks like you're probably using an older version of Modkit (before 0.6.0). In versions <0.6.0 the pileup command will emit records where reads have mismatched bases with respect to the reference. So if you have called base modifications on all 4 primary sequence bases, any reference A with at least 1 mismatch will have a base modification record reporting on that read. The good news is that modkit dmr pair will only use base modifications on the reference primary base, so in other words it will discard any m5C or PseU records at reference A bases. My recommendation is to update to Modkit v0.6.0 and refine your command as below:

# modkit v0.6.0
modkit pileup \
  --ref %s/dmel-all-chromosome-current.fasta \
  --preload-references \
  --modified-bases m6A inosine \
  --filter-threshold T:%s --filter-threshold C:%s --filter-threshold A:%s \
  --mod-threshold a:%s --mod-threshold 17596:${Inosine_mod_threshold} \
  --with-header %s.final.sorted.bam %s.%s.%s.pileup.final.sorted.bam.bed --log-filepath %s%s.%s.pileup.header.log. 

You will get Inosine records in the bedmethyl, but there will not be any records for primary sequence bases C, G, or T. This bedmethyl table should be suitable for modkit dmr. If you're only interested in m6A changes, you can filter the table for records that have a change in m6A greater than some threshold or discard records where the MAP-based p-value is less than some threshold but the m6A values between samples is the same. If it were me, I would keep the Inosine calls for downstream analysis - but if you give me more details on your experiment maybe I can advise otherwise.

If you must remove the Inosine calls, run

modkit modbam adjust-mods ${bam} ignore_inosine.bam --ignore 17596

first then use this BAM as input to pileup.