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