satijalab
Featureplot blend issue
Dear Seurat team,
thank you for developing an amazing package!
I had a question regarding the blend feature for FeaturePlot. When I ran this code around a year ago (Seurat V4.3.0), I was able to visualize the expression of two genes with intensities reflected as a color gradient (i.e. Dark red, dim red, dark blue, dim blue). However, I'm running the same code now (Seurat V5.1.0) for the same set of genes but the colors reflecting the expression intensities of the two genes are more discrete (i.e. one tone of red and blue), and do not reflect their expression levels. The script and code remained the same (even the blend threshold at 0.1) and I tried re-running the same code with V4.3.0 but didn't solve the issue - any ideas to how I can re-make the original blend feature plot reflecting expression intensities?
Before (i.e. expression intensities of two genes are reflected as a color gradient)
Current (i.e. no expression intensities are reflected)
Best, taimu
Hi,
I encounter the same exact issue! Have you been able to solve it?
Kind regards,
M.Z
Same observation. Still not solved.
First, I thought it was because I was trying to plot two "scores" calculated by AddModuleScore using blend = TRUE, but this is also the case for just "basic" 2 gene expressions: the gradient is not present for each gene...; instead, it's just discrete values that are plotted. Somehow, the discrete values are either 0, 9, 90, or 99; which correspond to double negatives, single positives (gene_1 or score_1), single positives (gene_2 or score_2), and double positives respectively. I have found those values by setting combine = FALSE in FeaturePlot(), and then putting each patchworked graph next to each other (that way the legends are displayed, and we can see 0, 9, 90, 99).
Minimal example:
library(Seurat)
library(SeuratData)
library(patchwork)
set.seed(123456789)
data("pbmc3k")
pbmc3k = UpdateSeuratObject(pbmc3k)
pbmc3k = NormalizeData(object = pbmc3k)
pbmc3k = FindVariableFeatures(object = pbmc3k)
pbmc3k = ScaleData(object = pbmc3k)
pbmc3k = RunPCA(object = pbmc3k)
pbmc3k = RunUMAP(object = pbmc3k, dims = 1:20)
my_plot = FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE, combine = FALSE)
( my_plot[[1]] | my_plot[[2]] ) / ( my_plot[[3]] | my_plot[[4]] )
Also, modifying blend.threshold does not change anything (except the color scheme/matrix) about the "positiveness" of each gene or score...
Modifying blend.threshold does not work too:
my_plot = FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE)
( my_plot[[1]] | my_plot[[2]] ) / ( my_plot[[3]] | my_plot[[4]] )
my_plot_0.1 = FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE, blend.threshold = 0.1)
( my_plot_0.1[[1]] | my_plot_0.1[[2]] ) / ( my_plot_0.1[[3]] | my_plot_0.1[[4]] )
my_plot_0.9 = FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE, blend.threshold = 0.9)
( my_plot_0.9[[1]] | my_plot_0.9[[2]] ) / ( my_plot_0.9[[3]] | my_plot_0.9[[4]] )
Of course, without using combine = FALSE, the results are the same, so the problem does not come from that either...
Same results without combine = FALSE:
FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"))
FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE)
FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE, blend.threshold = 0.1)
FeaturePlot(object = pbmc3k, features = c("GAPDH", "CD8A"), blend = TRUE, blend.threshold = 0.9)
Therefore, I suggest to implement new parameters to the functions. The first would be something like blend.type = "continuous" by default, but it could be also set to blend.type = "discrete". That way, the users could decide which type of blend graph type they want to see.
Then, if the user selected blend.type = "discrete", it could be good to allow (or force) user to decide additional parameters like blend.threshold.1 for gene_1 or score_1, and also blend.threshold.2 for gene_2 or score_2. Hopefully, the users would have found by themselves the right discrete thresholds to use for each gene or score based on RidgePlot() for example.
Then, if the user selected blend.type = "continuous", I don't really know what would be the best parameter to implement then. I feel I have not enough knowledge about that. But maybe keeping one unique blend.threshold parameter would be the best in that case, I don't know honestly.
I hope you have all the required information in order to solve that issue. If you need anything else, please feel free to ask. If you see any error/mistake in the codes above, please tell me how to correct them please. Thanks in advance for your response, Best regards,
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Monterey 12.7.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.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/Paris
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] patchwork_1.3.0 pbmcMultiome.SeuratData_0.1.4
[3] pbmc3k.SeuratData_3.1.4 ifnb.SeuratData_3.1.0
[5] SeuratData_0.2.2.9001 Seurat_5.1.0
[7] SeuratObject_5.0.2 sp_2.1-4
loaded via a namespace (and not attached):
[1] deldir_2.0-4 pbapply_1.7-2 gridExtra_2.3
[4] rlang_1.1.4 magrittr_2.0.3 RcppAnnoy_0.0.22
[7] matrixStats_1.4.1 ggridges_0.5.6 compiler_4.4.1
[10] spatstat.geom_3.3-4 png_0.1-8 vctrs_0.6.5
[13] reshape2_1.4.4 stringr_1.5.1 crayon_1.5.3
[16] pkgconfig_2.0.3 fastmap_1.2.0 labeling_0.4.3
[19] utf8_1.2.4 promises_1.3.2 purrr_1.0.2
[22] jsonlite_1.8.9 goftest_1.2-3 later_1.4.1
[25] spatstat.utils_3.1-1 irlba_2.3.5.1 parallel_4.4.1
[28] cluster_2.1.6 R6_2.5.1 ica_1.0-3
[31] stringi_1.8.4 RColorBrewer_1.1-3 spatstat.data_3.1-4
[34] reticulate_1.40.0 parallelly_1.39.0 spatstat.univar_3.1-1
[37] lmtest_0.9-40 scattermore_1.2 Rcpp_1.0.13-1
[40] tensor_1.5 future.apply_1.11.3 zoo_1.8-12
[43] sctransform_0.4.1 httpuv_1.6.15 Matrix_1.7-1
[46] splines_4.4.1 igraph_2.1.1 tidyselect_1.2.1
[49] abind_1.4-8 spatstat.random_3.3-2 codetools_0.2-20
[52] miniUI_0.1.1.1 spatstat.explore_3.3-3 listenv_0.9.1
[55] lattice_0.22-6 tibble_3.2.1 plyr_1.8.9
[58] withr_3.0.2 shiny_1.9.1 ROCR_1.0-11
[61] Rtsne_0.17 future_1.34.0 fastDummies_1.7.4
[64] survival_3.7-0 polyclip_1.10-7 fitdistrplus_1.2-1
[67] pillar_1.9.0 KernSmooth_2.23-24 plotly_4.10.4
[70] generics_0.1.3 RcppHNSW_0.6.0 ggplot2_3.5.1
[73] munsell_0.5.1 scales_1.3.0 globals_0.16.3
[76] xtable_1.8-4 glue_1.8.0 lazyeval_0.2.2
[79] tools_4.4.1 data.table_1.16.2 RSpectra_0.16-2
[82] RANN_2.6.2 leiden_0.4.3.1 dotCall64_1.2
[85] cowplot_1.1.3 grid_4.4.1 tidyr_1.3.1
[88] colorspace_2.1-1 nlme_3.1-166 cli_3.6.3
[91] rappdirs_0.3.3 spatstat.sparse_3.1-0 spam_2.11-0
[94] fansi_1.0.6 viridisLite_0.4.2 dplyr_1.1.4
[97] uwot_0.2.2 gtable_0.3.6 digest_0.6.37
[100] progressr_0.15.1 ggrepel_0.9.6 htmlwidgets_1.6.4
[103] farver_2.1.2 htmltools_0.5.8.1 lifecycle_1.0.4
[106] httr_1.4.7 mime_0.12 MASS_7.3-61
I think the original color gradient implementation is probably a better one, I had a hard time to understand what the blend.threshold actually does, as it turns out to be a threshold used to compare with normalized expression values, more misleading is that changing this threshold will not do much to what you see on the plot most of time. The final cell color is not a direct comparison (e.g., "above threshold → pure red/green"). Instead, after transforming gene expressions into normalized values, the differences relative to blend.threshold are fed into a sigmoid function, which produces some weights. These weights will determine what color the cells gets, it is very complex and not intuitive. if you want to ask a follow up question, like how many cells are colored as red, green or yellow, you simply have no idea. So what is the actual use to have this complex plot?
Hello Seurat team, I wanted to follow up and see if there was any update to this issue. I see that this issue is far-reaching given that there was another issue raised with the same problem (https://github.com/satijalab/seurat/issues/8526), and that my lab is also experiencing this issue. Any update or insight would be very helpful. Thank you!
Removing the binarization in FeaturePlot() like below seems to restore the expected behavior. However, I am not sure when this was changed and if this change has any side-effects. ... data.feature[data.feature < min.use] <- min.use data.feature[data.feature > max.use] <- max.use if (brewer.gran != 2) { data.feature <- if (all(data.feature == 0)) { rep_len(x = 0, length.out = length(x = data.feature)) } else { #as.numeric(x = as.factor(x = cut(x = as.numeric(x = data.feature), # breaks = 2))) data.feature } } data[[f]] <- data.feature } ...
Interestingly, the help page of FeaturePlot indicates that RColorBrewer color set names can be used in the cols argument; however, this failed in all cases tested. This is relevant here because the brewer.gran variable in the code snippet above is derived from the number of colors in such a RColorBrewer set (code: brewer.gran <- ifelse(test = length(x = cols) == 1, yes = brewer.pal.info[cols, ]$maxcolors, no = length(x = cols)) ). In previous Seurat versions, the brewer.gran variable was used to set the number of breakpoints in the commented code (breaks = brewer.gran). This was then changed to breaks = 2 in newer versions.
I have the feeling that there is some old code or old concepts lingering around that were not fully cleaned up.
Removing the binarization in FeaturePlot() like below seems to restore the expected behavior. However, I am not sure when this was changed and if this change has any side-effects. ... data.feature[data.feature < min.use] <- min.use data.feature[data.feature > max.use] <- max.use if (brewer.gran != 2) { data.feature <- if (all(data.feature == 0)) { rep_len(x = 0, length.out = length(x = data.feature)) } else { #as.numeric(x = as.factor(x = cut(x = as.numeric(x = data.feature), # breaks = 2))) data.feature } } data[[f]] <- data.feature } ...
Interestingly, the help page of FeaturePlot indicates that RColorBrewer color set names can be used in the cols argument; however, this failed in all cases tested. This is relevant here because the brewer.gran variable in the code snippet above is derived from the number of colors in such a RColorBrewer set (code: brewer.gran <- ifelse(test = length(x = cols) == 1, yes = brewer.pal.info[cols, ]$maxcolors, no = length(x = cols)) ). In previous Seurat versions, the brewer.gran variable was used to set the number of breakpoints in the commented code (breaks = brewer.gran). This was then changed to breaks = 2 in newer versions. I have the feeling that there is some old code or old concepts lingering around that were not fully cleaned up.
@hberger Thanks for the discovery! I am now using the following snippet to circumvent the brewer.gran logic without having to modify the source code[^1]:
trace("FeaturePlot", quote(brewer.gran <- 2), at = 22L)
[^1]: Tested on Seurat 5.3.0 from CRAN, position to attach the tracer (at = 22L) may vary.
