-
Notifications
You must be signed in to change notification settings - Fork 6
Description
Hello!
I'm trying to use aldex.glm() to explore differential abundance for a complex design involving both Genotype (B6 vs. TCR) and stress condition (Stress vs. NoStress). My code is as follows (dat_pr is my filtered phyloseq object):
#create model matrix
meta = data.frame(dat_pr_OBIT_f@sam_data)
mm = with(meta, model.matrix(~GENOTYPE + Stress + GENOTYPE:Stress))
#extract OTU table
otus = data.frame(dat_pr_OBIT_f@otu_table)
# check alignment
all(meta$SAMPLE_ID == colnames(otus))
#CLR transform values with model matrix
glm <- aldex.clr(otus, mm, mc.samples=128, denom="all")
glm.test <- aldex.glm(glm)
glm.test = rownames_to_column(glm.test, "ASV")
As mentioned in a previous issue, originally it seems that the glm function had a "drop1" argument, which based on my knowledge means that differentially abundant taxa for the first factor would only be tested within the baseline group of the second factor. For example, differentially abundant taxa B6 vs. TCR (model.GENOTYPETCR) would only be tested in the "NoStress" condition and vice versa. Is this still the case? If so, how can I construct a model that tests for the global effects of genotype on taxon abundance while still accounting/controlling for stress, without removing sample groups, similar to something like DESeq2 or corncob?
Thanks so much in advance!
Miranda