Dear Alpesh, I've commited a change to edgeR today, so that you can now produce the qq-plot of the goodness of fit statistics by gof(fit, plot=TRUE) Best wishes Gordon --------------------------------------------- Professor Gordon K Smyth, Bioinformatics Division, Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, Vic 3052, Australia. smyth at wehi.edu.au http://www.wehi.edu.au http://www.statsci.org/smyth On Wed, 23 May 2012, Gordon K Smyth wrote: > Dear Alpesh, > > What you're doing looks correct. > > However, if you are doing this analysis (Figure 2 from > http://nar.oxfordjournals.org/content/40/10/4288 ) to decide whether you need > to use tagwise dispersion for your data, note that we have never observed a > real data set for which tagwise estimation is not required. Hence it seems > unnecessary to make these plots as a routine diagnostic. > > Yes, dispersion=0 is Poisson. > > To plot highlighted blue points, see ?points as well as > > xy <- qqnorm(pcom) > > etc. > > Best wishes > Gordon > >> Date: Mon, 21 May 2012 10:15:37 -0500 >> From: Alpesh Querer <alpeshq at="" gmail.com=""> >> To: Bioconductor mailing list <bioconductor at="" r-project.org=""> >> Subject: [BioC] edgeR, comparing models >> >> Hello edgeR patrons, >> >> I have RNA-seq data for multiple samples with biological replicates, I want >> to look at the goodness of fit for >> fitting Poisson and NB models used by edgeR for common, trended and >> tag-wise dispersion. Does setting dispersion=0 >> in glmFit use the Poisson model? Also I am using the following code to >> generate and compare qqplots for the models(figure 2 of >> the most recent edgeR manuscript),please let me know if I am using the >> appropriate method. Also how to plot the blue points (genes with poor fit)? >> >> y<-estimateGLMCommonDisp(y,design) >> >> fitcommon <- glmFit(y,design) >> >> y <- estimateGLMTrendedDisp(y,design) >> >> fittrended <- glmFit(y,design) >> >> y <- estimateGLMTagwiseDisp(y,design) >> >> >> fittagwise <- glmFit(y,design) >> >> gcom <- gof(fitcommon) >> gtren <- gof(fittrended) >> gtag <- gof(fittagwise) >> >> zcom <- zscoreGamma(gcom$gof.statistics,shape=gcom$df/2,scale=2) >> ztren <- zscoreGamma(gtren$gof.statistics,shape=gtren$df/2,scale=2) >> ztag <- zscoreGamma(gtag$gof.statistics,shape=gtag$df/2,scale=2) >> >> >> pcom=zcom[is.finite(zcom) ] >> ptren=ztren[is.finite(ztren) ] >> ptag=ztag[is.finite(ztag) ] >> >> par(mfrow=c(3,1)) >> >> qqnorm(pcom) >> qqline(pcom) >> >> qqnorm(ptren) >> qqline(ptren) >> >> qqnorm(ptag) >> qqline(ptag) >> >> >>> sessionInfo() >> R version 2.15.0 (2012-03-30) >> Platform: i386-pc-mingw32/i386 (32-bit) >> >> locale: >> [1] LC_COLLATE=English_United States.1252 >> [2] LC_CTYPE=English_United States.1252 >> [3] LC_MONETARY=English_United States.1252 >> [4] LC_NUMERIC=C >> [5] LC_TIME=English_United States.1252 >> >> attached base packages: >> [1] splines stats graphics grDevices utils datasets methods >> [8] base >> >> other attached packages: >> [1] edgeR_2.6.3 limma_3.12.0 BiocInstaller_1.4.4 >> >> loaded via a namespace (and not attached): >> [1] tools_2.15.0 >> >> >> Thanks, >> Al > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}