[BioC] limma fit interpretation doubt

[James W. MacDonald]

Hi David,

On 2/26/2014 11:04 AM, David Moriña Soler wrote:
> Hi Jim,
>
> thank you very much for your help! If I'm right, the solution you propose
>
> > contrast <- makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, 
> Treat2diff = Treat2.6h-Treat2.0h, levels = design)
>
> gives the difference in expression between times inside each treatment 
> group, while I'm more interested in the 'conventional' treatment main 
> effect,  as you suggest, a main effect for treatment, after 
> controlling for time. Is there a way to achieve that using limma?

Of course! You just use model.matrix():

 > treat <- factor(rep(paste0("Treat", 1:2), each = 2, times = 6))
 > treat
  [1] Treat1 Treat1 Treat2 Treat2 Treat1 Treat1 Treat2 Treat2 Treat1 Treat1
[11] Treat2 Treat2 Treat1 Treat1 Treat2 Treat2 Treat1 Treat1 Treat2 Treat2
[21] Treat1 Treat1 Treat2 Treat2
Levels: Treat1 Treat2
 > time <- factor(rep(paste0("Hr", c(0,6)), times = 12))
 > time
  [1] Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 
Hr0 Hr6 Hr0
[20] Hr6 Hr0 Hr6 Hr0 Hr6
Levels: Hr0 Hr6
 > design <- model.matrix(~treat + time)
 > colnames(design) <- gsub("treat|time", "", colnames(design))
## this step just necessary because R likes to append things to the 
column names...
 > head(design)
   (Intercept) Treat2 Hr6
1           1      0   0
2           1      0   1
3           1      1   0
4           1      1   1
5           1      0   0
6           1      0   1
 >


In this parameterization Treat2 is implicitly the difference between 
Treat2 and Treat1, and Hr6 is implicitly the difference between the 6 
hour and 0 hour times so you don't use a contrasts.fit() step:

fit <- lmFit(eset, design)
fit <- eBayes(fit)
topTable(fit, 2) ## treatment main effect
topTable(fit, 3) ## time main effect

But note that these main effects are nonsensical in the presence of a 
significant interaction term. So you could hypothetically do

 > design <- model.matrix(~treat*time)
 > colnames(design) <- gsub("treat|time", "", colnames(design))
 > head(design)
   (Intercept) Treat2 Hr6 Treat2:Hr6
1           1      0   0          0
2           1      0   1          0
3           1      1   0          0
4           1      1   1          1
5           1      0   0          0
6           1      0   1          0

And the Treat2:Hr6 coefficient is the interaction term. You would then 
want to first find the genes with significant interaction, and then 
ignore those genes when testing for significant main effects.

Best,

Jim




>
> Thanks!
>
> David
>
>
> On 26/02/14 16:55, James W. MacDonald wrote:
>> Hi David,
>>
>> On Wednesday, February 26, 2014 8:09:21 AM, David Moriña Soler wrote:
>>> Dear Bioconductor users,
>>>
>>> I just began to use limma package for the analysis of microarray data.
>>> I want to include in the analysis two factors (treatment and time) and
>>> the interaction between them. My design matrix looks like
>>> > head(design)
>>> Treat1.0h Treat1.6h Treat2.0h Treat2.6h
>>> 1           0           0          1          0
>>> 2           0           0          0          1
>>> 3           1           0          0          0
>>> 4           0           1          0          0
>>> 5           0           0          1          0
>>> 6           0           0          0          1
>>>
>>> Then, if I run
>>> > fit <-
>>> lmFit(data.norm,design,block=pData(data)$Subject,correlation=corfit$consensus), 
>>>
>>>
>>>
>>> I have this output for the coefficients:
>>> > head(fit$coefficients)
>>>                        Treat1.0h      Treat1.6h Treat2.0h Treat2.6h
>>> KCNE4              5.020670    4.981786   5.038805   4.924326
>>> IRG1               6.119265    6.015868   6.095171   6.027943
>>> SNAR-G2            8.242385    8.186429   8.144942   8.230391
>>> MBNL3             10.438644   10.417312  10.417042  10.358303
>>> HOXC4              8.985834    8.854698   8.969801   8.939682
>>> ENST00000319813    3.913602    4.102653   4.000681   3.960431
>>>
>>> How can I obtain for example the effect of the treatment, using
>>> contrasts and eBayes?
>>
>> It depends on what you mean by 'the effect of the treatment'. Are you 
>> looking for a main effect for treatment, after controlling for time 
>> (e.g., a conventional main effect)? Or are you looking for the 
>> difference in treatment at each time separately?
>>
>> The easiest way to create a contrasts matrix is to use 
>> makeContrasts(), which is as simple as deciding what comparisons you 
>> want to make. So for instance, let's say you want to know the 
>> difference in expression between the treatments at each time:
>>
>> contrast <- makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, 
>> Treat2diff = Treat2.6h-Treat2.0h, levels = design)
>>
>> If you want the interaction too, it's easy to do that as well
>>
>> contrast <-  makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, 
>> Treat2diff = Treat2.6h-Treat2.0h, Interaction = 
>> Treat1.6h-Treat1.0h-Treat2.6h+Treat2.0h, levels = design)
>>
>> Best,
>>
>> Jim
>>
>>
>>>
>>> Thank you very much,
>>>
>>> David
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives:
>>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>> -- 
>> James W. MacDonald, M.S.
>> Biostatistician
>> University of Washington
>> Environmental and Occupational Health Sciences
>> 4225 Roosevelt Way NE, # 100
>> Seattle WA 98105-6099
>
>

-- 
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099