`modkit pileup` --bgzf but output not compressed

[jkh00]

Hii,

As implied in the title, I noticed when running modkit pileup with the --bgzf option, the outout bedmethyl was, however not compressed.

This was the command i ran:

modkit \
    pileup \
    --motif CHH 0 --motif CHG 0 --motif CG 0 --modified-bases 6mA --combine-mods --log-filepath pileup.log \
    --bgzf \
    --bgzf-threads 16 \
    --threads 16 \
    --prefix test \
    --ref $ref \
    $inbam \
    test.bed.bgz

and I checked the output using file test.bed.bgz, it shows as ASCII text. Not sure if I have written the params wrong or was this a bug?

Thanks again in advance!

[jkh00]

update:

when I ran pileup for only a single base, the compression works:

modkit \
    pileup \
    --modified-bases 6mA --combine-mods --log-filepath pileup.log \
    --bgzf \
    --bgzf-threads 32 \
    --threads 32 \
    --prefix test \
    --ref $ref \
    $inbam \
    test.bed.bgz

in the log file:

[modkit-logging/src/lib.rs::69][2025-12-19 14:01:40][DEBUG] command line: modkit pileup --modified-bases 6mA --combine-mods --log-filepath 6mA_pileup.log --bgzf --bgzf-threads 32 --threads 32 --prefix test --ref merged_genomic.fna aligned.bam 6mA_pileup.bed.bgz
[modkit-core/src/pileup/subcommand.rs::743][2025-12-19 14:01:40][INFO] discarded 2 contigs with zero aligned reads
[modkit-core/src/pileup/subcommand.rs::507][2025-12-19 14:01:40][INFO] parsed 1 base modification(s), combining modification counts together for primary sequence bases 'A'
[modkit-core/src/pileup/subcommand.rs::591][2025-12-19 14:01:40][INFO] adding single-base motif: 'A 0'
[modkit-core/src/pileup/subcommand.rs::855][2025-12-19 14:01:40][INFO] using bgzf compression with 32 compression threads
[modkit-core/src/command_utils.rs::102][2025-12-19 14:01:40][INFO] attempting to sample 10042 reads

I also checked the bedmethyl output, it was compressed.
and when I ran for all C contexts together with 6mA (the command in the first msg), this is the log I got:

[modkit-logging/src/lib.rs::69][2025-12-19 13:21:40][DEBUG] command line: modkit pileup --motif CHH 0 --motif CHG 0 --motif CG 0 --modified-bases 6mA --combine-mods --log-filepath pileup.log --bgzf --bgzf-threads 16 --threads 16 --prefix test --ref merged_genomic.fna aligned.bam test.bed.bgz
[modkit-core/src/pileup/subcommand.rs::743][2025-12-19 13:21:40][INFO] discarded 2 contigs with zero aligned reads
[modkit-core/src/pileup/subcommand.rs::486][2025-12-19 13:21:40][INFO] more than one motif requires use of general pileup processor
[modkit-core/src/pileup/subcommand.rs::846][2025-12-19 13:21:40][DEBUG] using multiple-motif writer
[modkit-core/src/command_utils.rs::102][2025-12-19 13:21:40][INFO] attempting to sample 10042 reads

[ArtRand]

Hello @jkh00,

I think I know where the problem is, let me get a fix over to you. I'm curious about your first command:

modkit \
    pileup \
    --motif CHH 0 --motif CHG 0 --motif CG 0 --modified-bases 6mA --combine-mods --log-filepath pileup.log \
    --bgzf \
    --bgzf-threads 16 \
    --threads 16 \
    --prefix test \
    --ref $ref \
    $inbam \
    test.bed.bgz

In your workflow is it easier to use one command and get one output than to split the motifs into different commands/files? Something like this:

modkit \
    pileup \
    --modified-bases 6mA [--combine-mods] \
    --log-filepath pileup.log \
    --bgzf \
    --bgzf-threads 8 \
    --threads 8 \ # you probably don't see a big speedup with 32 threads
    --ref $ref \
    $inbam \
    test_6mA.bed.bgz
for motif in "CHH" "CHG" "CG"; do
  modkit \
      pileup \
      --motif ${motif} 0 --modified-bases C --combine-mods --log-filepath pileup.log \
      --bgzf \
      --bgzf-threads 8 \
      --threads 8 \
      --ref $ref \
      $inbam \
      test_${motif}.bed.bgz
done

I'm mostly curious since I imagine that you'd want to split up the records by base modification/motif anyways. Now that the pileup is quicker maybe running it for each motif isn't much or any slower.

As for why this isn't working: multiple motifs like in your first command don't use the optimized pileup routine, they fall back to something that's quicker than ~v0.5.1 but not quite as quick as the faster ones. It's this "fall-back" logic that isn't correct and is the reason your output isn't compressed. I'll get that fixed - thanks for the report.

[jkh00]

thanks @ArtRand for the thorough comment!

Regarding to the command running practice, I'm currently trying to incorporate modkit pileup into a nextflow pipeline, methylong.

My initial idea, and what's being implemented currently, is to run the pileup step only once (only CG/ all 5mC contexts/ with 6mA). This way, we end up with a single bedmethyl file, and then any downstream post processing can be more user specific. For example, creating a context-specific bedgraph with a min. base coverage filter. There's no strict rule on how this should be done, though, and thinking about it more, I agree that your idea of splitting motifs makes more sense, especially for 6mA and 5mC. That said, for the 5mC contexts, I would probably stick to piling them up all at once, as it’s easier to incorporate into the current pipeline structure.

I also ran into another issue while testing. when i run pileup with --motif CHH 0 --motif CHG 0 --motif CG 0 --modified-bases 6mA, I don't see any 6mA pileup in the resulting bedmethyl, which is kinda strange. However, if I run pileup with only --modified-bases 6mA, it works as expected. Not sure if this is another bug.

Looking forward to the fix!