[BioC] limma small vs large number of samples

Ryan rct at thompsonclan.org
Tue Apr 29 04:29:45 CEST 2014


Hi Giovanni,

The explanation is quite simple. When you include all the data, limma 
can arrive at a much more robust estimate of the variance, so it can 
more confidently call differentially expressed genes. In other words, 
regardless of which specific contrast you are testing, the variance is 
estimated from all of the data included in the model.

-Ryan

On Mon Apr 28 17:54:32 2014, Giovanni Bucci wrote:
> Hello everybody,
>
> I have 32 samples, 4 factors with 2 levels each. Each level has 2
> replicates.
>
>> str(gxexprs)
>   num [1:15584, 1:32] 7.94 6.67 9.93 9.62 12.19 ...
>
>> Group
>   [1] R52VQ R52VQ R52VE R52VE R52EQ R52EQ R52EE R52EE R95VQ R95VQ R95VQ R95VE
> [13] R95VE R95VE R95EQ R95EQ R95EQ R95EE R95EE R95EE R97VQ R97VQ R97VQ R97VE
> [25] R97VE R97VE R97EQ R97EQ R97EQ R97EE R97EE R97EE
> 16 Levels: R52VQ R52VE R52EQ R52EE R95VQ R95VE R95EQ R95EE R97VQ ... R97EE
>
>
> design <- model.matrix(~0+Group)
> fit <- lmFit(gxexprs, design)
>
> contrast.matrix <- makeContrasts(contrasts="R52VQ - R52VE",levels=design)
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)
>
> TTable = topTable(fit2)
>
> global_p_val = TTable$P.Val
>
>
> gxexprs = gxexprs[, 1:4]
>
>
> ## same code as above but the expression matrix has only the first 4
> columns which represent the contrast tested above
>
> design <- model.matrix(~0+Group)
> fit <- lmFit(gxexprs, design)
>
> contrast.matrix <- makeContrasts(contrasts="R52VQ - R52VE",levels=design)
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)
>
> TTable = topTable(fit2)
>
> local_p_val = TTable$P.Val
>
> local_p_val has much greater values than global_p_val even though they
> represent the same comparison.
>
> What is the explanation for this?
>
> Can you point to some diagnostic functions that will show the difference?
>
> Thank you,
>
> Giovanni
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor



More information about the Bioconductor mailing list