This vignette outlines the steps of inference, analysis and visualization of niche for a Visium HD spatial dataset using SpatialNiche. We showcase SpatialNiche’s diverse functionalities by applying it to a Visium HD data on cells from Lymphangioleiomyomatosis (LAM, diseased) human lung from a patient.

SpatialNiche requires raw Visium HD spaceranger output and optionally processed gene expression data of cells/bins from Seurat objects as the user input and models the niche based on a specific cell type or gene expression.

Upon infering the niche, SpatialNichet provides functionality for further data exploration, analysis, and visualization.

Load the required libraries

library(SpatialNiche)
library(ggplot2) #for visualization
library(dplyr) #for data wrangling
library(scattermore) #for spatial plotting

Part I: Data input & processing and initialization of SpatialNiche list containing barcodes and images tibble

SpatialNiche requires two user inputs: one is the path to outs directory of the spaceranger output of the Visium HD dataset, and the other is the bins size you want to load.

Starting from spaceranger output

localdir <- "../../../4_em/outs" # path to outs directory of the spaceranger output of the Visium HD dataset

DataHD<-Load_Spatial_Data(localdir, size = '008um') #loading the data using the path and the bin size

bcsHD<-DataHD$barcodes #pulling the barcodes tibble out for easier access

head(bcsHD) #look at the top entries
#> # A tibble: 6 × 12
#>   barcode   tissue   row   col imagerow imagecol imagerow_scaled imagecol_scaled
#>   <chr>     <fct>  <int> <int>    <dbl>    <dbl>           <dbl>           <dbl>
#> 1 s_008um_… 0          0     0   27281.   -69.7             620.          -1.58 
#> 2 s_008um_… 0          0     1   27282.   -38.0             620.          -0.863
#> 3 s_008um_… 0          0     2   27282.    -6.19            620.          -0.141
#> 4 s_008um_… 0          0     3   27282.    25.6             620.           0.582
#> 5 s_008um_… 0          0     4   27282.    57.4             620.           1.30 
#> 6 s_008um_… 0          0     5   27282.    89.1             620.           2.03 
#> # ℹ 4 more variables: imagerow_scaled_round <dbl>, imagecol_scaled_round <dbl>,
#> #   height <int>, width <int>

Part II: Adding dconvolution information to SpatialNiche barcodes tibble

SpatialNiche supports import of spacexr’s RCTD based deconvolution information. Thre RCTD result can be directly loaded as an rds object and imported directly into the barcodes tibble.

Starting from RCTD output

DeconvolutionHD<-Load_RCTD('../../../RCTD_d4_SCT.rds') #load rds object of RCTD results

bcsHD<-Add_Deconvolution_Info(bcsHD,DeconvolutionHD,addWeights=FALSE) #use in-built function to import RCTD deconvolution results

bcsHD<-bcsHD %>% na.omit() #remove blank barcodes with no deconvolution information

head(bcsHD[, c(1,13:15)]) #look at the top entries with deconvolution information
#> # A tibble: 6 × 4
#>   barcode             DeconvolutionClass DeconvolutionLabel1 DeconvolutionLabel2
#>   <chr>               <chr>              <chr>               <chr>              
#> 1 s_008um_00102_0021… singlet            IM                  pDC                
#> 2 s_008um_00102_0021… reject             IM                  LAMCORE-3          
#> 3 s_008um_00102_0021… reject             IM                  LAMCORE-1          
#> 4 s_008um_00102_0021… singlet            Mast-Basophil       Plasma             
#> 5 s_008um_00102_0022… reject             AF2                 AT2                
#> 6 s_008um_00102_0022… singlet            LAMCORE-1           Plasma

Part III: Creating niche information

SpatialNiche creates niche information based on a specific cell type or gene expression.The niche information is defined based on a sepcific distance away from the bins/cells of interest.

Defining niche based on a cell type

PeripheryBCs<-Get_Periphery(barcodes=bcsHD, celltype='LAMCORE-1',distance=50,data.dir = localdir) #find barcodes in the cell specific niche

bcsHD$Periphery<-NA #assign empty column to barcodes tibble

bcsHD$Periphery[bcsHD$barcode%in%PeripheryBCs]<-"50 micron" #assign name to bins in the niche 

bcsHD$Periphery[is.na(bcsHD$Periphery)]<-"Tissue" #assign name to bins outside the niche 

bcsHD$Periphery[bcsHD$Periphery=="Tissue" & bcsHD$DeconvolutionLabel1=='LAMCORE-1']<-'LAMCORE-1' #assign name to bins specifying the niche 

head(bcsHD[, c(1,14:16)]) #look at the top entries with cell based niche information
#> # A tibble: 6 × 4
#>   barcode               DeconvolutionLabel1 DeconvolutionLabel2 Periphery
#>   <chr>                 <chr>               <chr>               <chr>    
#> 1 s_008um_00102_00216-1 IM                  pDC                 Tissue   
#> 2 s_008um_00102_00217-1 IM                  LAMCORE-3           Tissue   
#> 3 s_008um_00102_00218-1 IM                  LAMCORE-1           Tissue   
#> 4 s_008um_00102_00219-1 Mast-Basophil       Plasma              Tissue   
#> 5 s_008um_00102_00220-1 AF2                 AT2                 Tissue   
#> 6 s_008um_00102_00221-1 LAMCORE-1           Plasma              LAMCORE-1

Defining niche based on gene expression

em_d4 <- readRDS("../../../em_d4_SCT.rds") #import seurat object

bcsHD<-Add_Expression(bcsHD,em_d4,c("PMEL")) #add gene expression data to barcodes tibble
#> Loading required package: Seurat
#> Loading required package: SeuratObject
#> Loading required package: sp
#> 'SeuratObject' was built with package 'Matrix' 1.7.0 but the current
#> version is 1.7.3; it is recomended that you reinstall 'SeuratObject' as
#> the ABI for 'Matrix' may have changed
#> 
#> Attaching package: 'SeuratObject'
#> The following objects are masked from 'package:base':
#> 
#>     intersect, t

PeripheryGCs<-Get_Gene_Periphery(barcodes=bcsHD, gene=c("PMEL"), seurat=em_d4,distance=50,data.dir = localdir) #find barcodes in the gene specific niche

bcsHD$Peripheryg<-NA #assign empty column to barcodes tibble

bcsHD$Peripheryg[bcsHD$barcode%in%PeripheryGCs]<-"50 micron" #assign name to bins in the niche 

bcsHD$Peripheryg[is.na(bcsHD$Peripheryg)]<-"Tissue" #assign name to bins outside the niche 

bcsHD$Peripheryg[bcsHD$Peripheryg=="Tissue" & bcsHD$PMEL>0]<-'PMEL' #assign name to bins specifying the niche 

head(bcsHD[, c(1,14:18)]) #look at the top entries with gene expression based niche information
#> # A tibble: 6 × 6
#>   barcode     DeconvolutionLabel1 DeconvolutionLabel2 Periphery  PMEL Peripheryg
#>   <chr>       <chr>               <chr>               <chr>     <dbl> <chr>     
#> 1 s_008um_00… IM                  pDC                 Tissue        0 Tissue    
#> 2 s_008um_00… IM                  LAMCORE-3           Tissue        0 50 micron 
#> 3 s_008um_00… IM                  LAMCORE-1           Tissue        0 50 micron 
#> 4 s_008um_00… Mast-Basophil       Plasma              Tissue        0 50 micron 
#> 5 s_008um_00… AF2                 AT2                 Tissue        0 50 micron 
#> 6 s_008um_00… LAMCORE-1           Plasma              LAMCORE-1     0 50 micron

Part IV: Visualizing niche results

Upon infering the niche information, SpatialNiche provides various functionality for further data exploration, analysis, and visualization. Specifically:

slicing the dataset into smaller ROIs

SliceA<-Get_Square('s_008um_00128_00230-1',250,bcsHD) #find barcodes in a Square ROI

head(SliceA) #look at the top barcodes in the Square ROI
#> [1] "s_008um_00112_00214-1" "s_008um_00112_00216-1" "s_008um_00112_00217-1"
#> [4] "s_008um_00112_00218-1" "s_008um_00112_00221-1" "s_008um_00112_00222-1"

SliceB<-Get_Rectangle('s_008um_00331_00466-1',300,0.65,1.1,bcsHD) #find barcodes in a Rectangle ROI, use xfactor and yfactor to control size of the rectangle

head(SliceB) #look at the top barcodes in the Rectangle ROI
#> [1] "s_008um_00311_00456-1" "s_008um_00311_00457-1" "s_008um_00311_00459-1"
#> [4] "s_008um_00311_00460-1" "s_008um_00311_00461-1" "s_008um_00311_00462-1"

SliceC<-Get_Circle('s_008um_00128_00230-1',125,bcsHD, localdir) #find barcodes in a Circle ROI

head(SliceC) #look at the top barcodes in the Circle ROI
#> [1] "s_008um_00113_00226-1" "s_008um_00113_00227-1" "s_008um_00113_00228-1"
#> [4] "s_008um_00113_00229-1" "s_008um_00113_00230-1" "s_008um_00113_00231-1"

Visualizing based on deconvolution information

bcsHD %>% filter(DeconvolutionClass=="singlet") %>%
  ggplot(aes(x = imagecol_scaled, y = -imagerow_scaled,color=DeconvolutionLabel1)) +
  Geom_Spatial(data = DataHD$images, aes(grob = grob), x = 0.5, y = 0.5)+
  geom_scattermore(pointsize = 2,pixels = rep(2000,2))+
  coord_cartesian(expand = FALSE) +
  xlab("") +
  ylab("") +
  theme_set(theme_bw(base_size = 10))+
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) #Visualize singlets by deconvolution label

bcsHD  %>%
  ggplot(aes(x = imagecol_scaled, y = -imagerow_scaled,color=DeconvolutionLabel1)) +
  Geom_Spatial(data = DataHD$images, aes(grob = grob), x = 0.5, y = 0.5) +
  geom_point(shape=16,size=2)+
  coord_cartesian(expand = FALSE) +
  xlab("") +
  ylab("") + 
  theme_set(theme_bw(base_size = 10))+
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) #Visualize all bins by deconvolution label

Visualizing based on niche information

bcA=bcsHD[bcsHD$Periphery!='Tissue',] #Create new tibble with only bins in the cell specific niche

bcA %>%
  ggplot(aes(x = imagecol_scaled, y = -imagerow_scaled,color=Periphery)) +
  Geom_Spatial(data = DataHD$images, aes(grob = grob), x = 0.5, y = 0.5) +
  geom_scattermore(pointsize = 2,pixels = rep(2000,2)) +
  scale_color_manual(values=c('green','red','grey')) +
  coord_cartesian(expand = FALSE) +
  xlab("") +
  ylab("") +
  theme_set(theme_bw(base_size = 10)) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) #Visualize bins in the niche

bcA=bcsHD[bcsHD$Peripheryg!='Tissue',] #Create new tibble with only bins in the gene specific niche

bcA %>%
  ggplot(aes(x = imagecol_scaled, y = -imagerow_scaled,color=Peripheryg)) +
  Geom_Spatial(data = DataHD$images, aes(grob = grob), x = 0.5, y = 0.5) +
  geom_scattermore(pointsize = 2,pixels = rep(2000,2)) +
  scale_color_manual(values=c('green','red','grey')) +
  coord_cartesian(expand = FALSE) +
  xlab("") +
  ylab("") +
  theme_set(theme_bw(base_size = 10)) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) #Visualize bins in the niche

Part V: Calculating distance information

SpatialNiche is also capable of creating distance information by translating the pixel sizes to the micron scale. The results are output as a distane matrix that can be converted to a dataframe for use.

Defining niche based on a cell type

em_d4G1 <- readRDS("../../../KWB/em_d4_G1.rds") #import seurat object

spatial.locs = Seurat::GetTissueCoordinates(em_d4G1, scale = NULL, cols = c("imagerow", "imagecol")) #get coordinates from seurat object

spatial.locs = spatial.locs[,1:2] #remove barcodes from coordinates table

spot.size = 8 #assigning spot size

conversion.factor = spot.size/em_d4G1@images$slice1.008um@scale.factors$spot #calculate conversion factor

spatial.factors = data.frame(ratio = conversion.factor, tol = spot.size/2) #creating spatial factor

d.spatial <- Compute_Cell_Distance(coordinates = spatial.locs, ratio = spatial.factors$ratio, tol = spatial.factors$tol) #creating distance matrix

min(d.spatial[d.spatial!=0]) # this value should approximately equal 8um for 10X VisiumHD data
#> [1] 7.999177

df=as.data.frame(as.matrix(d.spatial)) #convert distance matrix to dataframe

head(df[,1:6]) #look at the distance between the top barcodes
#>                       s_008um_00110_00228-1 s_008um_00122_00239-1
#> s_008um_00110_00228-1                0.0000             130.22140
#> s_008um_00122_00239-1              130.2214               0.00000
#> s_008um_00125_00234-1              129.2325              46.64416
#> s_008um_00124_00223-1              118.9164             128.98896
#> s_008um_00131_00230-1              168.7433             101.81457
#> s_008um_00133_00224-1              186.7426             148.79694
#>                       s_008um_00125_00234-1 s_008um_00124_00223-1
#> s_008um_00110_00228-1             129.23247             118.91640
#> s_008um_00122_00239-1              46.64416             128.98896
#> s_008um_00125_00234-1               0.00000              88.35838
#> s_008um_00124_00223-1              88.35838               0.00000
#> s_008um_00131_00230-1              57.68334              79.19059
#> s_008um_00133_00224-1             102.44176              72.43592
#>                       s_008um_00131_00230-1 s_008um_00133_00224-1
#> s_008um_00110_00228-1             168.74335             186.74257
#> s_008um_00122_00239-1             101.81457             148.79694
#> s_008um_00125_00234-1              57.68334             102.44176
#> s_008um_00124_00223-1              79.19059              72.43592
#> s_008um_00131_00230-1               0.00000              50.59329
#> s_008um_00133_00224-1              50.59329               0.00000

Part VI: Save the SpatialNiche files

saveRDS(DataHD, 'DataHD.rds') #save SpatialNiche list containing barcodes and images tibble

saveRDS(PeripheryBCs, 'PeripheryBCs.rds') #save barcodes in the cell specific niche

saveRDS(PeripheryGCs, 'PeripheryGCs.rds') #save barcodes in the gene specific niche
sessionInfo()
#> R version 4.4.0 (2024-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 22631)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] Seurat_5.2.1       SeuratObject_5.0.2 sp_2.2-0           scattermore_1.2   
#> [5] dplyr_1.1.4        ggplot2_3.5.2      SpatialNiche_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>   [1] deldir_2.0-4           pbapply_1.7-2          gridExtra_2.3         
#>   [4] rlang_1.1.5            magrittr_2.0.3         RcppAnnoy_0.0.22      
#>   [7] spatstat.geom_3.3-6    matrixStats_1.5.0      ggridges_0.5.6        
#>  [10] compiler_4.4.0         reshape2_1.4.4         png_0.1-8             
#>  [13] vctrs_0.6.5            stringr_1.5.1          pkgconfig_2.0.3       
#>  [16] fastmap_1.2.0          labeling_0.4.3         utf8_1.2.6            
#>  [19] promises_1.3.2         rmarkdown_2.29         tzdb_0.5.0            
#>  [22] purrr_1.0.4            bit_4.6.0              xfun_0.52             
#>  [25] cachem_1.1.0           jsonlite_2.0.0         goftest_1.2-3         
#>  [28] later_1.4.1            spatstat.utils_3.1-3   jpeg_0.1-11           
#>  [31] tiff_0.1-12            irlba_2.3.5.1          parallel_4.4.0        
#>  [34] cluster_2.1.8.1        R6_2.6.1               ica_1.0-3             
#>  [37] spatstat.data_3.1-6    stringi_1.8.7          bslib_0.9.0           
#>  [40] RColorBrewer_1.1-3     reticulate_1.42.0      spatstat.univar_3.1-2 
#>  [43] parallelly_1.43.0      lmtest_0.9-40          jquerylib_0.1.4       
#>  [46] Rcpp_1.0.14            assertthat_0.2.1       knitr_1.50            
#>  [49] tensor_1.5             future.apply_1.11.3    zoo_1.8-13            
#>  [52] sctransform_0.4.1      httpuv_1.6.15          Matrix_1.7-3          
#>  [55] splines_4.4.0          igraph_2.1.4           tidyselect_1.2.1      
#>  [58] abind_1.4-8            rstudioapi_0.17.1      dichromat_2.0-0.1     
#>  [61] yaml_2.3.10            spatstat.random_3.3-3  spatstat.explore_3.4-2
#>  [64] codetools_0.2-20       miniUI_0.1.2           listenv_0.9.1         
#>  [67] plyr_1.8.9             lattice_0.22-7         tibble_3.2.1          
#>  [70] shiny_1.10.0           withr_3.0.2            ROCR_1.0-11           
#>  [73] evaluate_1.0.3         Rtsne_0.17             future_1.40.0         
#>  [76] fastDummies_1.7.5      survival_3.8-3         polyclip_1.10-7       
#>  [79] fitdistrplus_1.2-2     pillar_1.10.2          KernSmooth_2.23-26    
#>  [82] plotly_4.10.4          generics_0.1.4         RcppHNSW_0.6.0        
#>  [85] scales_1.4.0           globals_0.17.0         xtable_1.8-4          
#>  [88] glue_1.8.0             lazyeval_0.2.2         tools_4.4.0           
#>  [91] data.table_1.17.0      RSpectra_0.16-2        RANN_2.6.2            
#>  [94] bmp_0.3                dotCall64_1.2          cowplot_1.1.3         
#>  [97] grid_4.4.0             tidyr_1.3.1            colorspace_2.1-1      
#> [100] readbitmap_0.1.5       nlme_3.1-168           patchwork_1.3.0       
#> [103] cli_3.6.4              spatstat.sparse_3.1-0  spam_2.11-1           
#> [106] viridisLite_0.4.2      arrow_19.0.1           uwot_0.2.3            
#> [109] gtable_0.3.6           sass_0.4.10            digest_0.6.37         
#> [112] progressr_0.15.1       ggrepel_0.9.6          rjson_0.2.23          
#> [115] htmlwidgets_1.6.4      farver_2.1.2           htmltools_0.5.8.1     
#> [118] lifecycle_1.0.4        httr_1.4.7             mime_0.13             
#> [121] bit64_4.6.0-1          MASS_7.3-65