problems with limma

Could you give a bit more detail about your experimental design?  You're using affy, so you're working with single channel data - so nzw, akr and bas all have six arrays?


-----Original Message-----
From:	r-help-bounces at stat.math.ethz.ch on behalf of r.ghezzo at staff.mcgill.ca
Sent:	Mon 12/20/2004 8:45 PM
To:	r-help at stat.math.ethz.ch
Cc:	
Subject:	[R] problems with limma
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

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html