I would like to ask your opinion about some results I have with Humann3 / MaAsLin2.
I have metatranscriptomic (and no metagenomic) data of 2 sample groups: Control (n=5) and Treatment (n=5). I would like to find out what pathways/functional modules are differentially abundant between groups.
With Humann output, I joined and normalized (relab) the tables. I also splitted the tables to have both stratified and unstratified results. With Maaslin2, I considered only one fixed effect (Diet), with the following code:
fit_data <- Maaslin2( input_data=data, input_metadata=metadata, output='./results', fixed_effects = c('Diet'), min_abundance = 0.0, min_prevalence = 0.0, min_variance = 0.0, reference = 'Diet,CTL')
It runs without errors (just a harmless warning about xtfrm data frames). My issues are related to the result itself. Most of the pathways are negatively modulated, according to “coef” negative values (which is equivalent to fold change, with exponential and log2 conversions, as I saw in another topic).
So, the questions are:
Am I doing something wrong to have such bias towards negative modulation? I have tested both gene families and pathways levels, but the apparent bias happens in both cases. Is there something I could do to test if this is really a bias?
Could this have some relation to community composition? I used a custom taxonomic profile with 16S-derived relative abundances of species. But as I understood, the profile is only used to get reference genomes that are above the abundance threshold, right?
Thank you for your attention!