Dear Lev,
One of the purposes of a mailing list is that many people may reply,
so please do not address your question specifically to me.
As explained in Section 10.1 of the limma User's Guide, the
B-statistic requires a prior estimate of the proportion of DE genes.
By default this is set to 1%.
Therefore, the B-statistic will tend to underestimate significance if
the true proportion of DE genes is actually more than 1% and
overestimate if the true proportion is less than 1%. In your case,
the proportion of DE genes appears to be massively more than 1%,
hence you'd expect the B-statistic to underestimate significance.
That is, you'd expect the B-statistics to be too small.
Since the moderated t does not require a prior estimate, you'd expect
the p-values to suggest more DE than the B-statistics whenever the
proportion of DE genes in your data is large.
Best wishes
Gordon
At 12:27 AM 14/07/2007, Lev Soinov wrote:
>Dear Gordon,
>
>We are analysing a dataset of 14920 genes obtained with the AB1700
platform.
>It has three treatments L1, L2, L1+L2 and control. The data is in
>the form of expression data matrix with the first column as pobe ID
>and 14 other columns correspond to 4 above conditions. Using the
>code below, we obtain a huge number of genes with adjusted p values
><0.05, about 5000 for the comparison between L1 and control for
>example. At the same time B values corresponding to these probes are
>very small, i.e. we are getting B<-4 in the bottom of the list of
>probes with adj.p<0.05.
>Could you please comment on possible causes for this? Is it normal?
>With kind regards,
>Lev.
>
> >s<-scan('Data.txt',what='character')
>Read 223800 items
> > sm<-matrix(s,byrow=TRUE,ncol=15)
> > dim(sm)
>[1] 14920 15
> > rownames(sm)<-sm[,1]
> > sm<-sm[,2:ncol(sm)]
> > snn<-apply(sm,2,as.numeric)
> > rownames(snn)<-rownames(sm)
> > signals<-snn
> > dim(signals)
>[1] 14920 14
> > temp<-normalizeBetweenArrays(log2(signals), method='quantile')
> > design <- model.matrix(~0 +factor(c(1,1,1,1,2,2,2,3,3,3,3,4,4,4)))
> > colnames(design) <- c("Control","L1","L2","L1L2")
> > contrast.matrix <-
> makeContrasts(L1-Control,L2-Control,L1L2-Control,levels=design)
> > fit <- lmFit(temp, design)
> > fit2 <- contrasts.fit(fit, contrast.matrix)
> > fit2 <- eBayes(fit2)
> > topTable(fit2, coef=1, adjust='BH')
> ID logFC t P.Value adj.P.Val B
>8790 182417 5.813459 38.16912 1.072876e-15 1.479057e-11 25.47446
>6945 165482 8.573261 35.59768 2.856130e-15 1.479057e-11 24.69238
>7132 167208 6.247484 35.49523 2.973975e-15 1.479057e-11 24.65950
>10941 202780 4.881978 33.98673 5.467499e-15 2.039377e-11 24.15906
>1102 109858 5.076380 33.01348 8.214210e-15 2.451120e-11 23.81910
>3785 135458 3.686867 32.09869 1.217373e-14 3.027202e-11 23.48654
>12355 215284 5.035617 30.68240 2.288454e-14 4.877676e-11 22.94512
>6515 161539 5.885744 29.64292 3.704033e-14 6.908021e-11 22.52582
>9789 191745 8.189347 28.65188 5.953817e-14 9.870106e-11 22.10752
>8568 180293 4.749725 27.88955 8.671664e-14 1.293812e-10 21.77270
>
>
>The bottom of the adj.p<0.05 list:
> ID logFC t P.Value
> adj.P.Val B
>207673 -0.293302 -2.70223 0.01703 0.042251 -4.2848
>213498 -0.675519 -2.70219 0.017033 0.042251 -4.2849
>186148 -0.419934 -2.70201 0.017039 0.042258 -4.2853
>233859 -0.422533 -2.70185 0.017044 0.042263 -4.2856
>188263 -0.330067 -2.70179 0.017046 0.042263 -4.2857
