16 Differential abundance
An important component of many microbiome studies is differential abundance analysis (DAA). It is mainly used to identify differences in the abundances of individual taxa (at any taxonomic level) between two or more groups, such as treatment versus control groups.
The goal of DAA is thus to identify biomarkers of a certain phenotype or condition, and gain understanding of a complex system by looking at its isolated components. For example, the identification of a bacterial taxon that is more abundant in healthy patients compared to diseased patients can lead to important insights into the underlying mechanisms of the disease. In other words, differentially abundant taxa can be involved in the dynamics of the disease, which in turn helps to understand the system as a whole.
16.1 Challenges in differential abundance analysis
As discussed in Section 11.1, microbiome data display unique properties, such as compositionality and skewed/zero-inflated distributions. Evidently, these properties need to be addressed properly also in DAA. Unfortunately, there is currently no consensus on how that should be done. Instead, there exist dozens of methods employing varying techniques and, consequently, providing highly discordant results.
The discordancy of results provided by different DAA methods was recently well illustrated by Nearing et al. (2022). In this study, multiple DAA methods were applied to 38 different datasets and their results were compared to one another. Similar findings were also made by Yang and Chen (2022) who evaluated DAA methods 106 real datasets (along with simulated data). They found especially that on many datasets almost any taxon was detected by at least one DAA method (see Figure 7 in their article). They thus raised concerns of potential cherry-picking, that is, a researcher would be likely to find the taxa in favor of their hypothesis by just trying out several different DAA methods.
16.2 Recommendations based on recent literature
Usually, when a new DAA method is introduced, the authors demonstrate it to perform better compared to some of the existing alternatives. That is often achieved by employing simulated data and possibly one or more real datasets. The authors may have, however, selected their simulation models, example datasets, evaluation metrics, compared methods, etc. specifically to show the introduced method in a favorable light. Assessing the performance of different methods based solely on the studies where they were originally introduced may therefore not be advisable.
A more neutral evaluation of the performance of different DAA methods can be found from several independent benchmarking studies (e.g. Weiss et al. (2017), Hawinkel et al. (2019), Calgaro et al. (2020), Wallen (2021), Nearing et al. (2022), Yang and Chen (2022), Cappellato, Baruzzo, and Camillo (2022), Roche and Mukherjee (2022), Swift et al. (2023), Wirbel et al. (2024), Pelto et al. (2025), Gamboa-Tuz et al. (2025)). Unfortunately, these studies have ended up with very discordant conclusions on the methods. However, if we look specifically at some benchmarking studies that have used only real datasets (instead of simulated data), we can see that analysis of relative abundances (TSS normalized counts) with a simple non-parametric method or linear regression has performed well.
Nearing et al. (2022) studied the performance of DAA methods using 38 datasets. They evaluated the number of detected taxa (ASVs), consistency of the detected taxa between methods, replication of results between studies and error rate on datasets with randomly assigned groups. They emphasized conservatism (low error rate) in their recommendations and ended up recommending ALDEx2 (Fernandes et al. 2014) and ANCOMBC (Mandal et al. 2015). Additionally, they recommended Maaslin2 (Mallick, Rahnavard, and McIver 2020) (i.e. analyzing arcsine square root transformed relative abundances with linear regression) as a more powerful option.
In addition to simulations, Yang and Chen (2022) estimated the FDR of different DAA methods on 106 real datasets with randomly created groups. They found that, the observed FDR was at the nominal level when, for instance, analyzing relative abundances either with the Wilcoxon test or with linear regression (Maaslin2).
Pelto et al. (2025) studied how the results from different DAA methods replicate between random partitions of datasets and between datasets from independent studies. They found that when considering consistency along with sensitivity, the best performance was obtained by analyzing relative abundances with a non-parametric method or (after log transformation) with linear regression/t-test (Maaslin2) or analyzing presence/absence of taxa with logistic regression.
Gamboa-Tuz et al. (2025) employed three real datasets to benchmark different DAA methods against biological ground truth. They found that methods employing the TSS normalization (with e.g. the Wilcoxon test) and RNA-seq-derived methods performed best.
We therefore recommend performing DAA by analyzing relative (TSS normalized) abundances by a non-parametric method (e.g., the Wilcoxon test or ordinal regression model) or, after a log transformation, with t-test/linear regression. In addition to their good performance on real data they have the benefit of being intuitive and easy to understand.
We note, however, that the performance of DAA methods likely depends on the characteristics of the analyzed data (Weiss et al. (2017), Calgaro et al. (2020)). Consequently, while the methods recommended above seem to be good general default methods, they may not necessarily be the optimal methods in many cases. If one wants to evaluate which method may perform best specifically on their data, one option to use the benchdamic package by Calgaro et al. (2022) (see vignette here).
Lastly, based on recent studies on classification accuracy of healthy vs. non-healthy patients, it seems that a high proportion of information in microbiome data is included in the mere presence/absence of taxa (Giliberti (2022) Karwowska et al. (2025)). A viable, and even simpler, alternative to DAA is thus differential prevalence analysis (see also Pelto et al. (2025)).
16.3 Performing differential abundance analysis
Following the recommendations above, we show in this section how to perform DAA by analyzing (log transformed) relative abundances with linear regression models. We first illustrate the analysis in the simplest case where two groups are compared, and no other covariates are included in the analysis. After that we illustrate more complex analyses, namely, the analysis with additional covariate(s) included and the analysis with interaction of two variables. Furthermore, we illustrate how to perform differential prevalence analysis and the differential analysis of positive (non-zero) abundances.
16.3.1 Preparing the data
We use the Tengeler2020 dataset from mia package. The dataset includes samples from people with ADHD (N = 13) and without ADHD (Control group, N = 14). The people are from three different cohorts.
| Cohort_1 | Cohort_2 | Cohort_3 | |
|---|---|---|---|
| ADHD | 4 | 5 | 4 |
| Control | 6 | 5 | 3 |
Before starting the analysis, we agglomerate the features by genus and filter them by a prevalence threshold of 10%. Furthermore, we create assays for relative abundances and presence/absence of taxa. Lastly, we transform variable patient_status to a factor with Control as the reference level. (It is a standard practice to use the control group as the reference group in many statistical analyses.)
# Agglomerate by genus and subset by prevalence
tse <- subsetByPrevalent(tse, rank = "Genus", prevalence = 0.10)
# Transform count data and create assay of relative abundances
tse <- transformAssay(tse, assay.type = "counts", method = "relabundance")
# Transform count data and create assay of presence/absence of taxa
tse <- transformAssay(tse, assay.type = "counts", method = "pa")
# Transform patient_status to a factor and set "Control" as the reference level
tse$patient_status <- factor(tse$patient_status, levels = c("Control", "ADHD"))
# We also transform cohort as a factor
tse$cohort <- factor(tse$cohort, levels = c("Cohort_1", "Cohort_2", "Cohort_3"))16.3.2 Simple two-group comparison
Here we perform the simplest version of DAA, namely, the comparison of two groups without any other covariates included in the analysis. In our case this means comparing the relative abundances of genera between the people with and without ADHD.
We employ the Maaslin2 function from the Maaslin2 package (Mallick, Rahnavard, and McIver 2020). The function
- transforms counts to relative abundances,
- replaces zeros by half of the smallest non-zero relative abundance, (separately for each taxon),
- applies log2 transformation to the relative abundances, and
- runs a linear regression model for each taxon using the transformed relative abundances as the response variable.
To run the function we need to provide it with the count matrix and metadata as inputs. We also need to provide it with the variable indicating the compared groups, i.e. patient_status in our case. Furthermore, we need to specify the folder for some output files. Lastly, we prevent Maaslin2() to create certain plots. (Note that Maaslin2() prints a lot of some output which is suppressed here for convenience.)
knitr::opts_chunk$set(eval = FALSE)library(Maaslin2)
obj_daa_basic <- Maaslin2(
input_data = assay(tse), # The count matrix
input_metadata = colData(tse) |> as.data.frame(), # The metadata
fixed_effects = "patient_status", # The predictor variable
output = "output", # The name of the folder for output files.
plot_heatmap = FALSE,
plot_scatter = FALSE
)Next, we extract the data.frame of results from the list object returned by Maaslin2(). We also calculate the 95% confidence intervals and transform the results to a slightly more convenient format. (The sample size N is used as an approximation for the residual degrees of freedom (df) in the calculation of confidence intervals.)
# Extract the results from the list object and calculate the 95% CIs
res_daa_basic <- obj_daa_basic$results |>
mutate(
ci_lwr = coef - qt(.975, df = N) * stderr,
ci_upr = coef + qt(.975, df = N) * stderr
) |>
select(feature, name, coef, pval, qval, ci_lwr, ci_upr)
# Print the results for taxa with the smallest p-values
res_daa_basic |>
arrange(pval) |>
head() |>
kable()The columns of res_daa_basic are as follows:
-
featureindicates the name of a taxon (here genus). -
nameindicates the predictor variable. Aspatient_statusis a categorical variable, it also tells the group used as the non-reference level (hereADHDasControlwas set as the reference level). -
coefis the log2-fold change estimate, a standard way to express effect size in DAA. Most importantly, positive values indicate abundance being higher in the non-reference group (hereADHD). More precisely,2 ^ coefwould be interpreted as the ratio of geometric means of the relative abundances between the ADHD and Control groups. -
pvalis the unadjusted p-value. -
qval(q-value) is the FDR adjusted p-value. This is generally used to define statistical significance. Typically, q-values < .05 are considered (statistically) significant. -
ci_lwrandci_uprare the lower and upper bounds of the 95% confidence intervals for the log2-fold change.
Now, if we use the standard criterion for statistical significance, namely, q-value < 0.05, we can see that the DAA result is statistically significant for the genera named as [Ruminococcus] gauvreauii group, Faecalibacterium, Catabacter, [Clostridium] innocuum group and Erysipelatoclostridium. (Maaslin2 has some problems with special characters such as [ or ].) Moreover, for instance, the positive log2-fold change estimate of Faecalibacterium tells that the average relative abundance of Faecalibacterium is higher in the ADHD group.
16.3.2.1 Visualizing the results
While the numbers above provide summarized information about the group differences for each taxon, it is often useful to further examine these associations in more detail. That can be done, for instance, using strip plots with boxplots. We show here how such plots can be created for the significant taxa by using plotBoxplot function from the miaViz package.
library(miaViz)
# Create strip plots with boxplots (and prevalence bars)
plotBoxplot(
tse,
features = c(
"[Ruminococcus]_gauvreauii_group",
"Faecalibacterium",
"Catabacter",
"Erysipelatoclostridium",
"[Clostridium]_innocuum_group"
),
assay.type = "relabundance", # To plot relative abundances
x = "patient_status", # The variable on the x (horizontal) axis
facet.by = "rownames", # To create a separate facet for each feature
scales = "free_y", # To allow different y-axis scales
add.proportion = TRUE, # Add bars indicating prevalence
point.offset = "quasirandom", # To plot points nicely
point.size = 2, # The size of the points
point.color = "black", # The color of the edges of the points
point.alpha = 1, # The transparency of the edges of the points
box.alpha = .20
) + # The transparency of the boxplots
theme(strip.text = element_text(size = 7)) # The size of the facet labelsIn each facet,
- the thicker horizontal line in each boxplot indicates the median,
- each box indicates the interquartile range (IQR), namely, the range between the 25th (Q1) and 75th percentile (Q3),
- the lower whisker indicates the smallest value above \(Q1 - 1.5 \times IQR\), and
- the lower whisker indicates the largest value below \(Q3 + 1.5 \times IQR\).
- The horizontal bars at the bottom indicate the prevalence of each genus in each group.
The strip plots and boxplots clearly show, for instance, how the relative abundance of Faecalibacterium is higher in the ADHD group, as was indicated by the positive the log2-fold estimate (coef) in the results table above. Moreover, for instance, the prevalence bars of Erysipelatoclostridium show how Erysipelatoclostridium is more prevalent in the control group.
Lastly, we show how the results for all genera can be summarized visually by plotting the log2-fold estimates and their 95% confidence intervals.
# First arrange the genera by log2-fold estimate
top_features <- res_daa_basic[order(res_daa_basic[["coef"]], decreasing = TRUE), "feature"]
res_daa_basic[["genus_ord"]] <- factor(res_daa_basic[["feature"]], levels = top_features)
# Plot the log2-fold changes and the 95% confidence intervals
ggplot(res_daa_basic, aes(x = coef, y = genus_ord)) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_point() +
geom_errorbarh(aes(xmin = ci_lwr, xmax = ci_upr), height = 0.2) +
labs(x = "Log2-fold change (95% CI)") +
theme_light() +
theme(
axis.title.y = element_blank(),
axis.text.y = element_text(size = 7),
plot.title = element_text(hjust = 0.5)
)16.3.3 DAA with additional covariates
Here we show how to perform DAA with an additional covariate included in the analysis. One usually wants to include a variable as a covariate if they suspect that it may act as a confounding variable, namely, a variable that causes a spurious association between the variable of main interest and the abundances.
We use here cohort as a covariate. We proceed as above except that we now set fixed_effects = c("patient_status", "cohort"). Furthermore, we now calculate the q-values manually so that they are calculated separately for each variable in the model (indicated by the variable metadata in obj_daa_cov$results). This is because Maaslin2() calculates the q-values over all variables in the model which is not what here as we are only interested in the variable patient_status.
# Run Maaslin2 with patient_status and cohort as predictor variables
obj_daa_cov <- Maaslin2(
input_data = assay(tse),
input_metadata = colData(tse) |> as.data.frame(),
fixed_effects = c("patient_status", "cohort"),
output = "output",
plot_heatmap = FALSE,
plot_scatter = FALSE
)
# Extract the results and calculate the 95% CIs and q-values
res_daa_cov <- obj_daa_cov$results |>
mutate(
ci_lwr = coef - qt(.975, df = N) * stderr,
ci_upr = coef + qt(.975, df = N) * stderr
) |>
group_by(metadata) |>
mutate(qval = p.adjust(pval, method = "fdr")) |>
ungroup() |>
select(feature, name, coef, pval, qval, ci_lwr, ci_upr)Note that there now appear results also for the variable cohort. As cohort consists of three groups, namely, Cohort_1 (the reference group), Cohort_2 and Cohort_3, there are results separately for the comparison between Cohort_2 vs. Cohort_1 (cohortCohort_2) and for the comparison between Cohort_3 vs. Cohort_1 (cohortCohort_3).
We can now examine how the inclusion of cohort as a covariate affects the results for patient_status.
When comparing these results to the results in Section 16.3.2, we can see that, for instance, the q-values of Faecalibacterium and Catabacter are now somewhat smaller. Nevertheless, the results look overall rather similar to the results without the covariate included. Therefore, cohort does not seem to be a confounder causing spurious associations between patient status and the abundances of any taxa.
16.3.4 DAA with an interaction
Here we show a more complex version of DAA, namely, a case where the interaction effect of two variables is analyzed. We give this example here for illustration purposes but note that such complex analysis may not be advisable in a small dataset such as our example dataset.
We analyze the interaction between patient_status and cohort. We proceed as in Section 16.3.3 above except that we now add also the interaction term into the formula. As Maaslin2() does not directly support the analysis with interactions we must first create the interaction dummy variables in to the tse object. The main effects patient_status and cohort as well as the interaction dummy variables are then included in the Maaslin2() function.
# Create dummy variables for the interaction term
tse$patient_statusADHD_cohortCohort_2 <-
ifelse(tse$patient_status == "ADHD" & tse$cohort == "Cohort_2", 1, 0)
tse$patient_statusADHD_cohortCohort_3 <-
ifelse(tse$patient_status == "ADHD" & tse$cohort == "Cohort_3", 1, 0)
# Run Maaslin2 with patient_status, cohort and their interaction terms
obj_daa_ia <- Maaslin2(
input_data = assay(tse),
input_metadata = colData(tse) |> as.data.frame(),
fixed_effects = c(
"patient_status",
"cohort",
"patient_statusADHD_cohortCohort_2",
"patient_statusADHD_cohortCohort_3"
),
output = "output",
plot_heatmap = FALSE,
plot_scatter = FALSE
)
# Extract the results and calculate the 95% CIs and q-values
res_daa_ia <- obj_daa_ia$results |>
mutate(
ci_lwr = coef - qt(.975, df = N) * stderr,
ci_upr = coef + qt(.975, df = N) * stderr,
metadata = ifelse(
metadata %in% c(
"patient_statusADHD_cohortCohort_2",
"patient_statusADHD_cohortCohort_3"
),
"patient_satus:cohort",
metadata
)
) |>
group_by(metadata) |>
mutate(qval = p.adjust(pval, method = "fdr")) |>
ungroup() |>
select(feature, name, coef, pval, qval, ci_lwr, ci_upr)The results object res_daa_ia now also includes the interaction terms patient_statusADHD_cohortCohort_2 and patient_statusADHD_cohortCohort_3.
We examine the term patient_status_cohortCohort_3. It indicates how the differences between the ADHD and Control groups vary between cohorts 1 and 3.
The smallest p-value is observed for Coprobacter. As the corresponding log2-fold change estimate is negative it means that the difference in the relative abundance of Coprobacter between the ADHD and Control groups is higher in Cohort 1 than in Cohort 3.
This interaction can be better illustrated by plotting the relative abundances of Coprobacter by patient_status and cohort.
# Plot the abundances of Coprobacter by patient_status and cohort
plotBoxplot(
tse,
features = c("Coprobacter"),
assay.type = "relabundance", # To plot relative abundances
x = "patient_status", # The variable on the x (horizontal) axis
facet.by = "cohort", # To create a separate plot for each cohort
scales = "fixed", # To have the same y-axis scale in each facet
add.proportion = TRUE, # Add bars indicating prevalence
point.offset = "quasirandom", # To plot points nicely
point.size = 2, # The size of the points
point.color = "black", # The color of the edges of the points
point.alpha = 1, # The transparency of the edges of the points
box.alpha = .20
) # The transparency of the boxplotsHere we can clearly see how the relative abundances of Coprobacter are higher in the ADHD group in Cohort 1 but they are higher in the Control group in Cohort 3. As the association is positive in Cohort 1 (the baseline cohort) and negative in Cohort 3, sign of interaction estimate is negative (the association becomes more negative when we move from the baseline cohort to a non-baseline cohort).
Note that a similar plot could also be used to illustrate a possible confounding effect.
knitr::opts_chunk$set(eval = TRUE)16.3.5 Differential prevalence analysis (DPA)
We show here how to perform DPA with logistic regression. We show the case where we compare subjects with and without ADHD while controlling for cohort. Due to the convenient implementation of logistic regression for DPA, we employ here the maaslin3() function from the maaslin3 package (Nickols et al. 2026).
To run the maaslin3() function, we need to provide it with
- a
tseobject, - a model formula, in this case
~ patient_status + cohort + total, and - a folder for output files
Note that we also control for the sequencing depth/library size/total read count (total) in the analysis. This is because the detection (observed presence) of taxa may depend on the sequencing depth.
library(scater)
library(maaslin3)
# Calculate the sequencing depth (the total number of reads/counts) for each
# sample
tse <- addPerCellQC(tse)
# Run maaslin3 with patient_status, cohort and reads as the predictor variables
obj_DPA_cov <- maaslin3(
tse,
formula = ~ patient_status + cohort + total, # Model formula
output = "output", # The name of the folder for some output files
plot_associations = FALSE, # To prevent some unnecessary plots
plot_summary_plot = FALSE,
verbosity = "ERROR"
) # To prevent some unnecessary outputNext, we extract the data.frame of results from the list object returned by maaslin3(). We remove the variable total from the output as it is merely a technical covariate that we are not interested in. We also calculate the 95% confidence intervals. Furthermore, as with Maaslin2, we calculate the q-values separately for patient_status and cohort. Lastly, we transform the results to a slightly more convenient format.
# Extract the results and calculate the 95% CIs and q-values
res_DPA_cov <- obj_DPA_cov$fit_data_prevalence$results |>
filter(name != "total") |>
mutate(
ci_lwr = coef - qnorm(.975) * stderr,
ci_upr = coef + qnorm(.975) * stderr
) |>
group_by(metadata) |>
mutate(qval = p.adjust(pval_individual, method = "fdr")) |>
ungroup() |>
select(feature, name, coef, pval = pval_individual, qval, ci_lwr, ci_upr)
# Print the results for taxa with the smallest p-values for patient_status
res_DPA_cov |>
filter(name == "patient_statusADHD") |>
arrange(pval) |>
head() |>
kable()The results are now in the same format as the DAA results in Section 16.3.3. The only difference is coef now indicates the log-odds ratio for the difference in prevalences. The odds ratios (OR) would thus be calculated as or = exp(coef). Note also that the results are not provided for taxa that were present in all samples (as it is not sensible to analyze presence/absence differences for such taxa).
By examining the results we see that the smallest p-value was observed for genus Erysipelatoclostridium. As the sign of the estimate is negative is means that the genus is more prevalent in the control group. (This is can be seen from the prevalence bars of Erysipelatoclostridium in the strip plots in Section 16.3.2.1.)
16.3.5.1 DPA with an interaction
Contrary to Maaslin2(), the maaslin3() function supports also the analysis with interaction terms. Therefore, in the analysis with interaction(s), one can just add the interaction term(s) in to the model formula by using the standard R syntax, e.g., formula = ~ patient_status + cohort + patient_status:cohort + reads or alternatively formula = ~ patient_status * cohort + reads
16.3.6 Analysis of positive (non-zero) abundances
It may sometimes be biologically relevant to analyze only the positive (non-zero) relative abundances of taxa. We thus briefly show here how to perform such analysis using the maaslin3() function. Purely for illustration purposes, we show here the analysis with an interaction term.
We proceed as in the previous sections. The difference compared to DPA (Section 16.3.5) is that we now have to extract fit_data_abundance (and not fit_data_prevalence) from the list object returned by maaslin3(). We also need to set median_comparison_abundance = FALSE in the function call as we want to estimate the differences of (positive) relative abundances and not make any compositional adjustments.
# Run maaslin3 function.
obj_daa_pos_ia <- maaslin3(
tse,
output = "output",
formula = ~ patient_status * cohort,
median_comparison_abundance = FALSE,
verbosity = "ERROR",
plot_associations = FALSE,
plot_summary_plot = FALSE
)
# Extract the results and calculate the 95% confidence intervals and q-values
res_daa_pos_ia <- obj_daa_pos_ia$fit_data_abundance$results |>
filter(!is.na(pval_individual) & N_not_zero > 6) |>
mutate(
ci_lwr = coef - qt(.975, df = N_not_zero - 6) * stderr,
ci_upr = coef + qt(.975, df = N_not_zero - 6) * stderr
) |>
group_by(metadata) |>
mutate(qval = p.adjust(pval_individual, method = "fdr")) |>
ungroup() |>
select(feature, name, coef,
pval = pval_individual, qval, ci_lwr, ci_upr,
N_not_zero
)
# The results for the interaction between patient_status and cohort (when
# comparing Cohort_3 to Cohort_1) for the taxa with the smallest p-values
res_daa_pos_ia |>
filter(name == "patient_statusADHD:cohortCohort_3") |>
arrange(pval) |>
head() |>
kable()The output is similar to the outputs in the previous sections. It, however, now includes also N_not_zero indicating the number of samples with positive (non-zero) relative abundances, that is, the sample size used in the analysis for each taxon.
The results are interpreted in a similar manner as in normal DAA (Section 16.3.4). For instance, the estimates (coef) indicate the log2-fold changes of positive relative abundances (as maaslin3() performs log2 transformation to the positive relative abundances).
16.3.7 Alternative options
16.3.7.1 Advanced compositionality adjustments
We have shown above (sections 16.3.2, 16.3.3, 16.3.4 and 16.3.6) how to perform DAA by analyzing relative abundances (TSS normalized counts) without any advanced adjustments for compositionality. As explained in Section 16.2, this choice was motivated by some recent benchmark studies that have employed only real data, especially the study by Gamboa-Tuz et al. (2025).
We note, however, that in some benchmark studies such simple versions of DAA have performed worse to some more advanced methods when compared on simulated data with strong compositional effects (e.g. Roche and Mukherjee (2022), Yang and Chen (2022)). We therefore give below some options for incorporating more advanced compositionality adjustments in DAA.
- In Maaslin2, the default normalization method (TSS) can be replaced by TMM normalization (Robinson and Oshlack 2010) by including
NORM = "TMM"in theMaaslin2()function. - In the analysis of positive abundances with maaslin3, the adjustment for compositional bias can be enabled by setting
median_comparison_abundance = TRUEinmaaslin3()function. - One can use DAA methods that are specifically designed to account for the compositionality of microbiome data. Such methods include, for instance,
- ANCOMBC (Lin and Peddada 2020), vignette. Note that we refer here to the older version of ANCOM-BC. There is also a newer version (ANCOM-BC2) available, but some problems have been observed in its performance (e.g., Pelto et al. (2025)).
-
LinDA(Zhou et al. 2022) from the MicrobiomeStat package, reference manual. -
ZicoSeq(Yang and Chen 2022) from the GUniFrac package, vignette.
16.3.7.2 A non-parametic alternative to MaAsLin2
As MaAsLin2 analyzes log transformed relative abundances, it needs to replace the zeros by some positive numbers before the transformation. If one wants to avoid this somewhat arbitrary imputation, we suggest performing DAA by using non/semi-parametric ordinal regression models (ORM). ORM can be seen as an extension of the familiar Wilcoxon test in the sense that it provides p-values similar to those of the Wilcoxon test in the simple two group comparison, but it can also incorporate covariates (and interaction terms) (Harrell (2001), Liu et al. (2017)). ORM has also performed well in one benchmark study (Pelto et al. 2025).
As there does not currently exist a good implementation of ORM for DAA one must use a custom function. Such custom function with application examples is given in Appendix C. We also note that, in contrast to Maaslin2() function, the custom function provides confidence intervals and q-values (calculated for the variable of interest) automatically in its output.
16.3.7.3 A more powerful alternative to MaAsLin3 in DPA
While the maaslin3 package provides a convenient implementation of logistic regression for microbial DPA, we note, that it may provide somewhat conservative results, that is, it rarely provides false positive results but it may fail to detect some true effects. To obtain higher statistical power, we recommend using a version of logistic regression introduced by Firth (1993). As there currently exists no convenient implementation of Firth logistic regression for DPA, one needs to use a custom function. Such custom function with application examples is given in the Appendix C.