Hi Gordon Smyth, I found a subtle error in voomLmFit()
. Here is a reproducible example and a fix.
Best, Gabriel
library(edgeR)
library(variancePartition)
library(limma)
packageVersion("limma")
# 3.56.2
# example data is from variancePartition
# but I could have used another example
data(varPartDEdata)
dge = DGEList(counts = countMatrix)
dge = calcNormFactors(dge)
form <- ~ Disease
dsgn = model.matrix(form, metadata)
# set prior.weights to be 1
w = rep(1, ncol(dge))
# Since all prior weights are 1, these should give
# identical coefficients
fit1 = voomLmFit( dge, dsgn)
fit2 = voomLmFit( dge, dsgn, prior.weights=w)
# But the coefs are different
coef(fit1)[6,]
#> (Intercept) Disease1
#> 4.030300 1.006082
coef(fit2)[6,]
#> (Intercept) Disease1
#> 4.0520653 0.9855356
# Since the weights from the first gene are used, the coefs for the first gene are identical, but all other are different
It looks like the weights from the first gene are being used here
In edgeR::voomLmFit(
), the function asMatrixWeights()
is used which sets
attr(weights, "arrayweights") <- TRUE
code here
When this is TRUE, this condition is also TRUE, and only the weights from the first gene are used code here
If the code is modified to set attr(weights, "arrayweights")
to FALSE, then this loop
is used to evaluate one gene at a time using the proper row of weights and the code gives the expected results
code here
> sessionInfo( )
R version 4.3.0 (2023-04-21)
Platform: x86_64-apple-darwin22.4.0 (64-bit)
Running under: macOS Ventura 13.5
Matrix products: default
BLAS: /Users/gabrielhoffman/prog/R-4.3.0/lib/libRblas.dylib
LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] variancePartition_1.31.21 BiocParallel_1.34.2
[3] ggplot2_3.4.2 RUnit_0.4.32
[5] edgeR_3.42.4 limma_3.56.2
loaded via a namespace (and not attached):
[1] gtable_0.3.3 caTools_1.18.2 Biobase_2.60.0
[4] lattice_0.21-8 numDeriv_2016.8-1.1 vctrs_0.6.3
[7] tools_4.3.0 Rdpack_2.4 bitops_1.0-7
[10] generics_0.1.3 parallel_4.3.0 pbkrtest_0.5.2
[13] tibble_3.2.1 fansi_1.0.4 pkgconfig_2.0.3
[16] Matrix_1.5-4.1 KernSmooth_2.23-21 lifecycle_1.0.3
[19] stringr_1.5.0 compiler_4.3.0 gplots_3.1.3
[22] munsell_0.5.0 RhpcBLASctl_0.23-42 codetools_0.2-19
[25] lmerTest_3.1-3 pillar_1.9.0 nloptr_2.0.3
[28] tidyr_1.3.0 MASS_7.3-60 aod_1.3.2
[31] iterators_1.0.14 boot_1.3-28.1 nlme_3.1-162
[34] gtools_3.9.4 tidyselect_1.2.0 locfit_1.5-9.7
[37] stringi_1.7.12 fANCOVA_0.6-1 mvtnorm_1.2-2
[40] reshape2_1.4.4 dplyr_1.1.2 purrr_1.0.1
[43] splines_4.3.0 grid_4.3.0 colorspace_2.1-0
[46] cli_3.6.1 magrittr_2.0.3 utf8_1.2.3
[49] broom_1.0.5 corpcor_1.6.10 withr_2.5.0
[52] scales_1.2.1 backports_1.4.1 remaCor_0.0.17
[55] lme4_1.1-33 rbibutils_2.2.13 EnvStats_2.7.0
[58] rlang_1.1.1 Rcpp_1.0.11 glue_1.6.2
[61] BiocGenerics_0.46.0 minqa_1.2.5 plyr_1.8.8
[64] R6_2.5.1