[R] problems with limma

Tue Dec 21 00:21:02 CET 2004

It probably bounces because the email address is incorrect (it should be a
'w', not a 'v'). /Henrik Bengtsson

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of 
> r.ghezzo at staff.mcgill.ca
> Sent: Monday, December 20, 2004 9:45 PM
> To: r-help at stat.math.ethz.ch
> 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
> 
>