Hello everyone,

The good point of self-isolation is that I can review my code used for an important and expensive experiment I performed in my lab. First of all, a bit of introduction here:

• My project is about determining all the genes changed following incubation of mast cells with Listeria.
• We incubated cells (incubated or not with Listeria) at 12 and 24 hours. We then extracted RNA and kept homogenate in the -80 freezer. Indeed, given the fact mast cells requires a lot of time to be differentiated in vitro (several weeks) we had to carry out the experiments in different batches. That was an unavoidable step.
• Finally, per group, we have 4 different combinations (untreat_12hrs, treat_12_hrs, entreat_24hrs, treat_24hrs)

My goal is to compare gene expression of samples incubated at the same time vs a control (i.e. untreated) condition.

to sum up please find enclosed the table matrix describing the whole experiment:

colData = DataFrame(icb=factor(rep(c("untreat_12hrs","treat_12hrs","untreat_24hrs","treat_24hrs"),times=5)),
                    btc=factor(rep(c("1","2","4","5","6"),each=4)),
                    tme=factor(rep(c("twelve","twentyfour"), each=2, times=5)))

So far, I have taken into account the batch effect only but I don't think this is the only variable I have to control, right? How would you proceed and how would you create the model matrix? I checked the DESeq2 vignette and more specifically the "Group-specific condition effects, individuals nested within groups" section.

colData$ind.n = factor(rep(c('1','2'), each=2, times=5))

and then

model.matrix(~ btc + btc:ind.n + btc:icb, colData)

I think I am doing something wrong. particularly, not entirely sure if the relationship created are in line with the above mentioned rationale.