I have conducted an analysis differential expression analysis using edgeR, but am receiving some criticism of the design of the statistical tests.  I am hoping someone can either verify that the approach I have taken is valid, or that I should be doing the type of analysis that has been suggested.  Some background on the experiment:

I am interested in gene expression changes when individuals from two different populations are fed on different food sources.  Specifically, each population typically specializes on a particular type of food that differs between the populations.  The experimental design involved feeding individuals from each population on their normal food source as well as the food source typically used by the other population.  In addition, each population was also fed an artificial food source.

Population: two levels

Food: three levels

Each combination had three biological replicates.

There are specific pre-planned comparisons I am interested in:

1. Compare gene expression between populations on each of the two typical food sources (I am not interested in comparing the populations on the artificial food source).
2. Compare gene expression on the food sources within each population (I am only interested in one comparison involving the artificial food source though)
3. Compare gene expression between the populations across all food sources 

Following the recommendations in the edgeR documentation (as I understood them in section 3.3.1), I set this up so that each treatment combination was a group and then made contrasts corresponding to the comparisons listed above.  One of the follow-up goals was to look at reaction norms for DE genes from comparison one, looking for genes following a few specific patterns (basically constructing a generic interaction plot, though some of the patterns I’m looking for would actually involve no interaction).  I reconstructed these by examining overlapping DE gene lists produced by comparison one and two.  So, for example, I used comparison one to plot the relative location of points for each food in each population (e.g. spaced apart if DE or on top of each other if not DE), and comparison two to connect the lines (e.g. upward sloping, downward sloping, or flat if the gene was not DE in comparison two).  Note: specific values are not of interest, but rather just the general pattern.

It has been suggested that there is a problem because I am not examining interactions formally, and using this as a starting point to then look at the reaction norms.  I interpret this to mean that I should analyze the data in a traditional full factorial format where a significant interaction term would give license to further dissect the interaction.  The first issue I have with this is that the artificial food source was only included for a few specific comparisons, and I wouldn’t be interested in some of the comparisons using this group that would naturally come out of a factorial analysis.  I suppose that I could look at those experiments separately and then analyze the data with the other two foods in a factorial design.  However, ultimately I think I would end up needing to construct the contrasts (one and two) that I set-up directly in the first place in order to investigate the reaction norms. Moreover, some of the comparisons would not come out as directly if I used the factorial method.  For example, to get all the genes that differ between populations on one food source I would have to pull some of the genes from the interaction list and also the population main effect list. My understanding from the manual (and more so from the limma manual) is that these two approaches would basically be equivalent. But, I want to make sure I am not missing something by not doing the full factorial analysis as obviously I would like to do the most appropriate analysis.  Is the approach I have taken valid?  Am I really “ignoring” the interaction by the approach I have taken as has been suggested?  Is there anything to be gained in this case by doing the full factorial and looking at interactions specifically?  If so, I may need to ask a follow-up regarding specifically how to set this up.  I apologize for the long post, and greatly appreciate any help anyone can give.    

Thanks,

Jeremy

Hi Jeremy,

I agree that you should statistically test the interaction between food type and population. If you set up each treatment combination as a coefficient in the model and used these to make your pairwise contrasts, then you can also use them to make specific interactions. For example, if you have the coefficients of:

P1.F1, P1.F2, P1.F3, P2.F1, P2.F2, P2.F3

I assume you made two contrasts like this to answer your first question of gene expression between populations on each of the two typical food sources:

P1.F1 - P2.F1  and P1.F2 - P2.F2

To make the formal interaction between just food types 1 & 2 and population, the contrast would be:

(P1.F1 - P2.F1) - (P1.F2 - P2.F2) 

HTH,

Jenny

Your question has both statistical and scientific components.

I will address the statistical issues first. No, you are not ignoring the interaction. The group-meany layout design matrix is exactly equivalent to fitting the full factorial model. It is just parametrized differently so that comparisons of interest are easier to extract.

There is nothing to be gained by using a factorial design matrix, because you can easily test any comparison you wish (including interactions) using the design matrix you have.

Now for scientific issues. It is scientific questions that drive what tests you should be doing, and we can't generally advise you on the specific scientific issues relating to your experiment. However I would guess from your post that it might be of interest to test for genes that have different reaction norms between the populations, and so far you are not doing that. If you are truly not interested in the artificial food source, then the test for differences in reaction norms is simply the test for the interaction contrast that Jenny gave you.