[BioC] What's the difference of the two approach for two groups in limma?

Wed Apr 6 22:57:52 CEST 2011

Hi Jianhong Ou,

On 4/6/2011 3:07 PM, Ou, Jianhong wrote:
> Hi all,
>
> In the usersguide.pdf of limma, there are two approach for two groups (8.4, Two Groups: Common Reference)
>
> At first, I thought they are same, but...
>
> [code]
>> sd<- 0.3*sqrt(4/rchisq(100,df=4))
>> y<- matrix(rnorm(100*6,sd=sd),100,6)
>> rownames(y)<- paste("Gene",1:100)
>> y[1:2,4:6]<- y[1:2,4:6] + 2
>> design<- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
>> fit<- lmFit(y,design)
>> fit<- eBayes(fit)
>> topTable(fit,coef=1,n=20,sort.by="logFC")

You made a mistake here. You want coef = 2, which is the contrast 
between groups 1 and 2. What your code is doing is testing to see if the 
average of the first group is equal to zero or not.

Best,

Jim

>          ID      logFC          t     P.Value adj.P.Val          B
> 33 Gene 33  1.2800739  4.9383702 0.001156486 0.1156486 -0.4485569
> 23 Gene 23 -0.9290832 -1.6732940 0.133021448 0.9913757 -4.9814041
> 27 Gene 27  0.8616516  3.3965657 0.009487765 0.4743882 -2.4838934
> 76 Gene 76  0.7509497  1.7891908 0.111590341 0.9913757 -4.8250883
> 54 Gene 54 -0.5448046 -1.5837617 0.152121135 0.9913757 -5.0988743
> 16 Gene 16 -0.4970755 -1.7635064 0.116040512 0.9913757 -4.8601096
> 94 Gene 94 -0.4921761 -2.2254793 0.056881214 0.9913757 -4.2056166
> 93 Gene 93 -0.4574475 -1.3155906 0.224961686 0.9913757 -5.4296133
> 87 Gene 87 -0.4030807 -2.1728549 0.061734408 0.9913757 -4.2822231
> 92 Gene 92 -0.3708210 -1.3967269 0.200232145 0.9913757 -5.3332685
> 38 Gene 38 -0.3693720 -1.6087183 0.146559311 0.9913757 -5.0664427
> 41 Gene 41  0.3296712  0.3589144 0.728998939 0.9913757 -6.2080737
> 51 Gene 51 -0.3236529 -1.0426366 0.327764650 0.9913757 -5.7250738
> 40 Gene 40  0.3003745  1.3450728 0.215687305 0.9913757 -5.3950122
> 83 Gene 83  0.2948265  1.5802816 0.152911783 0.9913757 -5.1033766
> 47 Gene 47  0.2914631  1.4141449 0.195242647 0.9913757 -5.3121375
> 39 Gene 39 -0.2828721 -1.5449138 0.161161348 0.9913757 -5.1488494
> 69 Gene 69  0.2757269  1.2502012 0.246750193 0.9913757 -5.5046099
> 19 Gene 19  0.2684281  1.2757686 0.238027563 0.9913757 -5.4755794
> 3   Gene 3 -0.2642721 -1.6381776 0.140233597 0.9913757 -5.0278441
>> design<- cbind(Grp1=c(1,1,1,0,0,0),Grp2=c(0,0,0,1,1,1))
>> cont.matrix<- makeContrasts(G1vsG2=Grp1-Grp2, levels=design)
>> fit1<- lmFit(y,design)
>> fit1<- contrasts.fit(fit1, cont.matrix)
>> fit1<- eBayes(fit1)
>> topTable(fit1,coef=1,n=20,sort.by="logFC")
>          ID      logFC          t      P.Value   adj.P.Val          B
> 2   Gene 2 -2.0639600 -7.2886732 8.728373e-05 0.008728373  1.9349063
> 1   Gene 1 -1.9752877 -5.0131716 1.053468e-03 0.035115593 -0.6626835
> 33 Gene 33  1.9059513  5.1993066 8.380943e-04 0.035115593 -0.4227207
> 23 Gene 23 -1.6921631 -2.1549881 6.347344e-02 0.793418043 -4.8748108
> 93 Gene 93 -1.5367983 -3.1252258 1.421720e-02 0.313735106 -3.3751312
> 27 Gene 27  1.0978401  3.0600772 1.568676e-02 0.313735106 -3.4759186
> 54 Gene 54 -0.7961582 -1.6365658 1.405731e-01 0.949255469 -5.6234924
> 37 Gene 37 -0.6926892 -2.2573983 5.412387e-02 0.773198083 -4.7195118
> 42 Gene 42 -0.6895975 -0.9536791 3.683231e-01 0.952227637 -6.4283974
> 76 Gene 76  0.6887321  1.1603291 2.795521e-01 0.952227637 -6.2154077
> 87 Gene 87 -0.6277337 -2.3927605 4.383747e-02 0.730624439 -4.5120531
> 51 Gene 51  0.6097573  1.3889781 2.024874e-01 0.952227637 -5.9464611
> 41 Gene 41 -0.5810686 -0.4473240 6.665707e-01 0.952227637 -6.7979655
> 15 Gene 15  0.5801261  1.8016224 1.094946e-01 0.898763677 -5.3938078
> 7   Gene 7 -0.5487280 -0.4418824 6.703406e-01 0.952227637 -6.8005963
> 80 Gene 80 -0.5459118 -1.5641098 1.566352e-01 0.949255469 -5.7209503
> 39 Gene 39 -0.5139565 -1.9848398 8.263317e-02 0.875815090 -5.1285443
> 16 Gene 16 -0.5085620 -1.2758028 2.380161e-01 0.952227637 -6.0835464
> 95 Gene 95  0.5040813  0.9693600 3.609117e-01 0.952227637 -6.4133578
> 69 Gene 69  0.4815805  1.5440263 1.613734e-01 0.949255469 -5.7475589
> [/code]
>
> I want to know why the output are different. Which one should I believe?
>
> Thanks a lot.
>
> Yours sincerely,
>
> Jianhong Ou
>
> jianhong.ou at umassmed.edu
>
> _______________________________________________
> 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
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues 