[BioC] Questiond about the descrepancy between edgeR and limma-voom
Yunshun Chen
yuchen at wehi.EDU.AU
Mon Mar 3 01:17:43 CET 2014
Hi Ming,
The edgeR likelihood ratio test (glmLRT) is slightly liberal in general,
hence gives much more DEGs.
You can try the quasi-likelihood F-test in edgeR instead, and then compare
the results with limma-voom.
Also, the genewise dispersions can be estimated in one step using
estimateDisp().
See the followings:
> y <- estimateDisp(y, design)
> qlf <- glmQLFTest(y, design, con.matrix)
> show(summary(de <- decideTestsDGE(qlf)))
Regards,
Yunshun
----------------------------------------------------------------------
Message: 24. Questiond about the descrepancy between edgeR and limma-voom
(Ming Yi)
Hi, Dear all:
I am looking for advice for what I observed in my analysis results using
the same data but two very close related methods: limma-voom and edgeR.
for purpose of consolidation and comparison, since I try to get as accurate
DEGs lists for my next step downstream analysis, for the same dataset, I
used both edgeR and limma-voom originally try consolidate each other.
Surprisingly, as shown below the numbers of DEGs for 4 constrasts I am
interested, there are consistency between the results of normal contrasts:
e.g., for contrasts RasOnly.Normal_RasPNot.Normal and
RasP.Normal_RasPNot.Normal (both involved comparison of normal samples of
different types) in terms of numbers of DEGs (in very close ranges),
however, there are large discrepancy between the results of tumor contrasts:
e.g., RasP.Tumor_RasPNot.Tumor and RasOnly.Tumor_RasPNot.Tumor. The edgeR
result has about 1k(416+420) DEGs for RasP.Tumor_RasPNot.Tumor and
890+552=1.3k DEGs for RasOnly.Tumor_RasPNot.Tumor at FDR 5%, whereas
limma-voom method only picked up a bit more than 100 or less DEGs for the
same contrasts (almost no DEGs at the default cutoff FDR or adjusted p-value
level) shown below. Of course, maybe cutoffs (e.g., FDR vs adjusted.p-value)
or model may have different impacts in the DEGs sets,however, the odd thing
is: for the same cutoffs, we did see quite consistency in normal contrasts,
but large dsicrepancy in the tumor contrasts. From almost no DEG (or a bit
more than 100) to more than 1k DEGs seem a big difference to me. Since that
would impact the downstream analysis such as pathway analysis etc. The other
main difference is for the same data, which are claimed raw counts of RSEM
for RNAseq data, were rounded to integers for edgeR and also we used the
original RSEM raw count values for limma-voom. Both methods used the
identical makeContrasts commands shown below. Any suggestions and advice
would be highly appreciated. Thanks a lot in advance!
Ming
ATRF
NCI-Frederick
Maryland, USA
EdgR result:
currLrt<-glmLRT(fit, contrast=con.matrix[,colnames(con.matrix)[ii]]);
show(summary(de <- decideTestsDGE(currLrt)));
[1] "RasP.Tumor_RasPNot.Tumor"
-1 416
0 17414
1 420
[1] "RasOnly.Tumor_RasPNot.Tumor"
[,1]
-1 890
0 16808
1 552
[1] "RasOnly.Normal_RasPNot.Normal"
[,1]
-1 4
0 18229
1 17
[1] "RasP.Normal_RasPNot.Normal"
[,1]
-1 481
0 16936
1 833
limma-voom result:
> show(summary(de <- decideTests(fit)));
RasP.Tumor_RasPNot.Tumor RasOnly.Tumor_RasPNot.Tumor
-1 29 28
0 18202 18125
1 19 97
RasOnly.Normal_RasPNot.Normal RasP.Normal_RasPNot.Normal
-1 0 682
0 18250 17184
1 0 384
critcal commands:
EdgR:
> y <- estimateGLMCommonDisp(y, design, verbose=TRUE);
Disp = 0.42594 , BCV = 0.6526
> y <- estimateGLMTrendedDisp(y, design);
Loading required package: splines
> y <- estimateGLMTagwiseDisp(y, design)
> fit <- glmFit(y, design);
> con.matrix<-makeContrasts(
+ RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,
+ RasOnly.Tumor_RasPNot.Tumor=RasOnly.Tumor-RasPNot.Tumor,
+ RasOnly.Normal_RasPNot.Normal=RasOnly.Normal-RasPNot.Normal,
+ RasP.Normal_RasPNot.Normal=RasP.Normal-RasPNot.Normal,
+ levels=design)
currLrt<-glmLRT(fit, contrast=con.matrix[,colnames(con.matrix)[ii]]);
show(summary(de <- decideTestsDGE(currLrt)));
limma-voom:
>design <- model.matrix(~0+RasTum);
> con.matrix<-makeContrasts(
+ RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor,
+ RasOnly.Tumor_RasPNot.Tumor=RasOnly.Tumor-RasPNot.Tumor,
+ RasOnly.Normal_RasPNot.Normal=RasOnly.Normal-RasPNot.Normal,
+ RasP.Normal_RasPNot.Normal=RasP.Normal-RasPNot.Normal,
+ levels=design)
> fit<- lmFit(v,design)
> fit<-contrasts.fit(fit, con.matrix)
> fit <- eBayes(fit)
> show(summary(de <- decideTests(fit)));
[[alternative HTML version deleted]]
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}
More information about the Bioconductor
mailing list