Entering edit mode
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}}
