moshpit pipeline differential abundance analysis

llenzi · June 26, 2025, 5:24pm

Hi all,

I am working on a metatrascriptomic project and for that I am testing moshpit (2025.4).

For now I followed the tutorial to assign taxonomy to reads, and then I used ancom-bc2 to identify differentially abundant features. However, talking with my boss he asked me to avoid the bracken step, because we are not sure on how it transforms the counts in the case of RNA sequences. So my question is, which is the easiest way to get the kraken report into a feature table? I notice there is a 'mosh annotate kraken2-to-mag-features' for the MAG pipeline, is there anything similar for the read assignment pipeline?
Thanks
Luca

Nicholas_Bokulich · July 1, 2025, 6:07am

Hi @llenzi ,
Yes there is the kraken2-to-features action, though this outputs a presence-absence table. Would that work for you?

Just curious: what is your boss's rationale for not trusting bracken with RNA sequences?

llenzi · July 1, 2025, 8:36am

Hi @Nicholas_Bokulich
thanks for your reply, unfortunately presence-absence won't work in my case because I need to perform differential abundance testing.

I was looking for some ref for potentially bracken problems on trancriptomics, I will ask t o my boss a better match. But my understanding is that is concern is linked to this thread:

github.com/jenniferlu717/Bracken

Bracken greatly inflating total bacterial reads

opened 05:51AM - 20 Jul 23 UTC

jpmcelderry

I've noticed an oddity in my samples when using Bracken v2.8. For some samples B…racken works exactly as expected, but for others it seems to randomly inflate the total number of assigned reads 1.5-2 fold, all the way from the domain to the species level. For instance in one sample the total number of reads assigned to Bacteria jumped from 291,400 with Kraken to 540,261 after Bracken, and only 225 viral reads to 356,440. - Kraken 98.01 8647798 541788 R 1 root 91.84 8103509 5971127 R1 131567 cellular organisms 20.86 1840382 355476 D 2759 Eukaryota 3.30 291400 38208 D 2 Bacteria 0.01 600 57 D 2157 Archaea 0.00 225 0 D 10239 Viruses - Bracken 100.00 8644372 0 R 1 root 95.87 8287368 0 R1 131567 cellular organisms 89.62 7746732 0 D 2759 Eukaryota 6.25 540261 0 D 2 Bacteria 0.00 374 0 D 2157 Archaea 4.12 356440 0 D 10239 Viruses In some samples, these reads are distributed so unevenly that a species with relatively low abundance suddenly becomes the most abundant species in the sample entirely. I believe what's happening is that reads from the "root" and "cellular organisms" levels are being distributed, but my understanding is that reads are meant to be distributed only from intermediate levels. Is this a bug, or intended?

The main point being that bracken normalize on kmers, which are obtained from the genome and assumed to be equally present in the samples. However this may not be true in the case of transcriptome.
Hope it makes sense, but any suggestion counter point is largely appreciated!
Cheers
Luca

Nicholas_Bokulich · July 1, 2025, 9:35am

Hi @llenzi ,
Thanks for clarifying! Unfortunately it looks like the bracken devs never replied to that issue — would be interesting to get their input.

kraken2-to-features does not output a frequency table because — at least in the context of metagenomics the raw frequencies are skewed by genome size and cannot be used as such. In the context of metatranscriptomics I would also expect the raw counts to be highly skewed. But it could be interesting to expose, to allow correction with other methods. We are always open for contributions if you want to contribute to Q2/moshpit to expand options for metatranscriptomics

system · August 1, 2025, 3:36pm

This topic was automatically closed 31 days after the last reply. New replies are no longer allowed.