modkit dmr pair fails with Error: “valid coverage is not equal to the sum of canonical and modified counts”

[NathaneilKnight]

Hi Modkit maintainers,

First of all, thank you very much for developing and maintaining modkit—it’s been extremely helpful for my analysis.

modkit pileup \
  ./mm10_sorted.bam \
  ./mm10.bed.gz \
  --modified-bases 5mC 5hmC 6mA \
  -t 8\
  --reference ./mm10.fa \
  --log ./mm10_log.txt \
  --bgzf \

Everything seems OK, but when I run the following command:

modkit dmr pair \
  -a ./A1_mm10.bed.gz \
  -a ./A2_mm10.bed.gz \
  -b ./B1_mm10.bed.gz \
  -b ./B2_mm10.bed.gz \
  -o ./B_vs_A_singlesite_dmr_C.bed \
  --segment ./B_vs_A_coarse-grained_segement_dmr_C.bed \
  --ref ./mm10.fa \
  --base C \
  --threads 64 \
  --log-filepath ./B_vs_A_singlesite_dmr_C.log

"valid coverage is not equal to the sum of canonical and modified counts" came up like this:

reading reference FASTA at "./mm10.fa"
66 common sequence(s) between FASTA and both samples
running single-site analysis
using default prior, Beta(α: 0.55, β: 0.55)
estimating max coverages from data
sampled 44386636 a records and 44395324 b records, calculating max coverages for 95th percentile
calculated max coverage for a: 21 and b: 21
running with replicates and matched samples
batch failed: invalid data, valid coverage (346) is not equal to the sum of canonical and modified counts (65882), [BedMethylLine { chrom: "chr9", interval: Interval { start: 3026160, stop: 3026161, val: () }, raw_mod_code: Code('h'), strand: Positive, count_methylated: 1, valid_coverage: 346, count_canonical: 65180, count_other: 701, count_delete: 18264, count_fail: 3550, count_diff: 2477, count_nocall: 0 }, BedMethylLine { chrom: "chr9", interval: Interval { start: 3026160, stop: 3026161, val: () }, raw_mod_code: Code('m'), strand: Positive, count_methylated: 701, valid_coverage: 346, count_canonical: 65180, count_other: 1, count_delete: 18264, count_fail: 3550, count_diff: 2477, count_nocall: 0 }] chrom chr9 starting at 3026160, ending at 3026161, stopping
failed to send on channel, sending on a disconnected channel
encountered error whilst processing, stopping
8695 batches processed
Error! invalid-bedmethyl-data

When I check the log file, I found the "invalid data" Error Only occured in chrUn_JH584304 and chr9.

B_vs_A_singlesite_dmr_C.log

Some additional information:

1. I haven't filtered bedMethyl files in any way.
2. I'm using the newest version (0.6.0) of Modkit.
3. There are 3 types of modification in my modbam files: 5mC, 5hmC, 6mA.
4. My modbam files were derived from Dorado align command with mm10 as reference,which is the same as the reference used in modkit.
   Thanks again for your time and for all the work you’ve put into Modkit!

Best regards,
Nathaneil

[ArtRand]

Hello @NathaneilKnight,

Sorry about the error. Do you know if this run has particularly high sequencing depth? My guess is that this issue is due to an existing bug in v0.6.0 when used on high-depth modBAMs. Could you try the build available in this comment?

Summary of how to use the build with high depth:

modkit pileup \
  ./mm10_sorted.bam \
  ./mm10.bed.gz \
  --modified-bases 5mC 5hmC 6mA \
  --high-depth \ # <-- need this with the test build
  -t 8 \
  --reference ./mm10.fa \
  --log ./mm10_log.txt \
  --bgzf \
  ...

If this doesn't sort out the problem, would it be possible to send me a BAM with just this region? If you can email it to art.rand[at]nanoporetech.com or attach it here that would be great.

samtools view -bh $bam chr9:3026160-3026161 > $file_to_send

Thanks.

[NathaneilKnight]

Hi, @ArtRand!

Thanks very much for your timely response.

Sorry about the error. Do you know if this run has particularly high sequencing depth? My guess is that this issue is due to an existing bug in v0.6.0 when used on high-depth modBAMs. Could you try the build available #547 (comment)?

The average sequencing depth is around 20–30×. After switching to the v0.6.1 build you provided and running pileup on chr9 only, the error disappeared when I subsequently ran dmr (commands below):

parallel -j 8 modkit pileup    ./{}_mm10_sorted.bam     ./{}_mm10.bed.gz    --modified-bases 5mC 5hmC    --cpg    -t 8    --reference ./mm10.fa    --log ./{}_mm10_log.txt --region chr9:1-124595110 --high-depth  --bgzf ::: $SAMPLES

#For dmr pair command, modkit is v6.0.0
modkit dmr pair \
  -a ./A1_mm10.bed.gz \
  -a ./A2_mm10.bed.gz \
  -b ./B1_mm10.bed.gz \
  -b ./B2_mm10.bed.gz \
  -o ./B_vs_A_singlesite_dmr_C.bed \
  --segment ./B_vs_A_coarse-grained_segement_dmr_C.bed \
  --ref ./mm10.fa \
  --base C \
  --threads 64 \
  --log-filepath ./B_vs_A_singlesite_dmr_C.log

In addition, I also tested modkit v0.5.1 with the following command, and modkit dmr pair still ran successfully without errors afterward:

parallel -j 8 /dist_modkit_v0.5.1_8fa79e3/modkit pileup    ./{}_mm10_sorted.bam     ./{}_mm10.bed.gz --max-depth 100000   --cpg    -t 8    --reference ./mm10.fa    --log-filepath ./log/0.5.1test/{}_mm10_log.txt --region chr9:1-124595110     ::: $SAMPLES

I then examined the previously problematic site chr9:3028084-3028085 and compared the outputs from the three versions:

#0.5.1
chr9    3028084 3028085 h       6202    -       3028084 3028085 255,0,0 6202    2.26    140     761     5301    299     1041    583     0
chr9    3028084 3028085 m       6202    -       3028084 3028085 255,0,0 6202    85.47   5301    761     140     299     1041    583     0

#0.6.1
chr9    3028084 3028085 h       6243    -       3028084 3028085 255,0,0 6243    2.19    137     682     5424    233     981     543     0
chr9    3028084 3028085 m       6243    -       3028084 3028085 255,0,0 6243    86.88   5424    682     137     233     981     543     0

#0.6.0
chr9    3028084 3028085 h       28      -       3028084 3028085 255,0,0 28      5728.57 1604    8303    55657   5039    10932   6284    0
chr9    3028084 3028085 m       28      -       3028084 3028085 255,0,0 28      198775.00       55657   8303    1604    5039    10932   6284       0

I noticed that the reported modification percentages differ slightly across versions but are generally similar. For v0.6.0, after “correcting” the apparent inconsistency, I would infer the ratio should be approximately $N_{mod}/N_{valid} = 55657/65564 = 84.89%$. At the moment, I cannot determine whether downsampling slightly affects the modified base proportions and/or sample comparability. In my tests, v0.6.1 and v0.5.1 both show the same C-base threshold in modkit pileup (since both were run on chr9 only), but v0.6.0 reports a different C-base threshold compared to the other two. I've not done some test in a more strict way.

--max-depth removed. Modkit v0.6.0 has a new algorithm that doesn’t use max-depth. If you have a modBAM with very high depth and you don’t want to tabulate counts for this depth it is currently recommended to subsample the modBAM before using Modkit.

Because I am not yet sure how downsampling might impact downstream analyses and between-sample comparability, I would prefer to tabulate counts for the full depth if possible. For sites that trigger the error, do you think it is feasible to manually repair $N_{valid}$ (for example, by setting it to $N_{mod} + N_{canonical}$), or by re-extracting the true depth at the site directly from the BAM (using some recommended approach)?

I would greatly appreciate your advice and suggestions.

Best regards,
Nathan

[NathaneilKnight]

Hi, @ArtRand,

Thanks again for the detailed explanation and for the v0.6.1 build — that resolved the downstream dmr pair error.

I wanted to follow up because I think now I understand what was happening in v0.6.0. Based on the source, it looks like the pileup counters are stored as u16. In high-depth loci, that means not only Nvalid_cov but also the component counts (Nmodified, Ncanonical, Nother_modified, etc.) can overflow/wrap, which explains why Nvalid_cov can become inconsistent and also why “manually repairing” Nvalid_cov in the bedMethyl output is not reliable (since the other counts may already be incorrect as well).

So I understand that extremely high depth is not currently supported in modkit without downsampling, and it’s not something that can safely be fixed by editing the output file.

Could I ask two related questions?

1. Is there any plan (or would you be open) to switching these counters from u16 to u32 (e.g., in CountsMatrix.inner / PileupFeatureCounts2), at least behind a feature flag or a CLI option?

2. If subsampling is the recommended approach for now, could you share some details about how modkit downsampling the data since I cannot access the source code of 0.5.1 or 0.6.1.

My main goal is to better understand the Modkit's behavior and to preserve between-sample comparability when subsampling is necessary.

Best regards,
Nathan

[ArtRand]

Hello @NathaneilKnight,

You are correct, the bug at high depth is due to unhandled wrapping additions. As a result you cannot manually "repair" the counts.

To your questions:

Is there any plan (or would you be open) to switching these counters from u16 to u32 (e.g., in CountsMatrix.inner / PileupFeatureCounts2), at least behind a feature flag or a CLI option?

Maybe, I need to do some more testing to check the performance impact of this. With v0.6.0 onwards pileup needs to be essentially as lightweight as possible, so I need to weigh the benefits of supporting coverage $\ge 65,535$ over any performance cost.

If subsampling is the recommended approach for now, could you share some details about how modkit downsampling the data since I cannot access the source code of 0.5.1 or 0.6.1.

Sure, in Modkit v0.5.1 pileup uses the htslib "pileup" routines that are also (I believe) used in samtools mpileup. In Modkit v0.6.0 (and onwards) it no longer uses the pileup algorithm from htslib. The algorithm in Modkit v0.6.1 (the build you've tested) downsamples the reads by first checking the number of reads that have been used at a given genomic position. If that number is $\ge$ the maximum depth (specified on the command line), no further counts are added to that position. I'd have to go and double check exactly how htslib does the downsampling. The difference you see is because they aren't exactly equivalent. For what it's worth, Modkit v0.6.1 also agrees with a simple pysam script where you iterate over all of the reads and keep a counter. I'll get the branch code pushed so you can take a look.

I appreciate that you don't want to add a bias unnecessarily. However, if your analysis goal is to look for differentially methylated cytosines or regions then very high depth already comes with some sharp edges. For example, at that depth statistical tests become overpowered so you have to use Cohen's h or other effect size measures. As you increase depth beyond that, I don't know what benefit you really get. I'm not saying there isn't a benefit, maybe you have an idea. What you could do, is separate your modBAM into non-overlapping sets of reads, run v0.6.1 pileup on both, then modkit bedmethyl merge to sum the counts together (the type in merge is wide enough to handle it). I'll continue to think about how to support this use case, thanks for your inquiry and attention to detail.