[R] problems with limma

Mon Dec 20 21:45:11 CET 2004

I try to send this message To Gordon Smyth at smyth at vehi,edu.au but it bounced
back, so here it is to r-help

I am trying to use limma, just downloaded it from CRAN. I use R 2.0.1 on Win XP
see the following:
> library(RODBC)
> chan1 <- odbcConnectExcel("D:/Data/mgc/Chips/Chips4.xls")
> dd <- sqlFetch(chan1,"Raw")   # all data  12000
> #
> nzw <- cbind(dd$NZW1C,dd$NZW2C,dd$NZW3C,dd$NZW1T,dd$NZW2T,dd$NZW3T)
> akr <- cbind(dd$AKR1C,dd$AKR2C,dd$AKR3C,dd$AKR1T,dd$AKR2T,dd$AKR3T)
> bas <- cbind(dd$NZW1C,dd$NZW2C,dd$NZW3C,dd$AKR1C,dd$AKR2C,dd$AKR3C)
> #
>  design<-matrix(c(1,1,1,1,1,1,0,0,0,1,1,1),ncol=2)
>  fit1 <- lmFit(nzw,design)
>  fit1 <- eBayes(fit1)
>  topTable(fit1,adjust="fdr",number=5)
              M         t      P.Value         B
12222  3679.480 121.24612 7.828493e-06 -4.508864
1903   3012.405 118.32859 7.828493e-06 -4.508866
9068   1850.232  92.70893 1.178902e-05 -4.508889
10635  2843.534  91.99336 1.178902e-05 -4.508890
561   18727.858  90.17085 1.178902e-05 -4.508893
> #
>  fit2 <- lmFit(akr,design)
>  fit2 <- eBayes(fit2)
>  topTable(fit2,adjust="fdr",number=5)
              M        t      P.Value         B
88     1426.738 80.48058 5.839462e-05 -4.510845
1964  36774.167 73.05580 5.839462e-05 -4.510861
5854   7422.578 68.60316 5.839462e-05 -4.510874
11890  1975.316 66.54480 5.839462e-05 -4.510880
9088   2696.952 64.16343 5.839462e-05 -4.510889
> #
>  fit3 <- lmFit(bas,design)
>  fit3 <- eBayes(fit3)
>  topTable(fit3,adjust="fdr",number=5)
             M         t      P.Value         B
6262  1415.088 100.78933 2.109822e-05 -4.521016
5660  1913.479  96.40903 2.109822e-05 -4.521020
11900 4458.489  94.30738 2.109822e-05 -4.521022
9358  1522.330  80.46641 3.346749e-05 -4.521041
11773 1784.483  73.76620 3.346749e-05 -4.521053
> #    Now lets do all together in Anova
> #
>  all <- cbind(nzw,akr)
>  ts <- c(1,1,1,2,2,2,3,3,3,4,4,4)
>  ts <- as.factor(ts)
>  levels(ts) <- c("nzwC","nzwT","akrC","akrT")
>  design <- model.matrix(~0+ts)
>  colnames(design) <- levels(ts)
>  fit4 <- lmFit(all,design)
>  cont.matrix <- makeContrasts(
+      Baseline = akrC - nzwC,
+      NZW_Smk = nzwT - nzwC,
+      AKR_Smk = akrT - akrC,
+      Diff = (akrT - akrC) - (nzwT - nzwC),
+      levels=design)
>   fit42 <- contrasts.fit(fit4,cont.matrix)
>   fit42 <- eBayes(fit42)
> #
>   topTable(fit42,coef="Baseline",adjust="fdr",number=5)
               M         t     P.Value         B
3189    942.0993  13.57485 0.004062283 -4.528799
8607   2634.1826  11.23476 0.006913442 -4.530338
10242  -942.2860 -10.99253 0.006913442 -4.530551
283    -609.0831 -10.79354 0.006913442 -4.530735
3224  -1564.2572 -10.19429 0.008089034 -4.531351
----------------------------------------------------
------------- Shouldn't this be equal to fit1 above?
----------------------------------------------------
>   topTable(fit42,coef="NZW_Smk",adjust="fdr",number=5)
             M         t   P.Value         B
7724 -246.5956 -8.687324 0.1615395 -4.591133
1403 -307.8660 -7.063312 0.4066814 -4.591363
3865 -253.4899 -6.585582 0.4598217 -4.591457
3032 -509.2413 -5.841901 0.8294166 -4.591640
2490 -240.3259 -5.338679 0.9997975 -4.591795
----------------------------------------------------
------------- Shouldn't this be equal to fit2 above?
------------- The P.Value are unreal!!
----------------------------------------------------
>   topTable(fit42,coef="AKR_Smk",adjust="fdr",number=5)
             M        t  P.Value         B
11547 151.6622 6.380978 0.917470 -4.595085
12064 324.0851 6.337235 0.917470 -4.595085
6752  964.5478 5.858994 0.952782 -4.595086
10251 152.7587 5.339843 0.952782 -4.595087
1440  189.6056 4.933151 0.952782 -4.595089
----------------------------------------------------
------------- Shouldn't this be equal to fit3 above?
------------- The P.Value are unreal!!
----------------------------------------------------
>   topTable(fit42,coef="Diff",adjust="fdr",number=5)
              M         t   P.Value         B
7724   302.6892  7.540195 0.4102211 -4.593201
1403   419.4962  6.805495 0.4102211 -4.593265
10251  270.5269  6.686796 0.4102211 -4.593277
3270   409.8391  6.414966 0.4192042 -4.593307
10960 -511.4711 -5.469247 0.9652171 -4.593435
> #
>
So the results I get from just pairwise comparisons are very significant, but
when I try the Anova way, the significance completely dissapears.
Am I doing something completely wrong?
This is data from Affimetrix mouse chips.
Thanks for any help
Heberto Ghezzo
Ph.D.
Meakins-Christie Labs
McGill University
Montreal - Canada