I recently started using the cqn package and have found a small error in the vignette. In code chunk 17, the glm.offset is put into the d.mont
object and then estimateGLMCommonDisp
is run to create a new object d.mont.cqn
:
###################################################
### code chunk number 17: edgeRdisp
###################################################
design <- model.matrix(~ d.mont$sample$group)
d.mont$offset <- cqn.subset$glm.offset
d.mont.cqn <- estimateGLMCommonDisp(d.mont, design = design)
All this is fine and the results are shown in chunk 19:
> ###################################################
> ### code chunk number 19: cqn.Rnw:264-265
> ###################################################
> summary(decideTestsDGE(elrt.cqn))
d.mont$sample$groupgrp2
Down 146
NotSig 22971
Up 435
However, in code chunk 20 the d.mont
object is re-used, but since it still contains the glm.offset, the results are exactly the same:
> ###################################################
> ### code chunk number 20: edgeRstd
> ###################################################
> d.mont.std <- estimateGLMCommonDisp(d.mont, design = design)
> efit.std <- glmFit(d.mont.std, design = design)
> elrt.std <- glmLRT(efit.std, coef = 2)
> summary(decideTestsDGE(elrt.std))
d.mont$sample$groupgrp2
Down 146
NotSig 22971
Up 435
The quickest fix would be to just set the offsets back to NULL:
> d.mont$offset <- NULL
> d.mont.std <- estimateGLMCommonDisp(d.mont, design = design)
> efit.std <- glmFit(d.mont.std, design = design)
> elrt.std <- glmLRT(efit.std, coef = 2)
> summary(decideTestsDGE(elrt.std))
d.mont$sample$groupgrp2
Down 211
NotSig 23086
Up 255
I might also suggest updating to the edgeR-quasi method as it's the new preferred way. Here are the codes I've been using:
> ###################################################
> ### code chunk number 17: edgeRdisp
> ###################################################
> design <- model.matrix(~ d.mont$sample$group)
> d.mont$offset <- cqn.subset$glm.offset
> d.mont.cqn <- estimateDisp(d.mont, design = design)
>
>
> ###################################################
> ### code chunk number 18: edgeRfit
> ###################################################
> efit.cqn <- glmQLFit(d.mont.cqn, design = design)
> eqlf.cqn <- glmQLFTest(efit.cqn, coef = 2)
> topTags(eqlf.cqn, n = 2)
Coefficient: d.mont$sample$groupgrp2
length gccontent logFC logCPM F PValue FDR
ENSG00000253701 359 0.6100279 7.504784 3.608416 55.69780 6.437842e-08 0.001516241
ENSG00000211642 365 0.5835616 -10.252487 6.362877 28.67522 1.329693e-05 0.156584594
>
>
> ###################################################
> ### code chunk number 19: cqn.Rnw:264-265
> ###################################################
> summary(decideTestsDGE(eqlf.cqn))
d.mont$sample$groupgrp2
Down 0
NotSig 23551
Up 1
>
> ###################################################
> ### code chunk number 20: edgeRstd
> ###################################################
> d.mont$offset <- NULL
> d.mont.std <- estimateDisp(d.mont, design = design)
> efit.std <- glmQLFit(d.mont.std, design = design)
> eqlf.std <- glmQLFTest(efit.std, coef = 2)
> summary(decideTestsDGE(eqlf.std))
d.mont$sample$groupgrp2
Down 10
NotSig 23534
Up 8
Although with the edgeR-quasi method, very few genes reach the FDR < 0.05 threshold in this example.
Thanks for the cqn package!
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] splines stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.28.0 limma_3.42.0 scales_1.1.0 cqn_1.32.0
[5] quantreg_5.52 SparseM_1.77 preprocessCore_1.48.0 nor1mix_1.3-0
[9] mclust_5.4.5
loaded via a namespace (and not attached):
[1] Rcpp_1.0.3 locfit_1.5-9.1 lattice_0.20-38 grid_3.6.1 R6_2.4.1
[6] lifecycle_0.1.0 MatrixModels_0.4-1 rlang_0.4.2 farver_2.0.1 Matrix_1.2-17
[11] tools_3.6.1 munsell_0.5.0 compiler_3.6.1 colorspace_1.4-1