FeaturePlot_scCustom labels are incorrect when calling metadata and gene(feature) at the same time with split.by argument

[as-keller]

First of all, thanks for the great package! I have observed a bug that did't always appear, so I wan't sure what the issue was at first. But I have a reproducible example here now. When using FeaturePlot_scCustom along with the split.by argument, the resulting plot is labelled incorrectly, but only when calling a feature from the metadata. When calling a gene and something from the metadata along with a gene, the plots remain in alphabetical order (I think), but the names are in the order listed in the code.

library(SeuratData)
library(Seurat)

InstallData("pbmc3k")
data('pbmc3k')
pbmc3k <- UpdateSeuratObject(pbmc3k)
pbmc3k <- SetIdent(pbmc3k, value = pbmc3k@meta.data$seurat_annotations)
pbmc3k <- subset(pbmc3k, idents = c('Naive CD4 T', 'CD14+ Mono')) 

pbmc3k <- NormalizeData(pbmc3k)
pbmc3k <- FindVariableFeatures(pbmc3k, selection.method = "vst")
all.genes <- rownames(pbmc3k)
pbmc3k <- ScaleData(pbmc3k, features = all.genes)
pbmc3k <- RunPCA(pbmc3k, features = VariableFeatures(pbmc3k))
pbmc3k <- JackStraw(pbmc3k, num.replicate = 100)
pbmc3k <- ScoreJackStraw(pbmc3k, dims = 1:20)
pbmc3k <- FindNeighbors(pbmc3k, dims = 1:20)
pbmc3k <- FindClusters(pbmc3k, resolution = 0.5) 
pbmc3k <- RunUMAP(pbmc3k, dims = 1:20)
pbmc3k <- RunTSNE(pbmc3k, dims = 1:20,  check_duplicates = FALSE)

library(presto)
pbmc3k <- SetIdent(pbmc3k, value = pbmc3k@meta.data$seurat_annotations)
markers <- FindMarkers(pbmc3k, ident.1 = 'Naive CD4 T')
features <- row.names(markers)[1:30]
pbmc3k <- AddModuleScore(pbmc3k, features = features, name = "mod")

library(scCustomize)

FeaturePlot_scCustom(pbmc3k, c('CD74', 'mod1'))
FeaturePlot_scCustom(pbmc3k, c('CD74', 'mod1'), split.by = 'seurat_annotations')
FeaturePlot_scCustom(pbmc3k, c('mod1', 'CD74'), split.by = 'seurat_annotations')

FeaturePlot_scCustom(pbmc3k, c('CD74', 'nCount_RNA'), split.by = 'seurat_annotations')
FeaturePlot_scCustom(pbmc3k, c('nCount_RNA', 'CD74'), split.by = 'seurat_annotations')

FeaturePlot_scCustom(pbmc3k, c('CD74', 'TYROBP'), split.by = 'seurat_annotations')
FeaturePlot_scCustom(pbmc3k, c('TYROBP', 'CD74'), split.by = 'seurat_annotations')

sessionInfo() output

R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.9

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

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

other attached packages:
[1] scCustomize_2.0.1       Seurat_5.0.1            SeuratObject_5.0.1      sp_2.1-2                pbmc3k.SeuratData_3.1.4
[6] panc8.SeuratData_3.0.2  ifnb.SeuratData_3.0.0   SeuratData_0.2.2.9001  

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3     shape_1.4.6            rstudioapi_0.15.0      jsonlite_1.8.7         magrittr_2.0.3        
  [6] spatstat.utils_3.0-4   ggbeeswarm_0.7.2       farver_2.1.1           GlobalOptions_0.1.2    vctrs_0.6.5           
 [11] ROCR_1.0-11            spatstat.explore_3.2-5 paletteer_1.5.0        janitor_2.2.0          forcats_1.0.0         
 [16] htmltools_0.5.7        sctransform_0.4.1      parallelly_1.36.0      KernSmooth_2.23-22     htmlwidgets_1.6.3     
 [21] ica_1.0-3              plyr_1.8.9             lubridate_1.9.3        plotly_4.10.3          zoo_1.8-12            
 [26] igraph_1.5.1           mime_0.12              lifecycle_1.0.4        pkgconfig_2.0.3        Matrix_1.6-4          
 [31] R6_2.5.1               fastmap_1.1.1          snakecase_0.11.1       fitdistrplus_1.1-11    future_1.33.0         
 [36] shiny_1.8.0            digest_0.6.33          colorspace_2.1-0       rematch2_2.1.2         patchwork_1.1.3       
 [41] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1          labeling_0.4.3         progressr_0.14.0      
 [46] timechange_0.2.0       fansi_1.0.6            spatstat.sparse_3.0-3  httr_1.4.7             polyclip_1.10-6       
 [51] abind_1.4-5            compiler_4.3.2         withr_2.5.2            fastDummies_1.7.3      MASS_7.3-60           
 [56] rappdirs_0.3.3         tools_4.3.2            vipor_0.4.5            lmtest_0.9-40          beeswarm_0.4.0        
 [61] httpuv_1.6.12          future.apply_1.11.0    goftest_1.2-3          glue_1.6.2             nlme_3.1-164          
 [66] promises_1.2.1         grid_4.3.2             Rtsne_0.16             cluster_2.1.6          reshape2_1.4.4        
 [71] generics_0.1.3         gtable_0.3.4           spatstat.data_3.0-3    tidyr_1.3.0            data.table_1.14.8     
 [76] utf8_1.2.4             spatstat.geom_3.2-7    RcppAnnoy_0.0.21       ggrepel_0.9.4          RANN_2.6.1            
 [81] pillar_1.9.0           stringr_1.5.1          ggprism_1.0.4          spam_2.10-0            RcppHNSW_0.5.0        
 [86] limma_3.56.2           later_1.3.1            circlize_0.4.15        splines_4.3.2          dplyr_1.1.4           
 [91] lattice_0.22-5         survival_3.5-7         deldir_2.0-2           tidyselect_1.2.0       miniUI_0.1.1.1        
 [96] pbapply_1.7-2          gridExtra_2.3          scattermore_1.2        matrixStats_1.1.0      stringi_1.8.2         
[101] lazyeval_0.2.2         codetools_0.2-19       tibble_3.2.1           cli_3.6.1              uwot_0.1.16           
[106] xtable_1.8-4           reticulate_1.34.0      munsell_0.5.0          Rcpp_1.0.11            globals_0.16.2        
[111] spatstat.random_3.2-2  png_0.1-8              ggrastr_1.0.2          parallel_4.3.2         ellipsis_0.3.2        
[116] ggplot2_3.4.4          presto_1.0.0           dotCall64_1.1-1        listenv_0.9.0          viridisLite_0.4.2     
[121] scales_1.3.0           ggridges_0.5.4         leiden_0.4.3.1         purrr_1.0.2            crayon_1.5.2          
[126] rlang_1.1.2            cowplot_1.1.1  

[samuel-marsh]

Hi @as-keller,

Thanks so much for kind words and thanks very much for the reproducible example! You are right but also found I could reproducible just by general alphabetical issue which allowed me to find bug pretty quickly. This should be fixed now in the develop branch (v2.0.1.9012). If you still have issue after updating package please let me know here and I'll reopen the issue.

Thanks!
Sam

[samuel-marsh]

Also wanted to add just as a note but when adding a module score you need to provide the features as a list otherwise it computes score for each feature in vector and then adds each of those to the meta.data. If you run the code below it will properly create score:

markers <- FindMarkers(pbmc3k, ident.1 = 'Naive CD4 T')
features <- row.names(markers)[1:30]

pbmc3k <- AddModuleScore(pbmc3k, features = list(features), name = "mod")

[as-keller]

Also wanted to add just as a note but when adding a module score you need to provide the features as a list otherwise it computes score for each feature in vector and then adds each of those to the meta.data. If you run the code below it will properly create score:

markers <- FindMarkers(pbmc3k, ident.1 = 'Naive CD4 T')
features <- row.names(markers)[1:30]

pbmc3k <- AddModuleScore(pbmc3k, features = list(features), name = "mod")

Yes, indeed I forgot to call it as a list! . Thanks for the quick fix!