Q-values in Maaslin2

Downstream analysis and statistics MaAsLin
1

looking at the results file of my data analysis (GO terms molecular functions in copies-per-million [880 functions], 8 control individuals and 8 exposed ones, exposure is the only variable), I am not sure if p.adjust function performs as expected:

  1. for correction = “fdr” I get an error - “ERROR::Please select a correction method from the list of available options : BH, holm, hochberg, hommel, bonferroni, BY
    Error: Option not valid”;
  2. when I use correction = “BH” or do not use correction argument (default is BH), all q - values are 0.73 regardless of p-values that range from 0.004 to 0.999;
  3. when I use correction = “BY” option, all q-values = 1.

Is this a possible outcome of the p-value correction or something is wrong?

2

  1. Maaslin2 does a bit of pattern matching to allow flexible capitalization of the correction argument, but “fdr” doesn’t resolve to one of the allowed choices like it does with p.adjust(). “fdr” resolves to “BH” inside of p.adjust() will give you the same thing as what you originally intended.

  2. Depending on your distribution of p-values, it is possible to have the adjusted q-values collapse to a smaller number of unique values. You can see this independently via a one-line simulation:

> runif(1000) |> p.adjust(method = "BH") |> unique() |> length()
[1] 35

This is a natural consequence of how BH q-values are calculated, the complexity of which is often underappreciated, so I’ll point out the source from p.adjust() here:

lp <- length(p)
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o)
pmin(1, cummin(n/i * p[o]))[ro]

Specifically that rolling cumulative minimum can output the same value as it marches along the sorted p-value vector.

  1. A similar sort of situation as 2. You can see the source for BY adjustment by entering p.adjust (no parentheses) in the R console.