[BioC] 3-way interaction with limma

[Gordon K Smyth]

Hi Elina,

The F-test you are using tests whether any of the four contrasts are 
non-zero.  Since each of the four contrasts compares a mated group back to 
corresponding non-mated group, it follows that the F-test is testing 
whether there is a mating response for any (one or more) female/male type.

The F-test does not specifically test whether the mating responses are 
different.  Indeed the F-test will be highly significant if all four 
contrasts are equally large, i.e., when there is a large but consistent 
mating response.

If you want to test whether the mating responses are different, you need 
make contrasts between the mating effects, as I suggested in my previous 
email.  There are three possible degrees of freedom for this.  You could 
for example define:

cont.matrix3 <- makeContrasts(
   A.DvsA.C = Mating.A.D - Mating.A.C,
   B.DvsB.C = Mating.B.D - Mating.B.C,
   BvsA = (Mating.A.C + Mating.A.C - Mating.B.C - Mating.B.D)/2,
levels=colnames(fit2))
fit3 <- contrasts.fit(fit2, cont.matrix3)
fit3 <- eBayes(fit3)

etc.

Best wishes
Gordon

On Fri, 8 Jul 2011, Elina Immonen wrote:

> Hi Gordon,
>
> Thanks very much for your helpful reply.
>
> Earlier I have gone through contrasting each of those mating types to 
> their corresponding virgins, and I have also looked at the differences 
> between A and B females across male types, as well as the effect of male 
> types across female types for mated females only. I have also done 
> female x male type interaction for mated females only, but here I don't 
> get any sig. differences with FDR 5% cutoff.
>
> What we are after in the end, however, is to find out for which genes 
> the mating responses (i.e. contrast to virgins) between the different 
> mating types are different, and not only what those responses are for 
> each of them.
>
> When I tried that attempt of doing the 3-way interaction (as in my 
> question), by creating contrast matrix of 24 contrast combinations by 
> rearranging the terms in different orders (terrible, I know!), I applied 
> F test on these and got out a group of genes significant with 5% FDR. I 
> plotted out the profiles of these genes for each mating type vs virgin, 
> and they indeed look very different.
>
> I just realised that perhaps I could get the same answer in the 
> following way (rather than trying to do the actual interaction):
>
> As you suggested, doing the individual contrasts separately:
>
>> cont.matrix <- makeContrasts(
>>  Mating.A.C = A.C.yes - A.none.no,
>>  Mating.A.D = A.D.yes - A.none.no,
>>  Mating.B.C = B.C.yes - B.none.no,
>>  Mating.B.D = B.D.yes - A.none.no,
>> levels=design)
>
> fit2 <-contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
> ttF <- topTableF(fit2, genelist=fit2$genes, adjust.method="BH", p.value=0.05, number=nrow(eset), sort.by="F")
>
>
> And then calling for topTableF on those. This gave me exactly the same 
> answer as that clumsy matrix of 24 contrasts. As I understand, F test 
> here compares if expression differences to virgins are different from 
> one another for any of the mating groups? GO and KEGG terms enriched for 
> this group of genes are meaningful.
>
> I'm pretty encouraged by these results, but I'd be grateful to hear your 
> opinion on whether what I see is correctly done and means what I think 
> it does.
>
> Thank you very much for your help!
>
> Best wishes,
> Elina
>
>
>
> On 6 Jul 2011, at 23:42, Gordon K Smyth wrote:
>
>> Hi Elina,
>>
>> If I understand your aims correctly, I don't think that your contrast is quite right.
>>
>> It's probably helpful to back to first principles here, rather than framing the question in terms of interactions.  My understanding is that you are interested in differential expression (DE) between mated and non-mated females, and in whether this DE depends on the type of male or female.  Basically, you are able to compute four possible mating effects:
>>
>>  Mating.A.C = A.C.yes - A.none.no
>>  Mating.A.D = A.D.yes - A.none.no
>>  Mating.B.C = B.C.yes - B.none.no
>>  Mating.B.D = B.D.yes - A.none.no
>>
>> I would suggest that you make each of these contrasts and examine then individually.  The limma code is almost exactly as above:
>>
>> cont.matrix <- makeContrasts(
>>  Mating.A.C = A.C.yes - A.none.no,
>>  Mating.A.D = A.D.yes - A.none.no,
>>  Mating.B.C = B.C.yes - B.none.no,
>>  Mating.B.D = B.D.yes - A.none.no,
>> levels=design)
>>
>> You can also make contrasts between any of the mating effects.  For example, does the type of the male make any difference for A females?
>>
>>  A.DvsC = Mating.A.D - Mating.A.C = A.D.yes - A.C.yes
>>
>> Note that this is not really an interaction.  It is just a comparison between two groups, because the non-mated baseline cancels out. Similarly, does the type of male make any difference for B females?
>>
>>  B.DvsC = B.D.yes - B.C.yes
>>
>> Approached this way, I think you'll find all the contrasts easy to interpret.
>>
>> Best wishes
>> Gordon
>>
>>> Date: Tue, 5 Jul 2011 11:12:23 +0100
>>> From: Elina Immonen <ei10 at st-andrews.ac.uk>
>>> To: bioconductor at r-project.org
>>> Subject: [BioC] 3-way interaction with limma
>>>
>>> Dear list,
>>>
>>> I am trying to pull out a contrast for 3-way interaction using limma, but as I'm pretty new to microarrays I'm puzzling whether I got it right or not. So far I haven't been able to find much info for doing an interaction model involving three factors, so any help is much appreciated.
>>>
>>> Here I'm interested in whether the female response to mating depends on female x male type interaction. That is, whether the difference btw. mated and virgin females depends on combination of females and males.
>>>
>>> I have 36 one-color Agilent custom arrays with the following nested design:
>>>
>>> Factor 1: mating status, levels= yes, no
>>> Factor 2: female type, levels= A, B
>>> Factor 3: male type, levels= C, D, none (dummy)
>>>
>>> Replication: 6 for each
>>>
>>>
>>> MatingStatus		FemaleType		MaleType
>>> yes				A				C
>>> yes				A				D
>>> no				A				none
>>> yes				B				C
>>> yes				B				D
>>> no				B				none
>>>
>>>
>>>
>>> Here is my script:
>>>
>>> MatingS <- factor(eset$MatingStatus)
>>> FT <- factor(eset$FemaleType)
>>> MT <- factor(eset$MaleType)
>>>
>>> Mating <- paste (FT, MT, MatingS, sep=".")
>>> Mating <- factor(Mating, levels=c("A.C.yes",  "A.D.yes", "B.C.yes", "B.D.yes", "A.none.no",  "B.none.no"))
>>>
>>> design <- model.matrix(~0+ Mating)
>>> colnames(design)<-levels(Mating)
>>>
>>> fit <- lmFit(eset, design)
>>>
>>> cont.matrix <- makeContrasts((A.C.yes-A.none.no)*0.25-(A.D.yes-A.none.no)*0.25-(B.C.yes-B.none.no)*0.25-(B.D.yes-B.none.no)*0.25,
>>> levels=design)
>>>
>>> fit2 <-contrasts.fit(fit, cont.matrix)
>>> fit2 <- eBayes(fit2)
>>>
>>> My questions are:
>>>
>>> 1. Is this appropriate way of doing the contrast?
>>> 2. Why do I get a different number of sig. genes if I swap the order of terms (i.e. bits under parentheses)? In my second approach I've done the contrast in a similar way but with the terms in different orders (all the possible combinations), and then used F test to get a list of genes significant in any of them.
>>> 3. How do I interpret fold changes from 3-way interaction?
>>>
>>> Many thanks for your help.
>>>
>>>
>>> Best,
>>> Elina
>>>
>>>
>>> ***********************
>>> Elina Immonen
>>> School of Biology
>>> Centre for Evolution, Genes and Genomics
>>> Dyers Brae House/Sir Harold Mitchell Building
>>> University of St Andrews
>>> KY16 9TH, St Andrews, Fife
>>> Scotland

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}