Hello,
This might have been asked before, but I couldn't find the answer. So I try to ask here what the difference between voom and trend is for RNAseq analysis.
I have a study with four groups, and a batch effect. The lib sizes are between 8M to 11M reads. In the user's guide it says that when the lib sizes are not more than 3 fold different, trend is the method to choose. So when I do this I get significant genes. However, when I do the analysis with voom, I get many more significantly expressed genes. Why is that, and how reliable are these extra genes? Does this mean that the extra genes from limma voom are false positive?
trend
> DGE1 <- DGEList(counts, group=Group) > keep <- rowSums(cpm(DGE1)>1) >= 3 > DGE2 <- DGE1[keep,] > y1 <- calcNormFactors(DGE2) > logCPM <- cpm(y1, log=TRUE, prior.count=0.5) > design <- model.matrix(~0 + targets$Group + targets$Donor) > fit <- lmFit(logCPM, design) > fit2 <- eBayes(contrasts.fit(fit, contr)) > summary(decideTests(fit2, method="separate")) PPvsPN PPvsNN PPvsNP PNvsNN PNvsNP NPvsNN -1 54 140 86 27 433 137 0 15897 15654 15796 16039 15095 15803 1 154 311 223 39 577 165
voom
> DGE1 <- DGEList(counts_minus, group=Group) > keep <- rowSums(cpm(DGE1)>1) >= 3 > DGE2 <- DGE1[keep,] > y1 <- calcNormFactors(DGE2) > design <- model.matrix(~0 + targets$Group + targets$Donor) > v <- voom(y1, design, plot=TRUE) > fit <- lmFit(v, design) > fit2 <- eBayes(contrasts.fit(fit, contr)) > summary(decideTests(fit2, method="separate")) PPvsPN PPvsNN PPvsNP PNvsNN PNvsNP NPvsNN -1 363 579 428 277 873 465 0 15166 14647 15023 15469 14252 15087 1 576 879 654 359 980 553
Well, that's embarrassing... Forgetting the "trend=T" part in the limma trend analysis! Must have something to do with Mondays...
Thanks for your quick help, the figures are more alike now with the real trend included:
As you've discovered, voom and limma-trend tend to give about the same number of DE genes, although it will vary from one dataset to another.