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
No description has been provided for this image
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
)
No description has been provided for this image
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
)
No description has been provided for this image
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
)
No description has been provided for this image
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