Hello, I've found either a bug in estimateDisp() or a mistake in the documentation about when it uses the classic qCML vs. the GLM CR dispersion estimates. The help page for estimateDisp() says "If there is no design matrix, it calculates the quantile conditional likelihood for each gene (tag) and then maximize it. The method is same as in the function estimateCommonDisp and estimateTagwiseDisp". However, when you call estimateDisp() without specifying a design matrix but with a DGEList object with more than one level in samples$group, it internally calculates a design matrix and gives you common/trended/tagwise dispersion estimates. Actually, even when samples$group has only one level, estimateDisp still internally creates a unit vector design matrix and outputs common/trended/tagwise dispersion estimates. So is it a bug in the code or a mistake in the documentation and estimateDisp() can't be used to get the same results as estimateCommonDisp() and estimateTagwiseDisp()? (example below) It was my understanding that for a single factor experiment (2 or more groups), the classic qCML was preferred over the GLM dispersion. However, if you have more than 2 groups and you want to calculate a oneway ANOVA or specialized contrasts, it seems you have to use the GLM approach. Is the GLM approach therefore becoming the best method for all experimental designs, even those with only a single factor? Thanks, Jenny #example - data as indicated in edgeRUsersGuide() section 4.3 > library(edgeR) Loading required package: limma > shell.exec(system.file("doc","edgeRUsersGuide.pdf",package="edgeR")) > > targets <- readTargets() > targets Lane Treatment Label Con1 1 Control Con1 Con2 2 Control Con2 Con3 2 Control Con3 Con4 4 Control Con4 DHT1 5 DHT DHT1 DHT2 6 DHT DHT2 DHT3 8 DHT DHT3 > x <- read.delim("pnas_expression.txt", row.names=1, stringsAsFactors=FALSE) > y <- DGEList(counts=x[,1:7], group=targets$Treatment, genes=data.frame(Length=x[,8])) > colnames(y) <- targets$Label > keep <- rowSums(cpm(y)>1) >= 3 > y <- y[keep,] > y$samples$lib.size <- colSums(y$counts) > y <- calcNormFactors(y) > y$samples group lib.size norm.factors Con1 Control 976847 1.0296636 Con2 Control 1154746 1.0372521 Con3 Control 1439393 1.0362662 Con4 Control 1482652 1.0378383 DHT1 DHT 1820628 0.9537095 DHT2 DHT 1831553 0.9525624 DHT3 DHT 680798 0.9583181 > > y2 <- y > > y <- estimateCommonDisp(y, verbose=TRUE) Disp = 0.02002 , BCV = 0.1415 > y <- estimateTagwiseDisp(y) > y2 <- estimateDisp(y2) Loading required package: splines > > names(y) [1] "counts" "samples" "genes" [4] "common.dispersion" "pseudo.counts" "pseudo.lib.size" [7] "AveLogCPM" "prior.n" "tagwise.dispersion" > names(y2) [1] "counts" "samples" "genes" [4] "common.dispersion" "trended.dispersion" "trend.method" [7] "AveLogCPM" "span" "tagwise.dispersion" [10] "prior.df" "prior.n" > > all.equal(y$tagwise.dispersion, y2$tagwise.dispersion) [1] "Mean relative difference: 0.1156227" > > et <- exactTest(y) > summary(decideTestsDGE(et)) [,1] -1 2096 0 12059 1 2339 > > et2 <- exactTest(y2) > summary(decideTestsDGE(et2)) [,1] -1 2078 0 12084 1 2332 > > > sessionInfo() R version 3.1.1 (2014-07-10) Platform: x86_64-w64-mingw32/x64 (64-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] locfit_1.5-9.1 edgeR_3.6.6 limma_3.20.8 [4] BiocInstaller_1.14.2 loaded via a namespace (and not attached): [1] grid_3.1.1 lattice_0.20-29 tools_3.1.1 Jenny Drnevich, Ph.D. Functional Genomics Bioinformatics Specialist High Performance Biological Computing Program and The Roy J. Carver Biotechnology Center University of Illinois, Urbana-Champaign NOTE NEW OFFICE LOCATION 2112 IGB 1206 W. Gregory Dr. Urbana, IL 61801 USA ph: 217-300-6543 fax: 217-265-5066 e-mail: drnevich@illinois.edu<mailto:drnevich@illinois.edu> [[alternative HTML version deleted]]