Hi everyone,

I am finalizing the design of an RNASeq experiment and am having doubts about how I should plan to approach differential expression analysis. My design includes three factors of interest: Surgery group (2 levels-sham and surgery), time point (2 levels), and sex (2 levels). I have some experience with DESeq2 but am considering switching to edgeR or limma. The biological questions I am hoping to address are:

1. What genes are differentially expressed between surgery and sham groups?
   1. Of the genes affected by surgery, which are differentially expressed between males and females?
   2. Of the genes affected by surgery, which are differentially expressed between the two time points?
2. What genes are differentially expressed between males and females overall?

I’m looking for help with the following questions:

1. Is DESeq2, edgeR, or limma more straightforward for this type of analysis?
2. How should I design my differential expression model with these questions in mind? Generally, I am interested in which of the following two options is most correct:
   1. Make a model including interaction terms for SurgeryGroup:TimePoint and SurgeryGroup:Sex (such as Example 2 of the DESeq2 help for the Results function, but with an additional interaction term).
   2. Combine all of the condition combinations into a single vector and individually do all of the relevant contrasts, as described here: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions. If I do this, is it legitimate to first do a contrast to look at main effect of sex, then to look at which genes respond differently to surgery in male vs. female mice? My confusion may be stemming from the examples given in the DESeq2 documentation.
   3. Has anyone come across good examples of a 3-factor design with interactions?
3. In discussing my design with others, the concern was raised that I have 8 groups total and that by making comparisons between most of them with interactions, I will be reducing my statistical power to detect differentially expressed genes. To what extent is this a legitimate concern?   

I realize these are basic questions, so thank you for bearing with me. I appreciate any advice and suggestions you can provide!

Thanks!

Kendra

I'm the senior author of the limma and edgeR packages, so naturally I like them myself. But none of the three packages would have any difficulty with your experiment. They would each be straightforward in their own ways.

Three factor designs follow the same principles as for two factor experiments. You can essentially follow two factor examples but just extrapolate them to three factors.

You can fit a 3-way interaction model or else put all the condition combinations into one factor. These are statistically equivalent approaches and therefore both equally correct. However I have long argued (since before edgeR or DESeq existed) that the one factor approach is vastly easier, more interprettable and better suited to genomic applications. See section 9.5 of the limma User's guide: http://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

The factorial model is already leading you into an incorrect way of thinking. You cannot test for main effects and look at interactions later. You have to start the other way around. Looking for the "overall" or "average" treatment of only meaningful once you have established that the effect is consistent over both sexes or over both times. Otherwise you have to be specific about which time or sex the treatment is applied to. Averaging diverse effects to get a "main effect" is not a meaningful thing to do.

I would instead pick a particular time point of interest, then look for genes that are differentially expressed between treatments for either males or females at that time point. Then test for which genes have different treatment effects at that time point between male and females. This is easy to do in limma or edgeR and presumably in DESeq2 as well.

You really have no choice but to deal with 8 treatment combinations. That's what you have.