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Ā [2]:
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Ā [3]:
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Ā [4]:
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Ā [5]:
options(repr.plot.width = 7, repr.plot.height = 6)
SpatFeaturePlot(
visiumHD,
features = "nr_feats",
points_size = 6,
points_shape = 21
)
InĀ [6]:
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Ā [7]:
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Ā [11]:
options(repr.plot.width = 11, repr.plot.height = 5)
p1 <- CellDimPlot(visiumHD, group_by = "leiden_clus", label = TRUE)
p2 <- CellDimPlot(visiumHD, group_by = "leiden_clus", reduction = "umap", label = TRUE)
p1 + p2
InĀ [8]:
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Ā [9]:
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] sf_1.0-20 devtools_2.4.5 [7] httr_1.4.7 RColorBrewer_1.1-3 [9] repr_1.1.7 profvis_0.4.0 [11] tools_4.4.3 sctransform_0.4.2 [13] backports_1.5.0 R6_2.6.1 [15] uwot_0.2.3 lazyeval_0.2.2 [17] urlchecker_1.0.1 withr_3.0.2 [19] sp_2.2-0 gridExtra_2.3 [21] GiottoUtils_0.2.5 progressr_0.15.1 [23] quantreg_6.00 cli_3.6.5 [25] Biobase_2.62.0 spatstat.explore_3.4-3 [27] fastDummies_1.7.5 iNEXT_3.0.1 [29] labeling_0.4.3 Seurat_5.3.0 [31] arrow_19.0.1 spatstat.data_3.1-6 [33] proxy_0.4-27 ggridges_0.5.6 [35] pbapply_1.7-2 pbdZMQ_0.3-14 [37] dbscan_1.2.2 R.utils_2.13.0 [39] stringdist_0.9.15 parallelly_1.45.0 [41] sessioninfo_1.2.3 VGAM_1.1-13 [43] rstudioapi_0.17.1 FNN_1.1.4.1 [45] generics_0.1.4 shape_1.4.6.1 [47] gtools_3.9.5 ica_1.0-3 [49] spatstat.random_3.4-1 dplyr_1.1.4 [51] Matrix_1.7-3 S4Vectors_0.40.2 [53] abind_1.4-5 R.methodsS3_1.8.2 [55] terra_1.8-42 lifecycle_1.0.4 [57] SummarizedExperiment_1.32.0 SparseArray_1.2.4 [59] Rtsne_0.17 grid_4.4.3 [61] promises_1.3.2 crayon_1.5.3 [63] miniUI_0.1.2 lattice_0.22-7 [65] beachmat_2.18.1 cowplot_1.1.3 [67] pillar_1.10.2 GenomicRanges_1.54.1 [69] rjson_0.2.23 future.apply_1.20.0 [71] codetools_0.2-20 glue_1.8.0 [73] spatstat.univar_3.1-3 data.table_1.17.4 [75] remotes_2.5.0 vctrs_0.6.5 [77] png_0.1-8 spam_2.11-1 [79] gtable_0.3.6 assertthat_0.2.1 [81] cachem_1.1.0 S4Arrays_1.2.1 [83] mime_0.13 tidygraph_1.3.0 [85] survival_3.8-3 SingleCellExperiment_1.24.0 [87] units_0.8-5 ellipsis_0.3.2 [89] scRepertoire_2.2.1 fitdistrplus_1.2-2 [91] ROCR_1.0-11 nlme_3.1-168 [93] usethis_3.1.0 bit64_4.6.0-1 [95] RcppAnnoy_0.0.22 evd_2.3-7.1 [97] GenomeInfoDb_1.38.8 rprojroot_2.0.4 [99] irlba_2.3.5.1 KernSmooth_2.23-26 [101] DBI_1.2.3 plotthis_0.7.1 [103] colorspace_2.1-1 BiocGenerics_0.48.1 [105] tidyselect_1.2.1 bit_4.6.0 [107] compiler_4.4.3 SparseM_1.84-2 [109] xml2_1.3.8 desc_1.4.3 [111] ggdendro_0.2.0 DelayedArray_0.28.0 [113] plotly_4.10.4 checkmate_2.3.2 [115] scales_1.4.0 classInt_0.4-11 [117] lmtest_0.9-40 rappdirs_0.3.3 [119] stringr_1.5.1 digest_0.6.37 [121] goftest_1.2-3 spatstat.utils_3.1-4 [123] XVector_0.42.0 htmltools_0.5.8.1 [125] GiottoVisuals_0.2.12 pkgconfig_2.0.3 [127] base64enc_0.1-3 MatrixGenerics_1.14.0 [129] fastmap_1.2.0 rlang_1.1.6 [131] GlobalOptions_0.1.2 htmlwidgets_1.6.4 [133] shiny_1.10.0 farver_2.1.2 [135] zoo_1.8-14 jsonlite_2.0.0 [137] BiocParallel_1.36.0 R.oo_1.27.1 [139] BiocSingular_1.18.0 RCurl_1.98-1.17 [141] magrittr_2.0.3 GenomeInfoDbData_1.2.11 [143] dotCall64_1.2 patchwork_1.3.0 [145] IRkernel_1.3.2 Rcpp_1.0.14 [147] evmix_2.12 ggnewscale_0.5.1 [149] viridis_0.6.5 reticulate_1.42.0 [151] truncdist_1.0-2 stringi_1.8.7 [153] ggalluvial_0.12.5 ggraph_2.2.1 [155] zlibbioc_1.48.2 MASS_7.3-64 [157] plyr_1.8.9 pkgbuild_1.4.8 [159] parallel_4.4.3 listenv_0.9.1 [161] ggrepel_0.9.6 forcats_1.0.0 [163] deldir_2.0-4 graphlayouts_1.2.2 [165] IRdisplay_1.1 splines_4.4.3 [167] gridtext_0.1.5 tensor_1.5 [169] circlize_0.4.16 colorRamp2_0.1.0 [171] igraph_2.0.3 uuid_1.2-1 [173] spatstat.geom_3.4-1 cubature_2.1.4 [175] RcppHNSW_0.6.0 ScaledMatrix_1.10.0 [177] reshape2_1.4.4 stats4_4.4.3 [179] pkgload_1.4.0 evaluate_1.0.3 [181] SeuratObject_5.1.0 tweenr_2.0.3 [183] httpuv_1.6.15 MatrixModels_0.5-4 [185] RANN_2.6.2 tidyr_1.3.1 [187] purrr_1.0.4 polyclip_1.10-7 [189] future_1.58.0 scattermore_1.2 [191] ggplot2_3.5.2 ggforce_0.4.2 [193] rsvd_1.0.5 xtable_1.8-4 [195] e1071_1.7-16 RSpectra_0.16-2 [197] later_1.4.2 class_7.3-23 [199] viridisLite_0.4.2 gsl_2.1-8 [201] tibble_3.2.1 memoise_2.0.1 [203] IRanges_2.36.0 cluster_2.1.8.1 [205] globals_0.18.0