plot_density not working

[ftencaten]

I've used plot_density several times but today I ran into this error and couldn't find more details about it. Is it possible related to the recent Seurat update?

library(Nebulosa) Loading required package: ggplot2 Loading required package: patchwork data <- SeuratObject::pbmc_small plot_density(data,"CD3E") Error in as.vector(x, "character") : cannot coerce type 'environment' to vector of type 'character' sessionInfo() R version 4.3.2 (2023-10-31) Platform: x86_64-apple-darwin20 (64-bit) Running under: macOS Sonoma 14.3.1

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-x86_64/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: America/New_York tzcode source: internal

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

other attached packages: [1] Nebulosa_1.12.0 patchwork_1.2.0 ggplot2_3.5.0

loaded via a namespace (and not attached): [1] RColorBrewer_1.1-3 rstudioapi_0.15.0 jsonlite_1.8.8
[4] magrittr_2.0.3 spatstat.utils_3.0-4 farver_2.1.1
[7] zlibbioc_1.48.0 vctrs_0.6.5 ROCR_1.0-11
[10] spatstat.explore_3.2-6 RCurl_1.98-1.14 S4Arrays_1.2.0
[13] htmltools_0.5.7 SparseArray_1.2.3 pracma_2.4.4
[16] sctransform_0.4.1 parallelly_1.37.0 KernSmooth_2.23-22
[19] htmlwidgets_1.6.4 ica_1.0-3 plyr_1.8.9
[22] plotly_4.10.4 zoo_1.8-12 igraph_2.0.2
[25] mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
[28] Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
[31] GenomeInfoDbData_1.2.11 MatrixGenerics_1.14.0 fitdistrplus_1.1-11
[34] future_1.33.1 shiny_1.8.0 digest_0.6.34
[37] colorspace_2.1-0 S4Vectors_0.40.2 Seurat_5.0.1
[40] tensor_1.5 RSpectra_0.16-1 irlba_2.3.5.1
[43] GenomicRanges_1.54.1 labeling_0.4.3 progressr_0.14.0
[46] fansi_1.0.6 spatstat.sparse_3.0-3 httr_1.4.7
[49] polyclip_1.10-6 abind_1.4-5 compiler_4.3.2
[52] withr_3.0.0 fastDummies_1.7.3 MASS_7.3-60.0.1
[55] DelayedArray_0.28.0 tools_4.3.2 lmtest_0.9-40
[58] httpuv_1.6.14 future.apply_1.11.1 goftest_1.2-3
[61] glue_1.7.0 nlme_3.1-164 promises_1.2.1
[64] grid_4.3.2 Rtsne_0.17 cluster_2.1.6
[67] reshape2_1.4.4 generics_0.1.3 gtable_0.3.4
[70] spatstat.data_3.0-4 tidyr_1.3.1 data.table_1.15.0
[73] sp_2.1-3 utf8_1.2.4 XVector_0.42.0
[76] BiocGenerics_0.48.1 spatstat.geom_3.2-8 RcppAnnoy_0.0.22
[79] ggrepel_0.9.5 RANN_2.6.1 pillar_1.9.0
[82] stringr_1.5.1 spam_2.10-0 RcppHNSW_0.6.0
[85] later_1.3.2 splines_4.3.2 dplyr_1.1.4
[88] lattice_0.22-5 survival_3.5-8 deldir_2.0-2
[91] ks_1.14.2 tidyselect_1.2.0 SingleCellExperiment_1.24.0 [94] miniUI_0.1.1.1 pbapply_1.7-2 gridExtra_2.3
[97] IRanges_2.36.0 SummarizedExperiment_1.32.0 scattermore_1.2
[100] stats4_4.3.2 Biobase_2.62.0 matrixStats_1.2.0
[103] stringi_1.8.3 lazyeval_0.2.2 codetools_0.2-19
[106] tibble_3.2.1 cli_3.6.2 uwot_0.1.16
[109] xtable_1.8-4 reticulate_1.35.0 munsell_0.5.0
[112] Rcpp_1.0.12 GenomeInfoDb_1.38.5 globals_0.16.2
[115] spatstat.random_3.2-2 png_0.1-8 parallel_4.3.2
[118] ellipsis_0.3.2 mclust_6.1 dotCall64_1.1-1
[121] bitops_1.0-7 listenv_0.9.1 mvtnorm_1.2-4
[124] viridisLite_0.4.2 scales_1.3.0 ggridges_0.5.6
[127] crayon_1.5.2 SeuratObject_5.0.1 leiden_0.4.3.1
[130] purrr_1.0.2 rlang_1.1.3 cowplot_1.1.3

[pedriniedoardo]

I encountered precisely the same issue, and I managed to identify the source with a reproducible example. The problem seems to be coming from the use of guide_legend() I the plotting function.

# lirbaries ---------------------------------------------------------------
library(ggplot2)
library(ks)

# function definition -----------------------------------------------------
.extract_feature_data <- function(exp_data, features) {
  # Extract data for input features
  i <- colnames(exp_data) %in% features

  # Test existence of feature in gene expression data
  j <- !features %in% colnames(exp_data)
  if (any(j)) {
    stop(
      "'", paste(features[j], collapse = ", "),
      "' feature(s) not present in meta.data or expression data"
    )
  }
  vars <- exp_data[, i, drop = FALSE]
  vars <- vars[, features, drop = FALSE]
  vars
}

calculate_density <- function(w, x, method, adjust = 1, map = TRUE) {
  if (method == "ks") {
    dens <- kde(x[, c(1, 2)],
                w = w / sum(w) * length(w)
    )
  } else if (method == "wkde") {
    dens <- wkde2d(
      x = x[, 1],
      y = x[, 2],
      w = w / sum(w) * length(w),
      adjust = adjust
    )
  }

  if (map) {
    get_dens(x, dens, method)
  } else {
    dens
  }
}

get_dens <- function(data, dens, method) {
  if (method == "ks") {
    ix <- findInterval(data[, 1], dens$eval.points[[1]])
    iy <- findInterval(data[, 2], dens$eval.points[[2]])
    ii <- cbind(ix, iy)
    z <- dens$estimate[ii]
  } else if (method == "wkde") {
    ix <- findInterval(data[, 1], dens$x)
    iy <- findInterval(data[, 2], dens$y)
    ii <- cbind(ix, iy)
    z <- dens$z[ii]
  }
  z
}

# var definition ----------------------------------------------------------
joint = FALSE
dims = c(1, 2)
method = c("ks", "wkde")
adjust = 1
size = 1
shape = 16
combine = TRUE
pal = "viridis"
raster = F

# load a sample data ------------------------------------------------------
test <- readRDS("../../out/object/pbmc_final.rds")

# wrangling ---------------------------------------------------------------
cell_embeddings <- as.data.frame(SeuratObject::Embeddings(test[["umap"]]))
exp_data <- SeuratObject::FetchData(test, vars = "CD4", slot = "data")
vars <- .extract_feature_data(exp_data, "CD4")

dim_names <- colnames(cell_embeddings)

# library(ks)
z <- calculate_density(vars[, 1], cell_embeddings, method = "ks", adjust)

head(data.frame(cell_embeddings, feature = z))

# fails
ggplot(data.frame(cell_embeddings, feature = z)) +
  aes_string(dim_names[1], dim_names[2], color = "feature") +
  geom_point(shape = shape, size = size) +
  xlab(gsub("_", " ", dim_names[1])) +
  ylab(gsub("_", " ", dim_names[2])) +
  ggtitle("CD4") +
  labs(color = guide_legend("Density")) +
  theme(
    text = element_text(size = 14),
    panel.background = element_blank(),
    axis.text.x = element_text(color = "black"),
    axis.text.y = element_text(color = "black"),
    axis.line = element_line(size = 0.25),
    strip.background = element_rect(color = "black", fill = "#ffe5cc")
  )

ggplot(data.frame(cell_embeddings, feature = z)) +
  aes_string(dim_names[1], dim_names[2], color = "feature") +
  geom_point(shape = shape, size = size) +
  xlab(gsub("_", " ", dim_names[1])) +
  ylab(gsub("_", " ", dim_names[2])) +
  ggtitle("CD4") +
  labs(color = "Density") +
  theme(
    text = element_text(size = 14),
    panel.background = element_blank(),
    axis.text.x = element_text(color = "black"),
    axis.text.y = element_text(color = "black"),
    axis.line = element_line(size = 0.25),
    strip.background = element_rect(color = "black", fill = "#ffe5cc")
  )

TL;DR changing the implementation in plotting.R of the function plot_density_ at line 32 from labs(color = guide_legend(legend_title)) + to labs(color = legend_title) + seems to be solving the issue

library(Nebulosa)
test <- readRDS("../../out/object/pbmc_final.rds")
plot_density(test,"CD4")

This is the same object I have used in case you want to reproduce it.

This is the renv.lock file for my environment.

[beigelk]

I forked the repo and made the change to plotting.R as @pedriniedoardo suggested (labs(color = guide_legend(legend_title)) + to labs(color = legend_title) +), and plot_density is now working as expected for me.

[samuel-marsh]

Thanks @beigelk @pedriniedoardo,

I was just gonna post here as well. The issue is not Seurat but ggplot2 v3.5.0 update which just released on CRAN.

🤞for quick PR and push to Bioconductor!

[ftencaten]

Thanks @pedriniedoardo. Changing labs(color = guide_legend(legend_title)) + to labs(color = legend_title) + in plotting.R did the trick.

[SSandraL]

@pedriniedoardo code works, thank you. How would I change it so I can plot multiple features?

[beigelk]

@SSandraL The issue has been fixed in the GitHub version of Nebulosa. Assuming you have the Bioconductor install of Nebulosa: if you uninstall Nebulosa and re-install the current GitHub version, the plot_density function should work correctly.

# Uninstall current Nebulosa
remove.packages("Nebulosa")

# Install Nebulosa from GitHub
library("devtools")
devtools::install_github("powellgenomicslab/Nebulosa")
library("Nebulosa")

# Load test data
pbmc = readRDS("~/test_data/pbmc_processed.Rds")

# Plot density for two genes
plot_density(pbmc, features = c("LYZ", "CD14"), joint = TRUE)

[apal6]

Hi,

I have used this package before and it worked but it is no more working. I am not sure why as there is literally no error but I am also not able to get the plot. What do you suspect?

Code: plot_density(SeuOs, c("Il21r", "Il2rb"), joint = TRUE)

Thank you

[samuel-marsh]

Hi,

I think the issue is the attempt at merging dimensionality reductions. It seems likely your object does not have unifying dimensionality reduction to plot.

Best, Sam

[apal6]

Hi @samuel-marsh , Do you think join layers would solve the issue?

Best, Aastha

[samuel-marsh]

@apal6 The issue is not layers being separate but it seems like there is not a unifying dimensionality reduction across all cells in the object. That should resolve the issue I would think.

[apal6]

What do you mean by unifying dimensionality reduction?

[samuel-marsh]

The snapshot you posted has messages "merging tsne...". This would seem to imply that each of the objects being merged has it's own tsne. When the objects are merged there is not a tsne that represents the merged object and therefore it cannot be plotted. You need to run an analysis pipeline on merged object first.

[apal6]

Oh got it. Thank you. So, running the dimentionality reduction on the merged object would solve this issue hopefully!?

[samuel-marsh]

I'd hope so. and if it doesn't I'd suggest posting new issue for devs to help resolve the issue.

[apal6]

No, it didn't solve the issue.

I did the whole analysis and yet the same thing.

An object of class Seurat 31061 features across 51442 samples within 2 assays Active assay: RNA (31053 features, 2000 variable features) 3 layers present: counts, data, scale.data 1 other assay present: HTO 3 dimensional reductions calculated: pca, tsne, umap I will post a new issue. Thank you!