Hi all,

I am new to R/bioconductor and RNAseq but I'm an experienced Bioinformatician and have performed due diligence in reading the docs and other questions and can't find what I'm after.

I am trying to answer the question:

What is differentially regulated between cell lines infected with WT and mutant strains, normalised to non-infected cells.

I have 12 samples:

Sample condition

S1 control
S2 WT
S3 MUTANT
S4 control
S5 WT
S6 MUTANT
S7 control
S8 WT
S9 MUTANT
S10 control
S11 WT
S12 MUTANT

So I want to do DE between WT and MUTANT and remove genes that are differentially expressed between control and WT and genes that are differentially expressed between control and MUTANT. Is there a way to specify this with the deseq2 'design' function? Or do I have to stick with my current approach of doing WT/MUTANT, control/WT and control/MUTANT separately and then using a Python script to find genes that are DE between WT/MUTANT only?

Thanks,

Liam

First of all, your current approach is wrong (even though it is commonly done your way).

This is because DESeq2 (and, similarly, most other methods) performs hypothesis tests: it gives you a list of all genes for which there is sufficient evidence to reject the null hypothesis of them not being differentially expressed. So, if a gene is on the list, it is most likely affected by the treatment, but if it is not on the list, this does not at all mean that it is not affected by the treatment. ("Absence of evidence is not evidence of absence of an effect.")

You want genes for which you have evidence that (a) they differ between mutant infection and wt infection and (b) that they do not differ (or do not differ much) between control and wt infection. So, for (b) you have to reverse the "burden of proof". In DESeq2, this is done with the "altHypothesis" option to the "results" function. 

You could run results two times, first to obtain a list of genes that differ between mut and wt, then a run to find genes for which the there is evidence that the effect of wt infection compared to control is less than some reasonable threshold ("lfcThreshold" argument of "results"). (Note that this is not the same as filtering the results list by LFC; see our paper on DESeq2.)

Now, to find the intersection between the two lists, you can, of course, use Python. However, it may be worth your time to learn enough of R to do such trivial manipulations without having to resort to another language. (Hint: Use "intersect" on the gene IDs.)

The other question is: is this really what you want? Let's say a gene gets somewhat upregulated by infection by wt strains but strongly downregulated by infection with mutant strains, this might be a very interesting finding -- but your approach would miss it.

My recommendation would be to plot the (shrunken) log fold changes for the effect of infection (wt versus control) and for the effect of the mutation (mut versus wt) in a scatter plot and decide from this plot how to proceed.

Edit: You wrote originally that you are looking for genes that differ between mut and wt but neither between wt and control nor between mut and control. Note that this is a logical impossibility! (If wt=ctrl and mut=ctrl, then mut=wt.) So, if you found anything that way, then this were mainly the cases where you have absence of evidence but no evidence of absence of an effect of infection.