Search
Question: edgeR: setting up proper contrasts for interactions
0
3 months ago by
fcamus0
fcamus0 wrote:

Hi all,

I am curious on how to go about setting up proper contrasts to test for interactions. The experimental design I have is

• 2 diets (B & E)
• 2 sexes (Male vs Fem)
• 2 reproductive status (M vs V)

There are 3 biological replicates, making a total of 24 samples (libraries).

I have setup my design matrix:

des.full2 <- model.matrix(~0+full.mat2)

M.Fem.B M.Fem.E M.Male.B M.Male.E V.Fem.B V.Fem.E V.Male.B V.Male.E
1        0       0        0        0       1       0        0        0
2        0       0        0        0       1       0        0        0
3        0       0        0        0       1       0        0        0
4        1       0        0        0       0       0        0        0
5        1       0        0        0       0       0        0        0
6        1       0        0        0       0       0        0        0
...

So far what I've been doing is individual contrasts to find DE genes between specific treatments, which is working out quite well.
What I'd like to do is get a list of genes that interact significantly between the treatments.

diet
sex
status
diet:sex
[...]
diet:sex:status

For the single factors I've gone for the average across all treatments:

overall.diet <- makeContrasts((M.Fem.B+V.Fem.B+M.Male.B+V.Male.B)/4 -
(M.Fem.E+V.Fem.E+M.Male.E+V.Male.E)/4, levels=des.full2)

How would one set up the contrasts for 2-way and 3-way interactions?

Thanks very much for all your help!

modified 3 months ago by Aaron Lun21k • written 3 months ago by fcamus0
3
3 months ago by
Aaron Lun21k
Cambridge, United Kingdom
Aaron Lun21k wrote:

Two-way interactions are pretty easy to set up and test with makeContrasts. Let's say you want to know the interaction between diet and reproductive status for male mice. The contrast would simply involve testing whether the diet effect in status M is the same as the diet effect in status V:

con <- makeContrasts(
(M.Male.B - M.Male.E) # effect of diet in status "M"
- (V.Male.B - V.Male.E), # effect of diet in status "V"
levels=design)


Obviously, this is identical to:

con <- makeContrasts(
(M.Male.B - V.Male.B) # effect of status in diet "B"
- (M.Male.E - V.Male.E), # effect of status in diet "E"
levels=design)


... and it is fairly easy to extend this to all your interactions of interest.

Three-way interactions are much more difficult to interpret. It's much like differences of two-way interactions (in the same way that two-way interactions are differences of log-fold changes). Continuing from the example above, we would get the interaction for female mice:

# omitting makeContrasts() for simplicity
(M.Fem.B - V.Fem.B) # effect of status in diet "B"
- (M.Fem.E - V.Fem.E), # effect of status in diet "E"


The difference of the male- and female-specific interactions would then be:

(M.Fem.B - V.Fem.B) - (M.Fem.E - V.Fem.E)
- ((M.Male.B - V.Male.B) - (M.Male.E - V.Male.E))


... which is your three-way interaction. Or perhaps, in a more mathematically elegant form:

V.Fem.E # all three factors
- (M.Fem.E + V.Male.E + V.Fem.B) # any two factors
+ (M.Male.E + V.Male.B + M.Fem.B) # any one factor
- M.Male.B # "baseline"


... where male, M and B are treated as the baseline for the respective factors.

Testing the three-way interaction is equivalent to asking whether the diet-status interaction in males is the same as the diet-status interaction in females. (Which is equivalent to asking whether the sex-diet interaction is the same across status, or if the sex-status interaction is the same across diets, etc.) Whether or not this is a sensible question is up to you.

Many thanks Aaron for your input,

I guess what we are thinking about is clustering the genes we have based on their significance, and not really dwelling too much into fold change (for this analysis). For instance we want to say "X number of genes respond to diet, sex, status, sex:status ... diet:sex:status". Once we have a list of gens that respond to each treatment/interaction we could use something like MBcluster.Seq to look at their expression patterns across all treatments. Does this sound legit?