OSCA-source / OSCA.multisample

The multi-sample subbook of OSCA
5 stars 7 forks source link

Change of default `legacy` in `edgeR::glmQLFit()` with edgeR v4.0.7 changes results in cluster-abundance.Rmd #16

Open PeteHaitch opened 8 months ago

PeteHaitch commented 8 months ago

From edgeR news file:

Changes in version 4.0.0 (2023-10-25)

  • New statistical methods implemented in glmQLFit() to ensure accurate estimation of the quasi-dispersion for data with small counts. The new method computes adjusted residual deviances with adjusted degrees of freedom to improve the chisquare approximation to the residual deviance. The new methodology includes the new argument 'top.proportion' for glmQLFit() to specify the proportion of highly expressed genes used to estimate the common NB dispersion used in the new method. The output DGEGLM object contains new components 'leverage', 'unit.deviance.adj', 'unit.df.adj', 'deviance.adj', 'df.residual.adj' and 'working.dispersion'. The new method can be turned on 'legacy=FALSE'. By default, glmQLFit() will give the same results as in previous releases of edgeR.

Recently, the default changed from legacy=TRUE to legacy=FALSE in the release branch (BioC 3.18): https://code.bioconductor.org/browse/edgeR/commit/1f0de5e1fa24e436315e13fad517e1bdac502fdd @gksmyth: It's a bit unexpected for the default value to change in the release version, so I'd like to confirm that this change was intended for the release version (BioC 3.18) and not just the devel version (BioC 3.19) of edgeR?

@LTLA, @alanocallaghan: With regards to inst/book/cluster-abundance.Rmd in OSCA.multisample, we can revert to the previous behaviour just by adding legacy=TRUE in the appropriate places which seems the best solution for BioC 3.18 (I'll make a PR). For BioC 3.19, we can either adapt the text to the new results (i.e. use the default legacy=FALSE) or stick with the old results and text (i.e. also make legacy=TRUE).


Example based on inst/book/cluster-abundance.Rmd

Extract example and setup DGEList object

library(rebook)
suppressPackageStartupMessages(library(SingleCellExperiment))
extractCached(
  system.file("book/pijuan-embryo.Rmd", package = "OSCA.multisample"),
  chunk = "dimensionality-reduction",
  objects = "merged")
#> <button class="rebook-collapse">View set-up code (Chapter \@ref(chimeric-mouse-embryo-10x-genomics))</button>
#> <div class="rebook-content">
#> 
#> ```r
#> #--- loading ---#
#> library(MouseGastrulationData)
#> sce.chimera <- WTChimeraData(samples=5:10)
#> sce.chimera
#> 
#> #--- feature-annotation ---#
#> library(scater)
#> rownames(sce.chimera) <- uniquifyFeatureNames(
#>     rowData(sce.chimera)$ENSEMBL, rowData(sce.chimera)$SYMBOL)
#> 
#> #--- quality-control ---#
#> drop <- sce.chimera$celltype.mapped %in% c("stripped", "Doublet")
#> sce.chimera <- sce.chimera[,!drop]
#> 
#> #--- normalization ---#
#> sce.chimera <- logNormCounts(sce.chimera)
#> 
#> #--- variance-modelling ---#
#> library(scran)
#> dec.chimera <- modelGeneVar(sce.chimera, block=sce.chimera$sample)
#> chosen.hvgs <- dec.chimera$bio > 0
#> 
#> #--- merging ---#
#> library(batchelor)
#> set.seed(01001001)
#> merged <- correctExperiments(sce.chimera, 
#>     batch=sce.chimera$sample, 
#>     subset.row=chosen.hvgs,
#>     PARAM=FastMnnParam(
#>         merge.order=list(
#>             list(1,3,5), # WT (3 replicates)
#>             list(2,4,6)  # td-Tomato (3 replicates)
#>         )
#>     )
#> )
#> 
#> #--- clustering ---#
#> g <- buildSNNGraph(merged, use.dimred="corrected")
#> clusters <- igraph::cluster_louvain(g)
#> colLabels(merged) <- factor(clusters$membership)
#> 
#> #--- dimensionality-reduction ---#
#> merged <- runTSNE(merged, dimred="corrected", external_neighbors=TRUE)
#> merged <- runUMAP(merged, dimred="corrected", external_neighbors=TRUE)
#> ```
#> 
#> </div>
abundances <- table(merged$celltype.mapped, merged$sample)
abundances <- unclass(abundances)

Performing the DA analysis

suppressPackageStartupMessages(library(edgeR))
# Attaching some column metadata.
extra.info <- colData(merged)[match(colnames(abundances), merged$sample), ]
y.ab <- DGEList(abundances, samples = extra.info)
keep <- filterByExpr(y.ab, group = y.ab$samples$tomato)
y.ab <- y.ab[keep, ]
design <- model.matrix(~factor(pool) + factor(tomato), y.ab$samples)
y.ab <- estimateDisp(y.ab, design, trend="none")

# Change in behaviour
# edgeR v4.0.9 > v4.0.6
fit.ab <- glmQLFit(y.ab, design, robust = TRUE, abundance.trend = FALSE)
# edgeR v4.0.6 behaviour (legacy = TRUE)
fit.ab.legacy <- glmQLFit(y.ab, design, robust = TRUE, abundance.trend = FALSE, legacy = TRUE)
# Change in behaviour
fit.ab$var.prior
#> [1] 2.995648
fit.ab.legacy$var.prior
#> [1] 1.254041

# Consequential change in results
res <- glmQLFTest(fit.ab, coef=ncol(design))
res.legacy <- glmQLFTest(fit.ab.legacy, coef=ncol(design))
topTags(res, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM           F       PValue
#> ExE ectoderm                   -6.61386847 13.02497 44.00496047 2.682303e-08
#> Mesenchyme                      1.16998749 16.29382 15.11355411 3.104565e-04
#> Allantois                       0.77885715 15.50702  5.68283632 2.113409e-02
#> Erythroid3                     -0.64025988 17.28041  5.21921194 2.680090e-02
#> Cardiomyocytes                  0.78103396 14.86430  4.84112733 3.263443e-02
#> Neural crest                   -0.78224850 14.76462  4.61438627 3.678044e-02
#> Endothelium                     0.75511256 14.28905  3.81255096 5.671491e-02
#> Haematoendothelial progenitors  0.62212276 14.72323  3.05073961 8.709715e-02
#> ExE mesoderm                    0.35404299 15.67835  1.28015130 2.634921e-01
#> Pharyngeal mesoderm             0.35492435 15.72073  1.27978174 2.635601e-01
#> Forebrain/Midbrain/Hindbrain   -0.29051795 16.54919  1.01467831 3.188358e-01
#> Def. endoderm                   0.59254011 12.40076  0.97192278 3.291447e-01
#> Surface ectoderm                0.30115212 15.97647  0.95756471 3.327076e-01
#> Erythroid2                     -0.25947602 15.91448  0.67463758 4.155004e-01
#> Caudal Mesoderm                -0.46915460 12.08574  0.48677204 4.887368e-01
#> Paraxial mesoderm              -0.18710612 15.77605  0.36209038 5.501793e-01
#> Intermediate mesoderm          -0.21718837 14.31334  0.32958277 5.685855e-01
#> Somitic mesoderm               -0.19641966 14.26757  0.25071372 6.188613e-01
#> Erythroid1                      0.16020607 14.62373  0.19124084 6.638472e-01
#> Gut                             0.10147578 15.18933  0.09319683 7.614717e-01
#> Spinal cord                    -0.09975658 15.29656  0.09233149 7.625471e-01
#> NMP                            -0.09935100 15.14231  0.08772911 7.683620e-01
#> Blood progenitors 2            -0.09595636 13.57505  0.04718809 8.289512e-01
#> Rostral neurectoderm            0.03196157 13.37318  0.00457413 9.463593e-01
#>                                         FDR
#> ExE ectoderm                   6.437527e-07
#> Mesenchyme                     3.725478e-03
#> Allantois                      1.471217e-01
#> Erythroid3                     1.471217e-01
#> Cardiomyocytes                 1.471217e-01
#> Neural crest                   1.471217e-01
#> Endothelium                    1.944511e-01
#> Haematoendothelial progenitors 2.612914e-01
#> ExE mesoderm                   6.142294e-01
#> Pharyngeal mesoderm            6.142294e-01
#> Forebrain/Midbrain/Hindbrain   6.142294e-01
#> Def. endoderm                  6.142294e-01
#> Surface ectoderm               6.142294e-01
#> Erythroid2                     7.122864e-01
#> Caudal Mesoderm                7.819788e-01
#> Paraxial mesoderm              8.027090e-01
#> Intermediate mesoderm          8.027090e-01
#> Somitic mesoderm               8.251484e-01
#> Erythroid1                     8.382130e-01
#> Gut                            8.382130e-01
#> Spinal cord                    8.382130e-01
#> NMP                            8.382130e-01
#> Blood progenitors 2            8.649926e-01
#> Rostral neurectoderm           9.463593e-01
topTags(res.legacy, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM            F       PValue
#> ExE ectoderm                   -6.56633394 13.02497 66.266873794 1.352022e-10
#> Mesenchyme                      1.16515001 16.29382 11.290522676 1.534506e-03
#> Allantois                       0.83452938 15.50702  5.312163453 2.554581e-02
#> Cardiomyocytes                  0.84840077 14.86430  5.203962678 2.701306e-02
#> Neural crest                   -0.77063783 14.76462  4.106177390 4.830292e-02
#> Endothelium                     0.75186077 14.28905  3.911510756 5.371320e-02
#> Erythroid3                     -0.64308342 17.28041  3.603665884 6.367366e-02
#> Haematoendothelial progenitors  0.65811354 14.72323  3.123945760 8.350641e-02
#> ExE mesoderm                    0.38053094 15.67835  1.180626551 2.826550e-01
#> Pharyngeal mesoderm             0.37931600 15.72073  1.169279448 2.849523e-01
#> Def. endoderm                   0.52621397 12.40076  1.160204113 2.868070e-01
#> Surface ectoderm                0.29698249 15.97647  0.724196419 3.989963e-01
#> Forebrain/Midbrain/Hindbrain   -0.28585716 16.54919  0.712172283 4.029099e-01
#> Caudal Mesoderm                -0.39962542 12.08574  0.576435242 4.514256e-01
#> Erythroid2                     -0.24833470 15.91448  0.498398155 4.836174e-01
#> Paraxial mesoderm              -0.15433244 15.77605  0.193158498 6.622731e-01
#> Erythroid1                      0.16446277 14.62373  0.193071427 6.623444e-01
#> Somitic mesoderm               -0.16330322 14.26757  0.179493241 6.737016e-01
#> Intermediate mesoderm          -0.15822423 14.31334  0.173867968 6.785559e-01
#> Spinal cord                    -0.07592222 15.29656  0.045011021 8.328821e-01
#> NMP                            -0.07425803 15.14231  0.042897724 8.367941e-01
#> Blood progenitors 2             0.03781271 13.57505  0.008261248 9.279571e-01
#> Rostral neurectoderm            0.02531259 13.37318  0.003465095 9.533040e-01
#> Gut                             0.01618639 15.18933  0.002007604 9.644476e-01
#>                                         FDR
#> ExE ectoderm                   3.244852e-09
#> Mesenchyme                     1.841407e-02
#> Allantois                      1.620784e-01
#> Cardiomyocytes                 1.620784e-01
#> Neural crest                   2.148528e-01
#> Endothelium                    2.148528e-01
#> Erythroid3                     2.183097e-01
#> Haematoendothelial progenitors 2.505192e-01
#> ExE mesoderm                   6.257608e-01
#> Pharyngeal mesoderm            6.257608e-01
#> Def. endoderm                  6.257608e-01
#> Surface ectoderm               7.438336e-01
#> Forebrain/Midbrain/Hindbrain   7.438336e-01
#> Caudal Mesoderm                7.737878e-01
#> Erythroid2                     7.737878e-01
#> Paraxial mesoderm              8.571232e-01
#> Erythroid1                     8.571232e-01
#> Somitic mesoderm               8.571232e-01
#> Intermediate mesoderm          8.571232e-01
#> Spinal cord                    9.563362e-01
#> NMP                            9.563362e-01
#> Blood progenitors 2            9.644476e-01
#> Rostral neurectoderm           9.644476e-01
#> Gut                            9.644476e-01

Assuming most labels do not change

y.ab2 <- calcNormFactors(y.ab)
y.ab2 <- estimateDisp(y.ab2, design, trend="none")

# Change in behaviour
# edgeR v4.0.9 > v4.0.6
fit.ab2 <- glmQLFit(y.ab2, design, robust = TRUE, abundance.trend = FALSE)
# edgeR v4.0.6 behaviour (legacy = TRUE)
fit.ab2.legacy <- glmQLFit(y.ab2, design, robust = TRUE, abundance.trend = FALSE, legacy = TRUE)
# Change in behaviour
fit.ab2$var.prior
#> [1] 3.273901
fit.ab2.legacy$var.prior
#> [1] 1.285965

# Consequential change in results
res2 <- glmQLFTest(fit.ab2, coef=ncol(design))
res2.legacy <- glmQLFTest(fit.ab2.legacy, coef=ncol(design))
topTags(res2, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM           F       PValue
#> ExE ectoderm                   -6.99313525 13.17310 48.64910998 8.034546e-09
#> Mesenchyme                      0.94765279 16.26914  9.51490569 3.376641e-03
#> Erythroid3                     -0.85965408 17.34945  8.86993012 4.535282e-03
#> Neural crest                   -1.04110638 14.78020  7.80658441 7.458222e-03
#> Forebrain/Midbrain/Hindbrain   -0.49530467 16.54751  2.83808594 9.854616e-02
#> Erythroid2                     -0.54332440 16.00028  2.80434875 1.005134e-01
#> Cardiomyocytes                  0.57126263 14.83561  2.46147378 1.232387e-01
#> Allantois                       0.52024123 15.51064  2.42563734 1.259343e-01
#> Endothelium                     0.56026575 14.26928  1.97707695 1.661395e-01
#> Paraxial mesoderm              -0.39555970 15.76844  1.55559925 2.183621e-01
#> Intermediate mesoderm          -0.45555658 14.32710  1.37693876 2.464132e-01
#> Haematoendothelial progenitors  0.36547770 14.73329  1.00292439 3.216259e-01
#> Somitic mesoderm               -0.37073003 14.26438  0.84330677 3.630444e-01
#> Caudal Mesoderm                -0.61992557 12.08076  0.79280995 3.776931e-01
#> Spinal cord                    -0.28724638 15.27568  0.73213760 3.964424e-01
#> NMP                            -0.29161531 15.13137  0.71948898 4.005218e-01
#> Blood progenitors 2            -0.31180005 13.57577  0.46929784 4.966054e-01
#> Def. endoderm                   0.35465653 12.39317  0.32212075 5.729811e-01
#> Erythroid1                     -0.13642883 14.64797  0.13229466 7.176618e-01
#> Gut                            -0.11683607 15.17778  0.11763225 7.331144e-01
#> ExE mesoderm                    0.10904610 15.67924  0.11672417 7.341062e-01
#> Pharyngeal mesoderm             0.10553807 15.72341  0.10872880 7.430308e-01
#> Rostral neurectoderm           -0.09598667 13.36717  0.03893561 8.444095e-01
#> Surface ectoderm                0.06068961 15.96412  0.03759877 8.470693e-01
#>                                         FDR
#> ExE ectoderm                   1.928291e-07
#> Mesenchyme                     3.628226e-02
#> Erythroid3                     3.628226e-02
#> Neural crest                   4.474933e-02
#> Forebrain/Midbrain/Hindbrain   3.778028e-01
#> Erythroid2                     3.778028e-01
#> Cardiomyocytes                 3.778028e-01
#> Allantois                      3.778028e-01
#> Endothelium                    4.430386e-01
#> Paraxial mesoderm              5.240690e-01
#> Intermediate mesoderm          5.376288e-01
#> Haematoendothelial progenitors 6.007827e-01
#> Somitic mesoderm               6.007827e-01
#> Caudal Mesoderm                6.007827e-01
#> Spinal cord                    6.007827e-01
#> NMP                            6.007827e-01
#> Blood progenitors 2            7.010899e-01
#> Def. endoderm                  7.639748e-01
#> Erythroid1                     8.105790e-01
#> Gut                            8.105790e-01
#> ExE mesoderm                   8.105790e-01
#> Pharyngeal mesoderm            8.105790e-01
#> Rostral neurectoderm           8.470693e-01
#> Surface ectoderm               8.470693e-01
topTags(res2.legacy, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM           F       PValue
#> ExE ectoderm                   -6.92150525 13.17310 70.36444089 5.737569e-11
#> Mesenchyme                      0.95125144 16.26914  6.78656826 1.219244e-02
#> Neural crest                   -1.00319300 14.78020  6.46396976 1.428911e-02
#> Erythroid3                     -0.85038546 17.34945  5.51725422 2.299442e-02
#> Cardiomyocytes                  0.64004365 14.83561  2.73457457 1.047239e-01
#> Allantois                       0.60540495 15.51064  2.50255680 1.202293e-01
#> Forebrain/Midbrain/Hindbrain   -0.49431263 16.54751  1.92848420 1.713349e-01
#> Endothelium                     0.54819016 14.26928  1.91676555 1.726164e-01
#> Erythroid2                     -0.48181653 16.00028  1.67736283 2.014700e-01
#> Haematoendothelial progenitors  0.42624087 14.73329  1.18478517 2.818191e-01
#> Caudal Mesoderm                -0.55524673 12.08076  1.05155857 3.102887e-01
#> Paraxial mesoderm              -0.36141224 15.76844  0.97289061 3.289064e-01
#> Intermediate mesoderm          -0.38320228 14.32710  0.93668538 3.379832e-01
#> Somitic mesoderm               -0.35788403 14.26438  0.79192048 3.779586e-01
#> Spinal cord                    -0.27467496 15.27568  0.54766167 4.628792e-01
#> NMP                            -0.27553697 15.13137  0.54288458 4.648259e-01
#> Def. endoderm                   0.28482950 12.39317  0.32020658 5.741195e-01
#> Gut                            -0.20898284 15.17778  0.30710432 5.820354e-01
#> ExE mesoderm                    0.15695330 15.67924  0.18283839 6.708575e-01
#> Pharyngeal mesoderm             0.15325751 15.72341  0.17329545 6.790551e-01
#> Blood progenitors 2            -0.17694133 13.57577  0.16845734 6.833128e-01
#> Rostral neurectoderm           -0.12830816 13.36717  0.08347809 7.738824e-01
#> Erythroid1                     -0.08376588 14.64797  0.04611054 8.308847e-01
#> Surface ectoderm                0.07586382 15.96412  0.04355602 8.355649e-01
#>                                         FDR
#> ExE ectoderm                   1.377017e-09
#> Mesenchyme                     1.143129e-01
#> Neural crest                   1.143129e-01
#> Erythroid3                     1.379665e-01
#> Cardiomyocytes                 4.809171e-01
#> Allantois                      4.809171e-01
#> Forebrain/Midbrain/Hindbrain   5.178491e-01
#> Endothelium                    5.178491e-01
#> Erythroid2                     5.372534e-01
#> Haematoendothelial progenitors 6.239689e-01
#> Caudal Mesoderm                6.239689e-01
#> Paraxial mesoderm              6.239689e-01
#> Intermediate mesoderm          6.239689e-01
#> Somitic mesoderm               6.479291e-01
#> Spinal cord                    6.972389e-01
#> NMP                            6.972389e-01
#> Def. endoderm                  7.760471e-01
#> Gut                            7.760471e-01
#> ExE mesoderm                   7.809290e-01
#> Pharyngeal mesoderm            7.809290e-01
#> Blood progenitors 2            7.809290e-01
#> Rostral neurectoderm           8.355649e-01
#> Erythroid1                     8.355649e-01
#> Surface ectoderm               8.355649e-01
# In particular, there are now more cell types that are DA.
summary(decideTests(res2))
#>        factor(tomato)TRUE
#> Down                    3
#> NotSig                 20
#> Up                      1
summary(decideTests(res2.legacy))
#>        factor(tomato)TRUE
#> Down                    1
#> NotSig                 23
#> Up                      0

Removing the offending labels

offenders <- "ExE ectoderm"
y.ab3 <- y.ab[setdiff(rownames(y.ab), offenders),, keep.lib.sizes=FALSE]
y.ab3 <- estimateDisp(y.ab3, design, trend="none")

# Change in behaviour
# edgeR v4.0.9 > v4.0.6
fit.ab3 <- glmQLFit(y.ab3, design, robust = TRUE, abundance.trend = FALSE)
# edgeR v4.0.6 behaviour (legacy = TRUE)
fit.ab3.legacy <- glmQLFit(y.ab3, design, robust = TRUE, abundance.trend = FALSE, legacy = TRUE)
# Change in behaviour
fit.ab3$var.prior
#> [1] 2.662151
fit.ab3.legacy$var.prior
#> [1] 1.186167

# Consequential change in results
res3 <- glmQLFTest(fit.ab3, coef=ncol(design))
res3.legacy <- glmQLFTest(fit.ab3.legacy, coef=ncol(design))
topTags(res, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM           F       PValue
#> ExE ectoderm                   -6.61386847 13.02497 44.00496047 2.682303e-08
#> Mesenchyme                      1.16998749 16.29382 15.11355411 3.104565e-04
#> Allantois                       0.77885715 15.50702  5.68283632 2.113409e-02
#> Erythroid3                     -0.64025988 17.28041  5.21921194 2.680090e-02
#> Cardiomyocytes                  0.78103396 14.86430  4.84112733 3.263443e-02
#> Neural crest                   -0.78224850 14.76462  4.61438627 3.678044e-02
#> Endothelium                     0.75511256 14.28905  3.81255096 5.671491e-02
#> Haematoendothelial progenitors  0.62212276 14.72323  3.05073961 8.709715e-02
#> ExE mesoderm                    0.35404299 15.67835  1.28015130 2.634921e-01
#> Pharyngeal mesoderm             0.35492435 15.72073  1.27978174 2.635601e-01
#> Forebrain/Midbrain/Hindbrain   -0.29051795 16.54919  1.01467831 3.188358e-01
#> Def. endoderm                   0.59254011 12.40076  0.97192278 3.291447e-01
#> Surface ectoderm                0.30115212 15.97647  0.95756471 3.327076e-01
#> Erythroid2                     -0.25947602 15.91448  0.67463758 4.155004e-01
#> Caudal Mesoderm                -0.46915460 12.08574  0.48677204 4.887368e-01
#> Paraxial mesoderm              -0.18710612 15.77605  0.36209038 5.501793e-01
#> Intermediate mesoderm          -0.21718837 14.31334  0.32958277 5.685855e-01
#> Somitic mesoderm               -0.19641966 14.26757  0.25071372 6.188613e-01
#> Erythroid1                      0.16020607 14.62373  0.19124084 6.638472e-01
#> Gut                             0.10147578 15.18933  0.09319683 7.614717e-01
#> Spinal cord                    -0.09975658 15.29656  0.09233149 7.625471e-01
#> NMP                            -0.09935100 15.14231  0.08772911 7.683620e-01
#> Blood progenitors 2            -0.09595636 13.57505  0.04718809 8.289512e-01
#> Rostral neurectoderm            0.03196157 13.37318  0.00457413 9.463593e-01
#>                                         FDR
#> ExE ectoderm                   6.437527e-07
#> Mesenchyme                     3.725478e-03
#> Allantois                      1.471217e-01
#> Erythroid3                     1.471217e-01
#> Cardiomyocytes                 1.471217e-01
#> Neural crest                   1.471217e-01
#> Endothelium                    1.944511e-01
#> Haematoendothelial progenitors 2.612914e-01
#> ExE mesoderm                   6.142294e-01
#> Pharyngeal mesoderm            6.142294e-01
#> Forebrain/Midbrain/Hindbrain   6.142294e-01
#> Def. endoderm                  6.142294e-01
#> Surface ectoderm               6.142294e-01
#> Erythroid2                     7.122864e-01
#> Caudal Mesoderm                7.819788e-01
#> Paraxial mesoderm              8.027090e-01
#> Intermediate mesoderm          8.027090e-01
#> Somitic mesoderm               8.251484e-01
#> Erythroid1                     8.382130e-01
#> Gut                            8.382130e-01
#> Spinal cord                    8.382130e-01
#> NMP                            8.382130e-01
#> Blood progenitors 2            8.649926e-01
#> Rostral neurectoderm           9.463593e-01
topTags(res.legacy, n = Inf)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC   logCPM            F       PValue
#> ExE ectoderm                   -6.56633394 13.02497 66.266873794 1.352022e-10
#> Mesenchyme                      1.16515001 16.29382 11.290522676 1.534506e-03
#> Allantois                       0.83452938 15.50702  5.312163453 2.554581e-02
#> Cardiomyocytes                  0.84840077 14.86430  5.203962678 2.701306e-02
#> Neural crest                   -0.77063783 14.76462  4.106177390 4.830292e-02
#> Endothelium                     0.75186077 14.28905  3.911510756 5.371320e-02
#> Erythroid3                     -0.64308342 17.28041  3.603665884 6.367366e-02
#> Haematoendothelial progenitors  0.65811354 14.72323  3.123945760 8.350641e-02
#> ExE mesoderm                    0.38053094 15.67835  1.180626551 2.826550e-01
#> Pharyngeal mesoderm             0.37931600 15.72073  1.169279448 2.849523e-01
#> Def. endoderm                   0.52621397 12.40076  1.160204113 2.868070e-01
#> Surface ectoderm                0.29698249 15.97647  0.724196419 3.989963e-01
#> Forebrain/Midbrain/Hindbrain   -0.28585716 16.54919  0.712172283 4.029099e-01
#> Caudal Mesoderm                -0.39962542 12.08574  0.576435242 4.514256e-01
#> Erythroid2                     -0.24833470 15.91448  0.498398155 4.836174e-01
#> Paraxial mesoderm              -0.15433244 15.77605  0.193158498 6.622731e-01
#> Erythroid1                      0.16446277 14.62373  0.193071427 6.623444e-01
#> Somitic mesoderm               -0.16330322 14.26757  0.179493241 6.737016e-01
#> Intermediate mesoderm          -0.15822423 14.31334  0.173867968 6.785559e-01
#> Spinal cord                    -0.07592222 15.29656  0.045011021 8.328821e-01
#> NMP                            -0.07425803 15.14231  0.042897724 8.367941e-01
#> Blood progenitors 2             0.03781271 13.57505  0.008261248 9.279571e-01
#> Rostral neurectoderm            0.02531259 13.37318  0.003465095 9.533040e-01
#> Gut                             0.01618639 15.18933  0.002007604 9.644476e-01
#>                                         FDR
#> ExE ectoderm                   3.244852e-09
#> Mesenchyme                     1.841407e-02
#> Allantois                      1.620784e-01
#> Cardiomyocytes                 1.620784e-01
#> Neural crest                   2.148528e-01
#> Endothelium                    2.148528e-01
#> Erythroid3                     2.183097e-01
#> Haematoendothelial progenitors 2.505192e-01
#> ExE mesoderm                   6.257608e-01
#> Pharyngeal mesoderm            6.257608e-01
#> Def. endoderm                  6.257608e-01
#> Surface ectoderm               7.438336e-01
#> Forebrain/Midbrain/Hindbrain   7.438336e-01
#> Caudal Mesoderm                7.737878e-01
#> Erythroid2                     7.737878e-01
#> Paraxial mesoderm              8.571232e-01
#> Erythroid1                     8.571232e-01
#> Somitic mesoderm               8.571232e-01
#> Intermediate mesoderm          8.571232e-01
#> Spinal cord                    9.563362e-01
#> NMP                            9.563362e-01
#> Blood progenitors 2            9.644476e-01
#> Rostral neurectoderm           9.644476e-01
#> Gut                            9.644476e-01

Testing against a log-fold change threshold

# Change in behaviour (inherited from earlier call to glmQLFit())
# edgeR v4.0.9 > v4.0.6
res.lfc <- glmTreat(fit.ab, coef=ncol(design), lfc=1)
# edgeR v4.0.6 behaviour (legacy = TRUE from earlier call to glmQLFit())
res.lfc.legacy <- glmTreat(fit.ab.legacy, coef=ncol(design), lfc=1)

# Consequential change in results
topTags(res.lfc)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC unshrunk.logFC   logCPM       PValue
#> ExE ectoderm                   -6.61386847    -7.04526674 13.02497 3.792118e-08
#> Mesenchyme                      1.16998749     1.17062995 16.29382 1.910310e-02
#> Def. endoderm                   0.59254011     0.59772300 12.40076 5.036371e-01
#> Endothelium                     0.75511256     0.75685752 14.28905 5.365296e-01
#> Caudal Mesoderm                -0.46915460    -0.47427498 12.08574 5.885726e-01
#> Neural crest                   -0.78224850    -0.78344649 14.76462 7.451345e-01
#> Allantois                       0.77885715     0.77953175 15.50702 7.897714e-01
#> Haematoendothelial progenitors  0.62212276     0.62304597 14.72323 8.174108e-01
#> Cardiomyocytes                  0.78103396     0.78217722 14.86430 8.309032e-01
#> Blood progenitors 2            -0.09595636    -0.09641837 13.57505 9.207207e-01
#>                                         FDR
#> ExE ectoderm                   9.101082e-07
#> Mesenchyme                     2.292371e-01
#> Def. endoderm                  9.998683e-01
#> Endothelium                    9.998683e-01
#> Caudal Mesoderm                9.998683e-01
#> Neural crest                   9.998683e-01
#> Allantois                      9.998683e-01
#> Haematoendothelial progenitors 9.998683e-01
#> Cardiomyocytes                 9.998683e-01
#> Blood progenitors 2            9.998683e-01
topTags(res.lfc.legacy)
#> Coefficient:  factor(tomato)TRUE 
#>                                      logFC unshrunk.logFC   logCPM       PValue
#> ExE ectoderm                   -6.56633394    -7.00146617 13.02497 1.106487e-12
#> Mesenchyme                      1.16515001     1.16578344 16.29382 2.158484e-02
#> Def. endoderm                   0.52621397     0.53107528 12.40076 6.512339e-01
#> Endothelium                     0.75186077     0.75363171 14.28905 6.660777e-01
#> Caudal Mesoderm                -0.39962542    -0.40355207 12.08574 7.717705e-01
#> Cardiomyocytes                  0.84840077     0.84982090 14.86430 8.975877e-01
#> Allantois                       0.83452938     0.83536507 15.50702 9.277594e-01
#> Haematoendothelial progenitors  0.65811354     0.65914517 14.72323 9.330722e-01
#> Rostral neurectoderm            0.02531259     0.02538617 13.37318 9.935277e-01
#> Neural crest                   -0.77063783    -0.77191539 14.76462 9.946959e-01
#>                                         FDR
#> ExE ectoderm                   2.655569e-11
#> Mesenchyme                     2.590181e-01
#> Def. endoderm                  9.999981e-01
#> Endothelium                    9.999981e-01
#> Caudal Mesoderm                9.999981e-01
#> Cardiomyocytes                 9.999981e-01
#> Allantois                      9.999981e-01
#> Haematoendothelial progenitors 9.999981e-01
#> Rostral neurectoderm           9.999981e-01
#> Neural crest                   9.999981e-01
Session info ``` r sessionInfo() #> R version 4.3.2 (2023-10-31) #> Platform: x86_64-pc-linux-gnu (64-bit) #> Running under: CentOS Linux 7 (Core) #> #> Matrix products: default #> BLAS: /stornext/System/data/apps/R/R-4.3.2/lib64/R/lib/libRblas.so #> LAPACK: /stornext/System/data/apps/R/R-4.3.2/lib64/R/lib/libRlapack.so; LAPACK version 3.11.0 #> #> locale: #> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C #> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 #> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 #> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C #> [9] LC_ADDRESS=C LC_TELEPHONE=C #> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C #> #> time zone: Australia/Melbourne #> tzcode source: system (glibc) #> #> attached base packages: #> [1] stats4 stats graphics grDevices utils datasets methods #> [8] base #> #> other attached packages: #> [1] edgeR_4.0.9 limma_3.58.1 #> [3] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0 #> [5] Biobase_2.62.0 GenomicRanges_1.54.1 #> [7] GenomeInfoDb_1.38.5 IRanges_2.36.0 #> [9] S4Vectors_0.40.2 BiocGenerics_0.48.1 #> [11] MatrixGenerics_1.14.0 matrixStats_1.2.0 #> [13] rebook_1.12.0 #> #> loaded via a namespace (and not attached): #> [1] SparseArray_1.2.3 bitops_1.0-7 lattice_0.22-5 #> [4] digest_0.6.34 evaluate_0.23 grid_4.3.2 #> [7] fastmap_1.1.1 Matrix_1.6-5 graph_1.80.0 #> [10] BiocManager_1.30.22 XML_3.99-0.16 codetools_0.2-19 #> [13] abind_1.4-5 cli_3.6.2 CodeDepends_0.6.5 #> [16] rlang_1.1.3 crayon_1.5.2 BiocStyle_2.30.0 #> [19] XVector_0.42.0 splines_4.3.2 reprex_2.1.0 #> [22] withr_3.0.0 DelayedArray_0.28.0 yaml_2.3.8 #> [25] S4Arrays_1.2.0 tools_4.3.2 dir.expiry_1.10.0 #> [28] locfit_1.5-9.8 GenomeInfoDbData_1.2.11 filelock_1.0.3 #> [31] lifecycle_1.0.4 zlibbioc_1.48.0 fs_1.6.3 #> [34] glue_1.7.0 Rcpp_1.0.12 statmod_1.5.0 #> [37] xfun_0.41 rstudioapi_0.15.0 knitr_1.45 #> [40] htmltools_0.5.7 rmarkdown_2.25 compiler_4.3.2 #> [43] RCurl_1.98-1.14 ```
gksmyth commented 8 months ago

Ah, thanks for alerting me. Absolutely it was not my intention to change any edgeR defaults in the release version of edgeR. Six days ago I was fixing a bug in the legacy=FALSE glmQLFit.default pipeline and I accidentally overwrote the release version of glmQLFit with the devel version while doing that. Or rather, I forgot that I needed to keep them distinct. My apologies.

I have now repaired my mistake and reset the default back to legacy=TRUE in edgeR 4.0.10. You should be able to get it from git already, otherwise it should be reflected in BiocManager::install in a few days.

gksmyth commented 8 months ago

@LTLA @alanocallaghan If you're able to wait for the edgeR 4.0.10 to percolate through, the issue will be fixed and there is no need to make any changes to OSCA.multisample.

PeteHaitch commented 8 months ago

Thank you very much for the prompt reply, Gordon. I'll test with the git version and we can keep OSCA.multisample as-is in BioC 3.18.

We'll still need to decide what to do in BioC 3.19.

alanocallaghan commented 6 months ago

In https://github.com/OSCA-source/OSCA.multisample/commit/9a50f0d8359b21f586f5debc0dda3bad5cfebde9 I just stuck with legacy=TRUE, but it's probably best to update. Not sure if I'll have the bandwidth to do so before the next release cycle but will try my best