Reduced dimension plotting is one of the essential tools for the analysis of single cell data. However, as the number of cells/nuclei in these these plots increases, the usefulness of these plots decreases. Many cells are plotted on top of each other obscuring information, even when taking advantage of transparency settings. This package provides binning strategies of cells/nuclei into hexagon cells. Plotting summarized information of all cells/nuclei in their respective hexagon cells presents information without obstructions. The package seemlessly works with the two most common object classes for the storage of single cell data; SingleCellExperiment
from the SingleCellExperiment package and Seurat
from the Seurat package.
In order to demonstrate the capabilities of the schex package, I will use the a dataset of Peripheral Blood Mononuclear Cells (PBMC) freely available from 10x Genomics. There are 2,700 single cells that were sequenced on the Illumina NextSeq 500. This data is handly availabe in the TENxPBMCData
package.
tenx_pbmc3k <- TENxPBMCData(dataset = "pbmc3k")
#> see ?TENxPBMCData and browseVignettes('TENxPBMCData') for documentation
#> loading from cache
rownames(tenx_pbmc3k) <- uniquifyFeatureNames(
rowData(tenx_pbmc3k)$ENSEMBL_ID,
rowData(tenx_pbmc3k)$Symbol_TENx
)
In the next few sections, I will perform some simple quality control steps including filtering and normalization. I will then calculate various dimension reductions and cluster the data. These steps do by no means constitute comprehensive handling of the data. For a more detailed guide the reader is referred to the following guides:
I filter cells with high mitochondrial content as well as cells with low library size or feature count.
rowData(tenx_pbmc3k)$Mito <- grepl("^MT-", rownames(tenx_pbmc3k))
colData(tenx_pbmc3k) <- cbind(
colData(tenx_pbmc3k),
perCellQCMetrics(tenx_pbmc3k,
subsets = list(Mt = rowData(tenx_pbmc3k)$Mito)
)
)
rowData(tenx_pbmc3k) <- cbind(
rowData(tenx_pbmc3k),
perFeatureQCMetrics(tenx_pbmc3k)
)
tenx_pbmc3k <- tenx_pbmc3k[, !colData(tenx_pbmc3k)$subsets_Mt_percent > 50]
libsize_drop <- isOutlier(tenx_pbmc3k$total,
nmads = 3, type = "lower", log = TRUE
)
feature_drop <- isOutlier(tenx_pbmc3k$detected,
nmads = 3, type = "lower", log = TRUE
)
tenx_pbmc3k <- tenx_pbmc3k[, !(libsize_drop | feature_drop)]
I filter any genes that have 0 count for all cells.
I normalize the data by using a simple library size normalization.
I use both Principal Components Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) in order to obtain reduced dimension representations of the data. Since there is a random component in the UMAP, I will also set a seed.
At this stage in the workflow we usually would like to plot aspects of our data in one of the reduced dimension representations. Instead of plotting this in an ordinary fashion, I will demonstrate how schex can provide a better way of plotting this.
First, I will calculate the hexagon cell representation for each cell for a specified dimension reduction representation. I decide to use nbins=40
which specifies that I divide my x range into 40 bins. Note that this might be a parameter that you want to play around with depending on the number of cells/ nuclei in your dataset. Generally, for more cells/nuclei, nbins
should be increased.
First I plot how many cells are in each hexagon cell. This should be relatively even, otherwise change the nbins
parameter in the previous calculation.
Next I colour the hexagon cells by some meta information, such as the majority of cells cluster membership and the median total count in each hexagon cell.
While for plotting the cluster membership the outcome is not too different from the classic plot, it is much easier to observe differences in the total count.
For convenience there is also a function that allows the calculation of label positions for factor variables. These can be overlayed with the package ggrepel
.
Finally, I will visualize the gene expression of the POMGNT1 gene in the hexagon cell representation.
gene_id <- "POMGNT1"
plot_hexbin_feature(tenx_pbmc3k,
type = "logcounts", feature = gene_id,
action = "mean", xlab = "UMAP1", ylab = "UMAP2",
title = paste0("Mean of ", gene_id)
)
Again it is much easier to observe differences in gene expression using the hexagon cell representation than the classic representation.
We can overlay the gene expression data with the clusters for convenience.
schex
output as ggplot
objectsThe schex
packages renders ordinary ggplot
objects and thus these can be treated and manipulated using the ggplot
grammar. For example the non-data components of the plots can be changed using the function theme
.
gene_id <- "CD19"
gg <- schex::plot_hexbin_feature(tenx_pbmc3k,
type = "logcounts", feature = gene_id,
action = "mean", xlab = "UMAP1", ylab = "UMAP2",
title = paste0("Mean of ", gene_id)
)
gg + theme_void()
The fact that schex
renders ggplot
objects can also be used to save these plots. Simply use ggsave
in order to save any created plot.
To find the details of the session for reproducibility, use this:
sessionInfo()
#> R version 4.4.0 RC (2024-04-16 r86468)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggrepel_0.9.5 scran_1.33.0
#> [3] scater_1.33.0 scuttle_1.15.0
#> [5] TENxPBMCData_1.21.0 HDF5Array_1.33.0
#> [7] rhdf5_2.49.0 DelayedArray_0.31.0
#> [9] SparseArray_1.5.0 S4Arrays_1.5.0
#> [11] abind_1.4-5 Matrix_1.7-0
#> [13] igraph_2.0.3 Seurat_5.0.3
#> [15] SeuratObject_5.0.1 sp_2.1-4
#> [17] dplyr_1.1.4 schex_1.19.0
#> [19] SingleCellExperiment_1.27.0 SummarizedExperiment_1.35.0
#> [21] Biobase_2.65.0 GenomicRanges_1.57.0
#> [23] GenomeInfoDb_1.41.0 IRanges_2.39.0
#> [25] S4Vectors_0.43.0 BiocGenerics_0.51.0
#> [27] MatrixGenerics_1.17.0 matrixStats_1.3.0
#> [29] ggplot2_3.5.1
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.22 splines_4.4.0
#> [3] later_1.3.2 filelock_1.0.3
#> [5] tibble_3.2.1 polyclip_1.10-6
#> [7] fastDummies_1.7.3 lifecycle_1.0.4
#> [9] edgeR_4.3.0 globals_0.16.3
#> [11] lattice_0.22-6 MASS_7.3-60.2
#> [13] magrittr_2.0.3 limma_3.61.0
#> [15] plotly_4.10.4 sass_0.4.9
#> [17] rmarkdown_2.26 jquerylib_0.1.4
#> [19] yaml_2.3.8 metapod_1.13.0
#> [21] httpuv_1.6.15 sctransform_0.4.1
#> [23] spam_2.10-0 spatstat.sparse_3.0-3
#> [25] reticulate_1.36.1 cowplot_1.1.3
#> [27] pbapply_1.7-2 DBI_1.2.2
#> [29] RColorBrewer_1.1-3 zlibbioc_1.51.0
#> [31] Rtsne_0.17 purrr_1.0.2
#> [33] rappdirs_0.3.3 tweenr_2.0.3
#> [35] GenomeInfoDbData_1.2.12 irlba_2.3.5.1
#> [37] listenv_0.9.1 spatstat.utils_3.0-4
#> [39] goftest_1.2-3 RSpectra_0.16-1
#> [41] dqrng_0.3.2 spatstat.random_3.2-3
#> [43] fitdistrplus_1.1-11 parallelly_1.37.1
#> [45] DelayedMatrixStats_1.27.0 leiden_0.4.3.1
#> [47] codetools_0.2-20 ggforce_0.4.2
#> [49] tidyselect_1.2.1 UCSC.utils_1.1.0
#> [51] farver_2.1.1 viridis_0.6.5
#> [53] ScaledMatrix_1.13.0 BiocFileCache_2.13.0
#> [55] spatstat.explore_3.2-7 jsonlite_1.8.8
#> [57] BiocNeighbors_1.23.0 progressr_0.14.0
#> [59] ggridges_0.5.6 survival_3.6-4
#> [61] systemfonts_1.0.6 tools_4.4.0
#> [63] ica_1.0-3 Rcpp_1.0.12
#> [65] glue_1.7.0 gridExtra_2.3
#> [67] xfun_0.43 withr_3.0.0
#> [69] BiocManager_1.30.22 fastmap_1.1.1
#> [71] bluster_1.15.0 rhdf5filters_1.17.0
#> [73] fansi_1.0.6 rsvd_1.0.5
#> [75] entropy_1.3.1 digest_0.6.35
#> [77] R6_2.5.1 mime_0.12
#> [79] colorspace_2.1-0 scattermore_1.2
#> [81] tensor_1.5 spatstat.data_3.0-4
#> [83] RSQLite_2.3.6 utf8_1.2.4
#> [85] tidyr_1.3.1 generics_0.1.3
#> [87] hexbin_1.28.3 data.table_1.15.4
#> [89] FNN_1.1.4 httr_1.4.7
#> [91] htmlwidgets_1.6.4 uwot_0.2.2
#> [93] pkgconfig_2.0.3 gtable_0.3.5
#> [95] blob_1.2.4 lmtest_0.9-40
#> [97] XVector_0.45.0 htmltools_0.5.8.1
#> [99] dotCall64_1.1-1 scales_1.3.0
#> [101] png_0.1-8 knitr_1.46
#> [103] reshape2_1.4.4 curl_5.2.1
#> [105] nlme_3.1-164 cachem_1.0.8
#> [107] zoo_1.8-12 stringr_1.5.1
#> [109] BiocVersion_3.20.0 KernSmooth_2.23-22
#> [111] vipor_0.4.7 parallel_4.4.0
#> [113] miniUI_0.1.1.1 concaveman_1.1.0
#> [115] AnnotationDbi_1.67.0 pillar_1.9.0
#> [117] grid_4.4.0 vctrs_0.6.5
#> [119] RANN_2.6.1 promises_1.3.0
#> [121] BiocSingular_1.21.0 beachmat_2.21.0
#> [123] dbplyr_2.5.0 xtable_1.8-4
#> [125] cluster_2.1.6 beeswarm_0.4.0
#> [127] evaluate_0.23 locfit_1.5-9.9
#> [129] cli_3.6.2 compiler_4.4.0
#> [131] rlang_1.1.3 crayon_1.5.2
#> [133] future.apply_1.11.2 labeling_0.4.3
#> [135] ggbeeswarm_0.7.2 plyr_1.8.9
#> [137] stringi_1.8.3 BiocParallel_1.39.0
#> [139] viridisLite_0.4.2 deldir_2.0-4
#> [141] Biostrings_2.73.0 munsell_0.5.1
#> [143] lazyeval_0.2.2 spatstat.geom_3.2-9
#> [145] V8_4.4.2 ExperimentHub_2.13.0
#> [147] RcppHNSW_0.6.0 patchwork_1.2.0
#> [149] sparseMatrixStats_1.17.0 bit64_4.0.5
#> [151] future_1.33.2 Rhdf5lib_1.27.0
#> [153] statmod_1.5.0 KEGGREST_1.45.0
#> [155] shiny_1.8.1.1 highr_0.10
#> [157] AnnotationHub_3.13.0 ROCR_1.0-11
#> [159] memoise_2.0.1 bslib_0.7.0
#> [161] bit_4.0.5