Hello,

I'm relatively new to DESeq2 and having problem with multiple factors. I'm analyzing 16S data from patients with 3 different treatments (FOS, PN, ABX). Some patients were treated with FOS, some with PN. I've combined these two treatments into one:

levels(dds$FOS_pn)
[1] "noFOSnoPN"   "noFOSyesPN"  "yesFOSnoPN"  "yesFOSyesPN"

For ABX treatment there are 4 different options:

levels(dds$ABX)
[1] "anaerobe" "gmneg"    "gmpos"    "no"

I tried the following design but I am getting an error:

> design(dds) <- ~ ABX + FOS_pn + ABX:FOS_pn
> ddsMF <- DESeq(dds)
Error in designAndArgChecker(object, betaPrior) : 
  full model matrix is less than full rank

My dds looks like (abxgmposneg = ABX):

> colData(dds)
DataFrame with 39 rows and 3 columns
          abxgmposneg study_day     FOS_pn
             <factor>  <factor>   <factor>
ch10163         gmneg         7 yesFOSnoPN
ch10164         gmneg         7 yesFOSnoPN
ch10165         gmneg         7 yesFOSnoPN
ch10166         gmneg         0 yesFOSnoPN
ch10167         no            2 noFOSnoPN 
...               ...       ...        ...
sodium306    anaerobe         7 noFOSyesPN
sodium307    no               1 yesFOSnoPN
sodium308    no               2 yesFOSnoPN
sodium309    no               3 yesFOSnoPN
sodium310    no               4 yesFOSnoPN

When I tried:

> design(dds) <- ~ abxgmposneg + FOS_pn

> resultsNames(ddsMF)
[1] "Intercept"                       "abxgmposneg_gmneg_vs_anaerobe"   "abxgmposneg_gmpos_vs_anaerobe"   "abxgmposneg_no_vs_anaerobe"     
[5] "FOS_pn_noFOSyesPN_vs_noFOSnoPN"  "FOS_pn_yesFOSnoPN_vs_noFOSnoPN"  "FOS_pn_yesFOSyesPN_vs_noFOSnoPN"
>

I am getting no errors, however I think I should do:

design(dds) <- ~ abxgmposneg + FOS_pn + abxgmposneg:FOS_pn

My goal is to calculate microbial differential abundance between FOS_pn groups with ABX treatment as a covariate (adjusted for ABX as a covariate). I have found many threads about this error and how to solve the problem, but I'm lost with several variables in each category (treatment).

I would appreciate your input.

I don't think we can diagnose the problem without seeing your full colData. But in general, what that means is you've got some abx categories which have exactly one Fos value, or vice versa.

If you want abx as a co-variate, the way you'd have batch as a co-variate, you probably want the simpler design. The interactions are more for "which genes change differently in response to treatment depending on the abx value?"

this may help

  dds$group <- factor(paste0(dds$abxgmposneg, dds$FOS_pn))
 design(dds) <- ~ group
 dds <- DESeq(dds)
resultsNames(dds)

see more here http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions