On Wed, December 22, 2004 12:11 am, r.ghezzo at staff.mcgill.ca
said:> ----- Forwarded message from r.ghezzo at staff.mcgill.ca -----
> Date: Mon, 20 Dec 2004 15:45:11 -0500
> From: r.ghezzo at staff.mcgill.ca
> Reply-To: r.ghezzo at staff.mcgill.ca
> Subject: [R] problems with limma
> To: r-help at stat.math.ethz.ch
>
> 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
Questions about limma should be sent to the Bioconductor mailing list (the
address is given in the
introduction of the Limma User's Guide for example).
> 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?
Your commands look correct but your data looks crazy. The M-values are supposed
to be log-fold
changes, and you're getting values in the 10s of thousands. Perhaps you are
trying analyse data
on the raw intensity scale and have huge differences between arrays and groups
or very large
outliers. Note that the B-statistics are telling you that there is absolutely
no differential
expression throughout. Warning bells should go up when you see such large
t-statistics with such
small B-statistics. Please look at your data and do some quality assessment.
At very least you
probably need to log-transform.
Your 4-group approach should give the same M-values as the 2-group approach, but
the standard
errors and t-statistics will change because your error estimates change.
Gordon
> This is data from Affimetrix mouse chips.
> Thanks for any help
> Heberto Ghezzo
> Ph.D.
> Meakins-Christie Labs
> McGill University
> Montreal - Canada
