scplotter to work with 10x Visium HD data prepared by Giotto¶
See: https://drieslab.github.io/giotto_workshop_2024/visium-hd.html
Go back to scplotter documentation: https://pwwang.github.io/scplotter/
InĀ [1]:
library(Giotto)
# Ensure Giotto can access a python env
genv_exists <- suppressMessages(checkGiottoEnvironment())
print(genv_exists)
python_path <- file.path(Sys.getenv("HOME"), "miniconda3", "envs", "giotto_env", "bin", "python")
Sys.setenv(RETICULATE_PYTHON = python_path)
invisible(capture.output(suppressMessages(set_giotto_python_path(python_path = python_path))))
# library(scplotter)
devtools::load_all()
Loading required package: GiottoClass
Newer devel version of GiottoClass on GitHub: 0.4.8 Giotto Suite 4.2.1
[1] TRUE
ā¹ Loading scplotter
InĀ [9]:
library(Giotto)
## Set instructions
results_folder <- "data/Human_Colorectal_Cancer_workshop.results"
instructions <- createGiottoInstructions(
save_dir = results_folder,
save_plot = FALSE,
show_plot = TRUE,
return_plot = TRUE,
python_path = python_path
)
## Provide the path to the visium folder
data_path <- "data/Human_Colorectal_Cancer_workshop/square_002um"
expression_path <- file.path(data_path, 'raw_feature_bc_matrix')
expr_results <- get10Xmatrix(path_to_data = expression_path,
gene_column_index = 1)
tissue_positions_path <- file.path(data_path, 'spatial/tissue_positions.parquet')
tissue_positions <- data.table::as.data.table(arrow::read_parquet(tissue_positions_path))
matrix_tile_dt <- data.table::as.data.table(Matrix::summary(expr_results))
genes <- expr_results@Dimnames[[1]]
samples <- expr_results@Dimnames[[2]]
matrix_tile_dt[, gene := genes[i]]
matrix_tile_dt[, pixel := samples[j]]
expr_pos_data <- data.table::merge.data.table(matrix_tile_dt,
tissue_positions,
by.x = 'pixel',
by.y = 'barcode')
expr_pos_data <- expr_pos_data[,.(pixel, pxl_row_in_fullres, pxl_col_in_fullres, gene, x)]
colnames(expr_pos_data) = c('pixel', 'x', 'y', 'gene', 'count')
giotto_points = createGiottoPoints(x = expr_pos_data[,.(x, y, gene, pixel, count)])
hexbin400 <- tessellate(extent = ext(giotto_points),
shape = 'hexagon',
shape_size = 400,
name = 'hex400')
# gpoints provides spatial gene expression information
# gpolygons provides spatial unit information (here = hexagon tiles)
visiumHD = createGiottoObjectSubcellular(gpoints = list('rna' = giotto_points),
gpolygons = list('hex400' = hexbin400),
instructions = instructions)
# create spatial centroids for each spatial unit (hexagon)
visiumHD = addSpatialCentroidLocations(gobject = visiumHD,
poly_info = 'hex400')
visiumHD
python already initialized in this session active environment : 'giotto_env' python version : 3.10
Warning message:
āPotentially unsafe or invalid elements have been discarded from R metadata.
ā¹ Type: "externalptr"
ā If you trust the source, you can set `options(arrow.unsafe_metadata = TRUE)` to preserve them.ā
Selecting col "gene" as feat_ID column
Selecting cols "x" and "y" as x and y respectively
367 polygons generated
polygonlist is a list with names
[ hex400 ] Process polygon info...
pointslist is a named list
[ rna ] Process point info...
Start centroid calculation for polygon information
layer: hex400
An object of class giotto [SUBCELLULAR INFO] polygons : hex400 features : rna [AGGREGATE INFO] spatial locations ---------------- [hex400] raw Use objHistory() to see steps and params used
InĀ [25]:
options(repr.plot.width = 14, repr.plot.height = 6)
# devtools::load_all()
p1 <- SpatDimPlot(
visiumHD,
image = "black",
shapes = TRUE,
features = GiottoClass::featIDs(visiumHD, 'rna')[1:20],
spat_unit = 'hex400',
shapes_feat_type = 'hex400',
shapes_border_size = 0.1,
shapes_border_color = "white",
shapes_fill_by = "black",
shapes_alpha = 0.5,
points_size = 0.25
)
p2 <- SpatDimPlot(
visiumHD,
image = "black",
shapes = TRUE,
features = GiottoClass::featIDs(visiumHD, 'rna')[1:20],
spat_unit = 'hex400',
shapes_feat_type = 'hex400',
shapes_border_size = 0.1,
shapes_border_color = "white",
shapes_fill_by = "black",
shapes_alpha = 0.5,
points_size = 0.25,
raster = TRUE
)
p1 + p2
InĀ [26]:
visiumHD = calculateOverlap(visiumHD,
spatial_info = 'hex400',
feat_info = 'rna')
# convert overlap results to bin by gene matrix
visiumHD = overlapToMatrix(visiumHD,
poly_info = 'hex400',
feat_info = 'rna',
name = 'raw')
# this action will automatically create an active spatial unit, ie. hexbin 400
activeSpatUnit(visiumHD)
# filter on gene expression matrix
visiumHD <- filterGiotto(visiumHD,
expression_threshold = 1,
feat_det_in_min_cells = 5,
min_det_feats_per_cell = 25)
# normalize and scale gene expression data
visiumHD <- normalizeGiotto(visiumHD,
scalefactor = 1000,
verbose = T)
# add cell and gene statistics
visiumHD <- addStatistics(visiumHD)
1. convert polygon to raster 2. overlap raster and points 3. add polygon information 4. add points information 5. create overlap polygon information
'hex400'
completed 1: preparation completed 2: subset expression data completed 3: subset spatial locations completed 4: subset cell metadata completed 5: subset feature metadata completed 6: subset spatial network(s) completed 7: subsetted dimension reductions completed 8: subsetted nearest network(s) completed 9: subsetted spatial enrichment results for hex400 --> hex400 found back in polygon layer: hex400 completed 10: subsetted spatial information data subset feature info: rna completed 11: subsetted spatial feature data
Feature type: rna Number of cells removed: 0 out of 367 Number of feats removed: 349 out of 4967
first scale feats and then cells Setting expression [hex400][rna] normalized Setting expression [hex400][rna] scaled calculating statistics for "normalized" expression
InĀ [30]:
options(repr.plot.width = 7, repr.plot.height = 6)
SpatFeaturePlot(
visiumHD,
features = "nr_feats",
points_size = 6,
points_shape = 21
)
InĀ [34]:
options(repr.plot.width = 7, repr.plot.height = 6)
SpatFeaturePlot(
visiumHD,
image = "black",
shapes = TRUE,
spat_unit = 'hex400',
use_overlap = TRUE,
shapes_fill_by = "nr_feats",
shapes_feat_type = 'hex400',
points_size = 0.1
)
InĀ [35]:
visiumHD <- calculateHVF(visiumHD,
zscore_threshold = 1)
visiumHD <- runPCA(visiumHD,
expression_values = 'normalized',
feats_to_use = 'hvf')
visiumHD <- runUMAP(visiumHD,
dimensions_to_use = 1:14,
n_threads = 10)
# sNN network (default)
visiumHD <- createNearestNetwork(visiumHD,
dimensions_to_use = 1:14,
k = 5)
## leiden clustering ####
visiumHD <- doLeidenClusterIgraph(visiumHD, resolution = 0.5, n_iterations = 1000, spat_unit = 'hex400')
"hvf" column was found in the feats metadata information and will be used to select highly variable features Setting dimension reduction [hex400][rna] pca Setting dimension reduction [hex400][rna] umap
InĀ [37]:
options(repr.plot.width = 7, repr.plot.height = 6)
visiumHD@cell_metadata$hex400$rna$leiden_clus <- as.factor(visiumHD@cell_metadata$hex400$rna$leiden_clus)
SpatFeaturePlot(
visiumHD,
image = "black",
shapes = TRUE,
spat_unit = 'hex400',
shapes_fill_by = "leiden_clus",
shapes_feat_type = 'hex400',
points_size = 0.1
)
InĀ [2]:
x <- sessionInfo()
x <- capture.output(print(x))
# hide the BLAS/LAPACK paths
x <- x[!startsWith(x, "BLAS/LAPACK:")]
cat(paste(x, collapse = "\n"))
R version 4.4.3 (2025-02-28) Platform: x86_64-conda-linux-gnu Running under: Red Hat Enterprise Linux 8.10 (Ootpa) Matrix products: default locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 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/Chicago tzcode source: system (glibc) attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] scplotter_0.4.0 Giotto_4.2.1 GiottoClass_0.4.7 loaded via a namespace (and not attached): [1] fs_1.6.6 matrixStats_1.5.0 [3] spatstat.sparse_3.1-0 bitops_1.0-9 [5] devtools_2.4.5 httr_1.4.7 [7] RColorBrewer_1.1-3 repr_1.1.7 [9] profvis_0.4.0 tools_4.4.3 [11] sctransform_0.4.2 backports_1.5.0 [13] R6_2.6.1 uwot_0.2.3 [15] lazyeval_0.2.2 urlchecker_1.0.1 [17] withr_3.0.2 sp_2.2-0 [19] gridExtra_2.3 GiottoUtils_0.2.5 [21] progressr_0.15.1 quantreg_6.00 [23] cli_3.6.5 Biobase_2.62.0 [25] spatstat.explore_3.4-3 fastDummies_1.7.5 [27] iNEXT_3.0.1 Seurat_5.3.0 [29] spatstat.data_3.1-6 ggridges_0.5.6 [31] pbapply_1.7-2 pbdZMQ_0.3-14 [33] stringdist_0.9.15 parallelly_1.45.0 [35] sessioninfo_1.2.3 VGAM_1.1-13 [37] rstudioapi_0.17.1 generics_0.1.4 [39] shape_1.4.6.1 gtools_3.9.5 [41] ica_1.0-3 spatstat.random_3.4-1 [43] dplyr_1.1.4 Matrix_1.7-3 [45] S4Vectors_0.40.2 abind_1.4-5 [47] terra_1.8-42 lifecycle_1.0.4 [49] SummarizedExperiment_1.32.0 SparseArray_1.2.4 [51] Rtsne_0.17 grid_4.4.3 [53] promises_1.3.2 crayon_1.5.3 [55] miniUI_0.1.2 lattice_0.22-7 [57] cowplot_1.1.3 pillar_1.10.2 [59] GenomicRanges_1.54.1 rjson_0.2.23 [61] future.apply_1.20.0 codetools_0.2-20 [63] glue_1.8.0 spatstat.univar_3.1-3 [65] data.table_1.17.4 remotes_2.5.0 [67] vctrs_0.6.5 png_0.1-8 [69] spam_2.11-1 gtable_0.3.6 [71] assertthat_0.2.1 cachem_1.1.0 [73] S4Arrays_1.2.1 mime_0.13 [75] tidygraph_1.3.0 survival_3.8-3 [77] SingleCellExperiment_1.24.0 ellipsis_0.3.2 [79] scRepertoire_2.2.1 fitdistrplus_1.2-2 [81] ROCR_1.0-11 nlme_3.1-168 [83] usethis_3.1.0 RcppAnnoy_0.0.22 [85] evd_2.3-7.1 GenomeInfoDb_1.38.8 [87] rprojroot_2.0.4 irlba_2.3.5.1 [89] KernSmooth_2.23-26 plotthis_0.7.0 [91] colorspace_2.1-1 BiocGenerics_0.48.1 [93] tidyselect_1.2.1 compiler_4.4.3 [95] SparseM_1.84-2 xml2_1.3.8 [97] desc_1.4.3 ggdendro_0.2.0 [99] DelayedArray_0.28.0 plotly_4.10.4 [101] checkmate_2.3.2 scales_1.4.0 [103] lmtest_0.9-40 rappdirs_0.3.3 [105] stringr_1.5.1 digest_0.6.37 [107] goftest_1.2-3 spatstat.utils_3.1-4 [109] XVector_0.42.0 htmltools_0.5.8.1 [111] GiottoVisuals_0.2.12 pkgconfig_2.0.3 [113] base64enc_0.1-3 MatrixGenerics_1.14.0 [115] fastmap_1.2.0 rlang_1.1.6 [117] GlobalOptions_0.1.2 htmlwidgets_1.6.4 [119] shiny_1.10.0 farver_2.1.2 [121] zoo_1.8-14 jsonlite_2.0.0 [123] RCurl_1.98-1.17 magrittr_2.0.3 [125] GenomeInfoDbData_1.2.11 dotCall64_1.2 [127] patchwork_1.3.0 IRkernel_1.3.2 [129] Rcpp_1.0.14 evmix_2.12 [131] ggnewscale_0.5.1 viridis_0.6.5 [133] reticulate_1.42.0 truncdist_1.0-2 [135] stringi_1.8.7 ggalluvial_0.12.5 [137] ggraph_2.2.1 zlibbioc_1.48.2 [139] MASS_7.3-64 plyr_1.8.9 [141] pkgbuild_1.4.8 parallel_4.4.3 [143] listenv_0.9.1 ggrepel_0.9.6 [145] forcats_1.0.0 deldir_2.0-4 [147] graphlayouts_1.2.2 IRdisplay_1.1 [149] splines_4.4.3 gridtext_0.1.5 [151] tensor_1.5 circlize_0.4.16 [153] colorRamp2_0.1.0 igraph_2.0.3 [155] uuid_1.2-1 spatstat.geom_3.4-1 [157] cubature_2.1.4 RcppHNSW_0.6.0 [159] reshape2_1.4.4 stats4_4.4.3 [161] pkgload_1.4.0 evaluate_1.0.3 [163] SeuratObject_5.1.0 tweenr_2.0.3 [165] httpuv_1.6.15 MatrixModels_0.5-4 [167] RANN_2.6.2 tidyr_1.3.1 [169] purrr_1.0.4 polyclip_1.10-7 [171] future_1.58.0 scattermore_1.2 [173] ggplot2_3.5.2 ggforce_0.4.2 [175] xtable_1.8-4 RSpectra_0.16-2 [177] later_1.4.2 viridisLite_0.4.2 [179] gsl_2.1-8 tibble_3.2.1 [181] memoise_2.0.1 IRanges_2.36.0 [183] cluster_2.1.8.1 globals_0.18.0