I came across a problem with glimmaMA()
that took me a while to figure out. I've been using it successfully for a while but my latest experiment using almost the same code would just produce blank plots with no error message. The only difference in the code was it was a non-model organism so I pulled annotation information from NCBI's gtf file, which has the symbols in a column named "gene". This turns out to cause some internal problem in glimmaMA()
that tries to put the rownames into a column named "gene" in the bottom table. No error is given, just a blank plot. glimmaMDS()
is not affected. Below is a reproducible example. On a happy note, I first noticed this problem with R 4.0.3 so this forced me to update to 4.1.0/BioC 3.13 more quickly than I usually do!
library(Glimma)
library(limma)
library(edgeR)
dge <- readRDS(system.file("RNAseq123/dge.rds", package = "Glimma"))
#This works:
glimmaMDS(dge)
design <- readRDS(
system.file("RNAseq123/design.rds", package = "Glimma"))
contr.matrix <- readRDS(
system.file("RNAseq123/contr.matrix.rds", package = "Glimma"))
v <- voom(dge, design)
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
#This works:
glimmaMA(efit, dge = dge)
#Change SYMBOL column name to "gene"
names(dge$genes)[2] <- "gene"
names(efit$genes)[2] <- "gene"
#This still works:
glimmaMDS(dge)
#This does not:
glimmaMA(efit, dge = dge)
#Just a blank plot with no errors
sessionInfo( )
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.34.0 limma_3.48.0 Glimma_2.2.0
loaded via a namespace (and not attached):
[1] MatrixGenerics_1.4.0 Biobase_2.52.0 httr_1.4.2
[4] jsonlite_1.7.2 bit64_4.0.5 splines_4.1.0
[7] assertthat_0.2.1 stats4_4.1.0 blob_1.2.1
[10] GenomeInfoDbData_1.2.6 yaml_2.2.1 pillar_1.6.1
[13] RSQLite_2.2.7 lattice_0.20-44 glue_1.4.2
[16] digest_0.6.27 GenomicRanges_1.44.0 RColorBrewer_1.1-2
[19] XVector_0.32.0 colorspace_2.0-1 htmltools_0.5.1.1
[22] Matrix_1.3-3 DESeq2_1.32.0 XML_3.99-0.6
[25] pkgconfig_2.0.3 genefilter_1.74.0 zlibbioc_1.38.0
[28] purrr_0.3.4 xtable_1.8-4 scales_1.1.1
[31] BiocParallel_1.26.0 tibble_3.1.2 annotate_1.70.0
[34] KEGGREST_1.32.0 generics_0.1.0 IRanges_2.26.0
[37] ggplot2_3.3.3 ellipsis_0.3.2 cachem_1.0.5
[40] SummarizedExperiment_1.22.0 BiocGenerics_0.38.0 survival_3.2-11
[43] magrittr_2.0.1 crayon_1.4.1 memoise_2.0.0
[46] evaluate_0.14 fansi_0.4.2 tools_4.1.0
[49] lifecycle_1.0.0 matrixStats_0.58.0 S4Vectors_0.30.0
[52] locfit_1.5-9.4 munsell_0.5.0 DelayedArray_0.18.0
[55] AnnotationDbi_1.54.0 Biostrings_2.60.0 compiler_4.1.0
[58] GenomeInfoDb_1.28.0 tinytex_0.31 rlang_0.4.11
[61] grid_4.1.0 RCurl_1.98-1.3 rstudioapi_0.13
[64] htmlwidgets_1.5.3 bitops_1.0-7 rmarkdown_2.8
[67] gtable_0.3.0 DBI_1.1.1 R6_2.5.0
[70] knitr_1.33 dplyr_1.0.6 fastmap_1.1.0
[73] bit_4.0.4 utf8_1.2.1 parallel_4.1.0
[76] Rcpp_1.0.6 vctrs_0.3.8 geneplotter_1.70.0
[79] png_0.1-7 tidyselect_1.1.1 xfun_0.23
I have recently observed similar issue with glimmaXY() function when I was using a data.frame for gene annotation for 'anno' parameter with NULL value for 'display.columns' parameter. However, I did not have "gene" column in the anno table. The issue got resolved when I explicitly mentioned column names for 'display.columns'.
Thanks for the report and extra thanks to Jenny for diagnosing the problem and providing a reproducible example. I've put in a fix such that the "gene" column is taken from "anno" when available. Hopefully this doesn't cause further issues, currently the only problem is that "gene" from "anno" may contain empty entries, it looks a bit alarming but does not break anything.