[BioC] contrast.fit and different p-values for same comparison
Naomi Altman
naomi at stat.psu.edu
Fri Feb 12 21:47:52 CET 2010
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
More information about the Bioconductor
mailing list