Abdullahi Muhammad
posted 8 days ago

How to select differential abundant ASVs/ Keggs from ancombc2 for enrichment analysis

Hello all, I have been working on my 16S amplcon data for a while now and I have gotten to the last of the downstream analysis where I am stuck and I dont know hwo to move forward. I have data set that I woud say loks like a full factorial; Genotype (4 levels; G1, G2, G3, & G4), Day (3 levels; D1, D2 & D3), Treatment (6 levels; Control, Atrazin, PFOS, Diclo, Arsenic, wastewater) and Replicates (3 biolgical replicates of the genotypes across the time points and treatment). I have run a differential abundance analysis using the function "ancombc2" that uses the lmerTest in its model. This i think suites my kind of data because it will allow me look for interaction among the variabels and I will also have a nested model with replicates as random effect. Please see below my code:

output2 = ancombc2(data = ps, assay_name = "counts", tax_level = "Genus",
                  fix_formula = "Treatment * Genotype * Day", rand_formula = "(1|Replicates)",p_adj_method = "holm", pseudo = 0, pseudo_sens = FALSE,prv_cut = 0.10, lib_cut = 0, s0_perc = 0.05,group = "Treatment", struc_zero = FALSE, neg_lb = FALSE,alpha = 0.05, n_cl = 2, verbose = FALSE,global = TRUE, pairwise = TRUE, dunnet = TRUE, trend = FALSE,iter_control = list(tol = 1e-2, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
                  trend_control = list(contrast = NULL, node = NULL, solver = "ECOS", B = 100))
# ps = phyloseq object

I assume that the pairwise comparison will be agaisnt the base "Treatment", am not too famiiar with the meaning of the ancombc output. The "output" has several files: global, prim pairs, and Dunn test. I can see in the 'prim' output interactions but most are false in terms of p-val but the 'global' has a different table structure with diff_abun column, W, adj_pval and the taxon. I other to move forward with this analysis, my aim is to identify ASVs,/ kegg genes that are enriched and then visualise this. but at this point I dont know how to selct the diff_adun ASVs to create a list that will be use for enrichement analysis. To clarify, I am using the amcombe package to run differential abundance analysis on both picurst2 kegg output and phyloseq object for ASVs I would be grateful if anyone could share their thoughts on this. Thank you I

Login or Signup to leave a comment