Modkit pileup (v. 0.6.0) - modification percentage over 100%

[SimonBache]

Hi,

I tried running modkit pileup on a bam file using modkit v0.6.0. For some lines, the resulting bedmethyl file was inconsistent :

I 1097112 1097113 a 8 - 1097112 1097113 255,0,0 8 112.50 9 65535 0 9441 7192 257 1649
I 1097118 1097119 a 44 - 1097118 1097119 255,0,0 44 102.27 45 65535 0 4039 5614 167 224
I 1275566 1275567 a 14 + 1275566 1275567 255,0,0 14 107.14 15 65535 0 276 2387 90 157
I 1973681 1973682 a 51 - 1973681 1973682 255,0,0 51 101.96 52 65535 0 1453 10546 377 909
I 1973687 1973688 a 63 - 1973687 1973688 255,0,0 63 101.59 64 65535 0 226 9461 154 1403
I 1973728 1973729 a 10 - 1973728 1973729 255,0,0 10 110.00 11 65535 0 2436 6060 836 2666
I 1973815 1973816 a 24 - 1973815 1973816 255,0,0 24 104.17 25 65535 0 2414 8050 1917 1798
I 1973896 1973897 a 21 - 1973896 1973897 255,0,0 21 104.76 22 65535 0 2549 6701 1584 1251
I 1973939 1973940 a 2 - 1973939 1973940 255,0,0 2 150.00 3 65535 0 1589 4421 1738 674

As you can see, the percent of m6A modification is higher than 100%.

When using modkit v.0.4.1 on the same input bam file, I don't have the same error :

I 1097112 1097113 a 4937 - 1097112 1097113 255,0,0 4937 0.00 0 4937 0 344 346 9 57
I 1097118 1097119 a 5354 - 1097118 1097119 255,0,0 5354 0.06 3 5351 0 50 259 6 13
I 1275566 1275567 a 7746 + 1275566 1275567 255,0,0 7746 0.01 1 7745 0 27 266 10 13
I 1973681 1973682 a 7117 - 1973681 1973682 255,0,0 7117 0.04 3 7114 0 77 759 29 59
I 1973687 1973688 a 7200 - 1973687 1973688 255,0,0 7200 0.07 5 7195 0 16 717 13 100
I 1973728 1973729 a 7125 - 1973728 1973729 255,0,0 7125 0.00 0 7125 0 186 443 61 186
I 1973815 1973816 a 6727 - 1973815 1973816 255,0,0 6727 0.01 1 6726 0 168 634 147 175
I 1973896 1973897 a 6801 - 1973896 1973897 255,0,0 6801 0.03 2 6799 0 197 539 118 99
I 1973939 1973940 a 6971 - 1973939 1973940 255,0,0 6971 0.00 0 6971 0 119 344 150 47

It looks like there is something wrong with the calculation of the percentage of modification in v0.6.0.

I attach the 2 bedmethyl with the problematic rows coming from the 2 versions of modkit :

bedmethyl_041.bed.txt

bedmethyl_060.bed.txt

Thank you for help.

[ArtRand]

Hello @SimonBache,

Thank you for the report. Any chance you could show me the command you used? Also is it possible to send or attach a sample of the modBAM that produces the issue? (This really helped me debug a previous issue migrating to v0.6.0). You can attach it here or email me art.rand[at]nanoporetech.com

Thanks!

[ArtRand]

Hello @SimonBache,

Thank you for sending over the data for debugging. The issue is that Modkit v0.6.0 does not handle very high coverage (>60,000X) correctly. I'll have a patch out to you ASAP. In the mean time you (and anyone else with exceptionally high coverage samples) should probably stay on v0.5.1. Sorry about the inconvenience.

[SimonBache]

Hello @ArtRand,

No problem, I will stay on the previous version for now.

Thanks.

[ArtRand]

Hello @SimonBache,

Here is a build that adds a flag --high-depth and optionally --max-depth that will handle a modBAM with depth like yours similarly to how Modkit $\le$ v0.5.1 did. Looking at the script you sent, there are a couple things you could add to speed this up:

# previously you had something like this

modkit pileup -t 1 $input $output
bgzip -c $output > ${outputgz}
tabix -p bed ${outputgz}

# you can do this to skip the bgzip step
modkit pileup $input $outputgz \
  --bgzf \ # output a compressed file suitable for tabix
  --modified-bases 6mA \ # use the optimized routines
  --ref $ref \ # optimized routines require a reference
  --interval-size 10000 \ # with very high depth, decreasing the interval size can increase performance
  -t 4 # have a little parallelism, often times you don't need more than 4 worker threads

tabix -p bed ${outputgz}

Hope this helps, there are a few more things that need to get into v0.6.1, but I think this behavior will be included as it is in this build.

modkit_v0.6.1_ecdb2eb.zip

[SimonBache]

Hello @ArtRand ,

Great, thank you very much for the build and your advices on the script !