FrederickHuangLin / ANCOMBC

Differential abundance (DA) and correlation analyses for microbial absolute abundance data
https://www.nature.com/articles/s41467-020-17041-7
104 stars 28 forks source link

queries regarding pairwise test and outputs? #117

Closed marwa38 closed 1 year ago

marwa38 commented 1 year ago

Hi team

Thank you so much for ANCOMBC2 and the updates regarding pairwise testing.

My aim is to compare 6 groups together (M,V,MM,VM,MMV,VMV) in a pairwise form and see the differentially abundant taxa in each comparison along with the p and q values. According to my understanding, this means that I won't use a reference/intercept in my study. My main question is whether I need to only use the argument pairwise = TRUE and mdfdr_control = list(fwer_ctrl_method = "holm", B = 1)) for multiple comparisons (without the need to add dunnet = TRUE)? The expected results comparing every two groups should be in out.intes.phylum$res_pair? as I couldn't find my M group in out.intes.phylum$res_pair is it advisable to use the same test p_adj_method = BH as in mdfdr_control = list(fwer_ctrl_method = "BH" ...?

Kindly find attached my codes

## intestine
set.seed(123)
# phylum 
out.intes.phylum = ancombc2(
  data = ps.prev.intes, 
  tax_level = "Phylum",
  fix_formula = "Regime", 
  p_adj_method = "BH",
  prv_cut = 0,
  group = "Regime", 
  struc_zero = TRUE, 
  alpha = 0.05, 
  pairwise = TRUE,
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 1))
out.intes.phylum # attached below as zipped file of .rds (use readRDS to import in case needed)

[out.intes.phylum.zip](https://github.com/FrederickHuangLin/ANCOMBC/files/10047309/out.intes.phylum.zip)

One last thing if you don't mind, could you please let me know if this figure was created by ancombc2()? image

Thank you M

FrederickHuangLin commented 1 year ago

Hi @marwa38,

As for your questions:

My main question is whether I need to only use the argument pairwise = TRUE and mdfdr_control = list(fwer_ctrl_method = "holm", B = 1)) for multiple comparisons (without the need to add dunnet = TRUE )?

Yes, there is no need to specify dunnet = TRUE. However, there is one problem with your code, you shouldn't set B = 1 since it means you only specify one bootstrap sample. Keep it as the default value B = 100 or larger is highly recommended.

I couldn't find my M group in out.intes.phylum$res_pair

I am not sure if I get what you mean. Given there are 6 groups in your data, did you see 15 LFCs (SEs, Ws, p-values, etc.) in the output (out.intes.phylum$res_pair)?

is it advisable to use the same test p_adj_method = BH as in mdfdr_control = list(fwer_ctrl_method = "BH" ...?

Nope. fwer_ctrl_method is recommended to set to holm as it is specifically used for correcting family-wise error (FWER).

One last thing if you don't mind, could you please let me know if this figure was created by ancombc2()?

Nope, this figure was generated using the original (pre-release) ANCOM-BC codes.

Hope these help!

Best, Huang

marwa38 commented 1 year ago

Thanks very much for your detailed answer @FrederickHuangLin

Keep it as the default value B = 100 or larger is highly recommended.

I reran the analysis and left as it is. Based on what I need to choose a number above 100?

Yes, there is no need to specify dunnet = TRUE

and no need also for the trend = TRUE regarding getting pairwise without having the intercept (M group)?

this figure was generated using the original (pre-release) ANCOM-BC codes.

any tutorial to create these kind of figures?

I am not sure if I get what you mean. Given there are 6 groups in your data, did you see 15 LFCs (SEs, Ws, p-values, etc.) in the output (out.intes.phylum$res_pair)?

yeah, but I am not getting my M group in the comparisons, it is the first group so it is used as intercept.

set.seed(123)
# phylum 
out.intes.phylum = ancombc2(
  data = ps.prev.intes, 
  tax_level = "Phylum",
  fix_formula = "Regime", 
  p_adj_method = "BH",
  prv_cut = 0,
  group = "Regime", 
  struc_zero = TRUE, 
  alpha = 0.05, 
  pairwise = TRUE,
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100))

can't I control the comparisons to be between certain groups of the 6 groups rather than all groups?

  1. M vs V
  2. MM vs VM
  3. MMV vs VMV
  4. M vs MM
  5. V vs VM
  6. MM vs MMV
  7. VM vs VMV
  8. M vs MM
  9. V vs VM
  10. MMV vs M
  11. VMV vs V

Thank you so much Marwa

FrederickHuangLin commented 1 year ago

Hi @marwa38,

I reran the analysis and left as it is. Based on what I need to choose a number above 100?

Increasing B will give you more precise p-values at the price of longer computing time. Usually, the default value is sufficient based on our simulation benchmark.

and no need also for the trend = TRUE regarding getting pairwise without having the intercept (M group)?

You are right. It is not necessary if the trend test is not of your interest.

any tutorial to create these kind of figures?

You can see here for a reference. I will add an example for creating such plots in a vignette in the next iteration. Thanks for your suggestion!

yeah, but I am not getting my M group in the comparisons, it is the first group so it is used as intercept.

There is no intercept in the pairwise comparisons. It is still unclear to me what you meant by "not getting M group." Suppose you have three groups: A, B, and C. In the pairwise outputs of ancombc2, you will get three columns: "B", "C", and "B - C", which means "B - (minus) A", "C - A", and "B - C". Note that "A" is automatically selected as the reference group since it is the lowest level in alphabetical order.

can't I control the comparisons to be between certain groups of the 6 groups rather than all groups?

Unfortunately, the current iteration of ANCOMBC does not support specifying certain contrasts.

Best, Huang