Hello,

I need some help creating a design for my dataset in DESeq2.  I have RNA-Seq data from samples that have been treated with a Control, with LPS, with LPS followed by molecule A, or with LPS followed by molecule B. There are 5 subjects that have been treated with each condition.  So my coldata is as follows:

   Condition Subject
1    Control       1
2        LPS       1
3      LPS_A       1
4      LPS_B       1
5    Control       2
6        LPS       2
7      LPS_A       2
8      LPS_B       2
9    Control       3
10       LPS       3
11     LPS_A       3
12     LPS_B       3
13   Control       4
14       LPS       4
15     LPS_A       4
16     LPS_B       4
17   Control       5
18       LPS       5
19     LPS_A       5
20     LPS_B       5

I would like to identify genes in which the response to LPS is altered by Molecule A, but not by Molecule B.  Therefore, I need to identify genes that are DE in the LPS group compared to Controls, and then identify which of those genes are further DE by A, but not by B. Since the same 5 individuals have been treated with each of the conditions, I would like to test within subjects.  

Is it possible to test for this using a multifactorial model in DESeq2?  If so, could someone please help me with the design term, and how to specify the contrasts?

Hi asafari,

the design you will need is:

design ~ subject + condition

Deseq2 will then fit a coefficient (think of this as a mean per group)
for every level of the two factors AND an intercept. Inspect the fitted
coefficients via:

resultNames(dds)

This should give something like:

 "Intercept"  "Subject1"    ...  "Subject5" "conditionControl", "ConditionLPS", "ConditionLPS_A", "ConditionLPS_B"

Now want to you are interested in are the DE genes for LPS vs Controls, this is
the comparison "ConditionLPS"  vs "conditionControl" so you can specify

results(dds, contrast=c("condition","LPS","Control"))

Now it becomes more complicated, I assume here you also want to know which genes
have a stronger effect for A than for B, this you can get via the difference of
differences:

(LPS_A - Control) - (LPS_B - Control)  = LPS_A - LPS_B

so

 results(dds, contrast=c("condition","LPS_A","LPS_B"))

however, this does not exclude genes which are DE for

LPS_B - Control, so you might filter for them using another test of

 LPS_B - Control

before as in:

results(dds, contrast=c("condition","LPS_B","Control"))

So in summary we would have three steps

1.  compare LPS_A - LPS_B

 results(dds, contrast=c("condition","LPS_A","LPS_B"))

2. filter for LPS vs Control DE genes:
 results(dds, contrast=c("condition","LPS","Control"))

3. filter for genes  that are non-DE in  LPS_B vs Control
 results(dds, contrast=c("condition","LPS_B","Control"))

Maybe there is a more elegant way to analyze this experiment but it
does not seem to be factorial (so that molecules and LPS alway occur together) so I don't know whether you can fit
main effects for LPS and the molecules as well as interaction terms.

Best wishes,

Bernd