[BioC] Questiond about the descrepancy between edgeR and limma-voom

[Yunshun Chen]

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}}