I have this simple question: how a gene modifies the transcriptional response of an organism to a treatment?

I collected RNAseq data for a wildtype (WT) and a mutant (Mut) under treatment (t) or no-treatment (u).

Factor Level 1 Level 2 Genotype WT Mut Condition t u

In the 4/4/2019 vignette Michael wrote:

“Initial note: Many users begin to add interaction terms to the design formula, when in fact a much simpler approach would give all the results tables that are desired. “

Then I was trying to follow the logic of this “much simpler approach”

The overall strategy is described like:

dds$group <- factor(paste0(dds$genotype, dds$condition))
design(dds) <- ~ group
dds <- DESeq(dds)
resultsNames(dds)
results(dds, contrast=c("group", "IB", "IA"))

in my case I expect that the grouped factors will be:

WT.t; WT.u; Mut.t and Mut.u

Now:

res1 <- results(dds, contrast=c("group", "WT.t", "WT.u"))

will provide the list of genes that respond to the treatment in the WT background. Coming back to my original question: “how my favorite gene “myFGen” modifies the transcriptional response of an organism to a treatment.” I will say that this is my “baseline”

res2 <-results(dds, contrast=c("group", "Mut.t", "Mut.u"))

will provide the list of genes that still respond to the treatment even the fact the organism does not have myFGen.

Comparing the number of significant differentially expressed genes between res1 with res2 would provide an idea of the impact of the mutation in the global response to the treatment. However, I do not see how to find those genes whose respond to treatment is most affected by the mutation in myFGen.

QUESTION 1: once I run all possible contrasts how I select the list of genes showing the strongest interaction? I mean those genes whose response to the treatment is more affected in mutant background.

res3 <-results(dds, contrast=c("group", "WT.u", "Mut.u"))

This list represents the genes that under normal condition are affected by the mutation in myFGen. I can look if any of the top-ranked genes from res3 are also top-ranked genes in res1:

top-ranked res 1 INTERSECT top-ranked res3

This intersected group are genes that are controlled by myFGen and also strong responder to the treatment, but this will not tell me how the response is modified. For that I need res2 and res1 again. I can plot the response of each gene in the intersection list in the WT background and the mutant background and compare them and visually judge the impact. This approach does not identify conditional responder to the mutation in myFGen, I mean genes that are only differentially expressed between genotypes under the treatment. For that I need to identify genes that are significant in res4 but not significant in res3.

res4 <-results(dds, contrast=c("group", "WT.t", "Mut.t"))

Finally, I can also find the intersection of significant genes between res4 and res3 to identify genes with differential expression between backgrounds, regardless the tested condition (t or u): sig- res 3 INTERSEC sig- res4. Then I need to plot all of them and visually identify those with altered response among genotypes.

QUESTION 2: coming back to the beginning of this post and the interaction section in the 2019 Deseq2 vignette => why this approach is simpler than running DEseq2 with and interaction term in the design formula?

FOR THE OTHER HAND

IF I run

ds <- DESeqDataSet(se, design= ~genotype + genotype:condition)

QUESTION3: can I still get contrasts like the above res1 and res2? I mean the sole respond of each background to the treatment.

If you are interested in the interaction, it’s not easy to do the combining factors approach. In the vignette I recommend this approach when users aren’t interested in differences but just in doing pairwise comparisons. You should use an interaction design is you want to test for the differences in the treatment effect.