SPIAT worksheet

Welcome

Welcome to the SPIAT spatial analysis workshop! In this worksheet we will explore how to analyse spatial proteomics data using SPIAT and simulate spatial point pattern using spaSim. Let’s get started!

This worksheet was adapted from the SPIAT tutorials, but uses a real dataset instead. Please refer to the tutorials for detailed information.

Please make sure you have installed all the required packages prior to the workshop. We won’t have time for installation during the workshop.

## Install the packages: 
## (Only run the code if the packages were not installed yet)
## (latest verion of R required >=4.3.0)

# install.packages(c("elsa", "Rtsne", "umap", "alphahull", "plotly", "survminer", "survival", "vroom"))
# install.packages(c("rmarkdown", "knitr"))
# 
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install(version = "3.18")
# 
# BiocManager::install("SPIAT")
# BiocManager::install("spaSim")
# BiocManager::install("ComplexHeatmap")

# load the libraries
library(SPIAT)
library(spaSim)

In this workshop we will demonstrate the analysis workflow with a public TNBC dataset (Keren et al., 2018).

Reading and formatting

Reading and formatting

format_image_to_spe() is the main function to read in data to SPIAT. format_image_to_spe() creates a SpatialExperiment object which is used in all subsequent functions. The key data points of interest for SPIAT are cell coordinates, marker intensities and cell phenotypes for each cell.

format_image_to_spe() has several options designed specifically for importing data generated from the inForm, HALO, CODEX and cellprofiler platforms. However, we advise using the more flexible ‘general’ option. Note that before using ‘general’, we advise pre-formatting your data.

# use `?` to get the usage of a function in the "help" panel
?format_image_to_spe

Now we’ll use the general option to read in the TNBC MIBI data. This MIBI file is not the raw MIBI data, but was generated after cell segmentation (by Greg Bass, CSL). This shows SPIAT’s capability of reading data from a wide range of formats.

# read the csv
path_to_file <- "cellData_with-centroids.csv" # put the path to the csv in the quote
mibi <- vroom::vroom(path_to_file)

## Rows: 201656 Columns: 59
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (59): SampleID, cellLabelInImage, Xc, Yc, cellSize, C, Na, Si, P, Ca, Fe...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mibi[1:5,] # look at the first 5 rows

## # A tibble: 5 × 59
##   SampleID cellLabelInImage    Xc    Yc cellSize     C     Na      Si      P
##      <dbl>            <dbl> <dbl> <dbl>    <dbl> <dbl>  <dbl>   <dbl>  <dbl>
## 1        1                2    34   159      146     0 -0.591  0.875  -2.58 
## 2        1                3    32   192      102     0 -0.499  0.0175 -1.22 
## 3        1                4    31   213       43     0 -1.49  -0.630  -1.91 
## 4        1                5    36   271      211     0 -1.01  -0.532  -1.74 
## 5        1                6    36   382      177     0  0.158 -0.710   0.517
## # ℹ 50 more variables: Ca <dbl>, Fe <dbl>, dsDNA <dbl>, Vimentin <dbl>,
## #   SMA <dbl>, Background <dbl>, B7H3 <dbl>, FoxP3 <dbl>, Lag3 <dbl>,
## #   CD4 <dbl>, CD16 <dbl>, CD56 <dbl>, OX40 <dbl>, PD1 <dbl>, CD31 <dbl>,
## #   PD.L1 <dbl>, EGFR <dbl>, Ki67 <dbl>, CD209 <dbl>, CD11c <dbl>, CD138 <dbl>,
## #   CD163 <dbl>, CD68 <dbl>, CSF.1R <dbl>, CD8 <dbl>, CD3 <dbl>, IDO <dbl>,
## #   Keratin17 <dbl>, CD63 <dbl>, CD45RO <dbl>, CD20 <dbl>, p53 <dbl>,
## #   Beta.catenin <dbl>, HLA.DR <dbl>, CD11b <dbl>, CD45 <dbl>, H3K9ac <dbl>, …

unique(mibi$SampleID) # how many unique samples here?

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 31 32 33 34 35 36 37 38 39 40 41 42 43 44

# slice the data frame to keep only the sample with `sampleID==4`
mibi_4 <- mibi[mibi$SampleID == 4, ]
dim(mibi_4)

## [1] 6643   59

# add a column of cell IDs
mibi_4$CellID <- paste("Cell", mibi_4$cellLabelInImage, sep="_")

# construct the elements of the SPE object
## the marker intensity matrix
intensity <- t(mibi_4[,13:52])
intensity <- rbind(intensity, Tumor = mibi_4$tumorYN)
intensity <- intensity[rownames(intensity) %in% c("CD4", "CD3", "CD8",
                                                "CD20", "MPO", "CD68",
                                                "Tumor"),]
colnames(intensity) <- mibi_4$CellID # add cell IDs
intensity[, 1:5] # print the intensity of first 5 cells

##         Cell_2   Cell_3   Cell_4   Cell_5   Cell_6
## CD4   -0.37537 -0.37537 -0.37537 -0.37537 -0.37537
## CD68  -0.37415 -0.37415 -0.37415 -0.37415 -0.37415
## CD8   -0.37400 -0.37400 -0.37400 -0.37400 -0.37400
## CD3   -0.44579 -0.44579 -0.44579 -0.44579 -0.08331
## CD20  -0.35542 -0.35542 -0.35542 -0.35542 -0.35542
## MPO   -0.13903 -0.13903 -0.13903 -0.13903 -0.13903
## Tumor  1.00000  1.00000  1.00000  1.00000  0.00000

## the coordinates
coord_x <- mibi_4$Xc
coord_y <- mibi_4$Yc
## we don't have phenotype info in the file. Will predict the phenotypes later

# now format the image with function `format_image_to_spe()`
mibi_4_spe <- format_image_to_spe(format = "general", phenotypes = NA,
                                  intensity_matrix = intensity,
                                  coord_x = coord_x, coord_y = coord_y)

Now we get the SPE object of the image with sampleID==4.

Predicting phenotypes

You’ll notice there is no phenotype info in the MIBI .csv file. Let’s use SPIAT to predict phenotypes based on the marker intensity distribution of all cells.

# Print the markers to predict
rownames(assay(mibi_4_spe))

## [1] "CD4"   "CD68"  "CD8"   "CD3"   "CD20"  "MPO"   "Tumor"

# Predict (`reference_phenotypes = FALSE` because there are no default phenotypes)
mibi_4_spe <- predict_phenotypes(mibi_4_spe, tumour_marker = "Tumor",
                   baseline_markers = c("CD4", "CD68", "CD8", "CD3", "CD20", "MPO"),
                   reference_phenotypes = FALSE)

## [1] "CD4  threshold intensity:  0.663931806123381"

## [1] "CD68  threshold intensity:  0.382038332923384"

## [1] "CD8  threshold intensity:  0.651638852799246"

## [1] "CD3  threshold intensity:  0.489898334410501"

## [1] "CD20  threshold intensity:  -0.302089348538088"

## [1] "MPO  threshold intensity:  0.667462603716241"

## [1] "Tumor  threshold intensity:  0.501418155086776"

Have a look at the structure of this object. There are three main measurements of interest (coordinates, phenotype, intensities).

# View(mibi_4_spe)
# Access phenotype as a column (accessed with $)
# Print the unique phenotypes
unique(mibi_4_spe$Phenotype)

##  [1] "Tumor"                 "None"                  "CD3"                  
##  [4] "MPO"                   "CD20"                  "CD68"                 
##  [7] "CD68,Tumor"            "CD20,Tumor"            "CD3,Tumor"            
## [10] "CD8,CD20"              "CD68,MPO"              "CD4,CD3,CD20"         
## [13] "CD8,CD3"               "CD8"                   "CD4"                  
## [16] "MPO,Tumor"             "CD8,Tumor"             "CD4,CD3"              
## [19] "CD4,CD68,CD3"          "CD68,CD20"             "CD4,CD68"             
## [22] "CD68,CD8"              "CD3,CD20"              "CD4,CD68,CD3,CD20"    
## [25] "CD4,CD68,CD20"         "CD68,CD3"              "CD8,CD3,CD20"         
## [28] "CD4,CD20"              "CD4,CD68,CD8,CD3"      "CD4,CD8,CD3"          
## [31] "CD4,CD68,CD20,MPO"     "CD4,CD8,CD3,CD20"      "CD68,CD20,MPO"        
## [34] "CD68,CD8,CD3,CD20"     "CD4,CD68,CD8,CD3,CD20" "CD68,CD8,CD3"         
## [37] "CD3,MPO"               "CD4,Tumor"             "CD4,CD3,MPO"          
## [40] "CD4,CD68,CD8,MPO"      "CD20,MPO"              "CD4,CD68,CD3,MPO"     
## [43] "CD4,CD3,Tumor"         "CD4,MPO"               "CD3,CD20,MPO"         
## [46] "CD8,CD3,CD20,MPO"      "CD4,CD68,CD8,CD20"     "CD4,CD3,CD20,MPO"     
## [49] "CD4,CD8,CD3,MPO"       "CD4,CD8"               "CD4,CD8,CD20"         
## [52] "CD68,CD8,CD3,CD20,MPO" "CD68,CD3,CD20,MPO"     "CD68,CD8,CD20"        
## [55] "CD4,CD68,MPO"          "CD68,CD3,CD20"         "CD68,CD8,CD3,MPO"     
## [58] "CD4,CD68,CD8"          "CD8,CD3,Tumor"         "CD68,CD8,MPO"         
## [61] "CD68,CD3,MPO"          "CD4,CD68,CD3,CD20,MPO" "CD8,MPO"

# What is the dimension of the intensity matrix of this object?
dim(assay(mibi_4_spe))

## [1]    7 6643

# Access marker intensities with the helper function 'assay'
# Print the marker intensities of the first 5 cells
assay(mibi_4_spe)[, 1:5]

##         Cell_2   Cell_3   Cell_4   Cell_5   Cell_6
## CD4   -0.37537 -0.37537 -0.37537 -0.37537 -0.37537
## CD68  -0.37415 -0.37415 -0.37415 -0.37415 -0.37415
## CD8   -0.37400 -0.37400 -0.37400 -0.37400 -0.37400
## CD3   -0.44579 -0.44579 -0.44579 -0.44579 -0.08331
## CD20  -0.35542 -0.35542 -0.35542 -0.35542 -0.35542
## MPO   -0.13903 -0.13903 -0.13903 -0.13903 -0.13903
## Tumor  1.00000  1.00000  1.00000  1.00000  0.00000

# Access coordinates with the helper function 'spatialCoords'
spatialCoords(mibi_4_spe)[1:5, ]

##      Cell.X.Position Cell.Y.Position
## [1,]              40              36
## [2,]              37              62
## [3,]              37              87
## [4,]              52             111
## [5,]              44             154

Define cell types

Then we define the cell types.

phenotypes <- c("Tumor", "CD68", "CD20", "CD4,CD3", "CD3", "CD8,CD3", 
                "CD4", "MPO", "CD8")
cell_types <- c("Tumor", "Macrophage", "B_cells",  "Helper_T_cell", "T_cell",
                "Cyto_T_cell",  "Helper_T_cell", "Neutrophil", "Cyto_T_cell")

defined_mibi_4 <- define_celltypes(
    mibi_4_spe, 
    categories = phenotypes, 
    category_colname = "Phenotype", 
    names = cell_types,
    new_colname = "Cell.Type")
unique(defined_mibi_4$Cell.Type)

## [1] "Tumor"         "Undefined"     "T_cell"        "Neutrophil"   
## [5] "B_cells"       "Macrophage"    "Cyto_T_cell"   "Helper_T_cell"

You can write a for loop to read all of the patients’ images. We will look into that later.

Quality control and data visualisation

We’ll now use SPIAT’s quality control functions to check phenotyping quality, detect uneven staining, and test for other potential technical artefacts.

Visualise marker levels

Phenotyping of cells can be verified comparing marker intensities of cells labelled positive and negative for a marker. Cells positive for a marker should have high levels of the marker. An unclear separation of marker intensities between positive and negative cells would suggest phenotypes have not been accurately assigned. We can use marker_intensity_boxplot() to produce a boxplot for cells phenotyped as being positive or negative for a marker.

markers <- rownames(assay(mibi_4_spe))
for (marker in markers){
  print(marker_intensity_boxplot(mibi_4_spe, marker))
}

Uneven marker staining or high background intensity can be identified with plot_cell_marker_levels(). This produces a scatter plot of the intensity of a marker in each cell. This should be relatively even across the image and all phenotyped cells. Cells that were not phenotyped as being positive for the particular marker are excluded.

plot_cell_marker_levels(mibi_4_spe, "CD3")

# Try plot all markers with a for loop~
for (marker in markers){
  print(plot_cell_marker_levels(mibi_4_spe, marker))
}

## Warning in plot_cell_marker_levels(mibi_4_spe, marker): NaNs produced

For large images, there is also the option of ‘blurring’ the image with function plot_marker_level_heatmap(), where the image is split into multiple small areas, and marker intensities are averaged within each. The image is blurred based on the num_splits parameter.

# Try it yourself~
for (marker in markers){
  print(plot_marker_level_heatmap(mibi_4_spe, num_splits = 100, marker = marker))
}

Identifying incorrect phenotypes

We may see cells with biologically implausible combination of markers present in the input data when using:

unique(mibi_4_spe$Phenotype)

##  [1] "Tumor"                 "None"                  "CD3"                  
##  [4] "MPO"                   "CD20"                  "CD68"                 
##  [7] "CD68,Tumor"            "CD20,Tumor"            "CD3,Tumor"            
## [10] "CD8,CD20"              "CD68,MPO"              "CD4,CD3,CD20"         
## [13] "CD8,CD3"               "CD8"                   "CD4"                  
## [16] "MPO,Tumor"             "CD8,Tumor"             "CD4,CD3"              
## [19] "CD4,CD68,CD3"          "CD68,CD20"             "CD4,CD68"             
## [22] "CD68,CD8"              "CD3,CD20"              "CD4,CD68,CD3,CD20"    
## [25] "CD4,CD68,CD20"         "CD68,CD3"              "CD8,CD3,CD20"         
## [28] "CD4,CD20"              "CD4,CD68,CD8,CD3"      "CD4,CD8,CD3"          
## [31] "CD4,CD68,CD20,MPO"     "CD4,CD8,CD3,CD20"      "CD68,CD20,MPO"        
## [34] "CD68,CD8,CD3,CD20"     "CD4,CD68,CD8,CD3,CD20" "CD68,CD8,CD3"         
## [37] "CD3,MPO"               "CD4,Tumor"             "CD4,CD3,MPO"          
## [40] "CD4,CD68,CD8,MPO"      "CD20,MPO"              "CD4,CD68,CD3,MPO"     
## [43] "CD4,CD3,Tumor"         "CD4,MPO"               "CD3,CD20,MPO"         
## [46] "CD8,CD3,CD20,MPO"      "CD4,CD68,CD8,CD20"     "CD4,CD3,CD20,MPO"     
## [49] "CD4,CD8,CD3,MPO"       "CD4,CD8"               "CD4,CD8,CD20"         
## [52] "CD68,CD8,CD3,CD20,MPO" "CD68,CD3,CD20,MPO"     "CD68,CD8,CD20"        
## [55] "CD4,CD68,MPO"          "CD68,CD3,CD20"         "CD68,CD8,CD3,MPO"     
## [58] "CD4,CD68,CD8"          "CD8,CD3,Tumor"         "CD68,CD8,MPO"         
## [61] "CD68,CD3,MPO"          "CD4,CD68,CD3,CD20,MPO" "CD8,MPO"

Removing cells with incorrect phenotypes: If you identify incorrect phenotypes or have any properties (columns) that you want to exclude you can use select_celltypes(). Set celltypes as the values that you want to keep or exclude for a column (feature_colname). Set keep as TRUE to include these cells and FALSE to exclude them.

data_subset <- select_celltypes(
  mibi_4_spe, keep=TRUE,
  celltypes = c("Tumor", "CD68", "CD20", "CD4,CD3", 
                 "CD3", "CD8,CD3", "CD4", "MPO", "CD8"),
  feature_colname = "Phenotype")
# have a look at what phenotypes are present
unique(data_subset$Phenotype)

## [1] "Tumor"   "CD3"     "MPO"     "CD20"    "CD68"    "CD8,CD3" "CD8"    
## [8] "CD4"     "CD4,CD3"

# All the cells with the rest of the combinations of phenotypes were removed.
# define cell types
defined_mibi_subset <- define_celltypes(
    data_subset, 
    categories = phenotypes, 
    category_colname = "Phenotype", 
    names = cell_types,
    new_colname = "Cell.Type")

Dimensionality reduction to identify misclassified cells

We can also check for specific misclassified cells using dimensionality reduction. SPIAT offers tSNE and UMAPs based on marker intensities to visualise cells. Cells of distinct types should form clearly different clusters.

The generated dimensionality reduction plots are interactive, and users can hover over each cell to obtain the cell ID. Users can then remove specific misclassified cells.

g <- dimensionality_reduction_plot(defined_mibi_subset, plot_type = "TSNE", 
                                   feature_colname = "Cell.Type")

Note that dimensionality_reduction_plot() only prints a static version of the UMAP or tSNE plot. To interact with this plot, pass the result to the ggplotly() function from the plotly package.

plotly::ggplotly(g)

Visualising tissues

Let’s go back to the image without filtering out the undefined cells - mibi_4_spe.

We can see the location of all cell types (or any column in the data) in the tissue with plot_cell_categories(). Each dot in the plot corresponds to a cell and cells are coloured by cell type. Any cell types present in the data but not in the cell types of interest will be put in the category “OTHER” and coloured lightgrey.

my_colors <- c("red", "blue","orange", "purple", "pink", "darkcyan", "darkgreen")
categories_of_interest <- c("Tumor", "Macrophage", "B_cells",  "Helper_T_cell", 
                            "T_cell", "Cyto_T_cell", "Neutrophil")
  
plot_cell_categories(spe_object = defined_mibi_4, 
                     categories_of_interest = categories_of_interest, 
                     colour_vector = my_colors, feature_colname = "Cell.Type")

We can visualise a selected marker in 3D with marker_surface_plot(). The image is blurred based on the num_splits parameter.

marker_surface_plot(defined_mibi_4, num_splits=15, marker="CD68")

To visualise multiple markers in 3D in a single plot we can use marker_surface_plot_stack(). This shows normalised intensity level of specified markers and enables the identification of co-occurring and mutually exclusive markers.

marker_surface_plot_stack(defined_mibi_4, num_splits=15, markers_to_plot=c("CD68", "CD20"))

Basic analyses

Cell percentages

Obtain the number and proportion of each cell type. We can exclude any cell types that are not of interest e.g. “Undefined” with celltypes_to_exclude.

p_cells <- calculate_cell_proportions(defined_mibi_4, 
                                      feature_colname ="Cell.Type",
                                      celltypes_to_exclude = c("Others","Undefined"),
                                      plot.image = TRUE)

plot_cell_percentages(cell_proportions = p_cells, 
                      cells_to_exclude = "Tumor", cellprop_colname="Proportion_name")

Cell distances

We can calculate the pairwise distances between two cell types (cell type A and cell type B) with calculate_pairwise_distances_between_cell_types(). This function calculates the distances of all cells of type A against all cells of type B.

distances <- calculate_pairwise_distances_between_celltypes(
  spe_object = defined_mibi_4, 
  cell_types_of_interest = c("Tumor", "B_cells", "T_cell"),
  feature_colname = "Cell.Type")

# Visualise the distances
plot_cell_distances_violin(distances)

We can also calculate summary statistics for the distances between each combination of cell types, the mean, median, min, max and standard deviation, with calculate_summary_distances_between_celltypes().

summary_distances <- calculate_summary_distances_between_celltypes(distances)
summary_distances

##              Pair      Mean       Min      Max    Median  Std.Dev Reference
## 1 B_cells/B_cells  921.4534 12.649111 2482.944  903.4517 486.0542   B_cells
## 2  B_cells/T_cell  929.3509  8.062258 2647.739  892.9138 470.1658   B_cells
## 3   B_cells/Tumor 1107.9081 13.038405 2776.874 1106.3453 460.9599   B_cells
## 4  T_cell/B_cells  929.3509  8.062258 2647.739  892.9138 470.1658    T_cell
## 5   T_cell/T_cell  927.6880  5.830952 2576.001  897.8469 467.7482    T_cell
## 6    T_cell/Tumor 1085.7099 12.041595 2706.849 1087.3452 473.7943    T_cell
## 7   Tumor/B_cells 1107.9081 13.038405 2776.874 1106.3453 460.9599     Tumor
## 8    Tumor/T_cell 1085.7099 12.041595 2706.849 1087.3452 473.7943     Tumor
## 9     Tumor/Tumor 1079.8868  6.082763 2735.171 1092.9318 558.4303     Tumor
##    Target
## 1 B_cells
## 2  T_cell
## 3   Tumor
## 4 B_cells
## 5  T_cell
## 6   Tumor
## 7 B_cells
## 8  T_cell
## 9   Tumor

An example of the interpretation of this result is: “average pairwise distance between cells of Tumor and B cells is 1107.9081.

These pairwise cell distances can then be visualised as a heatmap with plot_distance_heatmap(). This example shows the average pairwise distances between cell types. Note that the pairwise distances are symmetrical (the average distance between cell type A and cell type B is the same as the average distance between cell Type B and cell Type A).

plot_distance_heatmap(phenotype_distances_result = summary_distances, metric = "mean")

The same workflow also applies to minimum distance between cell types. Give it a go with the following functions: calculate_minimum_distances_between_celltypes(), plot_cell_distances_violin(), calculate_summary_distances_between_celltypes(), plot_distance_heatmap().

distances <- calculate_minimum_distances_between_celltypes(
  spe_object = defined_mibi_4, 
  cell_types_of_interest = c("Tumor", "B_cells", "T_cell"),
  feature_colname = "Cell.Type")

## [1] "Markers had been selected in minimum distance calculation: "
## [1] "Tumor"   "T_cell"  "B_cells"

plot_cell_distances_violin(distances)

summary_distances <- calculate_summary_distances_between_celltypes(distances)
summary_distances

##             Pair      Mean       Min      Max    Median   Std.Dev Reference
## 1 B_cells/T_cell  68.68646  8.062258 189.8420  61.20457  35.77786   B_cells
## 2  B_cells/Tumor  69.06331 13.038405 285.3507  56.58622  48.82640   B_cells
## 3 T_cell/B_cells  57.30047  8.062258 662.5300  40.92617  69.56874    T_cell
## 4   T_cell/Tumor  57.00635 12.041595 247.1295  41.69221  42.70335    T_cell
## 5  Tumor/B_cells 196.84322 13.038405 769.3458 165.60495 139.10386     Tumor
## 6   Tumor/T_cell 159.54208 12.041595 467.8996 134.12308 102.09552     Tumor
##    Target
## 1  T_cell
## 2   Tumor
## 3 B_cells
## 4   Tumor
## 5 B_cells
## 6  T_cell

plot_distance_heatmap(phenotype_distances_result = summary_distances, metric = "mean")

Cell colocalisation

With SPIAT we can quantify cell colocalisation, which refers to how much two cell types are colocalising and thus potentially interacting.

Cells in Neighbourhood (CIN)

We can calculate the average percentage of cells of one cell type (target) within a radius of another cell type (reference) using average_percentage_of_cells_within_radius().

average_percentage_of_cells_within_radius(spe_object = defined_mibi_4, 
                                          reference_celltype = "B_cells", 
                                          target_celltype = "Macrophage", 
                                          radius=100, feature_colname="Cell.Type")

## [1] 9.434629

Alternatively, this analysis can also be performed based on marker intensities rather than cell types. Here, we use average_marker_intensity_within_radius() to calculate the average intensity of the target_marker within a radius from the cells positive for the reference marker. Note that it pools all cells with the target marker that are within the specific radius of any reference cell. Results represent the average intensities within a radius.

average_marker_intensity_within_radius(spe_object = defined_mibi_4,
                                       reference_marker ="CD20",
                                       target_marker = "CD68",
                                       radius=100)

## [1] 1.57434

plot_average_intensity() calculates the average intensity of a target marker for a number of user-supplied radii values, and plots the intensity level at each specified radius as a line graph.

plot_average_intensity(spe_object=defined_mibi_4, 
                       reference_marker="CD20", 
                       target_marker="CD68", radii=seq(20, 100, 5))

What can you interpret from this plot?

Mixing Score (MS) and Normalised Mixing Score (NMS)

Mixing score (MS) in SPIAT is defined as the number of target-reference interactions divided by the number of reference-reference interactions within a specified radius. The higher the score the greater the mixing of the two cell types. The normalised score is normalised for the number of target and reference cells in the image.

# What is the mixing score between Tumor cells and Macrophage in our sample?
mixing_score_summary(spe_object = defined_mibi_4, 
                     reference_celltype = "Tumor", 
                     target_celltype = "Macrophage", 
                     radius=20, feature_colname ="Cell.Type")

##   Reference     Target Number_of_reference_cells Number_of_target_cells
## 2     Tumor Macrophage                      2373                    585
##   Reference_target_interaction Reference_reference_interaction Mixing_score
## 2                           38                             893   0.04255319
##   Normalised_mixing_score
## 2              0.08627023

Cross K function

Cross K function calculates the number of target cell types across a range of radii from a reference cell type, and compares the behaviour of the input image with an image of randomly distributed points using a Poisson point process. There are four patterns that can be distinguished from K-cross function - independent, mixed clusters, mixed one cluster, separate clusters.

Try different cell types with different radii using calculate_cross_functions() . What patterns do you identify?

df_cross <- calculate_cross_functions(defined_mibi_4, method = "Kcross", 
                                      cell_types_of_interest = c("Tumor", "B_cells"), 
                                      feature_colname ="Cell.Type",
                                      dist = 500)

We can calculate the area under the curve (AUC) of the cross K-function. In general, this tells us the two types of cells are:

• negative values: separate clusters

• positive values: mixing of cell types

AUC_of_cross_function(df_cross)

## [1] -0.07022154

Cross-K Intersection (CKI)

There is another pattern in cross K curve which has not been previously appreciated, which is when there is a “ring” of one cell type, generally immune cells, surrounding the area of another cell type, generally tumour cells. For this pattern, the observed and expected curves in cross K function cross or intersect, such as the cross K plot above.

While this pattern does not appear in the example data, you can experiment with this pattern using the simulated data from the SPIAT package, or refer to the CKI section in the tutorial.

Aggregated entropy gradient

The function entropy_gradient_aggregated() proposes the aggregated entropy gradient as a self-contained metric to define the attraction and repulsion of cell types. You don’t have to set arbitrary threshold, but only use the shape of the curve to determine the colocalisation. The details of this algorithm were discussed thoroughly in Fig. 4 and Supplementary Fig. 4 of the paper.

Let’s have a look what results this function gives us.

table(defined_mibi_4$Cell.Type)

## 
##       B_cells   Cyto_T_cell Helper_T_cell    Macrophage    Neutrophil 
##           561           168           438           585            83 
##        T_cell         Tumor     Undefined 
##           182          2373          2253

gradient_pos <- c(50, 75, 100, 125, 150, 175, 200, 250, 300, 350, 400, 450, 500, 550, 600)
grad <- entropy_gradient_aggregated(defined_mibi_4, 
                                    cell_types_of_interest = c("Macrophage", "Tumor"),
                                    feature_colname = "Cell.Type",
                                    radii =  gradient_pos)

grad

## $gradient_df
##    Celltype1  Celltype2    Pos_50    Pos_75   Pos_100   Pos_125   Pos_150
## 1 Macrophage      Tumor 0.7712019 0.8864590 0.9345247 0.9635185 0.9824104
## 2      Tumor Macrophage 0.1780731 0.2114039 0.2322834 0.2534902 0.2739196
##     Pos_175   Pos_200   Pos_250   Pos_300   Pos_350   Pos_400   Pos_450
## 1 0.9920997 0.9971467 0.9999434 0.9967434 0.9913417 0.9831932 0.9736309
## 2 0.2904796 0.3070304 0.3398882 0.3717076 0.4017652 0.4341359 0.4666640
##    Pos_500   Pos_550   Pos_600
## 1 0.962555 0.9502988 0.9354429
## 2 0.498329 0.5292052 0.5602200
## 
## $peak
## [1] 8

Plot the gradients.

v <- as.numeric(grad$gradient_df[1,3:(length(gradient_pos)+2)])
plot(v, type = "b", lty = 2, pch = 16, cex = 1)

The trend indicates there is repulsion between tumor and macrophages. Let’s re-visit the cell category plot of the image and see the locations of these cells.

plot_cell_categories(spe_object = defined_mibi_4, 
                     categories_of_interest = categories_of_interest, 
                     colour_vector = my_colors, feature_colname = "Cell.Type")

Note that the shape of the curve is determined by several factors, including using which cell type as the reference cells. The discussion of this is included in the Supplementary Fig. 4 of the paper.

Spatial Heterogeneity

Cell colocalisation metrics allow capturing a dominant spatial pattern in an image. However, patterns are unlikely to be distributed evenly in a tissue, but rather there will be spatial heterogeneity of patterns. To measure this, SPIAT splits the image into smaller images (either using a grid or concentric circles around a reference cell population), followed by calculation of a spatial metric of a pattern of interest (e.g. cell colocalisation, entropy), and then measures the Prevalence and Distinctiveness of the pattern.

Entropy

Entropy in spatial analysis refers to the balance in the number of cells of distinct populations. An entropy score can be obtained for an entire image. However, the entropy of one image does not provide us spatial information of the image.

calculate_entropy(defined_mibi_4, cell_types_of_interest = c("Tumor","Macrophage"), 
                  feature_colname = "Cell.Type")

## [1] 0.7174431

We therefore propose the concept of Localised Entropy which calculates an entropy score for a predefined local region.

Fishnet grid

One approach to calculate localised metric is to split the image into fishnet grid squares. For each grid square, grid_metrics() calculates the metric for that square and visualise the raster image. You can choose any metric as the localised metric. Here we use entropy as an example.

For cases where the localised metric is not symmetrical (requires specifying a target and reference cell type), the first item in the vector used for cell_types_of_interest marks the reference cells and the second item the target cells. Here we are using Entropy, which is symmetrical, so we can use any order of cell types in the input.

grid <- grid_metrics(defined_mibi_4, FUN = calculate_entropy, n_split = 20,
                     cell_types_of_interest=c("Tumor","Macrophage"), 
                     feature_colname = "Cell.Type")

After calculating the localised entropy for each of the grid squares, we can apply metrics like percentages of grid squares with patterns (Prevalence) and Moran’s I (Distinctiveness).

For the Prevalence, we need to select a threshold over which grid squares are considered ‘positive’ for the pattern. The selection of threshold depends on the pattern and metric the user chooses to find the localised pattern. Here we chose 0.73 for entropy because 0.73 is roughly the entropy of two cell types when their ratio is 1:4 or 4:1.

You might be curious of how the entropy is calculated based on the ratio. Use the following helper function to calculate the entropy for any composition of cell counts.

entropy_helper <- function(v){
   s <- sum(v); l <- length(v); e <- 0
   for (i in 1:l) e <- e + -(v[i]/s)*log2(v[i]/s) 
   return(e)
 }

What is the entropy when the ratio of two cell types is 1:4?

entropy_helper(c(1,4))

## [1] 0.7219281

calculate_percentage_of_grids(grid, threshold = 0.73, above = TRUE)

## [1] 21.5

calculate_spatial_autocorrelation(grid, metric = "globalmoran")

## [1] 0.2730985

Characterise tissue structure

In certain analysis the focus is to understand the spatial distribution of a certain type of cell populations relative to the tissue regions.

One example of this functionality is to characterise the immune population in tumour structures. The following analysis will focus on the tumour/immune example, including determining whether there is a clear tumour margin, automatically identifying the tumour margin, and finally quantifying the proportion of immune populations relative to the margin. However, these analyses can also be generalised to other tissue and cell types.

Determining whether there is a clear tumour margin

In some instances tumour cells are distributed in such a way that there are no clear tumour margins. While this can be derived intuitively in most cases, SPIAT offers a way of quantifying the ‘quality’ of the margin for downstream analyses. This is meant to be used to help flag images with relatively poor margins, and therefore we do not offer a cutoff value.

To determine if there is a clear tumour margin, SPIAT can calculate the ratio of tumour bordering cells to tumour total cells (R-BC). This ratio is high when there is a disproportional high number of tumour margin cells compared to internal tumour cells.

R_BC(defined_mibi_4, cell_type_of_interest = "Tumor", "Cell.Type")

## [1] 0.1508639

The result is 0.1508639. This low value means there are relatively low number of bordering cells compared to total tumour cells, meaning that this image has clear tumour margins.

Automatic identification of the tumour margin

We can identify margins with identify_bordering_cells(). This function leverages off the alpha hull method (Pateiro-Lopez, Rodriguez-Casal, and. 2019) from the alpha hull package. Here we use tumour cells (Tumour_marker) as the reference to identify the bordering cells but any cell type can be used.

formatted_border <- identify_bordering_cells(defined_mibi_4, 
                                             reference_cell = "Tumor", 
                                             feature_colname="Cell.Type",
                                             ahull_alpha = 30,
                                             n_to_exclude = 25)

# Try with different ahull_alpha values and n_to_exclude values

# Then get the number of tumour clusters
attr(formatted_border, "n_of_clusters")

## [1] 3

There are 3 tumour clusters in the image.

Classification of cells based on their locations relative to the margin

First, we calculate the distance of cells to the tumour margin.

formatted_distance <- calculate_distance_to_margin(formatted_border)

## [1] "Markers had been selected in minimum distance calculation: "
## [1] "Non-border" "Border"

Next, we classify cells based on their location. As a distance cutoff, we use a distance of 5 cells from the tumour margin. The function first calculates the average minimum distance between all pairs of nearest cells and then multiples this number by 5. Users can change the number of cell layers to increase/decrease the margin width or use another argument margin_dist.

names_of_immune_cells <- c("B_cells",  "Helper_T_cell", "T_cell",
                "Cyto_T_cell",  "Helper_T_cell", "Neutrophil", "Cyto_T_cell")

formatted_structure <- define_structure(
  formatted_distance, cell_types_of_interest = names_of_immune_cells, 
  feature_colname = "Cell.Type", margin_dist = 90)

categories <- unique(formatted_structure$Structure)

We can plot and colour these structure categories.

plot_cell_categories(formatted_structure, feature_colname = "Structure")

Then calculate the proportions of immune cells in each of the locations.

immune_proportions <- calculate_proportions_of_cells_in_structure(
  spe_object = formatted_structure, 
  cell_types_of_interest = names_of_immune_cells, feature_colname ="Cell.Type")

immune_proportions

##                Cell.Type                            Relative_to
## 1                B_cells             All_cells_in_the_structure
## 2          Helper_T_cell             All_cells_in_the_structure
## 3                 T_cell             All_cells_in_the_structure
## 4            Cyto_T_cell             All_cells_in_the_structure
## 5          Helper_T_cell             All_cells_in_the_structure
## 6             Neutrophil             All_cells_in_the_structure
## 7            Cyto_T_cell             All_cells_in_the_structure
## 8                B_cells All_cells_of_interest_in_the_structure
## 9          Helper_T_cell All_cells_of_interest_in_the_structure
## 10                T_cell All_cells_of_interest_in_the_structure
## 11           Cyto_T_cell All_cells_of_interest_in_the_structure
## 12         Helper_T_cell All_cells_of_interest_in_the_structure
## 13            Neutrophil All_cells_of_interest_in_the_structure
## 14           Cyto_T_cell All_cells_of_interest_in_the_structure
## 15               B_cells  The_same_cell_type_in_the_whole_image
## 16         Helper_T_cell  The_same_cell_type_in_the_whole_image
## 17                T_cell  The_same_cell_type_in_the_whole_image
## 18           Cyto_T_cell  The_same_cell_type_in_the_whole_image
## 19         Helper_T_cell  The_same_cell_type_in_the_whole_image
## 20            Neutrophil  The_same_cell_type_in_the_whole_image
## 21           Cyto_T_cell  The_same_cell_type_in_the_whole_image
## 22 All_cells_of_interest             All_cells_in_the_structure
##    P.Infiltrated.CoI P.Internal.Margin.CoI P.External.Margin.CoI P.Stromal.CoI
## 1        0.000000000           0.001834862            0.09242619   0.132987439
## 2        0.006337136           0.003669725            0.16688062   0.081649372
## 3        0.006337136           0.004587156            0.04364570   0.037684326
## 4        0.021546261           0.003669725            0.05391528   0.028672856
## 5        0.006337136           0.003669725            0.16688062   0.081649372
## 6        0.027883397           0.005504587            0.02439024   0.009830694
## 7        0.021546261           0.003669725            0.05391528   0.028672856
## 8        0.000000000           0.095238095            0.24242424   0.457276995
## 9        0.102040816           0.190476190            0.43771044   0.280751174
## 10       0.102040816           0.238095238            0.11447811   0.129577465
## 11       0.346938776           0.190476190            0.14141414   0.098591549
## 12       0.102040816           0.190476190            0.43771044   0.280751174
## 13       0.448979592           0.285714286            0.06397306   0.033802817
## 14       0.346938776           0.190476190            0.14141414   0.098591549
## 15       0.000000000           0.003565062            0.12834225   0.868092692
## 16       0.011415525           0.009132420            0.29680365   0.682648402
## 17       0.027472527           0.027472527            0.18681319   0.758241758
## 18       0.101190476           0.023809524            0.25000000   0.625000000
## 19       0.011415525           0.009132420            0.29680365   0.682648402
## 20       0.265060241           0.072289157            0.22891566   0.433734940
## 21       0.101190476           0.023809524            0.25000000   0.625000000
## 22       0.062103929           0.019266055            0.38125802   0.290824686

Finally calculate summaries of the distances for immune cells in the tumour structure.

immune_distances <- calculate_summary_distances_of_cells_to_borders(
  spe_object = formatted_structure, 
  cell_types_of_interest = names_of_immune_cells, feature_colname = "Cell.Type")

immune_distances

##                     Cell.Type               Area    Min_d    Max_d    Mean_d
## 1  All_cell_types_of_interest Within_border_area 17.00000 248.8534 108.54431
## 2  All_cell_types_of_interest             Stroma 14.76482 685.5983 230.71498
## 3                     B_cells Within_border_area 71.47027  82.8734  77.17184
## 4                     B_cells             Stroma 15.81139 659.7287 259.69499
## 5               Helper_T_cell Within_border_area 17.00000 145.8801  86.04687
## 6               Helper_T_cell             Stroma 16.00000 647.0023 210.97708
## 7                      T_cell Within_border_area 20.80865 144.9000  83.31473
## 8                      T_cell             Stroma 18.86796 637.4080 221.76403
## 9                 Cyto_T_cell Within_border_area 28.30194 178.4657 108.36877
## 10                Cyto_T_cell             Stroma 14.76482 652.8599 206.74196
## 11              Helper_T_cell Within_border_area 17.00000 145.8801  86.04687
## 12              Helper_T_cell             Stroma 16.00000 647.0023 210.97708
## 13                 Neutrophil Within_border_area 22.67157 248.8534 127.15875
## 14                 Neutrophil             Stroma 16.12452 685.5983 182.19381
## 15                Cyto_T_cell Within_border_area 28.30194 178.4657 108.36877
## 16                Cyto_T_cell             Stroma 14.76482 652.8599 206.74196
##     Median_d   St.dev_d
## 1  107.06307  52.557389
## 2  205.70488 152.477111
## 3   77.17184   8.063226
## 4  246.41631 145.901123
## 5   91.35097  45.132337
## 6  158.00316 158.573639
## 7   87.41484  40.520997
## 8  197.03778 143.424998
## 9  108.74741  37.114966
## 10 167.17057 155.395117
## 11  91.35097  45.132337
## 12 158.00316 158.573639
## 13 111.28321  63.630089
## 14 165.19383 142.348933
## 15 108.74741  37.114966
## 16 167.17057 155.395117

Cellular neighbourhood

The aggregation of cells can result in cellular neighbourhoods. A neighbourhood is defined as a group of cells that cluster together. These can be homotypic, containing cells of a single class (e.g.immune cells), or heterotypic (e.g.a mixture of tumour and immune cells).

Function identify_neighborhoods() identifies cellular neighbourhoods. Users can select a subset of cell types of interest if desired. The algorithm is Hierarchical Clustering algorithm - Euclidean distances between cells are calculated, and pairs of cells with a distance less than a specified radius are considered to be interacting, with the rest being non-interacting. Hierarchical clustering is then used to separate the clusters. Larger radii will result in the merging of individual clusters.

You need to specify a radius that defines the distance for an interaction. We suggest testing different radii and select the one that generates intuitive clusters upon visualisation. Cells not assigned to clusters are assigned as Cluster_NA in the output table. The argument min_neighborhood_size specifies the threshold of a neighborhood size to be considered as a neighborhood. Smaller neighbourhoods will be outputted, but will not be assigned a number.

Please test out different radii and then visualise the clustering results. To aid in this process, users can use the average_minimum_distance() function, which calculates the average minimum distance between all cells in an image, and can be used as a starting point.

average_minimum_distance(defined_mibi_4)

## [1] 19.15148

We then identify the cellular neighbourhoods using our hierarchical algorithm with a radius of 50, and with a minimum neighbourhood size of 100. Cells assigned to neighbourhoods smaller than 100 will be assigned to the Cluster_NA neighbourhood.

clusters <- identify_neighborhoods(
  defined_mibi_4, method = "hierarchical", min_neighborhood_size = 100,
  cell_types_of_interest = c("B_cells",  "Helper_T_cell", "T_cell",
                "Cyto_T_cell",  "Helper_T_cell", "Neutrophil", "Cyto_T_cell"), radius = 50, 
  feature_colname = "Cell.Type")

Black cells correspond to ‘free’, un-clustered cells.

We can visualise the cell composition of neighborhoods. To do this, we can use composition_of_neighborhoods() to obtain the percentages of cells with a specific marker within each neighborhood and the number of cells in the neighborhood.

In this example we select cellular neighbourhoods with at least 5 cells.

neighorhoods_vis <- 
  composition_of_neighborhoods(clusters, feature_colname = "Cell.Type")
neighorhoods_vis <- 
  neighorhoods_vis[neighorhoods_vis$Total_number_of_cells >=5,]

Finally, we can use plot_composition_heatmap() to produce a heatmap showing the marker percentages within each cluster, which can be used to classify the derived neighbourhoods.

plot_composition_heatmap(neighorhoods_vis, feature_colname="Cell.Type")

Challenge

So far we have walked through the analysis workflow. Now it’s your turn to do some analyses!

The TNBC MIBI data set provides the clinical information of 39 patients. Try applying some spatial metrics to these images and associate them with the clinical data!

You can try formatting the patients’ images yourself with the following commented code.

# # format all patients' images
# # add a column of the patient ID
# mibi$patients_id <- paste("Patient", mibi$SampleID, sep="_")
# ##Removing patients without coordinates. 
# data <- data[!(data$patients_id %in% 
#                  c("Patient_42", "Patient_43", "Patient_44")),]
# ## the map is for defining the cell types
# phenotype_map <- data.frame(Markers = c("Tumor", "CD68", "CD20", "CD4,CD3", 
#                                         "CD3", "CD8,CD3", "CD4", "MPO",
#                                         "CD8"),
#                             Name = c("Tumor", "Macrophage", "B_cells", 
#                                      "Helper_T_cell", "T_cell", "Cyto_T_cell", 
#                                      "Helper_T_cell", "Neutrophil", 
#                                      "Cyto_T_cell"))
# ## format each patient info into an spe object, predict the phenotypes, 
# ## define the cells and save
# MIBI_spe <- list()
# p_id <- unique(mibi$patients_id)
# for (patient in p_id){
#     # get the current patient info
#     local_patient <- subset(data, patients_id == patient)
#     local_patient$CellID <- paste("Cell", local_patient$cellLabelInImage, sep="_")
# 
#     # get the marker intensities
#     intensity <- t(local_patient[,13:52])
#     intensity <- rbind(intensity, Tumor = local_patient$tumorYN)
#     intensity <- intensity[rownames(intensity) %in% c("CD4", "CD3", "CD8",
#                                                       "CD20", "MPO", "CD68",
#                                                       "Tumor"),]
#     # Added Cell IDs to the intensity matrix ----
#     colnames(intensity) <- local_patient$CellID
# 
#     # format spe object
#     local_patient_image <- format_image_to_spe(format = "general",
#                                                intensity_matrix = intensity,
#                                                phenotypes = NA,
#                                                coord_x = local_patient$Xc,
#                                                coord_y = local_patient$Yc)
# 
#     # predict phenotypes
#     local_patient_image <- predict_phenotypes(local_patient_image,
#                                               thresholds = NULL,
#                                               tumour_marker = "Tumor",
#                                               baseline_markers=c("CD4", "CD3",
#                                                                  "CD8", "CD20",
#                                                                  "MPO", "CD68"),
#                                               nuclear_marker = NULL,
#                                               reference_phenotypes = FALSE,
#                                               markers_to_phenotype = NULL,
#                                               plot_distribution=FALSE)
#     pheno <- unique(colData(local_patient_image)$Phenotype)
# 
#     # define cell types
#     other_marker_combinations <- pheno[!(pheno %in% phenotype_map$Markers)]
#     phenotypes_map <- rbind(phenotype_map,
#                             cbind(Markers=other_marker_combinations, Name="None"))
# 
#     local_patient_image <- define_celltypes(spe_object = local_patient_image,
#                                             categories=unique(pheno),
#                                             category_colname = "Phenotype",
#                                             names = phenotypes_map$
#                                                 Name[match(unique(pheno),
#                                                            phenotypes_map$Markers)],
#                                             new_colname = "Cell.Type2")
# 
#     attr(local_patient_image, "name") <- patient
#     MIBI_spe[[patient]] <- local_patient_image
#     print(patient)
# }

However, you can also load the objects that I have formatted for you.

load("MIBI_spe.Rda") # load the MIBI_spe objects

Here is only an example of what you can do with the SPIAT functions.

# calculate the mixing score of each patient
mixing_score <- list()
normalised_mixing_score <- list()
for(patient in names(MIBI_spe)){
   ms <- mixing_score_summary(
    spe_object = MIBI_spe[[patient]], reference_celltype = "Tumor",
    target_celltype = "Macrophage", radius = 20, feature_colname = "Cell.Type2")
  mixing_score[[patient]] <- ms$Mixing_score
  normalised_mixing_score[[patient]] <- ms$Normalised_mixing_score
  print(patient)
}

## [1] "Patient_1"
## [1] "Patient_2"
## [1] "Patient_3"
## [1] "Patient_4"
## [1] "Patient_5"
## [1] "Patient_6"
## [1] "Patient_7"
## [1] "Patient_8"
## [1] "Patient_9"
## [1] "Patient_10"
## [1] "Patient_11"
## [1] "Patient_12"
## [1] "Patient_13"
## [1] "Patient_14"
## [1] "Patient_15"
## [1] "Patient_16"
## [1] "Patient_17"
## [1] "Patient_18"
## [1] "Patient_19"
## [1] "Patient_20"
## [1] "Patient_21"
## [1] "Patient_22"
## [1] "Patient_23"
## [1] "Patient_24"
## [1] "Patient_25"
## [1] "Patient_26"
## [1] "Patient_27"
## [1] "Patient_28"
## [1] "Patient_29"
## [1] "There are no unique reference cells of specified celltype Tumor for target cell Macrophage"
## [1] "Patient_31"
## [1] "Patient_32"
## [1] "Patient_33"
## [1] "Patient_34"
## [1] "Patient_35"
## [1] "Patient_36"
## [1] "Patient_37"
## [1] "Patient_38"
## [1] "Patient_39"
## [1] "Patient_40"
## [1] "Patient_41"

# group the patients based on their normalised mixing score
nms <- na.omit(unlist(normalised_mixing_score))

low_nms <- normalised_mixing_score < 0.25 # this is an arbitrary threshold, adopted from the spaSim benchmarking
group_repulsion <- names(which(low_nms == TRUE))
group_attraction <- names(which(low_nms == FALSE))

# read in metadata
metadata <- read.delim("Formatted_clinical_data.csv", sep = ",")
metadata$Patient <- paste("Patient", metadata$InternalId, sep="_")
rownames(metadata) <- metadata$Patient

# join the data
metadata$Mixing_group <- NA
metadata$Mixing_group[rownames(metadata) %in% group_repulsion] <- "Group_Repulsion"
metadata$Mixing_group[rownames(metadata) %in% group_attraction] <- "Group_Attraction"

# survival
library(survminer)
library(survival)
survival <- survfit(Surv(Survival_days_capped., Censored) ~ Mixing_group, data=metadata)
  
g <- ggsurvplot(survival, 
             xlab = "Days", 
             ylab = "Overall survival probability", pval=TRUE)
print(g)

Introduction to spaSim

spaSim (spatial Simulator) is a simulator of tumour immune microenvironment spatial data. It includes a family of functions to simulate a diverse set of cell localisation patterns in tissues. Patterns include background cells (one cell type or multiple cell types of different proportions), tumour/immune clusters, immune rings and double immune rings and stripes (blood/lymphatic vessels).

Simulations from spaSim can be applied to test and benchmark spatial tools and metrics (like the ones in SPIAT). The output of spaSim are images in SpatialExperiment object format.

In this workshop, we will only cover how to simulate individual images. For more information about how to simulate a set of images, please refer to the spaSim tutorial.

# simulate background cells
set.seed(610)
bg <- simulate_background_cells(n_cells = 5000,
                                width = 2000,
                                height = 2000,
                                method = "Hardcore",
                                min_d = 10,
                                oversampling_rate = 1.6,
                                Cell.Type = "Others")

# simulate mixed background
mix_bg <- simulate_mixing(bg_sample = bg,
                          idents = c("Tumour", "Immune", "Others"),
                          props = c(0.2, 0.3, 0.5), 
                          plot_image = TRUE,
                          plot_colours = c("red","darkgreen","lightgray"))

# simulate clusters
cluster_properties <- list(
  C1 =list(name_of_cluster_cell = "Tumour", size = 500, shape = "Oval", 
           centre_loc = data.frame(x = 600, y = 600),infiltration_types = c("Immune1", "Others"), 
           infiltration_proportions = c(0.1, 0.05)), 
  C2 = list(name_of_cluster_cell = "Immune1", size = 600,  shape = "Irregular", 
            centre_loc = data.frame(x = 1500, y = 500), infiltration_types = c("Immune", "Others"),
            infiltration_proportions = c(0.1, 0.05)))
# can use any defined image as background image, here we use mix_bg defined in the previous section
clusters <- simulate_clusters(bg_sample = mix_bg,
                              n_clusters = 2,
                              bg_type = "Others",
                              cluster_properties = cluster_properties,
                              plot_image = TRUE,
                              plot_categories = c("Tumour" , "Immune", "Immune1", "Others"),
                              plot_colours = c("red", "darkgreen", "darkblue", "lightgray"))

# simulate immune rings
immune_ring_properties <- list(
  I1 = list(name_of_cluster_cell = "Tumour", size = 500, 
            shape = "Circle", centre_loc = data.frame(x = 930, y = 1000), 
            infiltration_types = c("Immune1", "Immune2", "Others"), 
            infiltration_proportions = c(0.15, 0.05, 0.05),
            name_of_ring_cell = "Immune1", immune_ring_width = 150,
            immune_ring_infiltration_types = c("Immune2", "Others"), 
            immune_ring_infiltration_proportions = c(0.1, 0.15)))
rings <- simulate_immune_rings(
  bg_sample = bg,
  bg_type = "Others",
  n_ir = 1,
  ir_properties = immune_ring_properties,
  plot_image = TRUE,
  plot_categories = c("Tumour", "Immune1", "Immune2", "Others"),
  plot_colours = c("red", "darkgreen", "darkblue", "lightgray"))

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.3.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Australia/Melbourne
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] survival_3.5-7              survminer_0.4.9            
##  [3] ggpubr_0.6.0                ggplot2_3.4.4              
##  [5] spaSim_1.4.0                SPIAT_1.4.1                
##  [7] SpatialExperiment_1.12.0    SingleCellExperiment_1.24.0
##  [9] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [11] GenomicRanges_1.54.1        GenomeInfoDb_1.38.0        
## [13] IRanges_2.36.0              S4Vectors_0.40.1           
## [15] BiocGenerics_0.48.1         MatrixGenerics_1.14.0      
## [17] matrixStats_1.0.0          
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2           bitops_1.0-7            tibble_3.2.1           
##   [4] R.oo_1.25.0             polyclip_1.10-6         lifecycle_1.0.3        
##   [7] rstatix_0.7.2           doParallel_1.0.17       lattice_0.22-5         
##  [10] vroom_1.6.4             crosstalk_1.2.0         backports_1.4.1        
##  [13] magrittr_2.0.3          plotly_4.10.3           sass_0.4.7             
##  [16] rmarkdown_2.25          jquerylib_0.1.4         yaml_2.3.7             
##  [19] sp_2.1-1                spatstat.sparse_3.0-3   cowplot_1.1.1          
##  [22] RColorBrewer_1.1-3      abind_1.4-5             zlibbioc_1.48.0        
##  [25] Rtsne_0.16              purrr_1.0.2             R.utils_2.12.2         
##  [28] RCurl_1.98-1.13         pracma_2.4.2            circlize_0.4.15        
##  [31] GenomeInfoDbData_1.2.11 KMsurv_0.1-5            ggrepel_0.9.4          
##  [34] spatstat.utils_3.0-4    terra_1.7-55            pheatmap_1.0.12        
##  [37] goftest_1.2-3           spatstat.random_3.2-1   codetools_0.2-19       
##  [40] DelayedArray_0.28.0     alphahull_2.5           tidyselect_1.2.0       
##  [43] shape_1.4.6             raster_3.6-26           elsa_1.1-28            
##  [46] farver_2.1.1            sgeostat_1.0-27         spatstat.explore_3.2-5 
##  [49] jsonlite_1.8.7          GetoptLong_1.0.5        ellipsis_0.3.2         
##  [52] ggridges_0.5.4          iterators_1.0.14        foreach_1.5.2          
##  [55] dbscan_1.1-11           tools_4.3.2             Rcpp_1.0.11            
##  [58] glue_1.6.2              gridExtra_2.3           SparseArray_1.2.0      
##  [61] xfun_0.41               dplyr_1.1.3             withr_2.5.2            
##  [64] fastmap_1.1.1           fansi_1.0.5             digest_0.6.33          
##  [67] R6_2.5.1                colorspace_2.1-0        gtools_3.9.4           
##  [70] tensor_1.5              spatstat.data_3.0-3     R.methodsS3_1.8.2      
##  [73] utf8_1.2.4              tidyr_1.3.0             generics_0.1.3         
##  [76] data.table_1.14.8       httr_1.4.7              htmlwidgets_1.6.2      
##  [79] S4Arrays_1.2.0          pkgconfig_2.0.3         gtable_0.3.4           
##  [82] ComplexHeatmap_2.18.0   XVector_0.42.0          survMisc_0.5.6         
##  [85] htmltools_0.5.7         carData_3.0-5           clue_0.3-65            
##  [88] scales_1.2.1            png_0.1-8               splancs_2.01-44        
##  [91] knitr_1.45              km.ci_0.5-6             rstudioapi_0.15.0      
##  [94] tzdb_0.4.0              reshape2_1.4.4          rjson_0.2.21           
##  [97] nlme_3.1-163            zoo_1.8-12              cachem_1.0.8           
## [100] GlobalOptions_0.1.2     stringr_1.5.0           parallel_4.3.2         
## [103] pillar_1.9.0            grid_4.3.2              vctrs_0.6.4            
## [106] RANN_2.6.1              car_3.1-2               dittoSeq_1.14.0        
## [109] xtable_1.8-4            cluster_2.1.4           evaluate_0.23          
## [112] magick_2.8.1            cli_3.6.1               compiler_4.3.2         
## [115] rlang_1.1.2             crayon_1.5.2            ggsignif_0.6.4         
## [118] labeling_0.4.3          interp_1.1-4            plyr_1.8.9             
## [121] stringi_1.7.12          viridisLite_0.4.2       deldir_1.0-9           
## [124] munsell_0.5.0           lazyeval_0.2.2          spatstat.geom_3.2-7    
## [127] Matrix_1.6-1.1          bit64_4.0.5             highr_0.10             
## [130] apcluster_1.4.11        mmand_1.6.3             broom_1.0.5            
## [133] bslib_0.5.1             bit_4.0.5