Hi,

I've looked at a number of posts on multiple treatment comparisons, but couldn't find an example of this exact situation. I'm looking at RNA expression in two cell types WT and KO. For each cell type I have 5 different treatment conditions: Mock, A, B, C, D.

Cell    Treatment
WT    mock
WT    mock
WT    mock
WT    A
WT    A
WT    A
WT    B
WT    B
WT    B
WT    C
WT    C
WT    C
WT    D
WT    D
WT    D
KO    mock
KO    mock
KO    mock
KO    A
KO    A
KO    A
KO    B
KO    B
KO    B
KO    C
KO    C
KO    C
KO    D
KO    D
KO    D​

First I'd like to ask is there a difference in expression between any of the treatment conditions (mock, A, B, C, D) within a cell type. Looking at some other forum posts, it looks like people recommend just grouping them and then using contrasts to compare everything to mock. Is this the best method?

dds$group <- factor(paste0(dds$Cell, dds$Treatment))
design(dds) <- ~ group
results(dds, contrast=c('group', 'AWT', 'MOCKWT')
results(dds, contrast=c('group', 'BWT', 'MOCKWT')
results(dds, contrast=c('group', 'CWT', 'MOCKWT')
results(dds, contrast=c('group', 'DWT', 'MOCKWT')

However, next I would like to ask whether there is a difference in treatment vs mock comparing the WT and KO cells. Essentially, is there a difference in the genes that change between A vs mock in WT and A vs mock in KO, B vs mock in WT and B vs mock in KO, etc. I think I would need to rewrite the design model with an interaction term, but I'm not quite sure what it would look like. Any suggestions would be greatly appreciated.

Thanks

Yes, grouping them is the best method. You could set up a factorial model with interaction terms, but that's unnecessarily complicated and you'll end up with the same results anyway.

For the contrasts you want to do, let's lay down some groundwork first. Assuming you're using limma, let's do:

design <- model.matrix(~group)
colnames(design) <- levels(group)

You can then test for a differential A vs. mock effect between WT and KO by running:

con <- makeContrasts((WTA - WTmock) - (KOA-KOmock), levels=design)

... and feeding the resulting contrast into contrasts.fit. As the code suggests, you're literally comparing the log-fold change in WT A vs. mock against the corresponding log-fold change in KO. Rejection of the null corresponds to significant differences in treatment effects between WT and KO. This is equivalent to an interaction term but is a lot easier to interpret.

Obviously, this strategy can be extended to test for WT/KO-specific DE between any pairs of treatments, e.g., B vs. mock, A vs. B.

I'm assuming you're using the latest version of DESeq2, v1.10. If not, see here how to upgrade Bioconductor to the latest version:

http://bioconductor.org/install

I'm also assuming mock is the reference level of treatment (see the vignette on how to set this).

There are a couple ways to do this for DESeq2, if you want to compare treatment within cell type and also treatment differences across cell type. Maybe the easiest is a design of ~ Cell + Cell:Treatment.

Then the Treatment effects (A vs mock within a cell type) are:

results(dds, name="CellWT.TreatmentA")

And to compare A vs mock across cell type:

results(dds, contrast=list("CellKO.TreatmentA","CellWT.TreatmentA"))

Sorry to hijack this thread, but I have a very similar experimental design (2 genotypes and 4 treatments) and use a ~ Genotype + Genotype:Treatment 

Using the naming of the origin post, I would like to identify differentially expressed genes that satisfy the following criteria:

"same behaviour in genotypes WT and KO for treatments (mock and A and B) but different behaviour in KO compared to WT in treatment C"

I do not see exactly how I could get to these.