[BioC] contrast.fit and different p-values for same comparison
Lakshmanan Iyer
laxvid at gmail.com
Fri Feb 12 20:17:54 CET 2010
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
>
More information about the Bioconductor
mailing list