Hej,

I am analysing data were bacterial communities were exposed to a toxin. I have paired data (before / after) for each community. However, we applied the toxin in a concentration gradient (with no replicates per level but for the control). Now I am interested in finding those OTU's where the difference (before/after exposure within each community) correlate with exposure concentration.

I tried ~BOT + Sampling + Tox but I get the error message the model matrix is not full rank, so the model cannot be fit as specified. One or more variables or interaction terms in the design formula are linear combinations of the others and must be removed.

I probably can make it work by dividing the Tox column into two or three groups (e.g. control, low, high) but this grouping would be somewhat arbitrary. Any pointers would be highly appreciated!

(ps.: I know that the data structure is not ideal but the design was not chosen primarily for DE)

The data look like:

This is covered in the vignette here.

This design will help you to find increases or decreases in the LFCs related to toxicity > 0 compared to toxicity = 0. Because you don't have replicates at the time points, you are a bit limited in what you can find. You have to specify a functional form of what you are looking for. If you want to do something more complicated than finding increases or decreases, I'd suggest working with a statistician.

Make columns of colData that look like this:

sample <- factor(rep(c(1,2,3,1),c(2,2,2,12)))
time <- factor(rep(0:1,9))
tox <- c(rep(0,6),rep(1:6,each=2)) # note: a numeric, not factor

Then you can use a design of ~sample + tox + time + tox:time, and simply run DESeq() and then you can test the interaction term using results(dds). The interaction term will be non-zero if there is an increase of decrease in the LFC due to end/mid as you have above, going along with toxicity and relative to toxicity=0. The null hypothesis for this results table includes those genes where LFCs for toxicity>0 are the same as toxicity=0.