Hi Ming,
Sorry, the decideTestsDGE() in edgeR won't give you the same output format
(one column for each contrast) as the decideTests() in limma-voom.
I think you have to take each column of your contrast matrix and do the
quasi-likelihood F-test separately, and then combine the results from each
test (like what you did with glmLRT).
Regards,
Yunshun
From: Ming Yi [mailto:yi02@hotmail.com]
Sent: Wednesday, 5 March 2014 4:18 AM
To: yuchen@wehi.EDU.AU
Cc: Bioconductor mailing list
Subject: RE: [BioC] Questiond about the descrepancy between edgeR and
limma-voom
Hi, YunShun:
I did try the quasi test using your code and did change some of my original
codes a bit (use the code you suggested and also I took out the three
estimate commands for GLM): here is what I did
>y <- estimateDisp(y, design)
> qlf <- glmQLFTest(y, design, con.matrix)
Error in glmFit.default(y = y$counts, design = design, dispersion =
dispersion, :
Length of dispersion vector incompatible with count matrix. Dispersion
argument must be either of length 1 (i.e. common dispersion) or length equal
to the number of rows of y (i.e. individual dispersion value for each
tag/gene).
I checked the manual of glmQLFTest and so I used the qlf <- glmQLFTest(y,
design,contrast= con.matrix), which seem working.
Then I tried the command you suggested: show(summary(de <-
decideTestsDGE(qlf))) and I got error message again, the qlf object seems
not work with the command show(summary(de <- decideTestsDGE(qlf))):
> show(summary(de <- decideTestsDGE(qlf)))
Error in show(summary(de <- decideTestsDGE(qlf))) :
error in evaluating the argument 'object' in selecting a method for
function 'show': Error in array(x, c(length(x), 1L), if (!is.null(names(x)))
list(names(x), :
'data' must be of a vector type, was 'NULL'
I also tried what I usually do in limma-voom and also got error:
> show(summary(de <- decideTests(qlf)))
Error in show(summary(de <- decideTests(qlf))) :
error in evaluating the argument 'object' in selecting a method for
function 'show': Error in decideTests(qlf) : Need MArrayLM object
Any idea? Thx and best
Ming
> From: yuchen@wehi.EDU.AU
> To: yi02@hotmail.com
> CC: bioconductor@r-project.org
> Subject: RE: [BioC] Questiond about the descrepancy between edgeR and
limma-voom
> Date: Mon, 3 Mar 2014 11:17:43 +1100
>
> 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 inte...{{dropped:15}}