Dear All, I am facing a problem that has been discussed previously between Bjorn Usadel and Gordon Smyth on the 9-11th Nov 2005. It is about producing two different topTables using the same initial eset, but trying with 2 different design lay outs. However my problem seems to be exaxctly the same as the one stated by Bjorn, it isn't fixed if I follow Gordon's indications to overcome the "probesets with zero variance" problem: I have 4 chips, 2 biological replicates of condition A and condition B on a custom affymetrix array. R version 2.2.0, 2005-10-06, i386-pc-mingw32 limma "2.2.0" Biobase "1.8.0" Here is what I've done: > design2 <- cbind(A=1, BvsA=c(0,0,1,1)) > fit2 <- lmFit(eset, design2) > fit2eB <- eBayes(fit2) > toptable(fit2eB, n=10) M t P.Value B 853 12.40696 143.0466 5.357015e-06 8.649517 974 10.69453 133.1012 5.357015e-06 8.580885 4813 10.49636 129.6253 5.357015e-06 8.553634 4718 11.01015 128.9608 5.357015e-06 8.548209 3587 10.24772 126.1652 5.357015e-06 8.524578 3812 12.35814 126.1029 5.357015e-06 8.524036 2170 11.39107 125.5051 5.357015e-06 8.518801 4265 12.62847 123.7714 5.357015e-06 8.503256 3372 11.47227 123.5476 5.357015e-06 8.501209 345 9.98423 123.4283 5.357015e-06 8.500113 however if I do the following: > design <- model.matrix(~ -1+factor(c(1,1,2,2))) > colnames(design) <- c("A", "B") > fit1 <- lmFit(eset, design) > contrast.matrix <- makeContrasts(A-B, levels=design) > fit12 <- contrasts.fit(fit1, contrast.matrix) > fit1eB <- eBayes(fit12) > toptable(fit1eB, n=10) M t P.Value B 311 -3.997113 -13.961019 0.5794648 -0.5112598 1327 -1.461801 -11.334987 0.5794648 -0.6889117 113 -1.690073 -10.880814 0.5794648 -0.7308602 4408 -3.066882 -10.019535 0.5794648 -0.8232774 1825 -3.576034 -9.781223 0.5794648 -0.8523026 1224 -1.099785 -9.445264 0.5794648 -0.8961448 289 -2.800736 -9.306995 0.5794648 -0.9152559 288 -1.689312 -8.759499 0.5794648 -0.9977157 2892 -2.513426 -8.675603 0.5794648 -1.0113879 3311 3.005392 8.377018 0.5794648 -1.0625077 Following the instructions given by Gordon, then I looked at: > i <- (fit2eB$sigma==0) # same result with fit1eB > sum(i) [1] 0 so the probe-sets with zero variance doesn't seem to be the reason here... I, of course, would be tempted to believe the 1st option (giving differentially expressed genes with B > 8) but it turns out that 96% of the genes are differentially expressed in this 1st option, which is quite unlikely! I can not understand why is it so. Any suggestions and/or indication of what I may have done wrong would be gratefully appreciated. All the best, Celine -- Celine Carret PhD Pathogen Microarrays group The Wellcome Trust Sanger Institute Hinxton, Cambridge CB10 1SA, UK. tel. +44 (0)1223 834 244 ext.7123 fax. +44 (0)1223 494 919 email: ckc at sanger.ac.uk http://www.sanger.ac.uk/PostGenomics/PathogenArrays/

> Date: Mon, 05 Dec 2005 19:23:30 +0000 > From: Celine Carret <ckc at="" sanger.ac.uk=""> > Subject: [BioC] (Limma) different toptable results for the same > dataset using 2 different designs > To: bioconductor at stat.math.ethz.ch > > Dear All, > > I am facing a problem that has been discussed previously between Bjorn > Usadel and Gordon Smyth on the 9-11th Nov 2005. > It is about producing two different topTables using the same initial > eset, but trying with 2 different design lay outs. However my problem > seems to be exaxctly the same as the one stated by Bjorn, it isn't fixed > if I follow Gordon's indications to overcome the "probesets with zero > variance" problem: > I have 4 chips, 2 biological replicates of condition A and condition B > on a custom affymetrix array. > R version 2.2.0, 2005-10-06, i386-pc-mingw32 > limma "2.2.0" > Biobase "1.8.0" > > Here is what I've done: > > design2 <- cbind(A=1, BvsA=c(0,0,1,1)) > > fit2 <- lmFit(eset, design2) > > fit2eB <- eBayes(fit2) > > toptable(fit2eB, n=10) topTable(fit2eB, coef=2) Gordon > M t P.Value B > 853 12.40696 143.0466 5.357015e-06 8.649517 > 974 10.69453 133.1012 5.357015e-06 8.580885 > 4813 10.49636 129.6253 5.357015e-06 8.553634 > 4718 11.01015 128.9608 5.357015e-06 8.548209 > 3587 10.24772 126.1652 5.357015e-06 8.524578 > 3812 12.35814 126.1029 5.357015e-06 8.524036 > 2170 11.39107 125.5051 5.357015e-06 8.518801 > 4265 12.62847 123.7714 5.357015e-06 8.503256 > 3372 11.47227 123.5476 5.357015e-06 8.501209 > 345 9.98423 123.4283 5.357015e-06 8.500113 > > however if I do the following: > > design <- model.matrix(~ -1+factor(c(1,1,2,2))) > > colnames(design) <- c("A", "B") > > fit1 <- lmFit(eset, design) > > contrast.matrix <- makeContrasts(A-B, levels=design) > > fit12 <- contrasts.fit(fit1, contrast.matrix) > > fit1eB <- eBayes(fit12) > > toptable(fit1eB, n=10) > M t P.Value B > 311 -3.997113 -13.961019 0.5794648 -0.5112598 > 1327 -1.461801 -11.334987 0.5794648 -0.6889117 > 113 -1.690073 -10.880814 0.5794648 -0.7308602 > 4408 -3.066882 -10.019535 0.5794648 -0.8232774 > 1825 -3.576034 -9.781223 0.5794648 -0.8523026 > 1224 -1.099785 -9.445264 0.5794648 -0.8961448 > 289 -2.800736 -9.306995 0.5794648 -0.9152559 > 288 -1.689312 -8.759499 0.5794648 -0.9977157 > 2892 -2.513426 -8.675603 0.5794648 -1.0113879 > 3311 3.005392 8.377018 0.5794648 -1.0625077 > > Following the instructions given by Gordon, then I looked at: > > i <- (fit2eB$sigma==0) # same result > with fit1eB > > sum(i) > [1] 0 > > so the probe-sets with zero variance doesn't seem to be the reason here... > I, of course, would be tempted to believe the 1st option (giving > differentially expressed genes with B > 8) but it turns out that 96% of > the genes are differentially expressed in this 1st option, which is > quite unlikely! > I can not understand why is it so. > Any suggestions and/or indication of what I may have done wrong would > be gratefully appreciated. > > All the best, > Celine > > > > > -- > Celine Carret PhD > Pathogen Microarrays group > The Wellcome Trust Sanger Institute > Hinxton, Cambridge CB10 1SA, UK. > tel. +44 (0)1223 834 244 ext.7123 > fax. +44 (0)1223 494 919 > email: ckc at sanger.ac.uk > http://www.sanger.ac.uk/PostGenomics/PathogenArrays/