scplotter to work with Multi-omics Spatial CITE-Seq data prepared by Giotto¶
See: https://drieslab.github.io/Giotto_website/articles/spatial_citeseq_human_skin.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/Giotto_Spatial_CITE-Seq.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/Giotto_Spatial_CITE-Seq"
x <- data.table::fread(file.path(data_path, "GSE213264_RAW/GSM6578065_humanskin_RNA.tsv.gz"))
spatial_coords <- data.frame(cell_ID = x$X)
spatial_coords <- cbind(spatial_coords,
stringr::str_split_fixed(spatial_coords$cell_ID,
pattern = "x",
n = 2))
colnames(spatial_coords)[2:3] = c("sdimx", "sdimy")
spatial_coords$sdimx <- as.integer(spatial_coords$sdimx)
spatial_coords$sdimy <- as.integer(spatial_coords$sdimy)
spatial_coords$sdimy <- spatial_coords$sdimy*(-1)
rna_matrix <- data.table::fread(file.path(data_path, "GSE213264_RAW/GSM6578065_humanskin_RNA.tsv.gz"))
rna_matrix <- rna_matrix[rna_matrix$X %in% spatial_coords$cell_ID,]
rna_matrix <- rna_matrix[match(spatial_coords$cell_ID, rna_matrix$X),]
rna_matrix <- t(rna_matrix[,-1])
colnames(rna_matrix) <- spatial_coords$cell_ID
protein_matrix <- data.table::fread(file.path(data_path, "GSE213264_RAW/GSM6578074_humanskin_protein.tsv.gz"))
protein_matrix <- protein_matrix[protein_matrix$X %in% spatial_coords$cell_ID,]
protein_matrix <- protein_matrix[match(spatial_coords$cell_ID, protein_matrix$X),]
protein_matrix <- t(protein_matrix[,-1])
colnames(protein_matrix) <- spatial_coords$cell_ID
my_giotto_object <- createGiottoObject(expression = list(rna = list(raw = rna_matrix),
protein = list(raw = protein_matrix)),
expression_feat = list("rna", "protein"),
spatial_locs = spatial_coords,
instructions = instructions)
my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
do_manual_adj = TRUE,
scale_factor = 0.5,
mg_object = file.path(data_path, "skin.jpg"),
negative_y = TRUE)
my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
images = list(my_giotto_image),
spat_loc_name = "raw")
force(my_giotto_object)
python already initialized in this session active environment : 'giotto_env' python version : 3.10 do_manual_adj == TRUE Boundaries will be adjusted by given values.
An object of class giotto >Active spat_unit: cell >Active feat_type: rna dimensions : 15486, 1691 (features, cells) [SUBCELLULAR INFO] [AGGREGATE INFO] expression ----------------------- [cell][rna] raw [cell][protein] raw spatial locations ---------------- [cell] raw attached images ------------------ images : image Use objHistory() to see steps and params used
In [3]:
options(repr.plot.width = 6, repr.plot.height = 6)
# devtools::load_all()
SpatDimPlot(my_giotto_object, points_size = 2, points_color_by = "lightblue")
In [ ]:
# RNA
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)
# Protein
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)
# RNA
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
norm_methods = "standard",
scalefactor = 10000,
verbose = TRUE)
# Protein
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
scalefactor = 6000,
verbose = TRUE)
# RNA
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna")
# Protein
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized")
# RNA
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
name = "rna.pca")
# Protein
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
scale_unit = TRUE,
center = FALSE,
method = "factominer")
# RNA
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
dimensions_to_use = 1:10,
dim_reduction_name = "rna.pca")
# Protein
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
dimensions_to_use = 1:10)
# RNA
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
type = "sNN",
dim_reduction_to_use = "pca",
dim_reduction_name = "rna.pca",
dimensions_to_use = 1:10,
k = 20)
# Protein
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
type = "sNN",
name = "protein_sNN.pca",
dimensions_to_use = 1:10,
k = 20)
# RNA
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "sNN",
name = "leiden_clus",
resolution = 1)
# Protein
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
nn_network_to_use = "sNN",
network_name = "protein_sNN.pca",
name = "leiden_clus",
resolution = 1)
my_giotto_object@cell_metadata$cell$rna$leiden_clus <- as.character(my_giotto_object@cell_metadata$cell$rna$leiden_clus)
my_giotto_object@cell_metadata$cell$protein$leiden_clus <- as.character(my_giotto_object@cell_metadata$cell$protein$leiden_clus)
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
type = "kNN",
dim_reduction_name = "rna.pca",
name = "rna_kNN.pca",
dimensions_to_use = 1:10,
k = 20)
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
type = "kNN",
name = "protein_kNN.pca",
dimensions_to_use = 1:10,
k = 20)
my_giotto_object <- runWNN(my_giotto_object,
# modality_1 = "rna",
# modality_2 = "protein",
# pca_name_modality_1 = "rna.pca",
# pca_name_modality_2 = "protein.pca",
k = 20)
my_giotto_object <- runIntegratedUMAP(my_giotto_object)
# modality1 = "rna",
# modality2 = "protein")
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "kNN",
network_name = "integrated_kNN",
name = "integrated_leiden_clus",
resolution = 0.7)
my_giotto_object <- createSpatialNetwork(gobject = my_giotto_object,
method = "kNN",
k = 6,
maximum_distance_knn = 5,
name = "spatial_network")
rank_spatialfeats <- binSpect(my_giotto_object,
bin_method = "rank",
calc_hub = TRUE,
hub_min_int = 5,
spatial_network_name = "spatial_network")
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
Feature type: rna Number of cells removed: 0 out of 1691 Number of feats removed: 363 out of 15486
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
Feature type: protein Number of cells removed: 0 out of 1691 Number of feats removed: 0 out of 283
first scale feats and then cells
Setting expression [cell][rna] normalized Setting expression [cell][rna] scaled Warning message in .rna_standard_normalization(gobject = gobject, raw_expr = raw_expr, : “Caution: Standard normalization was developed for RNA data ” first scale feats and then cells Setting expression [cell][protein] normalized Setting expression [cell][protein] scaled calculating statistics for "normalized" expression calculating statistics for "normalized" expression "hvf" was not found in the gene metadata information. all genes will be used. Setting dimension reduction [cell][rna] rna.pca "hvf" was not found in the gene metadata information. all genes will be used. Setting dimension reduction [cell][protein] protein.pca Setting dimension reduction [cell][rna] umap Setting dimension reduction [cell][protein] protein.umap The NN network name was not specified, default to the first: "rna_kNN.pca" The NN network name was not specified, default to the first: "protein_kNN.pca" Calculating rna-rna distance Calculating protein-protein distance Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics' Also defined by ‘spam’ Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics' Also defined by ‘spam’ Calculating integrated Nearest Neighbors Creating igraph Setting nearest neighbor network [cell][rna] integrated_kNN Calculating integrated UMAP Setting dimension reduction [cell][rna] integrated.umap Setting dimension reduction [cell][protein] integrated.umap Setting spatial network [cell] spatial_network This is the single parameter version of binSpect 1. matrix binarization complete 2. spatial enrichment test completed 3. (optional) average expression of high expressing cells calculated 4. (optional) number of high expressing cells calculated
In [45]:
im <- terra::rast(as.numeric(my_giotto_object@images$image@mg_object[[1]]), extent = c(1, 50, -50, 1))
im
class : SpatRaster dimensions : 2700, 2704, 3 (nrow, ncol, nlyr) resolution : 0.0181213, 0.01888889 (x, y) extent : 1, 50, -50, 1 (xmin, xmax, ymin, ymax) coord. ref. : source(s) : memory names : lyr.1, lyr.2, lyr.3 min values : 0, 0, 0 max values : 1, 1, 1
In [33]:
terra::res(im)
- 1
- 1
In [ ]:
options(repr.plot.width = 14, repr.plot.height = 7)
# devtools::load_all()
p1 <- SpatDimPlot(my_giotto_object, image = TRUE, group_by = "leiden_clus",
title = "RNA Leiden clustering")
p2 <- SpatDimPlot(my_giotto_object, image = TRUE, group_by = "leiden_clus",
spat_unit = "cell", feat_type = "protein", title = "Protein Leiden clustering")
p1 + p2
In [55]:
options(repr.plot.width = 7, repr.plot.height = 6)
# devtools::load_all()
my_giotto_object@cell_metadata$cell$rna$integrated_leiden_clus <- as.character(my_giotto_object@cell_metadata$cell$rna$integrated_leiden_clus)
# change the order
my_giotto_object@cell_metadata$cell$rna$integrated_leiden_clus <- factor(my_giotto_object@cell_metadata$cell$rna$integrated_leiden_clus,
levels = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9"))
p1 <- SpatDimPlot(my_giotto_object, image = TRUE, group_by = "integrated_leiden_clus",
title = "Integrated Leiden clustering", points_shape = 21, points_size = 2,
points_border_color = "grey20")
p1
In [60]:
options(repr.plot.width = 10, repr.plot.height = 14)
SpatFeaturePlot(my_giotto_object,
layer = "scaled",
features = rank_spatialfeats$feats[1:6],
points_palette = "Reds",
facet_ncol = 2)
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