`pileup` command is giving BED files that `dmr` command says `invalid-bedmethyl-data`.

[ArnavBharti]

pileup command is giving BED files that dmr command says invalid-bedmethyl-data. I have attached the output of DMR and the commands I have run.

OUTPUT (DMR)

> reading reference FASTA at "<genome_fasta_name>.fasta"
> 1521 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 1094523 a records and 905404 b records, calculating max coverages for 95th percentile
> calculated max coverage for a: 6 and b: 21
> running with replicates, but not matched samples
> batch failed: invalid data, valid coverage (1) is not equal to the sum of canonical and modified counts (0), [BedMethylLine { chrom: "AAKM01000036", interval: Interval { start: 7056, stop: 7057, val: () }, raw_mod_code: Code('h'), strand: Positive, count_methylated: 0, valid_coverage: 1, count_canonical: 0, count_other: 1, count_delete: 0, count_fail: 0, count_diff: 0, count_nocall: 0 }] chrom: AAKM01000036 starting at 7056, stopping
> 12 batches processed                                                                                                                                                                 > Error! invalid-bedmethyl-data

PILEUP LOGS

[modkit-logging/src/lib.rs::69][2026-02-07 10:41:20][DEBUG] command line: modkit pileup J_CC_5hmC.bam J_CC_5hmC.bed.gz --modified-bases 5hmC --reference /home/<genome_fasta_name>.fasta --log J_CC_5hmC.log --bgzf
[modkit-core/src/pileup/subcommand.rs::802][2026-02-07 10:41:20][INFO] discarded 1107 contigs with zero aligned reads
[modkit-core/src/pileup/subcommand.rs::576][2026-02-07 10:41:20][INFO] parsed 1 base modification(s). Base modifications other than 'C:h' will be counted as 'N_other'.
[modkit-core/src/pileup/subcommand.rs::653][2026-02-07 10:41:20][INFO] adding single-base motif: 'C 0'
[modkit-core/src/pileup/subcommand.rs::944][2026-02-07 10:41:20][INFO] using bgzf compression with 4 compression threads
[modkit-core/src/command_utils.rs::102][2026-02-07 10:41:20][INFO] attempting to sample 10042 reads
[modkit-core/src/reads_sampler/mod.rs::49][2026-02-07 10:41:20][DEBUG] found BAM index, sampling reads in 1000000 base pair chunks
[modkit-core/src/reads_sampler/sampling_schedule.rs::261][2026-02-07 10:41:20][DEBUG] removed 0 contigs from schedule with <= 1 reads
[modkit-core/src/reads_sampler/sampling_schedule.rs::141][2026-02-07 10:41:20][DEBUG] derived sampling schedule, sampling total 11547 reads from 1641 contigs, 0 unmapped reads
[modkit-core/src/reads_sampler/sampling_schedule.rs::164][2026-02-07 10:41:20][DEBUG] schedule
 chrom  count/frac

DMR Command

modkit dmr pair \
  -a UNC_ED_5hmC.bed.gz \
  -a UNC_YR1_5hmC.bed.gz \
  -a UNC_YR2_5hmC.bed.gz \
  -a UNC_EL_5hmC.bed.gz \
  -a UNC_FY_5hmC.bed.gz \
  -b CM_EL_5hmC.bed.gz \
  -b CM_EJ_5hmC.bed.gz \
  -o DMR_UNC_vs_CM_5hmC.bed \
  --ref "${REF}" \
  --base C \
  --threads "${THREADS}" \
  --log-filepath DMR_UNC_vs_CM_5hmC.log

PILEUP Command. mod = 5hmC

pileup () {
    local bam=$1
    local mod=$2
    local prefix=$3
    log "Running modkit pileup (${mod}) for ${prefix}"
    modkit pileup \
        "${bam}" \
        "${prefix}_${mod}.bed.gz" \
        --modified-bases "${mod}" \
        --reference "${REF}" \
        --log "${prefix}_${mod}.log" \
        --bgzf
}

I had run modified basecalling for 5mC 5hmC to create modified BAM files. @ArtRand

[Proy321]

Hello @ArtRand

It would be nice if you can help us with the same above query.

Thanks

[Proy321]

Hello @ArtRand

Gentle reminder for the above query.

Thanks

[msfalmeida]

Hi.
I am having a similar issue when running dmr pair to compare two haplotype bedmethyl files generated using pileup with --phased option using the most current release of the software.

In my case I think either a differently modified base or a deletion may be causing confusion as my log shows:

encountered invalid bedMethyl record(s), invalid data, valid coverage (14) is not equal to the sum of canonical and modified counts (13), [BedMethylLine { chrom: "3", interval: Interval { start: 3144586
, stop: 3144587, val: () }, raw_mod_code: Code('m'), strand: Both, count_methylated: 3, valid_coverage: 14, count_canonical: 10, count_other: 1, count_delete: 1, count_fail: 0, count_diff: 0, count_nocall
: 1 }] chrom: 3 starting at 3144586, stopping

Let me know if more details would be useful or how I can solve this issue.

Thanks

[ArtRand]

Hello @ArnavBharti and @msfalmeida,

Sorry about the delay.

The issue is the for DMR you need have a bedMethyl with all of the modifications explicitly recorded. So if you ran basecalling with 5mC and 5hmC you need to generate the bedMethyl with --modified-bases 5mC 5hmC for it to work with the current version of Modkit DMR.

@msfalmeida could you paste the exact command you're using?

[msfalmeida]

Dear @ArtRand,
Thank you for your quick reply.

I have basecalled both 5mC and 5hmC but then run modkit pileup for 5mC alone as follows:

/home/bioc1179/Software/dist_modkit_v0.6.1_481e3c9/modkit pileup $alignment_bam $output \ --cpg \ --modified-bases 5mC \ --combine-strands \ --reference $reference \ --phased \ --log modkit_pileup_log.txt \ --include-bed /home/bioc1179/Berrens/Resources/chr3_chrXenrichment.bed \ --threads 8

I have sorted, bgzip and tabix the output bed files.

I then attempted dmr pair like this:

/home/bioc1179/Software/dist_modkit_v0.6.1_481e3c9/modkit dmr pair \ -a C7H8_NPCd28_mm39.pileup.hp1.sort.bed.gz \ -b C7H8_NPCd28_mm39.pileup.hp2.sort.bed.gz \ -o ${dmr_result} \ -r $regions \ --ref $reference \ --base C \ --threads 8 \ --min-valid-coverage 5\ --log-filepath dmr.log

I have also tried as seen in another issue entry to start from the separated bam for the 2 haplotypes but it return the same error for the same position on chromosome 3 as I shown in my previous message.

So if I want to try to re-run pileup with both 5mC and 5hmC, how can I specify in dmr pair that I am only interested in 5mC?

Thanks.

[ArtRand]

@msfalmeida

Correct, you need to run Modkit pileup with --modified-bases 5mC 5hmC. There is no way in Modkit DMR to specify that you're only interested in 5mC. However, the MAP-based p-value, effect size, and Cohen's statistics are calculated assuming a "modified vs un-modified" model. I'm actively working on an improvement so you can specify that you're looking for changes in only one modification (or within modification types). One thing you could try is using dmr normally, then filter out regions/sites where the 5mC percentages/counts between conditions are similar, I appreciate this requires a bunch of awkward parsing etc. thus why I'd rather have a better solution.

[msfalmeida]

@ArtRand thank you for the suggestion and for your efforts in developing modkit.

From a biology perspective, I think it would make sense to look at the modifications separately. Until you come up with a solution for specifying it in DMR I wonder if it's possible to basecall 5mC only, but I don't remember that being an option. Do you know if that is possible?

[ArtRand]

@msfalmeida

Ok let me give it a think. An easy thing would be to combine the 5hmC counts with the unmodified counts. Then a site/region which changes 5hmC concentration between two conditions would not be different as far as the metrics are concerned. Obviously, I haven't given this as much consideration as is likely required. When you run with the normal settings, do you find regions/sites that are different but not due to changes in 5mC?