bootstrap_enrichment_test throws "Error in exp_mats[[cc]][s, ] <- sort(expD[, cc]) : replacement has length zero"

[Kfarls]

1. Bug description

When running EWCE with my own gene lists and single cell data, I get this error:
"Error in exp_mats[[cc]][s, ] <- sort(expD[, cc]) :
replacement has length zero"
traceback is:

4. | generate_bootstrap_plots_exp_mats(exp_mats = exp_mats, sct_data = sct_data, annotLevel = annotLevel, bootstrap_list = bootstrap_list, reps = reps, combinedGenes = combinedGenes, hits = hits, verbose = verbose)

5. | compute_gene_scores(sct_data = sct_data, annotLevel = annotLevel, bootstrap_list = bootstrap_list, hits = hits, combinedGenes = combinedGenes, verbose = verbose)

6. | get_summed_proportions(hits = hits, sct_data = sct_data, annotLevel = annotLevel, reps = reps, no_cores = no_cores, geneSizeControl = geneSizeControl, controlledCT = controlledCT, control_network = control_network, verbose = verbose)

7. EWCE::bootstrap_enrichment_test(sct_data = ctd, sctSpecies = "human",
   genelistSpecies = "human", hits = hits, reps = reps, annotLevel = annotLevel)

2. Reproducible example

Code

reps <- 100
# Use level 1 annotations)
annotLevel <- 1
hits <- GWAS_psych

psych_results <- EWCE::bootstrap_enrichment_test(sct_data = ctd,
                                                hits = hits, 
                                                sctSpecies_origin = "human",
                                                sctSpecies = "human",
                                                genelistSpecies = "human",
                                                reps = reps,
                                                annotLevel = annotLevel,
                                                geneSizeControl = TRUE,
                                                verbose = TRUE)

Data

(If possible, upload a small sample of your data so that we can reproduce the bug on our end. If that's not possible, please at least include a screenshot of your data and other relevant details.)
I have tried running EWCE with several of my own gene lists, the gene lists work fine when running with the example ctds however i get the above error if I run on my own ctd.
str(GWAS_neurol)
chr [1:45] "ANO3" "APOE" "DYNC1I1" "FHL5" "GALNT2" "PVRL2" "SCN2A" "ADSS" "ASH1L" "ASTN2" "ATG13" "BCL7C" "BTN2A2" "C14orf37" "CAMK1D" ...
ctd generated like so:
counts <- GetAssayData(seurat, slot = "counts")
temp <- list() # initalise empty list
celltype <- unique(metadata$major_clust)
temp <- sapply(1:length(celltype), function(x){
cellname <- celltype[x]
keep <- which(metadata$major_clust == cellname)
counts_subset <- counts[, keep]
keep <- which(tabulate(counts_subset@i + 1) > (0.05 * ncol(counts_subset)))
genes <- rownames(counts_subset)[keep]
updatedlist <- append(temp, list(genes))
})
a <- bitr(geneID = rownames(counts.filtered),
fromType = "ENTREZID",
toType = "SYMBOL",
OrgDb = "org.Hs.eg.db",
drop = FALSE) %>%
mutate(SYMBOL = if_else(is.na(SYMBOL), as.character(ENTREZID),
SYMBOL))
rownames(counts.filtered) <- a$SYMBOL
ctd <- generate_celltype_data(exp = counts.filtered,
annotLevels = annotLevels,
groupName = "pfc_dev",
input_species = "human",
output_species = "human",
savePath = "../results/single_cell/EWCE/")
ctd <- load_rdata("../results/single_cell/EWCE/ctd_pfc_dev.rda")

3. Session info

R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/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: Australia/Perth
tzcode source: internal

attached base packages:
[1] stats4 stats graphics grDevices datasets utils methods base

other attached packages:
[1] devtools_2.4.5 usethis_2.2.2 ewceData_1.8.0 ExperimentHub_2.8.1 AnnotationHub_3.8.0
[6] BiocFileCache_2.8.0 dbplyr_2.3.4 ggdendro_0.1.23 tidyverse_2.0.0 hdWGCNA_0.2.23
[11] EWCE_1.9.3 RNOmni_1.0.1.2 dplyr_1.1.3 tibble_3.2.1 tidyr_1.3.0
[16] readr_2.1.4 purrr_1.0.2 stringr_1.5.0 lubridate_1.9.2 forcats_1.0.0
[21] reshape2_1.4.4 ggplot2_3.4.2 patchwork_1.1.3 gplots_3.1.3 clusterProfiler_4.8.1
[26] rtracklayer_1.60.0 GenomicRanges_1.52.0 GenomeInfoDb_1.36.3 SeuratObject_4.1.3 Seurat_4.3.0.1
[31] org.Hs.eg.db_3.17.0 AnnotationDbi_1.62.2 IRanges_2.34.1 S4Vectors_0.38.1 Biobase_2.60.0
[36] BiocGenerics_0.46.0 TCseq_1.23.0 pacman_0.5.1

[Al-Murphy]

Hey, would it be possible to share your ctd? It would nbot be possible for us to debug otherwise.

Thanks,
Alan.

[Lina0125]

I got the same error, when length of hits are small, this error wont happen, like 15 hits in background with 1350 genes. I tried several times with the same ctd file, only got this error when hits length are big, I'm saying like 500 hits in a background with 1350 genes

[Al-Murphy]

Can you share reproducible code demonstrating the issue? As I said, we can't debug effectively otherwise. Thanks!

[Lina0125]

Hi, Unfortunately I can not share any data, the code I used were

full_results <- EWCE::bootstrap_enrichment_test(sct_data = ctd_mtg_Ensemble2_final, 
                                                bg = background,
                                                sctSpecies = "human",
                                                genelistSpecies = "human",
                                                hits = background[1:500],
                                                reps = 10000,
                                                annotLevel = 1)

I'm not super familiar with this method, I thought it might because some hits not appear in ctd? Or some package problem. It's very strange I can get it run when I use background[1:15] as hits

[Al-Murphy]

EWCE should filter to the genes in the background or at least print a statement to say there was little/no overlap. I can't really help with example data - I would suggest using the vignette ctd and test your ideas of what the problem is on that to get the same error message. If you manage to do that and share the code, I will be able to help.

Thanks,
Alan.

[MathieuBo]

Hi all,
I am facing the same error. I think the error comes from the single cell object I built or the fact that the algorithm uses previously sampled genes.

I have created a new single cell data from an h5ad file I got from cellxgene

Here is the code I used to assemble the object

HB_linarsson <- scKirby::ingest_data(obj = './ccabf939-39df-4ff0-8f7f-2d9f90fb366c.h5ad')

# This dataset uses ENSEMBL IDs as gene name so convert to HGNC
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

gene_conversion <- getBM(
  attributes = c("ensembl_gene_id", "hgnc_symbol", "external_gene_name"),
  filters = "ensembl_gene_id",
  values = row.names(HB_linarsson),
  mart = ensembl
)

# Create named vector for mapping
id_mapping <- setNames(gene_conversion$hgnc_symbol, gene_conversion$ensembl_gene_id)

# Replace rownames
new_gene_names <-  id_mapping[rownames(HB_linarsson)]
rownames(HB_linarsson) <- new_gene_names

# Remove missing gene names
HB_linarsson <- HB_linarsson[!is.na(rownames(HB_linarsson)), ]

# Build EWCE object
annotLevels <- list(level1class=HB_linarsson$supercluster_term,
                    level2class=HB_linarsson$cell_type)

HB_ready <- EWCE::generate_celltype_data(exp = HB_linarsson,
                                          annotLevels = annotLevels, 
                                          groupName = 'HB', )

hb <- EWCE::load_rdata(HB_ready)

try({
  plt <- EWCE::plot_ctd(ctd = hb,
                        level = 1,
                        genes = c("SLC17A7", "GFAP", "GAD2", 'CUX2', 'RORB'),
                        metric = "mean_exp")
})

Now let's use the default hit list - this works well

hits_ewce <- ewceData::example_genelist()

le_results_linarsson <- EWCE::bootstrap_enrichment_test(sct_data = hb,
                                                     hits = hits_ewce, 
                                                     reps = reps,
                                                     sctSpecies = 'human',
                                                     sctSpecies_origin = "human",
                                                     genelistSpecies = "human",
                                                     annotLevel = 1,
                                                     geneSizeControl = TRUE, verbose = TRUE )


try({
  plot_list <- EWCE::ewce_plot(total_res = le_results_linarsson$results,
                               mtc_method = "BH")
  print(plot_list$plain)
})

But now if I use my own list or a random sample for the human genome, I run into the error. Here is the code I used:

human.bg <- ewceData::mouse_to_human_homologs()$HGNC.symbol

gene.list.2 <- sample(human.bg,size = 30)

second_results <- EWCE::bootstrap_enrichment_test(sct_data = hb,
                                                  sctSpecies = "human",
                                                  hits = gene.list.2, 
                                                  reps = reps,
                                                  annotLevel = 1,)

Here is EWCE output:

1 core(s) assigned as workers (9 reserved).
Retrieving all genes using: homologene.
Retrieving all organisms available in homologene.
Mapping species name: human
Common name mapping found for human
1 organism identified from search: 9606
Gene table with 19,129 rows retrieved.
Returning all 19,129 genes from human.
Returning 19,129 unique genes from entire human genome.
Standardising CellTypeDataset
Checking gene list inputs.
Running with gene size control.
Warning: sctSpecies_origin not provided. Setting to 'mouse' by default.
Retrieving all genes using: gprofiler
Retrieving all organisms available in gprofiler.
Using stored gprofiler_orgs.
Mapping species name: human
Common name mapping found for human
1 organism identified from search: hsapiens
Gene table with 63,049 rows retrieved.
Returning all 63,049 genes from human.
63,049 human Ensembl IDs and 41,001 human Gene Symbols imported.
see ?ewceData and browseVignettes('ewceData') for documentation
loading from cache
Controlled bootstrapping network generated.
27 hit gene(s) remain after filtering.
Computing gene scores.
Using previously sampled genes.
Error in exp_mats[[cc]][s, ] <- sort(expD[, cc]) :
number of items to replace is not a multiple of replacement length

If I ran the code above 10 times, maybe 1 or 2 times I will get a result.
It seems to get worse if the hit list gets bigger.

Thanks for your help!!

Session Info

R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1

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: Europe/London
tzcode source: internal

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] ewceData_1.16.0 ExperimentHub_2.16.0 AnnotationHub_3.16.0 BiocFileCache_2.16.0 dbplyr_2.5.0 BiocGenerics_0.54.0 generics_0.1.4
[8] biomaRt_2.64.0

loaded via a namespace (and not attached):
[1] bitops_1.0-9 fs_1.6.6 matrixStats_1.5.0 spatstat.sparse_3.1-0 httr_1.4.7
[6] RColorBrewer_1.1-3 tools_4.5.0 sctransform_0.4.2 backports_1.5.0 R6_2.6.1
[11] DT_0.33 lazyeval_0.2.2 uwot_0.2.3 anndata_0.8.0 withr_3.0.2
[16] sp_2.2-0 prettyunits_1.2.0 gridExtra_2.3 progressr_0.15.1 cli_3.6.5
[21] Biobase_2.68.0 scKirby_0.1.4 spatstat.explore_3.4-3 fastDummies_1.7.5 labeling_0.4.3
[26] Seurat_5.3.0 spatstat.data_3.1-6 readr_2.1.5 ggridges_0.5.6 pbapply_1.7-2
[31] yulab.utils_0.2.0 R.utils_2.13.0 parallelly_1.45.0 orthogene_1.14.0 limma_3.64.1
[36] rstudioapi_0.17.1 RSQLite_2.4.1 gridGraphics_0.5-1 echoconda_0.99.10 ica_1.0-3
[41] spatstat.random_3.4-1 car_3.1-3 dplyr_1.1.4 zip_2.3.3 homologene_1.4.68.19.3.27
[46] Matrix_1.7-3 S4Vectors_0.46.0 abind_1.4-8 R.methodsS3_1.8.2 lifecycle_1.0.4
[51] yaml_2.3.10 carData_3.0-5 SummarizedExperiment_1.38.1 SparseArray_1.8.0 Rtsne_0.17
[56] grid_4.5.0 blob_1.2.4 promises_1.3.3 crayon_1.5.3 dir.expiry_1.16.0
[61] miniUI_0.1.2 lattice_0.22-7 cowplot_1.1.3 KEGGREST_1.48.0 pillar_1.10.2
[66] GenomicRanges_1.60.0 future.apply_1.20.0 codetools_0.2-20 glue_1.8.0 ggfun_0.1.8
[71] spatstat.univar_3.1-3 data.table_1.17.6 vctrs_0.6.5 png_0.1-8 treeio_1.32.0
[76] spam_2.11-1 gtable_0.3.6 assertthat_0.2.1 cachem_1.1.0 openxlsx_4.2.8
[81] S4Arrays_1.8.1 mime_0.13 survival_3.8-3 SingleCellExperiment_1.30.1 RNOmni_1.0.1.2
[86] statmod_1.5.0 fitdistrplus_1.2-2 ROCR_1.0-11 nlme_3.1-168 ggtree_3.16.0
[91] EWCE_1.16.0 bit64_4.6.0-1 progress_1.2.3 filelock_1.0.3 RcppAnnoy_0.0.22
[96] GenomeInfoDb_1.44.0 irlba_2.3.5.1 KernSmooth_2.23-26 splitstackshape_1.4.8 DBI_1.2.3
[101] tidyselect_1.2.1 bit_4.6.0 compiler_4.5.0 curl_6.3.0 httr2_1.1.2
[106] basilisk.utils_1.20.0 xml2_1.3.8 ggdendro_0.2.0 DelayedArray_0.34.1 plotly_4.10.4
[111] scales_1.4.0 lmtest_0.9-40 rappdirs_0.3.3 stringr_1.5.1 digest_0.6.37
[116] goftest_1.2-3 piggyback_0.1.5 spatstat.utils_3.1-4 basilisk_1.20.0 XVector_0.48.0
[121] htmltools_0.5.8.1 pkgconfig_2.0.3 MatrixGenerics_1.20.0 echodata_0.99.17 fastmap_1.2.0
[126] rlang_1.1.6 htmlwidgets_1.6.4 UCSC.utils_1.4.0 shiny_1.10.0 farver_2.1.2
[131] zoo_1.8-14 jsonlite_2.0.0 BiocParallel_1.42.1 R.oo_1.27.1 RCurl_1.98-1.17
[136] magrittr_2.0.3 Formula_1.2-5 GenomeInfoDbData_1.2.14 ggplotify_0.1.2 dotCall64_1.2
[141] patchwork_1.3.0 Rcpp_1.0.14 ape_5.8-1 babelgene_22.9 reticulate_1.42.0
[146] stringi_1.8.7 MASS_7.3-65 plyr_1.8.9 parallel_4.5.0 listenv_0.9.1
[151] ggrepel_0.9.6 deldir_2.0-4 Biostrings_2.76.0 splines_4.5.0 tensor_1.5.1
[156] hms_1.1.3 igraph_2.1.4 ggpubr_0.6.0 spatstat.geom_3.4-1 ggsignif_0.6.4
[161] RcppHNSW_0.6.0 reshape2_1.4.4 stats4_4.5.0 gprofiler2_0.2.3 BiocVersion_3.21.1
[166] SeuratObject_5.1.0 BiocManager_1.30.26 tzdb_0.5.0 httpuv_1.6.16 grr_0.9.5
[171] RANN_2.6.2 tidyr_1.3.1 purrr_1.0.4 polyclip_1.10-7 future_1.58.0
[176] scattermore_1.2 ggplot2_3.5.2 broom_1.0.8 xtable_1.8-4 RSpectra_0.16-2
[181] tidytree_0.4.6 rstatix_0.7.2 later_1.4.2 viridisLite_0.4.2 tibble_3.3.0
[186] aplot_0.2.6 memoise_2.0.1 AnnotationDbi_1.70.0 IRanges_2.42.0 cluster_2.1.8.1
[191] HGNChelper_0.8.15 globals_0.18.0

[Lina0125]

Hi, Unfortunately I can not share any data, the code I used were

full_results <- EWCE::bootstrap_enrichment_test(sct_data = ctd_mtg_Ensemble2_final, 
                                                bg = background,
                                                sctSpecies = "human",
                                                genelistSpecies = "human",
                                                hits = background[1:500],
                                                reps = 10000,
                                                annotLevel = 1)

I'm not super familiar with this method, I thought it might because some hits not appear in ctd? Or some package problem. It's very strange I can get it run when I use background[1:15] as hits

Not sure if it will work for your case, but in my case, I removed HES1 gene from my background genes and gene set, it works fine after that. My colleague also met the same problem and we targeted the error gene and it works fine after removing this error gene.

[MathieuBo]

Good to know, thank you!!
How did you define the background list in the example above?

[Lina0125]

Good to know, thank you!! How did you define the background list in the example above?

Hi, if you are doing differential expression or abundance analysis, the bg gene list are all genes you used for your previous differential expression or abundance analysis, the gene list is differentially expressed genes (pFDR<0.05 or other cutoff). For methodology you could check their paper.

[NathanSkene]

Lina, was there something unusual about the expression of HES1? Did it have NAs or something strange? What led you to delete that one gene?

[Lina0125]

Lina, was there something unusual about the expression of HES1? Did it have NAs or something strange? What led you to delete that one gene?

Could be, in the ctd file, or in the gene expression data I used for generating ctd, but I did not check though.