[BUG] Error during extract and reformat enrichment results when using `registration_wrapper`

[boyiguo1]

Describe the bug

Thanks for creating and maintaining the awesome spatialLIBD package. Recently, I was trying to conduct spatial registration by running registration_wrapper function on a SpatialExperiment object. I encountered the error message Error in combn(colnames(registration_model)[regis_cols], 2) : n < m (please find the reproducible example below), which I have trouble interpreting the error message and fix it.

The spatialLIBD is at v1.14.1 installed via Bioconductor version '3.18'

Provide a minimally reproducible example (reprex)

# Load packages -----------------------------------------------------------
library(here)
#> here() starts at /private/var/folders/ff/m8yvwsq51sjg43xscd3x1tmr0000gn/T/RtmpWA9JNL/reprex-774f49f557da-stony-moth
library(SpatialExperiment)
#> Loading required package: SingleCellExperiment
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
library(spatialLIBD)
library(SingleR)
# library(tidyverse)

# Load spatialDLPFC data ---------------------------------------------------

spe <- fetch_data(type = "spatialDLPFC_Visium")
#> 2024-01-02 13:33:34.255759 loading file /Users/bguo6/Library/Caches/org.R-project.R/R/BiocFileCache/5449288a26e9_spe_filtered_final_with_clusters_and_deconvolution_results.rds%3Fdl%3D1

set.seed(1)
scr_idx <- sample(unique(spe$sample_id), 15)
tar_idx <- setdiff(unique(spe$sample_id), scr_idx)

src_spe <- spe[,spe$sample_id %in% scr_idx]
tar_spe <- spe[,spe$sample_id %in% tar_idx]

# Spatial Registration ----------------------------------------------------

## Get enrichment statistics from source ---------------------------------
sce_modeling_results <- registration_wrapper(
  sce = src_spe,
  # sce = src_spe,
  var_registration = "BayesSpace_harmony_09",
  var_sample_id = "sample_id",
  gene_ensembl = "gene_id",
  gene_name = "gene_name"
)
#> 2024-01-02 13:34:59.526861 make pseudobulk object
#> 2024-01-02 13:35:07.619735 dropping 0 pseudo-bulked samples that are below 'min_ncells'.
#> 2024-01-02 13:35:07.650823 drop lowly expressed genes
#> 2024-01-02 13:35:07.781492 normalize expression
#> 2024-01-02 13:35:08.545616 create model matrix
#> 2024-01-02 13:35:08.696194 run duplicateCorrelation()
#> 2024-01-02 13:35:26.808509 The estimated correlation is: 0.960849320333061
#> 2024-01-02 13:35:26.811491 computing enrichment statistics
#> 2024-01-02 13:35:28.62727 extract and reformat enrichment results
#> Error in combn(colnames(registration_model)[regis_cols], 2): n < m

^{Created on 2024-01-02 with reprex v2.0.2}

Expected behavior

I would expect that registration_wrapper would run successfully.

R Session Information

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.3.2 (2023-10-31) #> os macOS Sonoma 14.2.1 #> system aarch64, darwin20 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/New_York #> date 2024-01-02 #> pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 4.3.0) #> AnnotationDbi 1.64.1 2023-11-02 [1] Bioconductor #> AnnotationHub 3.10.0 2023-10-26 [1] Bioconductor #> attempt 0.3.1 2020-05-03 [1] CRAN (R 4.3.0) #> beachmat 2.18.0 2023-10-24 [1] Bioconductor #> beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.3.0) #> benchmarkme 1.0.8 2022-06-12 [1] CRAN (R 4.3.0) #> benchmarkmeData 1.0.4 2020-04-23 [1] CRAN (R 4.3.0) #> Biobase * 2.62.0 2023-10-26 [1] Bioconductor #> BiocFileCache 2.10.1 2023-10-26 [1] Bioconductor #> BiocGenerics * 0.48.1 2023-11-02 [1] Bioconductor #> BiocIO 1.12.0 2023-10-26 [1] Bioconductor #> BiocManager 1.30.22 2023-08-08 [1] CRAN (R 4.3.0) #> BiocNeighbors 1.20.1 2023-12-18 [1] Bioconductor 3.18 (R 4.3.0) #> BiocParallel 1.36.0 2023-10-26 [1] Bioconductor #> BiocSingular 1.18.0 2023-10-24 [1] Bioconductor #> BiocVersion 3.18.1 2023-11-18 [1] Bioconductor 3.18 (R 4.3.2) #> Biostrings 2.70.1 2023-10-26 [1] Bioconductor #> bit 4.0.5 2022-11-15 [1] CRAN (R 4.3.0) #> bit64 4.0.5 2020-08-30 [1] CRAN (R 4.3.0) #> bitops 1.0-7 2021-04-24 [1] CRAN (R 4.3.0) #> blob 1.2.4 2023-03-17 [1] CRAN (R 4.3.0) #> bslib 0.6.1 2023-11-28 [1] CRAN (R 4.3.1) #> cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.0) #> cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.1) #> codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.2) #> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0) #> config 0.3.2 2023-08-30 [1] CRAN (R 4.3.0) #> cowplot 1.1.2 2023-12-15 [1] CRAN (R 4.3.1) #> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0) #> curl 5.2.0 2023-12-08 [1] CRAN (R 4.3.1) #> data.table 1.14.10 2023-12-08 [1] CRAN (R 4.3.1) #> DBI 1.2.0 2023-12-21 [1] CRAN (R 4.3.1) #> dbplyr 2.4.0 2023-10-26 [1] CRAN (R 4.3.1) #> DelayedArray 0.28.0 2023-10-24 [1] Bioconductor #> DelayedMatrixStats 1.24.0 2023-10-24 [1] Bioconductor #> digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.0) #> doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.3.0) #> dotCall64 1.1-1 2023-11-28 [1] CRAN (R 4.3.1) #> dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.3.1) #> DT 0.31 2023-12-09 [1] CRAN (R 4.3.1) #> edgeR 4.0.3 2023-12-09 [1] Bioconductor 3.18 (R 4.3.2) #> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0) #> evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.1) #> ExperimentHub 2.10.0 2023-10-26 [1] Bioconductor #> fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.1) #> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0) #> fields 15.2 2023-08-17 [1] CRAN (R 4.3.0) #> filelock 1.0.3 2023-12-11 [1] CRAN (R 4.3.1) #> foreach 1.5.2 2022-02-02 [1] CRAN (R 4.3.0) #> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.0) #> generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0) #> GenomeInfoDb * 1.38.1 2023-11-11 [1] Bioconductor #> GenomeInfoDbData 1.2.11 2023-11-08 [1] Bioconductor #> GenomicAlignments 1.38.0 2023-10-26 [1] Bioconductor #> GenomicRanges * 1.54.1 2023-10-30 [1] Bioconductor #> ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.3.0) #> ggplot2 3.4.4 2023-10-12 [1] CRAN (R 4.3.1) #> ggrepel 0.9.4 2023-10-13 [1] CRAN (R 4.3.1) #> glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0) #> golem 0.4.1 2023-06-05 [1] CRAN (R 4.3.0) #> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0) #> gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.0) #> here * 1.0.1 2020-12-13 [1] CRAN (R 4.3.0) #> htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.1) #> htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.1) #> httpuv 1.6.13 2023-12-06 [1] CRAN (R 4.3.1) #> httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.0) #> interactiveDisplayBase 1.40.0 2023-10-26 [1] Bioconductor #> IRanges * 2.36.0 2023-10-26 [1] Bioconductor #> irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.3.0) #> iterators 1.0.14 2022-02-05 [1] CRAN (R 4.3.0) #> jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.0) #> jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.1) #> KEGGREST 1.42.0 2023-10-26 [1] Bioconductor #> knitr 1.45 2023-10-30 [1] CRAN (R 4.3.1) #> later 1.3.2 2023-12-06 [1] CRAN (R 4.3.1) #> lattice 0.21-9 2023-10-01 [1] CRAN (R 4.3.2) #> lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.0) #> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.0) #> limma 3.58.1 2023-11-02 [1] Bioconductor #> locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.3.0) #> magick 2.8.2 2023-12-20 [1] CRAN (R 4.3.1) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0) #> maps 3.4.2 2023-12-15 [1] CRAN (R 4.3.1) #> Matrix 1.6-1.1 2023-09-18 [1] CRAN (R 4.3.2) #> MatrixGenerics * 1.14.0 2023-10-26 [1] Bioconductor #> matrixStats * 1.2.0 2023-12-11 [1] CRAN (R 4.3.1) #> memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.0) #> mime 0.12 2021-09-28 [1] CRAN (R 4.3.0) #> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0) #> paletteer 1.5.0 2022-10-19 [1] CRAN (R 4.3.0) #> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0) #> plotly 4.10.3 2023-10-21 [1] CRAN (R 4.3.1) #> png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0) #> promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.0) #> purrr 1.0.2 2023-08-10 [1] CRAN (R 4.3.0) #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.3.0) #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0) #> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.3.0) #> R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.3.1) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0) #> rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.0) #> RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.3.0) #> Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.0) #> RCurl 1.98-1.13 2023-11-02 [1] CRAN (R 4.3.1) #> rematch2 2.1.2 2020-05-01 [1] CRAN (R 4.3.0) #> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.3.0) #> restfulr 0.0.15 2022-06-16 [1] CRAN (R 4.3.0) #> rjson 0.2.21 2022-01-09 [1] CRAN (R 4.3.0) #> rlang 1.1.2 2023-11-04 [1] CRAN (R 4.3.0) #> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1) #> rprojroot 2.0.4 2023-11-05 [1] CRAN (R 4.3.0) #> Rsamtools 2.18.0 2023-10-26 [1] Bioconductor #> RSQLite 2.3.4 2023-12-08 [1] CRAN (R 4.3.1) #> rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0) #> rsvd 1.0.5 2021-04-16 [1] CRAN (R 4.3.0) #> rtracklayer 1.62.0 2023-10-26 [1] Bioconductor #> S4Arrays 1.2.0 2023-10-26 [1] Bioconductor #> S4Vectors * 0.40.2 2023-11-25 [1] Bioconductor 3.18 (R 4.3.2) #> sass 0.4.8 2023-12-06 [1] CRAN (R 4.3.1) #> ScaledMatrix 1.10.0 2023-10-24 [1] Bioconductor #> scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.1) #> scater 1.30.1 2023-11-16 [1] Bioconductor #> scuttle 1.12.0 2023-10-24 [1] Bioconductor #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0) #> shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.1) #> shinyWidgets 0.8.0 2023-08-30 [1] CRAN (R 4.3.0) #> SingleCellExperiment * 1.24.0 2023-10-24 [1] Bioconductor #> SingleR * 2.4.0 2023-10-24 [1] Bioconductor #> spam 2.10-0 2023-10-23 [1] CRAN (R 4.3.1) #> SparseArray 1.2.2 2023-11-08 [1] Bioconductor #> sparseMatrixStats 1.14.0 2023-10-26 [1] Bioconductor #> SpatialExperiment * 1.12.0 2023-10-26 [1] Bioconductor #> spatialLIBD * 1.14.1 2023-11-30 [1] Bioconductor 3.18 (R 4.3.0) #> statmod 1.5.0 2023-01-06 [1] CRAN (R 4.3.0) #> stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.1) #> stringr 1.5.1 2023-11-14 [1] CRAN (R 4.3.1) #> styler 1.10.2 2023-08-29 [1] CRAN (R 4.3.0) #> SummarizedExperiment * 1.32.0 2023-10-24 [1] Bioconductor #> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0) #> tidyr 1.3.0 2023-01-24 [1] CRAN (R 4.3.0) #> tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0) #> utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.1) #> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.1) #> vipor 0.4.7 2023-12-18 [1] CRAN (R 4.3.1) #> viridis 0.6.4 2023-07-22 [1] CRAN (R 4.3.0) #> viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.0) #> withr 2.5.2 2023-10-30 [1] CRAN (R 4.3.1) #> xfun 0.41 2023-11-01 [1] CRAN (R 4.3.1) #> XML 3.99-0.16 2023-11-29 [1] CRAN (R 4.3.1) #> xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.0) #> XVector 0.42.0 2023-10-26 [1] Bioconductor #> yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.1) #> zlibbioc 1.48.0 2023-10-26 [1] Bioconductor #> #> [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library #> #> ────────────────────────────────────────────────────────────────────────────── ```

[lahuuki]

Hi Boyi, This error was the result of src_spe$BayesSpace_harmony_09 being an integer, the root of the error is here where the pairwise comparison is set up. However this probably also cause unexpected behaviour with registration_stats_anova and registration_stats_enrichment as well since the model is not categorical but continuous.

This can be solved by replacing the numeric categories with syntactically valid character instead like this: src_spe$BayesSpace_harmony_09 <- paste0("k09d",src_spe$BayesSpace_harmony_09)

Non-syntatic characters names also cause issues like: Error in limma::makeContrasts(contrasts = z, levels = registration_model) : The levels must by syntactically valid names in R, see help(make.names).

We should include a check that var_registration is categorical and syntactically valid.

Thanks for pointing out this Bug!

[boyiguo1]

I think @lahuuki's explanation is perfect. I just want to add another corner case if other's not paying attention with categorical and syntactically valid.

The variable has to be a string or factor whose labels follows the syntactically valid.

The error will remain if someone uses the following line of code

> spe$BayesSpace_harmony_09 <- factor(spe$BayesSpace_harmony_09)
> levels(spe$BayesSpace_harmony_09)
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9"

There's no error with the following code

> spe$BayesSpace_harmony_09 <- factor(
+   spe$BayesSpace_harmony_09,
+   labels = paste0("k09d", levels(spe$BayesSpace_harmony_09))
+ )
> levels(spe$BayesSpace_harmony_09 )
[1] "k09d1" "k09d2" "k09d3" "k09d4" "k09d5" "k09d6" "k09d7"
[8] "k09d8" "k09d9"

[lcolladotor]

@lahuuki can you use Boyi's two reproducible examples to build a unit tests for these functions?

Thanks!