contrast.fit and different p-values for same comparison
1
0
Entering edit mode
@lakshmanan-iyer-1829
Last seen 9.3 years ago
United States
Hi I used limma... topTable on the "same contrast" for a complete dataset (36 arrays, 12 conditions, 3 replicates per condition) and a subset of the same (18 arrays). I used lmFit, followed by contrasts.fit and eBayes in both cases. However, I get different Ave Expression and P-values in the two cases. Why is this happening? I cannot readily explain this. Thanks for any help/clues. The first row is full dataset second row is the subset > rbind (x[x[,1]=="ILMN_2772632",],x1[x1[,1]=="ILMN_2772632",]) ID logFC AveExpr t P.Value adj.P.Val 37561 ILMN_2772632 4.198562 10.40451 22.86259 2.007820e-21 9.091611e-17 375611 ILMN_2772632 4.198562 9.85615 24.89322 1.224698e-17 1.631045e-14 B 37561 37.44148 375611 29.95781 ---------------------------------------- -best -Lax Here is the snippet of code with some output. library (limma); eset <- read.delim("results/expressions.tsv",sep="\t",stringsAsFactors=F) rownames(eset) <- eset[,1] eset <- eset[,-1] names(eset) <- gsub("^X","",names(eset)) samples <- read.delim("mySamples.tsv",sep="\t",stringsAsFactors=F) TS <- factor(samples$Condition, levels=unique(samples$Condition)) design <- model.matrix(~0 + TS) colnames(design) <- levels(TS) cont.matrix <- makeContrasts(Day1AgedPnxVsCtrl = Day1AgedPnx - Day1AgedControl, Day3AgedPnxVsCtrl = Day3AgedPnx - Day3AgedControl, Day7AgedVsCtrl = Day7AgedPnx - Day7AgedControl, Day1YoungPnxVsCtrl = Day1YoungPnx - Day1YoungControl, Day3YoungPnxVsCtrl = Day3YoungPnx - Day3YoungControl, Day7YoungPnxVsCtrl =Day7YoungPnx - Day7YoungControl, levels=design) fit <- lmFit(eset, design) fit2 <- contrasts.fit(fit,cont.matrix) fit2 <- eBayes(fit2) x <- topTable(fit2, coef="Day1AgedPnxVsCtrl",num=dim(eset)[1]) # TS1 <- factor(samples$Condition[1:18], levels=unique(samples$Condition[1:18])) design1 <- model.matrix(~0 + TS1) colnames(design1) <- levels(TS1) cont.matrix1 <- makeContrasts(Day1AgedPnxVsCtrl = Day1AgedPnx - Day1AgedControl, Day1YoungPnxVsCtrl = Day1YoungPnx - Day1YoungControl, levels=design1) fit1 <- lmFit(eset[,1:18], design1) fit21 <- contrasts.fit(fit1,cont.matrix1) fit21 <- eBayes(fit21) topTable(fit21, coef="Day1AgedPnxVsCtrl") x1 <- topTable(fit21, coef="Day1AgedPnxVsCtrl",num=dim(eset)[1]) > dim (eset) [1] 45281 36 > dim (design) [1] 36 12 > dim (cont.matrix) [1] 12 6 > dim (design1) [1] 18 6 > dim (cont.matrix1) [1] 6 2 > rbind (x[x[,1]=="ILMN_2772632",],x1[x1[,1]=="ILMN_2772632",]) ID logFC AveExpr t P.Value adj.P.Val 37561 ILMN_2772632 4.198562 10.40451 22.86259 2.007820e-21 9.091611e-17 375611 ILMN_2772632 4.198562 9.85615 24.89322 1.224698e-17 1.631045e-14 B 37561 37.44148 375611 29.95781 ------------- > sessionInfo() R version 2.10.1 (2009-12-14) i486-pc-linux-gnu locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] limma_3.2.1 loaded via a namespace (and not attached): [1] tools_2.10.1 >
limma limma • 1.1k views
ADD COMMENT
0
Entering edit mode
Naomi Altman ★ 6.0k
@naomi-altman-380
Last seen 3.7 years ago
United States
The point of eBayes is to use all of the available data. If you use only 18 arrays, the available data are only these arrays, so your results differ. Even if you do not use eBayes, the MSE used for the statistical tests will depend on the arrays included in the analysis. Regards, Naomi At 02:17 PM 2/12/2010, Lakshmanan Iyer wrote: >Hi >I used limma... topTable on the "same contrast" for a complete dataset >(36 arrays, 12 conditions, 3 replicates per condition) and a subset of >the same (18 arrays). >I used lmFit, followed by contrasts.fit and eBayes in both cases. > >However, I get different Ave Expression and P-values in the two cases. > >Why is this happening? I cannot readily explain this. > >Thanks for any help/clues. > >The first row is full dataset >second row is the subset > > > rbind (x[x[,1]=="ILMN_2772632",],x1[x1[,1]=="ILMN_2772632",]) > ID logFC AveExpr t P.Value adj.P.Val >37561 ILMN_2772632 4.198562 10.40451 22.86259 2.007820e-21 9.091611e-17 >375611 ILMN_2772632 4.198562 9.85615 24.89322 1.224698e-17 1.631045e-14 > B >37561 37.44148 >375611 29.95781 >---------------------------------------- >-best >-Lax >Here is the snippet of code with some output. > >library (limma); >eset <- read.delim("results/expressions.tsv",sep="\t",stringsAsFactors=F) >rownames(eset) <- eset[,1] >eset <- eset[,-1] >names(eset) <- gsub("^X","",names(eset)) >samples <- read.delim("mySamples.tsv",sep="\t",stringsAsFactors=F) >TS <- factor(samples$Condition, levels=unique(samples$Condition)) >design <- model.matrix(~0 + TS) >colnames(design) <- levels(TS) >cont.matrix <- makeContrasts(Day1AgedPnxVsCtrl = Day1AgedPnx - >Day1AgedControl, Day3AgedPnxVsCtrl = Day3AgedPnx - Day3AgedControl, >Day7AgedVsCtrl = Day7AgedPnx - Day7AgedControl, Day1YoungPnxVsCtrl = >Day1YoungPnx - Day1YoungControl, Day3YoungPnxVsCtrl = Day3YoungPnx - >Day3YoungControl, Day7YoungPnxVsCtrl =Day7YoungPnx - Day7YoungControl, >levels=design) > >fit <- lmFit(eset, design) >fit2 <- contrasts.fit(fit,cont.matrix) >fit2 <- eBayes(fit2) >x <- topTable(fit2, coef="Day1AgedPnxVsCtrl",num=dim(eset)[1]) ># >TS1 <- factor(samples$Condition[1:18], levels=unique(samples$Condition[1:18])) >design1 <- model.matrix(~0 + TS1) >colnames(design1) <- levels(TS1) >cont.matrix1 <- makeContrasts(Day1AgedPnxVsCtrl = Day1AgedPnx - >Day1AgedControl, Day1YoungPnxVsCtrl = Day1YoungPnx - Day1YoungControl, >levels=design1) >fit1 <- lmFit(eset[,1:18], design1) >fit21 <- contrasts.fit(fit1,cont.matrix1) >fit21 <- eBayes(fit21) >topTable(fit21, coef="Day1AgedPnxVsCtrl") >x1 <- topTable(fit21, coef="Day1AgedPnxVsCtrl",num=dim(eset)[1]) > > > dim (eset) >[1] 45281 36 > > dim (design) >[1] 36 12 > > dim (cont.matrix) >[1] 12 6 > > dim (design1) >[1] 18 6 > > dim (cont.matrix1) >[1] 6 2 > > rbind (x[x[,1]=="ILMN_2772632",],x1[x1[,1]=="ILMN_2772632",]) > ID logFC AveExpr t P.Value adj.P.Val >37561 ILMN_2772632 4.198562 10.40451 22.86259 2.007820e-21 9.091611e-17 >375611 ILMN_2772632 4.198562 9.85615 24.89322 1.224698e-17 1.631045e-14 > B >37561 37.44148 >375611 29.95781 > > >------------- > > sessionInfo() >R version 2.10.1 (2009-12-14) >i486-pc-linux-gnu > >locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C >[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > >attached base packages: >[1] stats graphics grDevices utils datasets methods base > >other attached packages: >[1] limma_3.2.1 > >loaded via a namespace (and not attached): >[1] tools_2.10.1 > > > >_______________________________________________ >Bioconductor mailing list >Bioconductor at stat.math.ethz.ch >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: >http://news.gmane.org/gmane.science.biology.informatics.conductor Naomi S. Altman 814-865-3791 (voice) Associate Professor Dept. of Statistics 814-863-7114 (fax) Penn State University 814-865-1348 (Statistics) University Park, PA 16802-2111
ADD COMMENT

Login before adding your answer.

Traffic: 354 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6