+
+```r
+assay(spe, "exprs") <- asinh(counts(spe)/1)
+dittoRidgePlot(spe, var = "CD3", group.by = "patient_id", assay = "exprs") +
+ ggtitle("CD3 - after transformation")
+```
+
+
+
+**Define interesting channels**
+
+For downstream analysis such as visualization, dimensionality reduction and
+clustering, only a subset of markers should be used. As convenience, we can
+store an additional entry in the `rowData` slot that specifies the markers of
+interest. Here, we deselect the nuclear markers, which were primarily used for
+cell segmentation, and keep all other biological targets. However, more informed
+marker selection should be performed to exclude lowly expressed marker or
+markers with low signal-to-noise ratio.
+
+
+```r
+rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))
+```
+
+**Define color schemes**
+
+We will define color schemes for different metadata entries of the data and
+conveniently store them in the `metadata` slot of the `SpatialExperiment` which
+will be helpful for downstream data visualizations. We will use colors from the
+`RColorBrewer` and `dittoSeq` packages but any other coloring package will
+suffice.
+
+
+```r
+library(RColorBrewer)
+color_vectors <- list()
+
+ROI <- setNames(brewer.pal(length(unique(spe$ROI)), name = "BrBG"),
+ unique(spe$ROI))
+patient_id <- setNames(brewer.pal(length(unique(spe$patient_id)), name = "Set1"),
+ unique(spe$patient_id))
+sample_id <- setNames(c(brewer.pal(6, "YlOrRd")[3:5],
+ brewer.pal(6, "PuBu")[3:6],
+ brewer.pal(6, "YlGn")[3:5],
+ brewer.pal(6, "BuPu")[3:6]),
+ unique(spe$sample_id))
+indication <- setNames(brewer.pal(length(unique(spe$indication)), name = "Set2"),
+ unique(spe$indication))
+
+color_vectors$ROI <- ROI
+color_vectors$patient_id <- patient_id
+color_vectors$sample_id <- sample_id
+color_vectors$indication <- indication
+
+metadata(spe)$color_vectors <- color_vectors
+```
+
+## Read in images {#read-images}
+
+The `cytomapper` package allows multi-channel image handling and visualization
+within the Bioconductor framework. The most common data format for multi-channel
+images or segmentation masks is the TIFF file format, which is used by `steinbock`
+and the `IMC segementation pipeline` to save images.
+
+Here, we will read in multi-channel images and segmentation masks into a
+[CytoImageList](https://www.bioconductor.org/packages/release/bioc/vignettes/cytomapper/inst/doc/cytomapper.html#5_The_CytoImageList_object)
+data container. It allows storing multiple multi-channel images and requires
+matched channels across all images within the object.
+
+The `loadImages` function is used to read in processed multi-channel images and
+their corresponding segmentation masks. **Of note**: the multi-channel images
+generated by `steinbock` are saved as 32-bit images while the segmentation masks
+are saved as 16-bit images. To correctly scale pixel values of the segmentation
+masks when reading them in, we will need to set `as.is = TRUE`.
+
+
+```r
+images <- loadImages("data/steinbock/img/")
+```
+
+```
+## All files in the provided location will be read in.
+```
+
+```r
+masks <- loadImages("data/steinbock/masks_deepcell/", as.is = TRUE)
+```
+
+```
+## All files in the provided location will be read in.
+```
+
+In the case of multi-channel images, it is beneficial to set the `channelNames`
+for easy visualization. Using the `steinbock` framework, the channel order of
+the single-cell data matches the channel order of the multi-channel images.
+However, it is recommended to make sure that the channel order is identical
+between the single-cell data and the images.
+
+
+```r
+channelNames(images) <- rownames(spe)
+images
+```
+
+```
+## CytoImageList containing 14 image(s)
+## names(14): Patient1_001 Patient1_002 Patient1_003 Patient2_001 Patient2_002 Patient2_003 Patient2_004 Patient3_001 Patient3_002 Patient3_003 Patient4_005 Patient4_006 Patient4_007 Patient4_008
+## Each image contains 40 channel(s)
+## channelNames(40): MPO HistoneH3 SMA CD16 CD38 HLADR CD27 CD15 CD45RA CD163 B2M CD20 CD68 Ido1 CD3 LAG3 / LAG33 CD11c PD1 PDGFRb CD7 GrzB PDL1 TCF7 CD45RO FOXP3 ICOS CD8a CarbonicAnhydrase CD33 Ki67 VISTA CD40 CD4 CD14 Ecad CD303 CD206 cleavedPARP DNA1 DNA2
+```
+
+For visualization shown in Section \@ref(image-visualization) we will need to
+add additional metadata to the `elementMetadata` slot of the `CytoImageList`
+objects. This slot is easily accessible using the `mcols` function.
+
+Here, we will store the matched `sample_id`, `patient_id` and `indication`
+information within the `elementMetadata` slot of the multi-channel images and
+segmentation masks objects. It is crucial that the order of the images in
+both `CytoImageList` objects is the same.
+
+
+```r
+all.equal(names(images), names(masks))
+```
+
+```
+## [1] TRUE
+```
+
+```r
+# Extract patient id from image name
+patient_id <- str_extract(names(images), "Patient[1-4]")
+
+# Retrieve cancer type per patient from metadata file
+indication <- meta$Indication[match(patient_id, meta$`Sample ID`)]
+
+# Store patient and image level information in elementMetadata
+mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
+ patient_id = patient_id,
+ indication = indication)
+```
+
+## Generate single-cell data from images
+
+An alternative way of generating a `SingleCellExperiment` object directly
+from the multi-channel images and segmentation masks is supported by the
+[measureObjects](https://bodenmillergroup.github.io/cytomapper/reference/measureObjects.html)
+function of the `cytomapper` package. For each cell present in the `masks`
+object, the function computes the mean pixel intensity per channel as well as
+morphological features (area, radius, major axis length, eccentricity) and the
+location of cells:
+
+
+```r
+cytomapper_sce <- measureObjects(masks, image = images, img_id = "sample_id")
+
+cytomapper_sce
+```
+
+```
+## class: SingleCellExperiment
+## dim: 40 47859
+## metadata(0):
+## assays(1): counts
+## rownames(40): MPO HistoneH3 ... DNA1 DNA2
+## rowData names(0):
+## colnames: NULL
+## colData names(10): sample_id object_id ... patient_id indication
+## reducedDimNames(0):
+## mainExpName: NULL
+## altExpNames(0):
+```
+
+## Accessing publicly available IMC datasets
+
+The [imcdatasets](https://github.com/BodenmillerGroup/imcdatasets)
+R/Bioconductor package provides a number of publicly available IMC datasets. For
+a complete introduction to the package, please refer to the
+[documentation](https://bioconductor.org/packages/release/data/experiment/vignettes/imcdatasets/inst/doc/imcdatasets.html).
+Here, we can read in example data of [@Damond2019] taken from patients diagnosed
+with Type I Diabetes. The example here consists of a `CytoImageList` object of
+100 images, a `CytoImageList` object of 100 segmentation masks and a
+`SingleCellExperiment` object containing 252059 cells. **Of note:** downloading the
+images takes quite some time and uses 8GB of memory.
+
+
+```r
+library(imcdatasets)
+
+pancreasImages <- Damond_2019_Pancreas(data_type = "images")
+pancreasMasks <- Damond_2019_Pancreas(data_type = "masks")
+pancreasSCE <- Damond_2019_Pancreas(data_type = "sce")
+```
+
+## Save objects
+
+Finally, the generated data objects can be saved for further downstream
+processing and analysis.
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+saveRDS(images, "data/images.rds")
+saveRDS(masks, "data/masks.rds")
+```
+
+
+
+
+
+
+
+
+
+
+
+## Session Info
+
+
+
+```r
+# Thresholded on 200 pixel counts
+plotSpotHeatmap(sce, log = FALSE, threshold = 200)
+```
+
+
+
+As we can see, nearly all median pixel intensities are > 200 counts for each spot.
+We also observe acquired channels for which no spot was placed (e.g., Xe134, Ir191, Ir193).
+
+### Optional pixel binning {#pixel_binning}
+
+In cases where median pixel intensities are low (< 200 counts), consecutive
+pixels can be summed to increase the robustness of the spillover estimation.
+The `imcRtools` package provides the `binAcrossPixels` function,
+which performs aggregation for each channel across `bin_size` consecutive pixels
+per spotted metal.
+
+
+```r
+# Define grouping
+bin_size = 10
+
+sce2 <- binAcrossPixels(sce, bin_size = bin_size)
+
+# Log10 median pixel counts per spot and channel
+plotSpotHeatmap(sce2)
+```
+
+
+
+```r
+# Thresholded on 200 pixel counts
+plotSpotHeatmap(sce2, log = FALSE, threshold = 200)
+```
+
+
+
+Here, we can see an increase in the median pixel intensities and accumulation of
+off-diagonal signal. Due to already high original pixel intensities, we will
+refrain from aggregating across consecutive pixels for this demonstration.
+
+### Filtering incorrectly assigned pixels
+
+The following step uses functions provided by the `CATALYST` package to
+"debarcode" the pixels. Based on the intensity distribution of all channels,
+pixels are assigned to their corresponding barcode; here this is the already
+known metal spot. This procedure serves the purpose to identify pixels that
+cannot be robustly assigned to the spotted metal. Pixels of such kind can be
+regarded as "noisy", "background" or "artefacts" that should be removed prior to
+spillover estimation.
+
+We will also need to specify which channels were spotted (argument `bc_key`).
+This information is directly contained in the `colData(sce)` slot.
+To facilitate visualization, we will order the `bc_key` by mass.
+
+The general workflow for pixel debarcoding is as follows:
+
+1. assign a preliminary metal mass to each pixel
+2. for each pixel, estimate a cutoff parameter for the distance between
+positive and negative pixel sets
+3. apply the estimated cutoffs to identify truly positive pixels
+
+
+```r
+library(CATALYST)
+
+bc_key <- as.numeric(unique(sce$sample_mass))
+bc_key <- bc_key[order(bc_key)]
+
+sce <- assignPrelim(sce, bc_key = bc_key)
+sce <- estCutoffs(sce)
+sce <- applyCutoffs(sce)
+```
+
+The obtained `SingleCellExperiment` now contains the additional `bc_id` entry.
+For each pixel, this vector indicates the assigned mass (e.g. `161`) or
+`0`, meaning unassigned.
+
+This information can be visualized in form of a heatmap:
+
+
+```r
+library(pheatmap)
+cur_table <- table(sce$bc_id, sce$sample_mass)
+
+# Visualize the correctly and incorrectly assigned pixels
+pheatmap(log10(cur_table + 1),
+ cluster_rows = FALSE,
+ cluster_cols = FALSE)
+```
+
+
+
+```r
+# Compute the fraction of unassigned pixels per spot
+cur_table["0",] / colSums(cur_table)
+```
+
+```
+## 113 115 141 142 143 144 145 146 147 148 149
+## 0.1985 0.1060 0.2575 0.3195 0.3190 0.3825 0.3545 0.4280 0.3570 0.4770 0.4200
+## 150 151 152 153 154 155 156 158 159 160 161
+## 0.4120 0.4025 0.4050 0.4630 0.4190 0.4610 0.3525 0.4020 0.4655 0.4250 0.5595
+## 162 163 164 165 166 167 168 169 170 171 172
+## 0.4340 0.4230 0.4390 0.4055 0.5210 0.3900 0.3285 0.3680 0.5015 0.4900 0.5650
+## 173 174 175 176 89
+## 0.3125 0.4605 0.4710 0.2845 0.3015
+```
+
+We can see here, that all pixels were assigned to the right mass and that all
+pixel sets are made up of > 800 pixels.
+
+However, in cases where incorrect assignment occurred or where few pixels were
+measured for some spots, the `imcRtools` package exports a simple helper
+function to exclude pixels based on these criteria:
+
+
+```r
+sce <- filterPixels(sce, minevents = 40, correct_pixels = TRUE)
+```
+
+In the `filterPixels` function, the `minevents` parameter specifies the threshold
+under which correctly assigned pixel sets are excluded from spillover
+estimation. The `correct_pixels` parameter indicates whether pixels that were
+assigned to masses other than the spotted mass should be excluded from spillover
+estimation. The default values often result in sufficient pixel filtering;
+however, if very few pixels (~100) are measured per spot, the `minevents`
+parameter value needs to be lowered.
+
+### Compute spillover matrix
+
+Based on the single-positive pixels, we use the `CATALYST::computeSpillmat()`
+function to compute the spillover matrix and `CATALYST::plotSpillmat()` to
+visualize it. The `plotSpillmat` function checks the spotted and acquired
+metal isotopes against a pre-defined `CATALYST::isotope_list()`. In this data,
+the `Ar80` channel was additionally acquired to check for deviations in signal
+intensity. `Ar80` needs to be added to a custom `isotope_list` object for
+visualization.
+
+
+```r
+sce <- computeSpillmat(sce)
+
+isotope_list <- CATALYST::isotope_list
+isotope_list$Ar <- 80
+
+plotSpillmat(sce, isotope_list = isotope_list)
+```
+
+```
+## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
+## ggplot2 3.3.4.
+## ℹ Please use "none" instead.
+## ℹ The deprecated feature was likely used in the CATALYST package.
+## Please report the issue at 
+
+```r
+# Save spillover matrix in variable
+sm <- metadata(sce)$spillover_matrix
+```
+
+**Of note: the visualization of the spillover matrix using CATALYST does currently
+not visualize spillover between the larger channels.** In this case, the
+spillover matrix is clipped at Yb171.
+
+As we can see, the largest spillover appears in `In113 --> In115` and we also
+observe the `+16` oxide impurities for e.g. `Nd148 --> Dy164`.
+
+We can save the spillover matrix for external use.
+
+
+```r
+write.csv(sm, "data/sm.csv")
+```
+
+## Single-cell data compensation
+
+The `CATALYST` package can be used to perform spillover compensation on the
+**single-cell mean intensities**. Here, the `SpatialExperiment` object generated
+in Section \@ref(read-data) is read in. The `CATALYST` package requires an entry
+to `rowData(spe)$channel_name` for the `compCytof` function to run. This entry
+should contain the metal isotopes in the form (mt)(mass)Di (e.g., `Sm152Di` for
+Samarium isotope with the atomic mass 152).
+
+The `compCytof` function performs channel spillover compensation on the mean
+pixel intensities per channel and cell. Here, we will not overwrite the assays
+in the `SpatialExperiment` object to later highlight the effect of compensation.
+As shown in Section \@ref(read-data), also the compensated counts are
+asinh-transformed using a cofactor of 1.
+
+
+```r
+spe <- readRDS("data/spe.rds")
+rowData(spe)$channel_name <- paste0(rowData(spe)$channel, "Di")
+
+spe <- compCytof(spe, sm,
+ transform = TRUE, cofactor = 1,
+ isotope_list = isotope_list,
+ overwrite = FALSE)
+```
+
+To check the effect of channel spillover compensation, the expression of markers
+that are affected by spillover (e.g., E-cadherin in channel Yb173 and CD303 in
+channel Yb174) can be visualized in form of scatter plots before and after
+compensation.
+
+
+```r
+library(dittoSeq)
+library(patchwork)
+before <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
+ assay.x = "exprs", assay.y = "exprs") +
+ ggtitle("Before compensation")
+
+after <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
+ assay.x = "compexprs", assay.y = "compexprs") +
+ ggtitle("After compensation")
+before + after
+```
+
+
+
+We observe that the spillover Yb173 --> Yb174 was successfully corrected.
+To facilitate further downstream analysis, the non-compensated assays can now be
+replaced by their compensated counterparts:
+
+
+```r
+assay(spe, "counts") <- assay(spe, "compcounts")
+assay(spe, "exprs") <- assay(spe, "compexprs")
+assay(spe, "compcounts") <- assay(spe, "compexprs") <- NULL
+```
+
+## Image compensation
+
+The [cytomapper](https://github.com/BodenmillerGroup/cytomapper) package allows channel
+spillover compensation directly on **multi-channel images**.
+The `compImage` function takes a `CytoImageList` object and the estimated
+spillover matrix as input. More info on how to work with `CytoImageList`
+objects can be seen in Section \@ref(image-visualization).
+
+At this point, we can read in the `CytoImageList` object containing multi-channel
+images as generated in Section \@ref(read-data).
+The `channelNames` need to be set according to their metal isotope in the form
+(mt)(mass)Di and therefore match `colnames(sm)`.
+
+
+```r
+library(cytomapper)
+
+images <- readRDS("data/images.rds")
+channelNames(images) <- rowData(spe)$channel_name
+```
+
+The CATALYST package provides the `adaptSpillmat` function that corrects the
+spillover matrix in a way that rows and columns match a predefined set of
+metals. Please refer to `?compCytof` for more information how metals in the
+spillover matrix are matched to acquired channels in the `SingleCellExperiment`
+object.
+
+The spillover matrix can now be adapted to exclude channels that were not kept
+for downstream analysis.
+
+
+```r
+adapted_sm <- adaptSpillmat(sm, channelNames(images),
+ isotope_list = isotope_list)
+```
+
+```
+## Compensation is likely to be inaccurate.
+## Spill values for the following interactions
+## have not been estimated:
+```
+
+```
+## Ir191Di -> Ir193Di
+```
+
+```
+## Ir193Di -> Ir191Di
+```
+
+The adapted spillover matrix now matches the `channelNames` of the
+`CytoImageList` object and can be used to perform pixel-level spillover
+compensation. Here, we parallelise the image compensation on all available minus 2 cores. When
+working on Windows, you will need to use the `SnowParam` function instead of
+`MultiCoreParam`.
+
+
+```r
+library(BiocParallel)
+
+images_comp <- compImage(images, adapted_sm,
+ BPPARAM = MulticoreParam())
+```
+
+As a sanity check, we will visualize the image before and after compensation:
+
+
+```r
+# Before compensation
+plotPixels(images[5], colour_by = "Yb173Di",
+ image_title = list(text = "Yb173 (Ecad) - before", position = "topleft"),
+ legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
+```
+
+
+
+```r
+plotPixels(images[5], colour_by = "Yb174Di",
+ image_title = list(text = "Yb174 (CD303) - before", position = "topleft"),
+ legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
+```
+
+
+
+```r
+# After compensation
+plotPixels(images_comp[5], colour_by = "Yb173Di",
+ image_title = list(text = "Yb173 (Ecad) - after", position = "topleft"),
+ legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
+```
+
+
+
+```r
+plotPixels(images_comp[5], colour_by = "Yb174Di",
+ image_title = list(text = "Yb174 (CD303) - after", position = "topleft"),
+ legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
+```
+
+
+
+For convenience, we will re-set the `channelNames` to their biological targtes:
+
+
+```r
+channelNames(images_comp) <- rownames(spe)
+```
+
+
+
+## Write out compensated images
+
+In the final step, the compensated images are written out as 16-bit TIFF
+files:
+
+
+```r
+library(tiff)
+dir.create("data/comp_img")
+lapply(names(images_comp), function(x){
+ writeImage(as.array(images_comp[[x]])/(2^16 - 1),
+ paste0("data/comp_img/", x, ".tiff"),
+ bits.per.sample = 16)
+})
+```
+
+## Save objects
+
+For further downstream analysis, the compensated `SpatialExperiment` and
+`CytoImageList` objects are saved replacing the former objects:
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+saveRDS(images_comp, "data/images.rds")
+```
+
+
+
+## Session Info
+
+
+
+We can see that nuclei are centered within the segmentation masks and all cell
+types are correctly segmented (note: to zoom into the image you can right click
+and select `Open Image in New Tab`). A common challenge here is to segment large (e.g.,
+epithelial cells - in cyan) _versus_ small (e.g., B cells - in red). However, the
+segmentation approach here appears to correctly segment cells across different
+sizes.
+
+An easier and interactive way of observing segmentation quality is to use the
+interactive image viewer provided by the
+[cytoviewer](https://github.com/BodenmillerGroup/cytoviewer) R/Bioconductor
+package [@Meyer2023]. Under "Image-level" > "Basic controls", up to six markers
+can be selected for visualization. The contrast of each marker can be adjusted.
+Under "Image-level" > "Advanced controls", click the "Show cell outlines" box
+to outline segmented cells on the images.
+
+
+```r
+library(cytoviewer)
+
+app <- cytoviewer(image = images,
+ mask = masks,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id")
+
+if (interactive()) {
+ shiny::runApp(app, launch.browser = TRUE)
+}
+```
+
+An additional approach to observe cell segmentation quality and potentially also
+antibody specificity issues is to visualize single-cell expression in form of a
+heatmap. Here, we sub-sample the dataset to 2000 cells for visualization
+purposes and overlay the cancer type from which the cells were extracted.
+
+
+```r
+library(dittoSeq)
+library(viridis)
+cur_cells <- sample(seq_len(ncol(spe)), 2000)
+
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = "indication",
+ annotation_colors = list(indication = metadata(spe)$color_vectors$indication))
+```
+
+
+
+We can differentiate between epithelial cells (Ecad+) and immune cells
+(CD45RO+). Some of the markers are detected in specific cells (e.g., Ki67, CD20,
+Ecad) while others are more broadly expressed across cells (e.g., HLADR, B2M,
+CD4).
+
+## Image-level quality control {#image-quality}
+
+Image-level quality control is often performed using tools that offer a
+graphical user interface such as [QuPath](https://qupath.github.io/),
+[FIJI](https://imagej.net/software/fiji/) and the previously mentioned
+[cytoviewer](https://github.com/BodenmillerGroup/cytoviewer) package. Viewers
+that were specifically developed for IMC data can be seen
+[here](https://bodenmillergroup.github.io/IMCWorkflow/viewers.html). In this
+section, we will specifically focus on quantitative metrics to assess image
+quality.
+
+It is often of interest to calculate the signal-to-noise ratio (SNR) for
+individual channels and markers. Here, we define the SNR as:
+
+$$SNR = I_s/I_n$$
+
+where $I_s$ is the intensity of the signal (mean intensity of pixels with true
+signal) and $I_n$ is the intensity of the noise (mean intensity of pixels
+containing noise). This definition of the SNR is just one of many and other
+measures can be applied. Finding a threshold that separates pixels containing
+signal and pixels containing noise is not trivial and different approaches can
+be chosen. Here, we use the `otsu` thresholding approach to find pixels of the
+"foreground" (i.e., signal) and "background" (i.e., noise). The SNR is then
+defined as the mean intensity of foreground pixels divided by the mean intensity
+of background pixels. We compute this measure as well as the mean signal
+intensity per image. The plot below shows the average SNR _versus_ the average
+signal intensity across all images.
+
+
+```r
+library(tidyverse)
+library(ggrepel)
+library(EBImage)
+
+cur_snr <- lapply(names(images), function(x){
+ img <- images[[x]]
+ mat <- apply(img, 3, function(ch){
+ # Otsu threshold
+ thres <- otsu(ch, range = c(min(ch), max(ch)), levels = 65536)
+ # Signal-to-noise ratio
+ snr <- mean(ch[ch > thres]) / mean(ch[ch <= thres])
+ # Signal intensity
+ ps <- mean(ch[ch > thres])
+
+ return(c(snr = snr, ps = ps))
+ })
+ t(mat) %>% as.data.frame() %>%
+ mutate(image = x,
+ marker = colnames(mat)) %>%
+ pivot_longer(cols = c(snr, ps))
+})
+
+cur_snr <- do.call(rbind, cur_snr)
+
+cur_snr %>%
+ group_by(marker, name) %>%
+ summarize(log_mean = log2(mean(value))) %>%
+ pivot_wider(names_from = name, values_from = log_mean) %>%
+ ggplot() +
+ geom_point(aes(ps, snr)) +
+ geom_label_repel(aes(ps, snr, label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
+```
+
+
+
+We observe PD1, LAG3 and cleaved PARP to have high SNR but low signal intensity
+meaning that in general these markers are not abundantly expressed. The Iridium
+intercalator (here marked as DNA1 and DNA2) has the highest signal intensity
+but low SNR. This might be due to staining differences between individual nuclei
+where some nuclei are considered as background. We do however observe high
+SNR and sufficient signal intensity for the majority of markers.
+
+Otsu thesholding and SNR calculation does not perform well if the markers are
+lowly abundant. In the next code chunk, we will remove markers that have
+a positive signal of below 2 per image.
+
+
+```r
+cur_snr <- cur_snr %>%
+ pivot_wider(names_from = name, values_from = value) %>%
+ filter(ps > 2) %>%
+ pivot_longer(cols = c(snr, ps))
+
+cur_snr %>%
+ group_by(marker, name) %>%
+ summarize(log_mean = log2(mean(value))) %>%
+ pivot_wider(names_from = name, values_from = log_mean) %>%
+ ggplot() +
+ geom_point(aes(ps, snr)) +
+ geom_label_repel(aes(ps, snr, label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
+```
+
+
+
+This visualization shows a reduces SNR for PD1, LAG3 and cleaved PARP which was
+previously inflated due to low signal.
+
+Another quality indicator is the image area covered by cells (or biological
+tissue). This metric identifies ROIs where little cells are present, possibly
+hinting at incorrect selection of the ROI. We can compute the percentage of
+covered image area using the metadata contained in the `SpatialExperiment`
+object:
+
+
+```r
+cell_density <- colData(spe) %>%
+ as.data.frame() %>%
+ group_by(sample_id) %>%
+ # Compute the number of pixels covered by cells and
+ # the total number of pixels
+ summarize(cell_area = sum(area),
+ no_pixels = mean(width_px) * mean(height_px)) %>%
+ # Divide the total number of pixels
+ # by the number of pixels covered by cells
+ mutate(covered_area = cell_area / no_pixels)
+
+# Visualize the image area covered by cells per image
+ggplot(cell_density) +
+ geom_point(aes(reorder(sample_id,covered_area), covered_area)) +
+ theme_minimal(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
+ ylim(c(0, 1)) +
+ ylab("% covered area") + xlab("")
+```
+
+
+
+We observe that two of the 14 images show unusually low cell coverage. These
+two images can now be visualized using `cytomapper`.
+
+
+```r
+# Normalize and clip images
+cur_images <- images[c("Patient4_005", "Patient4_007")]
+cur_images <- cytomapper::normalize(cur_images, separateImages = TRUE)
+cur_images <- cytomapper::normalize(cur_images, inputRange = c(0, 0.2))
+
+plotPixels(cur_images,
+ mask = masks[c("Patient4_005", "Patient4_007")],
+ img_id = "sample_id",
+ missing_colour = "white",
+ colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
+ colour = list(CD163 = c("black", "yellow"),
+ CD20 = c("black", "red"),
+ CD3 = c("black", "green"),
+ Ecad = c("black", "cyan"),
+ DNA1 = c("black", "blue")),
+ legend = list(colour_by.title.cex = 0.7,
+ colour_by.labels.cex = 0.7))
+```
+
+
+
+These two images display less dense tissue structure but overall the images are
+intact and appear to be segmented correctly.
+
+Finally, it can be beneficial to visualize the mean marker expression per image
+to identify images with outlying marker expression. This check does not
+indicate image quality _per se_ but can highlight biological differences. Here,
+we will use the `aggregateAcrossCells` function of the
+*[scuttle](https://bioconductor.org/packages/3.17/scuttle)* package to compute the mean expression per
+image. For visualization purposes, we again `asinh` transform the mean expression
+values.
+
+
+```r
+library(scuttle)
+
+image_mean <- aggregateAcrossCells(spe,
+ ids = spe$sample_id,
+ statistics="mean",
+ use.assay.type = "counts")
+assay(image_mean, "exprs") <- asinh(counts(image_mean))
+
+dittoHeatmap(image_mean, genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", cluster_cols = TRUE, scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("indication", "patient_id", "ROI"),
+ annotation_colors = list(indication = metadata(spe)$color_vectors$indication,
+ patient_id = metadata(spe)$color_vectors$patient_id,
+ ROI = metadata(spe)$color_vectors$ROI),
+ show_colnames = TRUE)
+```
+
+
+
+We observe extensive biological variation across the 14 images specifically for
+some of the cell phenotype markers including the macrophage marker CD206, the B
+cell marker CD20, the neutrophil marker CD15, and the proliferation marker Ki67.
+These differences will be further studied in the following chapters.
+
+## Cell-level quality control {#cell-quality}
+
+In the following paragraphs we will look at different metrics and visualization
+approaches to assess data quality (as well as biological differences) on the
+single-cell level.
+
+Related to the signal-to-noise ratio (SNR) calculated above on the pixel-level,
+a similar measure can be derived on the single-cell level. Here, we will use
+a two component Gaussian mixture model for each marker to find cells
+with positive and negative expression. The SNR is defined as:
+
+$$SNR = I_s/I_n$$
+
+where $I_s$ is the intensity of the signal (mean intensity of cells with
+positive signal) and $I_n$ is the intensity of the noise (mean intensity of
+cells lacking expression). To define cells with positive and negative marker
+expression, we fit the mixture model across the transformed counts of all cells
+contained in the `SpatialExperiment` object. Next, for each marker we calculate
+the mean of the non-transformed counts for the positive and the negative cells.
+The SNR is then the ratio between the mean of the positive signal and the mean
+of the negative signal.
+
+
+```r
+library(mclust)
+
+set.seed(220224)
+mat <- sapply(seq_len(nrow(spe)), function(x){
+ cur_exprs <- assay(spe, "exprs")[x,]
+ cur_counts <- assay(spe, "counts")[x,]
+
+ cur_model <- Mclust(cur_exprs, G = 2)
+ mean1 <- mean(cur_counts[cur_model$classification == 1])
+ mean2 <- mean(cur_counts[cur_model$classification == 2])
+
+ signal <- ifelse(mean1 > mean2, mean1, mean2)
+ noise <- ifelse(mean1 > mean2, mean2, mean1)
+
+ return(c(snr = signal/noise, ps = signal))
+})
+
+cur_snr <- t(mat) %>% as.data.frame() %>%
+ mutate(marker = rownames(spe))
+
+cur_snr %>% ggplot() +
+ geom_point(aes(log2(ps), log2(snr))) +
+ geom_label_repel(aes(log2(ps), log2(snr), label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
+```
+
+
+
+Next, we observe the distributions of cell size across the individual images.
+Differences in cell size distributions can indicate segmentation biases due to
+differences in cell density or can indicate biological differences due to cell
+type compositions (tumor cells tend to be larger than immune cells).
+
+
+```r
+dittoPlot(spe, var = "area",
+ group.by = "sample_id",
+ plots = "boxplot") +
+ ylab("Cell area") + xlab("")
+```
+
+
+
+```r
+summary(spe$area)
+```
+
+```
+## Min. 1st Qu. Median Mean 3rd Qu. Max.
+## 3.00 47.00 70.00 76.38 98.00 466.00
+```
+
+The median cell size is 70 pixels with a median major axis
+length of 11.3. The largest cell
+has an area of 466 pixels which relates to a diameter of
+21.6 pixels assuming a circular shape.
+Overall, the distribution of cell sizes is similar across images with images from
+`Patient4_005` and `Patient4_007` showing a reduced average cell size. These
+images contain fewer tumor cells which can explain the smaller average cell size.
+
+We detect very small cells in the dataset and will remove them.
+The chosen threshold is arbitrary and needs to be adjusted per dataset.
+
+
+```r
+sum(spe$area < 5)
+```
+
+```
+## [1] 65
+```
+
+```r
+spe <- spe[,spe$area >= 5]
+```
+
+Another quality indicator can be an absolute measure of cell density often
+reported in cells per mm$^2$.
+
+
+```r
+cell_density <- colData(spe) %>%
+ as.data.frame() %>%
+ group_by(sample_id) %>%
+ summarize(cell_count = n(),
+ no_pixels = mean(width_px) * mean(height_px)) %>%
+ mutate(cells_per_mm2 = cell_count/(no_pixels/1000000))
+
+ggplot(cell_density) +
+ geom_point(aes(reorder(sample_id,cells_per_mm2), cells_per_mm2)) +
+ theme_minimal(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
+ ylab("Cells per mm2") + xlab("")
+```
+
+
+
+The number of cells per mm$^2$ varies across images which also depends on the
+number of tumor/non-tumor cells. As we can see in the following sections, some
+immune cells appear in cell dense regions while other stromal regions are less
+dense.
+
+The data presented here originate from samples from different locations with
+potential differences in pre-processing and each sample was stained individually.
+These (and other) technical aspects can induce staining differences between
+samples or batches of samples. Observing potential staining differences can be
+crucial to assess data quality. We will use ridgeline visualizations to check
+differences in staining patterns:
+
+
+```r
+multi_dittoPlot(spe, vars = rownames(spe)[rowData(spe)$use_channel],
+ group.by = "patient_id", plots = "ridgeplot",
+ assay = "exprs",
+ color.panel = metadata(spe)$color_vectors$patient_id)
+```
+
+
+
+We observe variations in the distributions of marker expression across patients.
+These variations may arise partly from different abundances of cells in
+different images (e.g., Patient3 may have higher numbers of CD11c+ and PD1+
+cells) as well as staining differences between samples. While most of the
+selected markers are specifically expressed in immune cell subtypes, we can see
+that E-Cadherin (a marker for epithelial (tumor) cells) shows a similar
+expression range across all patients.
+
+Finally, we will use non-linear dimensionality reduction methods to project
+cells from a high-dimensional (40) down to a low-dimensional (2) space. For this
+the *[scater](https://bioconductor.org/packages/3.17/scater)* package provides the `runUMAP` and
+`runTSNE` function. To ensure reproducibility, we will need to set a seed;
+however different seeds and different parameter settings (e.g., the `perplexity`
+parameter in the `runTSNE` function) need to be tested to avoid
+over-interpretation of visualization artefacts. For dimensionality reduction, we
+will use all channels that show biological variation across the dataset.
+However, marker selection can be performed with different biological questions
+in mind. Here, both the `runUMAP` and `runTSNE` function are not deterministic,
+meaning they produce different results across different runs. We therefore
+set a `seed` in this chunk for reproducibility purposes.
+
+
+```r
+library(scater)
+
+set.seed(220225)
+spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs")
+spe <- runTSNE(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs")
+```
+
+After dimensionality reduction, the low-dimensional embeddings are stored in the
+`reducedDim` slot.
+
+
+```r
+reducedDims(spe)
+```
+
+```
+## List of length 2
+## names(2): UMAP TSNE
+```
+
+```r
+head(reducedDim(spe, "UMAP"))
+```
+
+```
+## UMAP1 UMAP2
+## Patient1_001_1 -4.810167 -3.777362
+## Patient1_001_2 -4.397347 -3.456036
+## Patient1_001_3 -4.369883 -3.445561
+## Patient1_001_4 -4.081614 -3.162119
+## Patient1_001_5 -6.234012 -2.433976
+## Patient1_001_6 -5.666597 -3.428058
+```
+
+Visualization of the low-dimensional embedding facilitates assessment of
+potential "batch effects". The `dittoDimPlot`
+function allows flexible visualization. It returns `ggplot` objects which
+can be further modified.
+
+
+```r
+library(patchwork)
+
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP")
+p2 <- dittoDimPlot(spe, var = "patient_id", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on TSNE")
+
+# visualize region of interest id
+p3 <- dittoDimPlot(spe, var = "ROI", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$ROI) +
+ ggtitle("ROI ID on UMAP")
+p4 <- dittoDimPlot(spe, var = "ROI", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$ROI) +
+ ggtitle("ROI ID on TSNE")
+
+# visualize indication
+p5 <- dittoDimPlot(spe, var = "indication", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$indication) +
+ ggtitle("Indication on UMAP")
+p6 <- dittoDimPlot(spe, var = "indication", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$indication) +
+ ggtitle("Indication on TSNE")
+
+(p1 + p2) / (p3 + p4) / (p5 + p6)
+```
+
+
+
+
+```r
+# visualize marker expression
+p1 <- dittoDimPlot(spe, var = "Ecad", reduction.use = "UMAP",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "Ecad") +
+ ggtitle("E-Cadherin expression on UMAP")
+p2 <- dittoDimPlot(spe, var = "CD45RO", reduction.use = "UMAP",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "CD45RO") +
+ ggtitle("CD45RO expression on UMAP")
+p3 <- dittoDimPlot(spe, var = "Ecad", reduction.use = "TSNE",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "Ecad") +
+ ggtitle("Ecad expression on TSNE")
+p4 <- dittoDimPlot(spe, var = "CD45RO", reduction.use = "TSNE",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "CD45RO") +
+ ggtitle("CD45RO expression on TSNE")
+
+(p1 + p2) / (p3 + p4)
+```
+
+
+
+We observe a strong separation of tumor cells (Ecad+ cells) between the
+patients. Here, each patient was diagnosed with a different tumor type. The
+separation of tumor cells could be of biological origin since tumor cells tend
+to display differences in expression between patients and cancer types and/or of
+technical origin: the panel only contains a single tumor marker (E-Cadherin) and
+therefore slight technical differences in staining causes visible separation
+between cells of different patients. Nevertheless, the immune compartment
+(CD45RO+ cells) mix between patients and we can rule out systematic staining
+differences between patients.
+
+## Save objects
+
+The modified `SpatialExperiment` object is saved for further downstream analysis.
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+```
+
+
+
+## Session Info
+
+
+
+We observe an imperfect merging of Patient3 into all other samples. This
+was already seen when displaying the merging information above.
+We now also visualize the expression of selected markers across all cells
+before and after batch correction.
+
+
+```r
+markers <- c("Ecad", "CD45RO", "CD20", "CD3", "FOXP3", "CD206", "MPO", "SMA", "Ki67")
+
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+```r
+# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_mnnCorrected",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+We observe that immune cells across patients are merged after batch correction
+using `fastMNN`. However, the tumor cells of different patients still cluster
+separately.
+
+## harmony correction
+
+The `harmony` algorithm performs batch correction by iteratively clustering and
+correcting the positions of cells in PCA space [@Korsunsky2019]. We will first
+perform PCA on the asinh-transformed counts and then call the `RunHarmony`
+function to perform data integration.
+
+Paper: [Fast, sensitive and accurate integration of single-cell data with Harmony](https://www.nature.com/articles/s41592-019-0619-0)
+Documentation: [harmony](https://portals.broadinstitute.org/harmony/index.html)
+
+Similar to the `fastMNN` function, `harmony` returns the corrected
+low-dimensional coordinates for each cell. These can be transfered to the
+`reducedDim` slot of the original `SpatialExperiment` object.
+
+
+
+
+```r
+library(harmony)
+library(BiocSingular)
+
+spe <- runPCA(spe,
+ subset_row = rowData(spe)$use_channel,
+ exprs_values = "exprs",
+ ncomponents = 30,
+ BSPARAM = ExactParam())
+
+set.seed(230616)
+out <- RunHarmony(spe, group.by.vars = "patient_id")
+
+# Check that order of cells is the same
+stopifnot(all.equal(colnames(spe), colnames(out)))
+
+reducedDim(spe, "harmony") <- reducedDim(out, "HARMONY")
+```
+
+
+
+The computational time of the `HarmonyMatrix` function call is
+1.3 minutes.
+
+### Visualization
+
+We will now again visualize the cells in low dimensions after UMAP embedding.
+
+
+```r
+set.seed(220228)
+spe <- runUMAP(spe, dimred = "harmony", name = "UMAP_harmony")
+```
+
+
+```r
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP before correction")
+p2 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP_harmony", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP after correction")
+
+plot_grid(p1, p2)
+```
+
+
+
+And we visualize selected marker expression as defined above.
+
+
+```r
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+```r
+# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_harmony",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+We observe a more aggressive merging of cells from different patients compared
+to the results after `fastMNN` correction. Importantly, immune cell and epithelial
+markers are expressed in distinct regions of the UMAP.
+
+## Seurat correction
+
+The `Seurat` package provides a number of functionalities to analyze single-cell
+data. As such it also allows the integration of cells across different samples.
+Conceptually, `Seurat` performs batch correction similarly to `fastMNN` by
+finding mutual nearest neighbors (MNN) in low dimensional space before
+correcting the expression values of cells [@Stuart2019].
+
+Paper: [Comprehensive Integration of Single-Cell Data](https://www.cell.com/cell/fulltext/S0092-8674(19)30559-8)
+Documentation: [Seurat](https://satijalab.org/seurat/index.html)
+
+To use `Seurat`, we will first create a `Seurat` object from the `SpatialExperiment`
+object and add relevant metadata. The object also needs to be split by patient
+prior to integration.
+
+
+```r
+library(Seurat)
+library(SeuratObject)
+seurat_obj <- as.Seurat(spe, counts = "counts", data = "exprs")
+seurat_obj <- AddMetaData(seurat_obj, as.data.frame(colData(spe)))
+
+seurat.list <- SplitObject(seurat_obj, split.by = "patient_id")
+```
+
+To avoid long run times, we will use an approach that relies on reciprocal PCA
+instead of canonical correlation analysis for dimensionality reduction and
+initial alignment. For an extended tutorial on how to use `Seurat` for data
+integration, please refer to their
+[vignette](https://satijalab.org/seurat/articles/integration_rpca.html).
+
+We will first define the features used for integration and perform PCA on cells
+of each patient individually. The `FindIntegrationAnchors` function detects MNNs between
+cells of different patients and the `IntegrateData` function corrects the
+expression values of cells. We slightly increase the number of neighbors to be
+considered for MNN detection (the `k.anchor` parameter). This increases the integration
+strength.
+
+
+
+
+```r
+features <- rownames(spe)[rowData(spe)$use_channel]
+seurat.list <- lapply(X = seurat.list, FUN = function(x) {
+ x <- ScaleData(x, features = features, verbose = FALSE)
+ x <- RunPCA(x, features = features, verbose = FALSE, approx = FALSE)
+ return(x)
+})
+
+anchors <- FindIntegrationAnchors(object.list = seurat.list,
+ anchor.features = features,
+ reduction = "rpca",
+ k.anchor = 20)
+
+combined <- IntegrateData(anchorset = anchors)
+```
+
+
+
+We now select the `integrated` assay and perform PCA dimensionality reduction.
+The cell coordinates in PCA reduced space can then be transferred to the
+original `SpatialExperiment` object. **Of note:** by splitting the object into
+individual batch-specific objects, the ordering of cells in the integrated
+object might not match the ordering of cells in the input object. In this case,
+columns will need to be reordered. Here, we test if the ordering of cells in the
+integrated `Seurat` object matches the ordering of cells in the main
+`SpatialExperiment` object.
+
+
+```r
+DefaultAssay(combined) <- "integrated"
+
+combined <- ScaleData(combined, verbose = FALSE)
+combined <- RunPCA(combined, npcs = 30, verbose = FALSE, approx = FALSE)
+
+# Check that order of cells is the same
+stopifnot(all.equal(colnames(spe), colnames(combined)))
+
+reducedDim(spe, "seurat") <- Embeddings(combined, reduction = "pca")
+```
+
+The computational time of the `Seurat` function calls is
+4.29 minutes.
+
+### Visualization
+
+As above, we recompute the UMAP embeddings based on `Seurat` integrated results
+and visualize the embedding.
+
+
+```r
+set.seed(220228)
+spe <- runUMAP(spe, dimred = "seurat", name = "UMAP_seurat")
+```
+
+Visualize patient IDs.
+
+
+```r
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP before correction")
+p2 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP_seurat", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP after correction")
+
+plot_grid(p1, p2)
+```
+
+
+
+Visualization of marker expression.
+
+
+```r
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+```r
+# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_seurat",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
+```
+
+
+
+Similar to the methods presented above, `Seurat` integrates immune cells correctly.
+When visualizing the patient IDs, slight patient-to-patient differences within tumor
+cells can be detected.
+
+Choosing the correct integration approach is challenging without having ground truth
+cell labels available. It is recommended to compare different techniques and different
+parameter settings. Please refer to the documentation of the individual tools
+to become familiar with the possible parameter choices. Furthermore, in the following
+section, we will discuss clustering and classification approaches in light of
+expression differences between samples.
+
+In general, it appears that MNN-based approaches are less conservative in terms
+of merging compared to `harmony`. On the other hand, `harmony` could well merge
+cells in a way that regresses out biological signals.
+
+## Save objects
+
+The modified `SpatialExperiment` object is saved for further downstream analysis.
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+```
+
+
+
+## Session Info
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("pg_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$pg_clusters))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+
+
+The `Rphenograph` function call took
+2.31 minutes.
+
+We can observe that some of the clusters only contain cells of a single patient.
+This can often be observed in the tumor compartment. In the next step, we
+use the integrated cells (see Section \@ref(batch-effects)) in low dimensional
+embedding for clustering. Here, the low dimensional embedding can
+be directly accessed from the `reducedDim` slot.
+
+
+```r
+mat <- reducedDim(spe, "fastMNN")
+
+set.seed(230619)
+out <- Rphenograph(mat, k = 45)
+
+clusters <- factor(membership(out[[2]]))
+
+spe$pg_clusters_corrected <- clusters
+
+dittoDimPlot(spe, var = "pg_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("Phenograph clusters on UMAP, integrated cells")
+```
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("pg_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$pg_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+Clustering using the integrated embedding leads to clusters that contain cells
+of different patients. Cluster annotation can now be performed by manually
+labeling cells based on their marker expression (see Notes in Section
+\@ref(clustering-notes)).
+
+### Shared nearest neighbour graph {#snn-graph}
+
+The [bluster](https://www.bioconductor.org/packages/release/bioc/html/bluster.html)
+package provides a simple interface to cluster cells using a number of different
+[clustering approaches](https://www.bioconductor.org/packages/release/bioc/vignettes/bluster/inst/doc/clusterRows.html) and different metrics to [access cluster stability](https://www.bioconductor.org/packages/release/bioc/vignettes/bluster/inst/doc/diagnostics.html).
+
+For simplicity, we will focus on graph based clustering as this is the most
+popular and a fast method for single-cell clustering. The `bluster` package
+provides functionalities to build k-nearest neighbor (KNN) graphs and its weighted
+version, shared nearest neighbor (SNN) graphs where nodes represent cells.
+The user can chose the number of neighbors to consider (parameter `k`),
+the edge weighting method (parameter `type`) and the community detection
+function to use (parameter `cluster.fun`). As all parameters affect the clustering
+results, the `bluster` package provides the `clusterSweep` function to test
+a number of parameter settings in parallel. In the following code chunk,
+we select the asinh-transformed mean pixel intensities and subset the markers
+of interest. The resulting matrix is transposed to fit to the requirements of
+the bluster package (cells in rows).
+
+We test two different settings for `k`, two for `type` and fix the `cluster.fun`
+to `louvain` as this is one of the most common approaches for community detection.
+This function call is parallelized by setting the `BPPARAM` parameter.
+
+
+
+
+```r
+library(bluster)
+library(BiocParallel)
+library(ggplot2)
+
+mat <- t(assay(spe, "exprs")[rowData(spe)$use_channel,])
+
+combinations <- clusterSweep(mat,
+ BLUSPARAM=SNNGraphParam(),
+ k=c(10L, 20L),
+ type = c("rank", "jaccard"),
+ cluster.fun = "louvain",
+ BPPARAM = MulticoreParam(RNGseed = 220427))
+```
+
+We next calculate two metrics to estimate cluster stability: the average
+silhouette width and the neighborhood purity.
+
+We use the `approxSilhouette` function to compute the silhouette width for each
+cell and compute the average across all cells per parameter setting. Please see
+`?silhouette` for more information on how the silhouette width is computed for
+each cell. A large average silhouette width indicates a cluster parameter
+setting for which cells that are well clustered.
+
+The `neighborPurity` function computes the fraction of cells around each cell
+with the same cluster ID. Per parameter setting, we compute the average
+neighborhood purity across all cells. A large average neighborhood purity
+indicates a cluster parameter setting for which cells that are well clustered.
+
+
+```r
+sil <- vapply(as.list(combinations$clusters),
+ function(x) mean(approxSilhouette(mat, x)$width),
+ 0)
+
+ggplot(data.frame(method = names(sil),
+ sil = sil)) +
+ geom_point(aes(method, sil)) +
+ theme_classic(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+ xlab("Cluster parameter combination") +
+ ylab("Average silhouette width")
+```
+
+
+
+```r
+pur <- vapply(as.list(combinations$clusters),
+ function(x) mean(neighborPurity(mat, x)$purity),
+ 0)
+
+ggplot(data.frame(method = names(pur),
+ pur = pur)) +
+ geom_point(aes(method, pur)) +
+ theme_classic(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+ xlab("Cluster parameter combination") +
+ ylab("Average neighborhood purity")
+```
+
+
+
+
+
+The cluster parameter sweep took
+8.81 minutes.
+
+Performing a cluster sweep takes some time as multiple function calls are run in parallel.
+We do however recommend testing a number of different parameter settings to
+assess clustering performance.
+
+Once parameter settings are known, we can either use the `clusterRows` function
+of the `bluster` package to cluster cells or its convenient wrapper function
+exported by the
+[scran](https://bioconductor.org/packages/release/bioc/html/scran.html) package.
+The `scran::clusterCells` function accepts a `SpatialExperiment` (or
+`SingleCellExperiment`) object which stores cells in columns. By default, the
+function detects the 10 nearest neighbours for each cell, performs rank-based
+weighting of edges (see `?makeSNNGraph` for more information) and uses the
+`cluster_walktrap` function to detect communities in the graph.
+
+As we can see above, the clustering approach in this dataset with `k` being 20
+and rank-based edge weighting leads to the highest silhouette width and highest
+neighborhood purity.
+
+
+
+
+```r
+library(scran)
+
+set.seed(220620)
+clusters <- clusterCells(spe[rowData(spe)$use_channel,],
+ assay.type = "exprs",
+ BLUSPARAM = SNNGraphParam(k=20,
+ cluster.fun = "louvain",
+ type = "rank"))
+
+spe$nn_clusters <- clusters
+
+dittoDimPlot(spe, var = "nn_clusters",
+ reduction.use = "UMAP", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SNN clusters on UMAP")
+```
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("nn_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$nn_clusters))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+
+
+The shared nearest neighbor graph clustering approach took
+1.31 minutes.
+
+This function was used by [@Tietscher2022] to cluster cells obtained by IMC. Setting
+`type = "jaccard"` performs clustering similar to `Rphenograph` above and [Seurat](https://satijalab.org/seurat/articles/pbmc3k_tutorial.html#cluster-the-cells-1).
+
+Similar to the results obtained by `Rphenograph`, some of the clusters are
+patient-specific. We can now perform clustering of the integrated cells
+by directly specifying which low-dimensional embedding to use:
+
+
+```r
+set.seed(220621)
+clusters <- clusterCells(spe,
+ use.dimred = "fastMNN",
+ BLUSPARAM = SNNGraphParam(k = 20,
+ cluster.fun = "louvain",
+ type = "rank"))
+
+spe$nn_clusters_corrected <- clusters
+
+dittoDimPlot(spe, var = "nn_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SNN clusters on UMAP, integrated cells")
+```
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("nn_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$nn_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+### Self organizing maps
+
+An alternative to graph-based clustering is offered by the
+[CATALYST](https://bioconductor.org/packages/release/bioc/html/CATALYST.html)
+package. The `cluster` function internally uses the
+[FlowSOM](https://bioconductor.org/packages/release/bioc/html/FlowSOM.html)
+package to group cells into 100 (default) clusters based on self organizing maps
+(SOM). In the next step, the
+[ConsensusClusterPlus](https://bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html)
+package is used to perform hierarchical consensus clustering of the previously
+detected 100 SOM nodes into 2 to `maxK` clusters. Cluster stability for each `k`
+can be assessed by plotting the `delta_area(spe)`. The optimal number
+of clusters can be found by selecting the `k` at which a plateau is reached.
+In the example below, an optimal `k` lies somewhere around 13.
+
+
+
+
+```r
+library(CATALYST)
+
+# Run FlowSOM and ConsensusClusterPlus clustering
+set.seed(220410)
+spe <- cluster(spe,
+ features = rownames(spe)[rowData(spe)$use_channel],
+ maxK = 30)
+
+# Assess cluster stability
+delta_area(spe)
+```
+
+
+
+```r
+spe$som_clusters <- cluster_ids(spe, "meta13")
+
+dittoDimPlot(spe, var = "som_clusters",
+ reduction.use = "UMAP", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SOM clusters on UMAP")
+```
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("som_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$som_clusters))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+
+
+Running FlowSOM clustering took 0.22 minutes.
+
+The `CATALYST` package does not provide functionality to perform `FlowSOM` and
+`ConsensusClusterPlus` clustering directly on the batch-corrected, integrated cells. As an
+alternative to the `CATALYST` package, the `bluster` package provides SOM
+clustering when specifying the `SomParam()` parameter. Similar to the `CATALYST`
+approach, we will first cluster the dataset into 100 clusters (also called
+"codes"). These codes are then further clustered into a maximum of 30 clusters
+using `ConsensusClusterPlus` (using hierarchical clustering and euclidean
+distance). The delta area plot can be accessed using the (not exported)
+`.plot_delta_area` function from `CATALYST`. Here, it seems that the plateau is
+reached at a `k` of 16 and we will store the final cluster IDs within the
+`SpatialExperiment` object.
+
+
+```r
+library(kohonen)
+library(ConsensusClusterPlus)
+
+# Select integrated cells
+mat <- reducedDim(spe, "fastMNN")
+
+# Perform SOM clustering
+set.seed(220410)
+som.out <- clusterRows(mat, SomParam(100), full = TRUE)
+
+# Cluster the 100 SOM codes into larger clusters
+ccp <- ConsensusClusterPlus(t(som.out$objects$som$codes[[1]]),
+ maxK = 30,
+ reps = 100,
+ distance = "euclidean",
+ seed = 220410,
+ plot = NULL)
+```
+
+
+```r
+# Visualize delta area plot
+CATALYST:::.plot_delta_area(ccp)
+```
+
+
+
+```r
+# Link ConsensusClusterPlus clusters with SOM codes and save in object
+som.cluster <- ccp[[16]][["consensusClass"]][som.out$clusters]
+spe$som_clusters_corrected <- as.factor(som.cluster)
+
+dittoDimPlot(spe, var = "som_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SOM clusters on UMAP, integrated cells")
+```
+
+
+
+
+```r
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("som_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$som_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
+```
+
+
+
+The `FlowSOM` clustering approach has been used by [@Hoch2022] to sub-cluster tumor
+cells as measured by IMC.
+
+### Compare between clustering approaches
+
+Finally, we can compare the results of different clustering approaches. For
+this, we visualize the number of cells that are shared between different
+clustering results in a pairwise fashion. In the following heatmaps a high match
+between clustering results can be seen for those clusters that are uniquely
+detected in both approaches.
+
+First, we will visualize the match between the three different approaches
+applied to the asinh-transformed counts.
+
+
+```r
+library(patchwork)
+library(pheatmap)
+library(gridExtra)
+
+tab1 <- table(paste("Rphenograph", spe$pg_clusters),
+ paste("SNN", spe$nn_clusters))
+tab2 <- table(paste("Rphenograph", spe$pg_clusters),
+ paste("SOM", spe$som_clusters))
+tab3 <- table(paste("SNN", spe$nn_clusters),
+ paste("SOM", spe$som_clusters))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab2 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab3 + 10), color = viridis(100))
+```
+
+
+
+We observe that `Rphenograph` and the shared nearest neighbor (SNN) approach by
+`scran` show similar results (first heatmap above). For example, Rphenograph
+cluster 20 (a tumor cluster) is perfectly captured by SNN cluster 12. On the
+other hand, the Neutrophil cluster (SNN cluster 6) is split into Rphenograph
+cluster 2 and Rphenograph cluster 6. A common approach
+is to now merge clusters that contain similar cell types and annotate them by
+hand (see below).
+
+Below, a comparison between the clustering results of the integrated cells
+is shown.
+
+
+```r
+tab1 <- table(paste("Rphenograph", spe$pg_clusters_corrected),
+ paste("SNN", spe$nn_clusters_corrected))
+tab2 <- table(paste("Rphenograph", spe$pg_clusters_corrected),
+ paste("SOM", spe$som_clusters_corrected))
+tab3 <- table(paste("SNN", spe$nn_clusters_corrected),
+ paste("SOM", spe$som_clusters_corrected))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab2 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab3 + 10), color = viridis(100))
+```
+
+
+
+In comparison to clustering on the non-integrated cells, the clustering results
+of the integrated cells show higher overlap. The SNN approach resulted in fewer
+clusters and therefore matches better with the SOM clustering approach.
+
+### Further clustering notes {#clustering-notes}
+
+The `bluster` package provides a number of metrics to assess cluster stability
+[here](https://www.bioconductor.org/packages/release/bioc/vignettes/bluster/inst/doc/diagnostics.html).
+For brevity we only highlighted the use of the silhouette width and the
+neighborhood purity but different metrics should be tested to assess cluster
+stability.
+
+To assign cell types to clusters, we manually annotate clusters based on their
+marker expression. For example, SNN cluster 12 (clustering of the integrated
+cells) shows high, homogeneous expression of CD20 and we might therefore label
+this cluster as B cells. The next chapter \@ref(single-cell-visualization) will
+highlight single-cell visualization methods that can be helpful for manual
+cluster annotations.
+
+An example how to label clusters can be seen below:
+
+
+```r
+library(dplyr)
+cluster_celltype <- recode(spe$nn_clusters_corrected,
+ "1" = "Tumor_proliferating",
+ "2" = "Myeloid",
+ "3" = "Tumor",
+ "4" = "Tumor",
+ "5" = "Stroma",
+ "6" = "Proliferating",
+ "7" = "Myeloid",
+ "8" = "Plasma_cell",
+ "9" = "CD8",
+ "10" = "CD4",
+ "11" = "Neutrophil",
+ "12" = "Bcell",
+ "13" = "Stroma")
+
+spe$cluster_celltype <- cluster_celltype
+```
+
+## Classification approach {#classification}
+
+In this section, we will highlight a cell type classification approach based
+on ground truth labeling and random forest classification. The rational for
+this supervised cell phenotyping approach is to use the information contained
+in the pre-defined markers to detect cells of interest. This approach was
+used by Hoch _et al._ to classify cell types in a metastatic melanoma IMC
+dataset [@Hoch2022].
+
+The antibody panel used in the example data set mainly focuses on immune cell
+types and little on tumor cell phenotypes. Therefore we will label the following
+cell types:
+
+* Tumor (E-cadherin positive)
+* Stroma (SMA, PDGFRb positive)
+* Plasma cells (CD38 positive)
+* Neutrophil (MPO, CD15 positive)
+* Myeloid cells (HLADR positive)
+* B cells (CD20 positive)
+* B next to T cells (CD20, CD3 positive)
+* Regulatory T cells (FOXP3 positive)
+* CD8+ T cells (CD3, CD8 positive)
+* CD4+ T cells (CD3, CD4 positive)
+
+The "B next to T cell" phenotype (`BnTcell`) is commonly observed in immune
+infiltrated regions measured by IMC. We include this phenotype to account for B
+cell/T cell interactions where precise classification into B cells or T cells is
+not possible. The exact gating scheme can be seen at
+[img/Gating_scheme.pdf](img/Gating_scheme.pdf).
+
+As related approaches, [Astir](https://github.com/camlab-bioml/astir) and
+[Garnett](https://cole-trapnell-lab.github.io/garnett/) use pre-defined panel
+information to classify cell phenotypes based on their marker expression.
+
+### Manual labeling of cells
+
+The [cytomapper](https://www.bioconductor.org/packages/release/bioc/html/cytomapper.html)
+package provides the `cytomapperShiny` function that allows gating of cells
+based on their marker expression and visualization of selected cells directly
+on the images.
+
+
+```r
+library(cytomapper)
+if (interactive()) {
+
+ images <- readRDS("data/images.rds")
+ masks <- readRDS("data/masks.rds")
+
+ cytomapperShiny(object = spe, mask = masks, image = images,
+ cell_id = "ObjectNumber", img_id = "sample_id")
+}
+```
+
+The labeled cells for this data set can be accessed at
+[10.5281/zenodo.6554544](https://zenodo.org/record/6554544) and were downloaded
+in Section \@ref(prerequisites). Gating is performed per image and the
+`cytomapperShiny` function allows the export of gated cells in form of a
+`SingleCellExperiment` or `SpatialExperiment` object. The cell label is stored
+in `colData(object)$cytomapper_CellLabel` and the gates are stored in
+`metadata(object)`. In the next section, we will read in and consolidate the
+labeled data.
+
+### Define color vectors
+
+For consistent visualization of cell types, we will now pre-define their colors:
+
+
+```r
+celltype <- setNames(c("#3F1B03", "#F4AD31", "#894F36", "#1C750C", "#EF8ECC",
+ "#6471E2", "#4DB23B", "grey", "#F4800C", "#BF0A3D", "#066970"),
+ c("Tumor", "Stroma", "Myeloid", "CD8", "Plasma_cell",
+ "Treg", "CD4", "undefined", "BnTcell", "Bcell", "Neutrophil"))
+
+metadata(spe)$color_vectors$celltype <- celltype
+```
+
+### Read in and consolidate labeled data
+
+Here, we will read in the individual `SpatialExperiment` objects containing the
+labeled cells and concatenate them. In the process of concatenating the
+`SpatialExperiment` objects along their columns, the `sample_id` entry is
+appended by `.1, .2, .3, ...` due to replicated entries.
+
+
+```r
+library(SingleCellExperiment)
+label_files <- list.files("data/gated_cells",
+ full.names = TRUE, pattern = ".rds$")
+
+# Read in SPE objects
+spes <- lapply(label_files, readRDS)
+
+# Merge SPE objects
+concat_spe <- do.call("cbind", spes)
+```
+
+In the following code chunk we will identify cells that were labeled multiple
+times. This occurs when different cell phenotypes are gated per image and can
+affect immune cells that are located inside the tumor compartment.
+
+We will first identify those cells that were uniquely labeled. In the next step,
+we will identify those cells that were labeled twice AND were labeled as Tumor
+cells. These cells will be assigned their immune cell label. Finally, we will
+save the unique labels within the original `SpatialExperiment` object.
+
+**Of note:** this concatenation strategy is specific for cell phenotypes contained in this
+example dataset. The gated cell labels might need to be processed in a slightly
+different way when working with other samples.
+
+For these tasks, we will define a filter function:
+
+
+```r
+filter_labels <- function(object,
+ label = "cytomapper_CellLabel") {
+ cur_tab <- unclass(table(colnames(object), object[[label]]))
+
+ cur_labels <- colnames(cur_tab)[apply(cur_tab, 1, which.max)]
+ names(cur_labels) <- rownames(cur_tab)
+
+ cur_labels <- cur_labels[rowSums(cur_tab) == 1]
+
+ return(cur_labels)
+}
+```
+
+This function is now applied to all cells and then only non-tumor cells.
+
+
+```r
+labels <- filter_labels(concat_spe)
+
+cur_spe <- concat_spe[,concat_spe$cytomapper_CellLabel != "Tumor"]
+
+non_tumor_labels <- filter_labels(cur_spe)
+
+additional_cells <- setdiff(names(non_tumor_labels), names(labels))
+
+final_labels <- c(labels, non_tumor_labels[additional_cells])
+
+# Transfer labels to SPE object
+spe_labels <- rep("unlabeled", ncol(spe))
+names(spe_labels) <- colnames(spe)
+spe_labels[names(final_labels)] <- final_labels
+spe$cell_labels <- spe_labels
+
+# Number of cells labeled per patient
+table(spe$cell_labels, spe$patient_id)
+```
+
+```
+##
+## Patient1 Patient2 Patient3 Patient4
+## Bcell 152 131 234 263
+## BnTcell 396 37 240 1029
+## CD4 45 342 167 134
+## CD8 60 497 137 128
+## Myeloid 183 378 672 517
+## Neutrophil 97 4 17 16
+## Plasma_cell 34 536 87 59
+## Stroma 84 37 85 236
+## Treg 139 149 49 24
+## Tumor 2342 906 1618 1133
+## unlabeled 7214 9780 7826 9580
+```
+
+Based on these labels, we can now train a random forest classifier to classify
+all remaining, unlabeled cells.
+
+### Train classifier
+
+In this section, we will use the
+[caret](https://topepo.github.io/caret/index.html) framework for machine
+learning in R. This package provides an interface to train a number of
+regression and classification models in a coherent fashion. We use a random
+forest classifier due to low number of parameters, high speed and an observed
+high performance for cell type classification [@Hoch2022].
+
+In the following section, we will first split the `SpatialExperiment` object
+into labeled and unlabeled cells. Based on the labeled cells, we split
+the data into a train (75% of the data) and test (25% of the data) dataset.
+We currently do not provide an independently labeled validation dataset.
+
+The `caret` package provides the `trainControl` function, which specifies model
+training parameters and the `train` function, which performs the actual model
+training. While training the model, we also want to estimate the best model
+parameters. In the case of the chosen random forest model (`method = "rf"`), we
+only need to estimate a single parameters (`mtry`) which corresponds to the
+number of variables randomly sampled as candidates at each split. To estimate
+the best parameter, we will perform a 5-fold cross validation (set within
+`trainControl`) over a tune length of 5 entries to `mtry`. In the following
+code chunk, the `createDataPartition` and the `train` function are not deterministic,
+meaning they return different results across different runs. We therefore set
+a `seed` here for both functions.
+
+
+
+
+```r
+library(caret)
+
+# Split between labeled and unlabeled cells
+lab_spe <- spe[,spe$cell_labels != "unlabeled"]
+unlab_spe <- spe[,spe$cell_labels == "unlabeled"]
+
+# Randomly split into train and test data
+set.seed(221029)
+trainIndex <- createDataPartition(factor(lab_spe$cell_labels), p = 0.75)
+
+train_spe <- lab_spe[,trainIndex$Resample1]
+test_spe <- lab_spe[,-trainIndex$Resample1]
+
+# Define fit parameters for 5-fold cross validation
+fitControl <- trainControl(method = "cv",
+ number = 5)
+
+# Select the arsinh-transformed counts for training
+cur_mat <- t(assay(train_spe, "exprs")[rowData(train_spe)$use_channel,])
+
+# Train a random forest classifier
+rffit <- train(x = cur_mat,
+ y = factor(train_spe$cell_labels),
+ method = "rf", ntree = 1000,
+ tuneLength = 5,
+ trControl = fitControl)
+
+rffit
+```
+
+```
+## Random Forest
+##
+## 10049 samples
+## 37 predictor
+## 10 classes: 'Bcell', 'BnTcell', 'CD4', 'CD8', 'Myeloid', 'Neutrophil', 'Plasma_cell', 'Stroma', 'Treg', 'Tumor'
+##
+## No pre-processing
+## Resampling: Cross-Validated (5 fold)
+## Summary of sample sizes: 8040, 8039, 8038, 8038, 8041
+## Resampling results across tuning parameters:
+##
+## mtry Accuracy Kappa
+## 2 0.9643726 0.9524051
+## 10 0.9780071 0.9707483
+## 19 0.9801973 0.9736577
+## 28 0.9787052 0.9716635
+## 37 0.9779095 0.9705890
+##
+## Accuracy was used to select the optimal model using the largest value.
+## The final value used for the model was mtry = 19.
+```
+
+
+
+Training the classifier took
+11.77 minutes.
+
+### Classifier performance
+
+We next observe the accuracy of the classifer when predicting cell phenotypes
+across the cross-validation and when applying the classifier to the test
+dataset.
+
+First, we can visualize the classification accuracy during parameter
+tuning:
+
+
+```r
+ggplot(rffit) +
+ geom_errorbar(data = rffit$results,
+ aes(ymin = Accuracy - AccuracySD,
+ ymax = Accuracy + AccuracySD),
+ width = 0.4) +
+ theme_classic(base_size = 15)
+```
+
+
+
+The best value for `mtry` is 19 and is used when predicting new data.
+
+It is often recommended to visualize the variable importance of the
+classifier. The following plot specifies which variables (markers) are
+most important for classifying the data.
+
+
+```r
+plot(varImp(rffit))
+```
+
+
+
+As expected, the markers that were used for gating (Ecad, CD3, CD20, HLADR,
+CD8a, CD38, FOXP3) were important for classification.
+
+To assess the accuracy, sensitivity, specificity, among other quality measures of
+the classifier, we will now predict cell phenotypes in the test data.
+
+
+```r
+# Select the arsinh-transformed counts of the test data
+cur_mat <- t(assay(test_spe, "exprs")[rowData(test_spe)$use_channel,])
+
+# Predict the cell phenotype labels of the test data
+set.seed(231019)
+cur_pred <- predict(rffit, newdata = cur_mat)
+```
+
+While the overall classification accuracy can appear high, we also want
+to check if each cell phenotype class is correctly predicted.
+For this, we will calculate the confusion matrix between predicted and actual
+cell labels. This measure may highlight individual cell phenotype classes that
+were not correctly predicted by the classifier. When setting `mode = "everything"`,
+the `confusionMatrix` function returns all available prediction measures including
+sensitivity, specificity, precision, recall and the F1 score per cell
+phenotype class.
+
+
+```r
+cm <- confusionMatrix(data = cur_pred,
+ reference = factor(test_spe$cell_labels),
+ mode = "everything")
+
+cm
+```
+
+```
+## Confusion Matrix and Statistics
+##
+## Reference
+## Prediction Bcell BnTcell CD4 CD8 Myeloid Neutrophil Plasma_cell Stroma
+## Bcell 186 2 0 0 0 0 6 0
+## BnTcell 4 423 1 0 0 0 0 0
+## CD4 0 0 163 0 0 2 3 2
+## CD8 0 0 0 199 0 0 8 0
+## Myeloid 0 0 2 1 437 0 0 0
+## Neutrophil 0 0 0 0 0 30 0 0
+## Plasma_cell 1 0 3 2 0 0 158 0
+## Stroma 0 0 2 0 0 0 0 108
+## Treg 0 0 0 0 0 0 3 0
+## Tumor 4 0 1 3 0 1 1 0
+## Reference
+## Prediction Treg Tumor
+## Bcell 0 1
+## BnTcell 0 1
+## CD4 0 5
+## CD8 0 3
+## Myeloid 0 0
+## Neutrophil 0 0
+## Plasma_cell 1 0
+## Stroma 0 0
+## Treg 89 2
+## Tumor 0 1487
+##
+## Overall Statistics
+##
+## Accuracy : 0.9806
+## 95% CI : (0.9753, 0.985)
+## No Information Rate : 0.4481
+## P-Value [Acc > NIR] : < 2.2e-16
+##
+## Kappa : 0.9741
+##
+## Mcnemar's Test P-Value : NA
+##
+## Statistics by Class:
+##
+## Class: Bcell Class: BnTcell Class: CD4 Class: CD8
+## Sensitivity 0.95385 0.9953 0.94767 0.97073
+## Specificity 0.99714 0.9979 0.99622 0.99650
+## Pos Pred Value 0.95385 0.9860 0.93143 0.94762
+## Neg Pred Value 0.99714 0.9993 0.99716 0.99809
+## Precision 0.95385 0.9860 0.93143 0.94762
+## Recall 0.95385 0.9953 0.94767 0.97073
+## F1 0.95385 0.9906 0.93948 0.95904
+## Prevalence 0.05830 0.1271 0.05142 0.06129
+## Detection Rate 0.05561 0.1265 0.04873 0.05949
+## Detection Prevalence 0.05830 0.1283 0.05232 0.06278
+## Balanced Accuracy 0.97549 0.9966 0.97195 0.98361
+## Class: Myeloid Class: Neutrophil Class: Plasma_cell
+## Sensitivity 1.0000 0.909091 0.88268
+## Specificity 0.9990 1.000000 0.99779
+## Pos Pred Value 0.9932 1.000000 0.95758
+## Neg Pred Value 1.0000 0.999095 0.99340
+## Precision 0.9932 1.000000 0.95758
+## Recall 1.0000 0.909091 0.88268
+## F1 0.9966 0.952381 0.91860
+## Prevalence 0.1306 0.009865 0.05351
+## Detection Rate 0.1306 0.008969 0.04723
+## Detection Prevalence 0.1315 0.008969 0.04933
+## Balanced Accuracy 0.9995 0.954545 0.94024
+## Class: Stroma Class: Treg Class: Tumor
+## Sensitivity 0.98182 0.98889 0.9920
+## Specificity 0.99938 0.99846 0.9946
+## Pos Pred Value 0.98182 0.94681 0.9933
+## Neg Pred Value 0.99938 0.99969 0.9935
+## Precision 0.98182 0.94681 0.9933
+## Recall 0.98182 0.98889 0.9920
+## F1 0.98182 0.96739 0.9927
+## Prevalence 0.03288 0.02691 0.4481
+## Detection Rate 0.03229 0.02661 0.4445
+## Detection Prevalence 0.03288 0.02810 0.4475
+## Balanced Accuracy 0.99060 0.99368 0.9933
+```
+
+To easily visualize these results, we can now plot the true positive rate
+(sensitivity) versus the false positive rate (1 - specificity). The size of the
+point is determined by the number of true positives divided by the total number
+of cells.
+
+
+```r
+library(tidyverse)
+
+data.frame(cm$byClass) %>%
+ mutate(class = sub("Class: ", "", rownames(cm$byClass))) %>%
+ ggplot() +
+ geom_point(aes(1 - Specificity, Sensitivity,
+ size = Detection.Rate,
+ fill = class),
+ shape = 21) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme_classic(base_size = 15) +
+ ylab("Sensitivity (TPR)") +
+ xlab("1 - Specificity (FPR)")
+```
+
+
+
+We observe high sensitivity and specificity for most cell types. Plasma cells
+show the lowest true positive rate with 88% being sufficiently high.
+
+Finally, to observe which cell phenotypes were wrongly classified, we can visualize
+the distribution of classification probabilities per cell phenotype class:
+
+
+```r
+set.seed(231019)
+cur_pred <- predict(rffit,
+ newdata = cur_mat,
+ type = "prob")
+cur_pred$truth <- factor(test_spe$cell_labels)
+
+cur_pred %>%
+ pivot_longer(cols = Bcell:Tumor) %>%
+ ggplot() +
+ geom_boxplot(aes(x = name, y = value, fill = name), outlier.size = 0.5) +
+ facet_wrap(. ~ truth, ncol = 1) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme(panel.background = element_blank(),
+ axis.text.x = element_text(angle = 45, hjust = 1))
+```
+
+
+
+The boxplots indicate the classification probabilities per class. The classifier
+is well trained if classification probabilities are only high for the one
+specific class.
+
+### Classification of new data
+
+In the final section, we will now use the tuned and tested random forest
+classifier to predict the cell phenotypes of the unlabeled data.
+
+First, we predict the cell phenotypes and extract their classification
+probabilities.
+
+
+```r
+# Select the arsinh-transformed counts of the unlabeled data for prediction
+cur_mat <- t(assay(unlab_spe, "exprs")[rowData(unlab_spe)$use_channel,])
+
+# Predict the cell phenotype labels of the unlabeled data
+set.seed(231014)
+cell_class <- as.character(predict(rffit,
+ newdata = cur_mat,
+ type = "raw"))
+names(cell_class) <- rownames(cur_mat)
+
+table(cell_class)
+```
+
+```
+## cell_class
+## Bcell BnTcell CD4 CD8 Myeloid Neutrophil
+## 817 979 3620 2716 6302 559
+## Plasma_cell Stroma Treg Tumor
+## 2692 4904 1170 10641
+```
+
+```r
+# Extract prediction probabilities for each cell
+set.seed(231014)
+cell_prob <- predict(rffit,
+ newdata = cur_mat,
+ type = "prob")
+```
+
+Each cell is assigned to the class with highest probability. There are however
+cases, where the highest probability is low meaning the cell can not be uniquely
+assigned to a class. We next want to identify these cells and label them as
+"undefined". Here, we select a maximum classification probability threshold
+of 40% but this threshold needs to be adjusted for other datasets. The adjusted
+cell labels are then stored in the `SpatialExperiment` object.
+
+
+```r
+library(ggridges)
+
+# Distribution of maximum probabilities
+tibble(max_prob = rowMax(as.matrix(cell_prob)),
+ type = cell_class) %>%
+ ggplot() +
+ geom_density_ridges(aes(x = max_prob, y = cell_class, fill = cell_class)) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme_classic(base_size = 15) +
+ xlab("Maximum probability") +
+ ylab("Cell type") +
+ xlim(c(0,1.2))
+```
+
+```
+## Picking joint bandwidth of 0.0238
+```
+
+
+
+```r
+# Label undefined cells
+cell_class[rowMax(as.matrix(cell_prob)) < 0.4] <- "undefined"
+
+# Store labels in SpatialExperiment onject
+cell_labels <- spe$cell_labels
+cell_labels[colnames(unlab_spe)] <- cell_class
+spe$celltype <- cell_labels
+
+table(spe$celltype, spe$patient_id)
+```
+
+```
+##
+## Patient1 Patient2 Patient3 Patient4
+## Bcell 179 527 431 458
+## BnTcell 416 586 594 1078
+## CD4 391 1370 699 1385
+## CD8 518 1365 479 1142
+## Myeloid 1369 2197 1723 2731
+## Neutrophil 348 9 148 176
+## Plasma_cell 650 2122 351 274
+## Stroma 633 676 736 3261
+## Treg 553 409 243 310
+## Tumor 5560 3334 5648 2083
+## undefined 129 202 80 221
+```
+
+We can now compare the cell labels derived by classification to the different
+clustering strategies. The first comparison is against the clustering results
+using the asinh-transformed counts.
+
+
+```r
+tab1 <- table(spe$celltype,
+ paste("Rphenograph", spe$pg_clusters))
+tab2 <- table(spe$celltype,
+ paste("SNN", spe$nn_clusters))
+tab3 <- table(spe$celltype,
+ paste("SOM", spe$som_clusters))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab2 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab3 + 10), color = viridis(100))
+```
+
+
+
+We can see that Tumor and Myeloid cells span multiple clusters while
+Neutrophiles are detected as an individual cluster by all clustering approaches.
+
+We next compare the cell classification against clustering results using the
+integrated cells.
+
+
+```r
+tab1 <- table(spe$celltype,
+ paste("Rphenograph", spe$pg_clusters_corrected))
+tab2 <- table(spe$celltype,
+ paste("SNN", spe$nn_clusters_corrected))
+tab3 <- table(spe$celltype,
+ paste("SOM", spe$som_clusters_corrected))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab2 + 10), color = viridis(100))
+```
+
+
+
+```r
+pheatmap(log10(tab3 + 10), color = viridis(100))
+```
+
+
+
+We observe a high agreement between the shared nearest neighbor clustering
+approach using the integrated cells and the cell phenotypes derived by
+classification.
+
+In the next sections, we will highlight visualization strategies to verify the
+correctness of the phenotyping approach. Specifically, Section
+\@ref(outline-cells) shows how to outline identified cell phenotypes on
+composite images.
+
+
+
+Finally, we save the updated `SpatialExperiment` object.
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+```
+
+
+
+## Session Info
+
+
+
+The `plotReducedDim` function of the `scater` package provides an alternative
+way for visualizing cells in low dimensions. Here, we loop over all type
+markers, generate one plot per marker and plot the indivudual plots side-by-side.
+
+
+```r
+# UMAP colored by expression for all markers - plotReducedDim
+plot_list <- lapply(rownames(spe)[rowData(spe)$marker_class == "type"], function(x){
+ p <- plotReducedDim(spe,
+ dimred = "UMAP_mnnCorrected",
+ colour_by = x,
+ by_exprs_values = "exprs",
+ point_size = 0.2)
+ return(p)
+ })
+
+plot_grid(plotlist = plot_list)
+```
+
+
+
+### Heatmap visualization
+
+Next, it is often useful to visualize single-cell expression per cell
+type in form of a heatmap. For this, we will use the `dittoHeatmap`
+function from the
+[DittoSeq](https://bioconductor.org/packages/release/bioc/html/dittoSeq.html)
+package.
+
+We sub-sample the dataset to 4000 cells for ease of visualization and
+overlay the cancer type and patient ID from which the cells were
+extracted.
+
+
+```r
+set.seed(220818)
+cur_cells <- sample(seq_len(ncol(spe)), 4000)
+
+# Heatmap visualization - DittoHeatmap
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$marker_class == "type"],
+ assay = "exprs",
+ cluster_cols = FALSE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("celltype", "indication", "patient_id"),
+ annotation_colors = list(indication = metadata(spe)$color_vectors$indication,
+ patient_id = metadata(spe)$color_vectors$patient_id,
+ celltype = metadata(spe)$color_vectors$celltype))
+```
+
+
+
+Similarly, we can visualize the mean marker expression per cell type for all
+cells bu first calculating the mean marker expression per cell type using the
+`aggregateAcrossCells` function from the
+[scuttle](https://bioconductor.org/packages/release/bioc/html/scuttle.html)
+package and then use `dittoHeatmap`. We will annotate the heatmap with the
+number of cells per cell type and we will used different ways for feature
+scaling.
+
+
+```r
+library(scuttle)
+
+## aggregate by cell type
+celltype_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
+ ids = spe$celltype,
+ statistics = "mean",
+ use.assay.type = "exprs",
+ subset.row = rownames(spe)[rowData(spe)$marker_class == "type"])
+
+# No scaling
+dittoHeatmap(celltype_mean,
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("celltype", "ncells"),
+ annotation_colors = list(celltype = metadata(spe)$color_vectors$celltype,
+ ncells = plasma(100)))
+```
+
+
+
+```r
+# Scaled to max
+dittoHeatmap(celltype_mean,
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scaled.to.max = TRUE,
+ heatmap.colors.max.scaled = inferno(100),
+ annot.by = c("celltype", "ncells"),
+ annotation_colors = list(celltype = metadata(spe)$color_vectors$celltype,
+ ncells = plasma(100)))
+```
+
+
+
+```r
+# Z score scaled
+dittoHeatmap(celltype_mean,
+ assay = "exprs",
+ cluster_cols = TRUE,
+ annot.by = c("celltype", "ncells"),
+ annotation_colors = list(celltype = metadata(spe)$color_vectors$celltype,
+ ncells = plasma(100)))
+```
+
+
+
+As illustrated above for not-, max-, and Z score-scaled expression values,
+different ways of scaling can have strong effects on visualization
+output and we encourage the user to test multiple options.
+
+Overall, we can observe cell-type specific marker expression (e.g., Tumor
+= Ecad high and B cells = CD20 high) in agreement with the gating scheme
+of Section \@ref(classification).
+
+### Violin plot visualization
+
+The `plotExpression` function from the
+[scater](https://bioconductor.org/packages/release/bioc/html/scater.html) package
+allows to plot the distribution of expression values across cell types
+for a chosen set of proteins. The output is a `ggplot` object which can be
+modified further.
+
+
+```r
+# Violin Plot - plotExpression
+plotExpression(spe[,cur_cells],
+ features = rownames(spe)[rowData(spe)$marker_class == "type"],
+ x = "celltype",
+ exprs_values = "exprs",
+ colour_by = "celltype") +
+ theme(axis.text.x = element_text(angle = 90))+
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+### Scatter plot visualization
+
+Moreover, a protein expression based scatter plot can be generated with
+`dittoScatterPlot` (returns a `ggplot` object). We overlay the plot with
+the cell type information.
+
+
+```r
+# Scatter plot
+dittoScatterPlot(spe,
+ x.var = "CD3",
+ y.var="CD20",
+ assay.x = "exprs",
+ assay.y = "exprs",
+ color.var = "celltype") +
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype) +
+ ggtitle("Scatterplot for CD3/CD20 labelled by celltype")
+```
+
+
+
+We can nicely observe how the "B next to T cell" phenotype (`BnTcell`)
+has high expression values for both CD20 and CD3.
+
+**Of note**, in a setting where the user aims to assign labels to
+clusters based on marker genes/proteins, all of the above plots can be
+particularly helpful.
+
+### Barplot visualization
+
+In order to display frequencies of cell types per sample/patient, the
+`dittoBarPlot` function will be used. Data can be represented as
+percentages or counts and again `ggplot` objects are outputted.
+
+
+```r
+# by sample_id - percentage
+dittoBarPlot(spe,
+ var = "celltype",
+ group.by = "sample_id") +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+```r
+# by patient_id - percentage
+dittoBarPlot(spe,
+ var = "celltype",
+ group.by = "patient_id") +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+```r
+# by patient_id - count
+dittoBarPlot(spe,
+ scale = "count",
+ var = "celltype",
+ group.by = "patient_id") +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+We can see that cell type frequencies change between samples/patients
+and that the highest proportion/counts of plasma cells and stromal
+cells can be observed for Patient 2 and Patient 4, respectively.
+
+### CATALYST-based visualization
+
+In the following, we highlight some useful visualization
+functions from the
+[CATALYST](https://bioconductor.org/packages/release/bioc/html/CATALYST.html)
+package.
+
+To this end, we will first convert the `SpatialExperiment` object into a
+CATALYST-compatible format.
+
+
+```r
+library(CATALYST)
+
+# Save SPE in CATALYST-compatible object with renamed colData entries and
+# new metadata information
+spe_cat <- spe
+
+spe_cat$sample_id <- factor(spe$sample_id)
+spe_cat$condition <- factor(spe$indication)
+spe_cat$cluster_id <- factor(spe$celltype)
+
+# Add celltype information to metadata
+metadata(spe_cat)$cluster_codes <- data.frame(celltype = factor(spe_cat$celltype))
+```
+
+All of the `CATALYST` functions presented below return `ggplot` objects,
+which allow flexible downstream adjustment.
+
+#### Pseudobulk-level MDS plot
+
+Pseudobulk-level multi-dimensional scaling (MDS) plots can be rendered
+with the exported `pbMDS` function.
+
+Here, we will use `pbMDS` to highlight expression similarities between
+cell types and subsequently for each celltype-sample-combination.
+
+
+```r
+# MDS pseudobulk by cell type
+pbMDS(spe_cat,
+ by = "cluster_id",
+ features = rownames(spe_cat)[rowData(spe_cat)$marker_class == "type"],
+ label_by = "cluster_id",
+ k = "celltype") +
+ scale_color_manual(values = metadata(spe_cat)$color_vectors$celltype)
+```
+
+
+
+```r
+# MDS pseudobulk by cell type and sample_id
+pbMDS(spe_cat,
+ by = "both",
+ features = rownames(spe_cat)[rowData(spe_cat)$marker_class == "type"],
+ k = "celltype",
+ shape_by = "condition",
+ size_by = TRUE) +
+ scale_color_manual(values = metadata(spe_cat)$color_vectors$celltype)
+```
+
+
+
+We can see that the pseudobulk-expression profile of neutrophils seems
+markedly distinct from the other cell types, while comparable cell types
+such as the T cell subtypes group together. Furthermore, pseudobulk
+cell-type profiles from SCCHN appear different from the other
+indications.
+
+#### Reduced dimension plot on CLR of proportions
+
+The `clrDR` function produces dimensionality reduction plots on centered
+log-ratios (CLR) of sample/cell type proportions across cell
+type/samples.
+
+As with `pbMDS`, the output plots aim to illustrate the degree of
+similarity between cell types based on sample proportions.
+
+
+```r
+# CLR on cluster proportions across samples
+clrDR(spe_cat,
+ dr = "PCA",
+ by = "cluster_id",
+ k = "celltype",
+ label_by = "cluster_id",
+ arrow_col = "sample_id",
+ point_pal = metadata(spe_cat)$color_vectors$celltype)
+```
+
+
+
+We can again observe that neutrophils have a divergent profile also in
+terms of their sample proportions.
+
+#### Pseudobulk expression boxplot
+
+The `plotPbExprs` generates combined box- and jitter-plots of aggregated marker
+expression per cell type and sample (image). Here, we further split the data by
+cancer type.
+
+
+```r
+plotPbExprs(spe_cat,
+ k = "celltype",
+ facet_by = "cluster_id",
+ ncol = 2,
+ features = rownames(spe_cat)[rowData(spe_cat)$marker_class == "type"]) +
+ scale_color_manual(values = metadata(spe_cat)$color_vectors$indication)
+```
+
+
+
+Notably, CD15 levels are elevated in SCCHN in comparison to all other
+indications for most cell types.
+
+## Sample-level {#sample-level}
+
+In the next section, we will shift the grouping-level focus from the
+cell type to the sample-level. Sample-levels will be further divided
+into the sample-(image) and patient-level.
+
+Although we will mostly repeat the functions from the previous section
+\@ref(cell-type-level), sample- and patient-level centered visualization
+can provide additional quality control and biological interpretation.
+
+### Dimensionality reduction visualization
+
+Visualization of low-dimensional embeddings, here comparing non-corrected and
+fastMNN-corrected UMAPs, and coloring it by sample-levels is often used
+for "batch effect" assessment as mentioned in Section
+\@ref(cell-quality).
+
+We will again use `dittoDimPlot`.
+
+
+```r
+## UMAP colored by cell type and expression - dittoDimPlot
+p1 <- dittoDimPlot(spe,
+ var = "sample_id",
+ reduction.use = "UMAP",
+ size = 0.2,
+ colors = viridis(100),
+ do.label = FALSE) +
+ scale_color_manual(values = metadata(spe)$color_vectors$sample_id) +
+ theme(legend.title = element_blank()) +
+ ggtitle("Sample ID")
+
+p2 <- dittoDimPlot(spe,
+ var = "sample_id",
+ reduction.use = "UMAP_mnnCorrected",
+ size = 0.2,
+ colors = viridis(100),
+ do.label = FALSE) +
+ scale_color_manual(values = metadata(spe)$color_vectors$sample_id) +
+ theme(legend.title = element_blank()) +
+ ggtitle("Sample ID")
+
+p3 <- dittoDimPlot(spe,
+ var = "patient_id",
+ reduction.use = "UMAP",
+ size = 0.2,
+ do.label = FALSE) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ theme(legend.title = element_blank()) +
+ ggtitle("Patient ID")
+
+p4 <- dittoDimPlot(spe,
+ var = "patient_id",
+ reduction.use = "UMAP_mnnCorrected",
+ size = 0.2,
+ do.label = FALSE) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ theme(legend.title = element_blank()) +
+ ggtitle("Patient ID")
+
+(p1 + p2) / (p3 + p4)
+```
+
+
+
+As illustrated in Section \@ref(batch-effects), we see that the fastMNN
+approach (right side of the plot) leads to mixing of cells across
+samples/patients and thus batch effect correction.
+
+### Heatmap visualization
+
+It can be beneficial to use a heatmap to visualize single-cell
+expression per sample and patient. Such a plot, which we will create
+using `dittoHeatmap`, can highlight biological differences across
+samples/patients.
+
+
+```r
+# Heatmap visualization - DittoHeatmap
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$marker_class == "type"],
+ assay = "exprs",
+ order.by = c("patient_id","sample_id"),
+ cluster_cols = FALSE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("celltype", "indication", "patient_id", "sample_id"),
+ annotation_colors = list(celltype = metadata(spe)$color_vectors$celltype,
+ indication = metadata(spe)$color_vectors$indication,
+ patient_id = metadata(spe)$color_vectors$patient_id,
+ sample_id = metadata(spe)$color_vectors$sample_id))
+```
+
+
+
+As in Section \@ref(image-quality), aggregated mean marker expression
+per sample/patient allow identification of samples/patients with
+outlying expression patterns.
+
+Here, we will focus on the patient level and use `aggregateAcrossCells`
+and `dittoHeatmap`. The heatmap will be annotated with the number of
+cells per patient and cancer type and displayed using two scaling
+options.
+
+
+```r
+# mean expression by patient_id
+patient_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
+ ids = spe$patient_id,
+ statistics = "mean",
+ use.assay.type = "exprs",
+ subset.row = rownames(spe)[rowData(spe)$marker_class == "type"])
+
+# No scaling
+dittoHeatmap(patient_mean,
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("patient_id","indication","ncells"),
+ annotation_colors = list(patient_id = metadata(spe)$color_vectors$patient_id,
+ indication = metadata(spe)$color_vectors$indication,
+ ncells = plasma(100)))
+```
+
+
+
+```r
+# Max expression scaling
+dittoHeatmap(patient_mean,
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scaled.to.max = TRUE,
+ heatmap.colors.max.scaled = inferno(100),
+ annot.by = c("patient_id","indication","ncells"),
+ annotation_colors = list(patient_id = metadata(spe)$color_vectors$patient_id,
+ indication = metadata(spe)$color_vectors$indication,
+ ncells = plasma(100)))
+```
+
+
+
+As seen before, CD15 levels are elevated in Patient 4 (SCCHN), while SMA
+levels are highest for Patient 4 (CRC).
+
+### Barplot visualization
+
+Complementary to displaying cell type frequencies per sample/patient, we
+can use `dittoBarPlot` to display sample/patient frequencies per cell
+type.
+
+
+```r
+dittoBarPlot(spe,
+ var = "patient_id",
+ group.by = "celltype") +
+ scale_fill_manual(values = metadata(spe)$color_vectors$patient_id)
+```
+
+
+
+```r
+dittoBarPlot(spe,
+ var = "sample_id",
+ group.by = "celltype") +
+ scale_fill_manual(values = metadata(spe)$color_vectors$sample_id)
+```
+
+
+
+`Patient2` has the highest and lowest proportion of plasma cells and
+neutrophils, respectively.
+
+### CATALYST-based visualization
+
+#### Pseudobulk-level MDS plot
+
+Expression-based pseudobulks for each sample can be compared with the
+`pbMDS` function.
+
+
+```r
+# MDS pseudobulk by sample_id
+pbMDS(spe_cat,
+ by = "sample_id",
+ color_by = "sample_id",
+ features = rownames(spe_cat)[rowData(spe_cat)$marker_class == "type"]) +
+ scale_color_manual(values = metadata(spe_cat)$color_vectors$sample_id)
+```
+
+
+
+There are marked differences in pseudobulk-expression patterns between
+samples and across patients, which can be driven by biological
+differences and also technical aspects such as divergent region
+selection.
+
+#### Reduced dimension plot on CLR of proportions
+
+The `clrDR` function can also be used to analyze similarity of samples
+based on cell type proportions.
+
+
+```r
+# CLR on sample proportions across clusters
+clrDR(spe_cat,
+ dr = "PCA",
+ by = "sample_id",
+ point_col = "sample_id",
+ k = "celltype",
+ point_pal = metadata(spe_cat)$color_vectors$sample_id) +
+ scale_color_manual(values = metadata(spe_cat)$color_vectors$celltype)
+```
+
+```
+## Scale for colour is already present.
+## Adding another scale for colour, which will replace the existing scale.
+```
+
+
+
+There are notable differences between samples based on their cell type
+proportions.
+
+Interestingly, `Patient3_001`, `Patient1_003`, `Patient4_007` and
+`Patient4_006` group together and the PC loadings indicate a strong
+contribution of BnT and B cells, which could propose formation of
+tertiary lymphoid structures (TLS). In section \@ref(spatial-viz), we
+will be able to confirm this hypothesis visually on the images.
+
+## Further examples {#rich-example}
+
+In the last section of this chapter, we will use the popular
+[ComplexHeatmap](https://bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html)
+package to create a visualization example that combines various
+cell-type- and sample-level information.
+
+[ComplexHeatmap](https://bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html)
+is highly versatile and is originally inspired from the
+[pheatmap](https://cran.r-project.org/web/packages/pheatmap/index.html)
+package. Therefore, many arguments have the same/similar names.
+
+For more details, we would recommend to read the [reference
+book](https://jokergoo.github.io/ComplexHeatmap-reference/book/).
+
+### Publication-ready ComplexHeatmap
+
+For this example, we will concatenate heatmaps and annotations
+horizontally into one rich heatmap list. The grouping-level for the
+visualization will again be the cell type information from Section
+\@ref(classification)
+
+Initially, we will create two separate `Heatmap` objects for cell type
+and state markers.
+
+Then, metadata information, including the cancer type proportion and
+number of cells/patients per cell type, will be extracted into
+`HeatmapAnnotation` objects.
+
+Notably, we will add spatial features per cell type, here the number of
+neighbors extracted from `colPair(spe)` and cell area, in another
+`HeatmapAnnotation` object.
+
+Ultimately, all objects are combined in a `HeatmapList` and visualized.
+
+
+```r
+library(ComplexHeatmap)
+library(circlize)
+library(tidyverse)
+set.seed(22)
+
+### 1. Heatmap bodies ###
+
+# Heatmap body color
+col_exprs <- colorRamp2(c(0,1,2,3,4),
+ c("#440154FF","#3B518BFF","#20938CFF",
+ "#6ACD5AFF","#FDE725FF"))
+
+# Create Heatmap objects
+# By cell type markers
+celltype_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
+ ids = spe$celltype,
+ statistics = "mean",
+ use.assay.type = "exprs",
+ subset.row = rownames(spe)[rowData(spe)$marker_class == "type"])
+
+h_type <- Heatmap(t(assay(celltype_mean, "exprs")),
+ column_title = "type_markers",
+ col = col_exprs,
+ name= "mean exprs",
+ show_row_names = TRUE,
+ show_column_names = TRUE)
+
+# By cell state markers
+cellstate_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
+ ids = spe$celltype,
+ statistics = "mean",
+ use.assay.type = "exprs",
+ subset.row = rownames(spe)[rowData(spe)$marker_class == "state"])
+
+h_state <- Heatmap(t(assay(cellstate_mean, "exprs")),
+ column_title = "state_markers",
+ col = col_exprs,
+ name= "mean exprs",
+ show_row_names = TRUE,
+ show_column_names = TRUE)
+
+
+### 2. Heatmap annotation ###
+
+### 2.1 Metadata features
+
+anno <- colData(celltype_mean) %>% as.data.frame %>% select(celltype, ncells)
+
+# Proportion of indication per celltype
+indication <- unclass(prop.table(table(spe$celltype, spe$indication), margin = 1))
+
+# Number of contributing patients per celltype
+cluster_PID <- colData(spe) %>%
+ as.data.frame() %>%
+ select(celltype, patient_id) %>%
+ group_by(celltype) %>% table() %>%
+ as.data.frame()
+
+n_PID <- cluster_PID %>%
+ filter(Freq>0) %>%
+ group_by(celltype) %>%
+ count(name = "n_PID") %>%
+ column_to_rownames("celltype")
+
+# Create HeatmapAnnotation objects
+ha_anno <- HeatmapAnnotation(celltype = anno$celltype,
+ border = TRUE,
+ gap = unit(1,"mm"),
+ col = list(celltype = metadata(spe)$color_vectors$celltype),
+ which = "row")
+
+ha_meta <- HeatmapAnnotation(n_cells = anno_barplot(anno$ncells, width = unit(10, "mm")),
+ n_PID = anno_barplot(n_PID, width = unit(10, "mm")),
+ indication = anno_barplot(indication,width = unit(10, "mm"),
+ gp = gpar(fill = metadata(spe)$color_vectors$indication)),
+ border = TRUE,
+ annotation_name_rot = 90,
+ gap = unit(1,"mm"),
+ which = "row")
+
+### 2.2 Spatial features
+
+# Add number of neighbors to spe object (saved in colPair)
+spe$n_neighbors <- countLnodeHits(colPair(spe, "neighborhood"))
+
+# Select spatial features and average over celltypes
+spatial <- colData(spe) %>%
+ as.data.frame() %>%
+ select(area, celltype, n_neighbors)
+
+spatial <- spatial %>%
+ select(-celltype) %>%
+ aggregate(by = list(celltype = spatial$celltype), FUN = mean) %>%
+ column_to_rownames("celltype")
+
+# Create HeatmapAnnotation object
+ha_spatial <- HeatmapAnnotation(
+ area = spatial$area,
+ n_neighbors = spatial$n_neighbors,
+ border = TRUE,
+ gap = unit(1,"mm"),
+ which = "row")
+
+### 3. Plot rich heatmap ###
+
+# Create HeatmapList object
+h_list <- h_type +
+ h_state +
+ ha_anno +
+ ha_spatial +
+ ha_meta
+
+# Add customized legend for anno_barplot()
+lgd <- Legend(title = "indication",
+ at = colnames(indication),
+ legend_gp = gpar(fill = metadata(spe)$color_vectors$indication))
+
+# Plot
+draw(h_list,annotation_legend_list = list(lgd))
+```
+
+
+
+This plot summarizes most of the information we have seen in this
+chapter previously. In addition, we can observe that tumor cells have
+the largest mean cell area, high number of neighbors and elevated Ki67
+expression. BnT cells have the highest number of neighbors on average,
+which is biological sound given their predominant location in highly
+immune infiltrated regions (such as TLS).
+
+### Interactive visualization
+
+For interactive visualization of the single-cell data the
+[iSEE](https://www.bioconductor.org/packages/release/bioc/html/iSEE.html) shiny
+application can be used. For a comprehensive tutorial, please refer to the
+[iSEE vignette](https://www.bioconductor.org/packages/release/bioc/vignettes/iSEE/inst/doc/basic.html).
+
+
+```r
+if (interactive()) {
+ library(iSEE)
+
+ iSEE(spe)
+}
+```
+
+## Session Info
+
+
+
+The plot above shows the tissue expression of the epithelial tumor marker
+E-cadherin on the 3 selected images. The `bcg` parameter (default `c(0, 1, 1)`)
+stands for "background", "contrast", "gamma" and controls these attributes of
+the image. This parameter takes a named list where each entry specifies these
+attributes per channel. The first value of the numeric vector will be added to
+the pixel intensities (background); pixel intensities will be multiplied by the
+second entry of the vector (contrast); pixel intensities will be exponentiated
+by the third entry of the vector (gamma). In most cases, it is sufficient to
+adjust the second (contrast) entry of the vector.
+
+The following example highlights the visualization of 6 markers (maximum allowed
+number of markers) at once per image. The markers indicate the spatial
+distribution of tumor cells (E-cadherin), T cells (CD3), B cells (CD20), CD8+ T
+cells (CD8a), plasma cells (CD38) and proliferating cells (Ki67).
+
+
+```r
+plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1),
+ CD8a = c(0, 5, 1),
+ CD38 = c(0, 8, 1),
+ Ki67 = c(0, 5, 1)))
+```
+
+
+
+### Adjusting colors
+
+The default colors for visualization are chosen by the additive RGB (red, green,
+blue) color model. For six markers the default colors are: red, green, blue,
+cyan (green + blue), magenta (red + blue), yellow (green + red). These colors
+are the easiest to distinguish by eye. However, you can select other colors for
+each channel by setting the `colour` parameter:
+
+
+```r
+plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1)),
+ colour = list(Ecad = c("black", "burlywood1"),
+ CD3 = c("black", "cyan2"),
+ CD20 = c("black", "firebrick1")))
+```
+
+
+
+The `colour` parameter takes a named list in which each entry specifies the
+colors from which a color gradient is constructed via `colorRampPalette`. These
+are usually vectors of length 2 in which the first entry is `"black"` and the
+second entry specifies the color of choice. Although not recommended, you can
+also specify more than two colors to generate a more complex color gradient.
+
+### Image normalization
+
+As an alternative to setting the `bcg` parameter, images can first be
+normalized. Normalization here means to scale the pixel intensities per channel
+between 0 and 1 (or a range specified by the `ft` parameter in the `normalize`
+function). By default, the `normalize` function scales pixel intensities across
+**all** images contained in the `CytoImageList` object (`separateImages = FALSE`).
+Each individual channel is scaled independently (`separateChannels = TRUE`).
+
+After 0-1 normalization, maximum pixel intensities can be clipped to enhance the
+contrast of the image (setting the `inputRange` parameter). In the following
+example, the clipping to 0 and 0.2 is the same as multiplying the pixel
+intensities by a factor of 5.
+
+
+```r
+# 0 - 1 channel scaling across all images
+norm_images <- cytomapper::normalize(cur_images)
+
+# Clip channel at 0.2
+norm_images <- cytomapper::normalize(norm_images, inputRange = c(0, 0.2))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
+```
+
+
+
+The default setting of scaling pixel intensities across all images ensures
+comparable intensity levels across images. Pixel intensities can also be
+scaled **per image** therefore correcting for staining/expression differences
+between images:
+
+
+```r
+# 0 - 1 channel scaling per image
+norm_images <- cytomapper::normalize(cur_images, separateImages = TRUE)
+
+# Clip channel at 0.2
+norm_images <- cytomapper::normalize(norm_images, inputRange = c(0, 0.2))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
+```
+
+
+
+As we can see, the marker Ki67 appears brighter on image 2 and 3 in comparison
+to scaling the channel across all images.
+
+Finally, the `normalize` function also accepts a named list input for the
+`inputRange` argument. In this list, the clipping range per channel can be set
+individually:
+
+
+```r
+# 0 - 1 channel scaling per image
+norm_images <- cytomapper::normalize(cur_images,
+ separateImages = TRUE,
+ inputRange = list(Ecad = c(0, 50),
+ CD3 = c(0, 30),
+ CD20 = c(0, 40),
+ CD8a = c(0, 50),
+ CD38 = c(0, 10),
+ Ki67 = c(0, 70)))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
+```
+
+
+
+## Cell visualization {#mask-visualization}
+
+In the following section, we will show examples on how to visualize single
+cells either as segmentation masks or outlined on composite images. This type
+of visualization allows to observe the spatial distribution of cell phenotypes,
+the visual assessment of morphological features and quality control in terms
+of cell segmentation and phenotyping.
+
+### Visualzing metadata
+
+The `cytomapper` package provides the `plotCells` function that accepts a
+`CytoImageList` object containing segmentation masks. These are defined as
+single channel images where sets of pixels with the same integer ID identify
+individual cells. This integer ID can be found as an entry in the `colData(spe)`
+slot and as pixel information in the segmentation masks. The entry in
+`colData(spe)` needs to be specified via the `cell_id` argument to the
+`plotCells` function. In that way, data contained in the `SpatialExperiment`
+object can be mapped to segmentation masks. For the current dataset, the cell
+IDs are stored in `colData(spe)$ObjectNumber`.
+
+As cell IDs are only unique within a single image, `plotCells` also requires
+the `img_id` argument. This argument specifies the `colData(spe)` as well as the
+`mcols(masks)` entry that stores the unique image name from which each cell was
+extracted. In the current dataset the unique image names are stored in
+`colData(spe)$sample_id` and `mcols(masks)$sample_id`.
+
+Providing these two entries that allow mapping between the `SpatialExperiment`
+object and segmentation masks, we can now color individual cells based on their
+cell type:
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype")
+```
+
+
+
+For consistent visualization, the `plotCells` function takes a named list as
+`color` argument. The entry name must match the `colour_by` argument.
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype))
+```
+
+
+
+If only individual cell types should be visualized, the `SpatialExperiment`
+object can be subsetted (e.g., to only contain CD8+ T cells). In the following
+example CD8+ T cells are colored in red and all other cells that are not
+contained in the dataset are colored in white (as set by the `missing_color`
+argument).
+
+
+```r
+CD8 <- spe[,spe$celltype == "CD8"]
+
+plotCells(cur_masks,
+ object = CD8,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = c(CD8 = "red")),
+ missing_colour = "white")
+```
+
+
+
+In terms of visualizing metadata, any entry in the `colData(spe)` slot can be
+visualized. The `plotCells` function automatically detects if the entry
+is continuous or discrete. In this fashion, we can now visualize the area of each
+cell:
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "area")
+```
+
+
+
+### Visualizating expression
+
+Similar to visualizing single-cell metadata on segmentation masks, we can
+use the `plotCells` function to visualize the aggregated pixel intensities
+per cell. In the current dataset pixel intensities were aggregated by computing
+the mean pixel intensity per cell and per channel. The `plotCells` function
+accepts the `exprs_values` argument (default `counts`) that allows selecting
+the assay which stores the expression values that should be visualized.
+
+In the following example, we visualize the asinh-transformed mean pixel
+intensities of the epithelial marker E-cadherin on segmentation masks.
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "Ecad",
+ exprs_values = "exprs")
+```
+
+
+
+We will now visualize the maximum number of
+allowed markers as composites on the segmentation masks. As above the markers
+indicate the spatial distribution of tumor cells (E-cadherin), T cells (CD3), B
+cells (CD20), CD8+ T cells (CD8a), plasma cells (CD38) and proliferating cells
+(Ki67).
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ exprs_values = "exprs")
+```
+
+
+
+While visualizing 6 markers on the pixel-level may still allow the distinction
+of different tissue structures, observing single-cell expression levels is
+difficult when visualizing many markers simultaneously due to often overlapping
+expression.
+
+Similarly to adjusting marker colors when visualizing pixel intensities, we
+can change the color gradients per marker by setting the `color` argument:
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ exprs_values = "exprs",
+ colour = list(Ecad = c("black", "burlywood1"),
+ CD3 = c("black", "cyan2"),
+ CD20 = c("black", "firebrick1")))
+```
+
+
+
+### Outlining cells on images {#outline-cells}
+
+The following section highlights the combined visualization of pixel- and
+cell-level information at once. For this, besides the `SpatialExperiment` object,
+the `plotPixels` function accepts two `CytoImageList` objects. One for the
+multi-channel images and one for the segmentation masks. By specifying the
+`outline_by` parameter, the outlines of cells can now be colored based on their
+metadata.
+
+The following example first generates a 3-channel composite images displaying
+the expression of E-cadherin, CD3 and CD20 before coloring the cells' outlines
+by their cell phenotype.
+
+
+```r
+plotPixels(image = cur_images,
+ mask = cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ outline_by = "celltype",
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1)),
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ thick = TRUE)
+```
+
+
+
+Distinguishing individual cell phenotypes is nearly impossible in the images
+above.
+
+However, the `SpatialExperiment` object can be subsetted to only contain cells
+of a single or few phenotypes. This allows the selective visualization of cell
+outlines on composite images.
+
+Here, we select all CD8+ T cells from the dataset and outline them on a 2-channel
+composite image displaying the expression of CD3 and CD8a.
+
+
+```r
+CD8 <- spe[,spe$celltype == "CD8"]
+
+plotPixels(image = cur_images,
+ mask = cur_masks,
+ object = CD8,
+ cell_id = "ObjectNumber", img_id = "sample_id",
+ colour_by = c("CD3", "CD8a"),
+ outline_by = "celltype",
+ bcg = list(CD3 = c(0, 5, 1),
+ CD8a = c(0, 5, 1)),
+ colour = list(celltype = c("CD8" = "white")),
+ thick = TRUE)
+```
+
+
+
+This type of visualization allows the quality control of two things: 1.
+segmentation quality of individual cell types can be checked and 2. cell
+phenotyping accuracy can be visually assessed against expected marker expression.
+
+## Adjusting plot annotations
+
+The `cytomapper` package provides a number of function arguments to adjust the
+visual appearance of figures that are shared between the `plotPixels` and
+`plotCells` function.
+
+For a full overview of the arguments please refer to `?plotting-param`.
+
+We use the following example to highlight how to adjust the scale bar, the image
+title, the legend appearance and the margin between images.
+
+
+```r
+plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1),
+ CD8a = c(0, 5, 1),
+ CD38 = c(0, 8, 1),
+ Ki67 = c(0, 5, 1)),
+ scale_bar = list(length = 100,
+ label = expression("100 " ~ mu * "m"),
+ cex = 0.7,
+ lwidth = 10,
+ colour = "grey",
+ position = "bottomleft",
+ margin = c(5,5),
+ frame = 3),
+ image_title = list(text = mcols(cur_images)$indication,
+ position = "topright",
+ colour = "grey",
+ margin = c(5,5),
+ font = 2,
+ cex = 2),
+ legend = list(colour_by.title.cex = 0.7,
+ margin = 10),
+ margin = 40)
+```
+
+
+
+## Displaying individual images
+
+By default, all images are displayed on the same graphics device. This can be
+useful when saving all images at once (see next section) to zoom into the
+individual images instead of opening each image individually. However, when
+displaying images in a markdown document these are more accessible when
+visualized individually. For this, the `plotPixels` and `plotCells` function
+accepts the `display` parameter that when set to `"single"` displays each
+resulting image in its own graphics device:
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ display = "single",
+ legend = NULL)
+```
+
+
+
+## Saving and returning images
+
+The final section addresses how to save composite images and how to return them
+for integration with other plots.
+
+The `plotPixels` and `plotCells` functions accept the `save_plot` argument which
+takes a named list of the following entries: `filename` indicates the location
+and file type of the image saved to disk; `scale` adjusts the resolution of the
+saved image (this only needs to be adjusted for small images).
+
+
+```r
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ save_plot = list(filename = "data/celltype_image.png"))
+```
+
+The composite images (together with their annotation) can also be returned. In
+the following code chunk we save two example plots to variables (`out1` and
+`out2`).
+
+
+```r
+out1 <- plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ return_plot = TRUE)
+```
+
+```r
+out2 <- plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ exprs_values = "exprs",
+ return_plot = TRUE)
+```
+
+The composite images are stored in `out1$plot` and `out2$plot` and can be
+converted into a graph object recognized by the
+[cowplot](https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html)
+package.
+
+The final function call of the following chunk plots both object next to each
+other.
+
+
+```r
+library(cowplot)
+library(gridGraphics)
+p1 <- ggdraw(out1$plot, clip = "on")
+p2 <- ggdraw(out2$plot, clip = "on")
+
+plot_grid(p1, p2)
+```
+
+
+
+## Interactive image visualization
+
+The
+[cytoviewer](https://bioconductor.org/packages/release/bioc/html/cytoviewer.html)
+package allows the interactive visualization of multi-channel images and
+segmentation masks. It also allows to map cellular metadata onto segmentation
+masks and outlining of cells on composite images. For a full introduction to the
+package, please refer to [the vignette](https://bioconductor.org/packages/release/bioc/vignettes/cytoviewer/inst/doc/cytoviewer.html).
+
+
+```r
+library(cytoviewer)
+
+app <- cytoviewer(image = images,
+ mask = masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id")
+
+if (interactive()) {
+ shiny::runApp(app, launch.browser = TRUE)
+}
+```
+
+## Session Info
+
+
+
+```r
+# knn interaction graph
+plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "celltype",
+ img_id = "sample_id",
+ draw_edges = TRUE,
+ colPairName = "knn_interaction_graph",
+ nodes_first = FALSE,
+ edge_color_fix = "grey") +
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype) +
+ ggtitle("knn interaction graph")
+```
+
+
+
+```r
+# expansion interaction graph
+plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "celltype",
+ img_id = "sample_id",
+ draw_edges = TRUE,
+ colPairName = "expansion_interaction_graph",
+ nodes_first = FALSE,
+ edge_color_fix = "grey") +
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype) +
+ ggtitle("expansion interaction graph")
+```
+
+
+
+```r
+# delaunay interaction graph
+plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "celltype",
+ img_id = "sample_id",
+ draw_edges = TRUE,
+ colPairName = "delaunay_interaction_graph",
+ nodes_first = FALSE,
+ edge_color_fix = "grey") +
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype) +
+ ggtitle("delaunay interaction graph")
+```
+
+
+
+Nodes can also be colored based on the cells' expression levels (e.g.,
+E-cadherin expression) and their size can be adjusted (e.g., based on measured
+cell area).
+
+
+```r
+plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "Ecad",
+ assay_type = "exprs",
+ img_id = "sample_id",
+ draw_edges = TRUE,
+ colPairName = "expansion_interaction_graph",
+ nodes_first = FALSE,
+ node_size_by = "area",
+ directed = FALSE,
+ edge_color_fix = "grey") +
+ scale_size_continuous(range = c(0.1, 2)) +
+ ggtitle("E-cadherin expression")
+```
+
+
+
+Finally, the `plotSpatial` function allows displaying all images at once. This
+visualization can be useful to quickly detect larger structures of interest.
+
+
+```r
+plotSpatial(spe,
+ node_color_by = "celltype",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+For a full documentation on the `plotSpatial` function, please refer to
+`?plotSpatial`.
+
+## Spatial community analysis
+
+The detection of spatial communities was proposed by [@Jackson2020]. Here, cells
+are clustered solely based on their interactions as defined by the spatial
+object graph. We can perform spatial community detection across all cells as
+displayed in the next code chunk. Communities with less than 10 cells are
+excluded. **Of note:** we set the seed outside of the function call for
+reproducibility porposes as internally the `louvain` modularity optimization
+function is used which gives different results over different runs.
+
+
+```r
+set.seed(230621)
+spe <- detectCommunity(spe,
+ colPairName = "neighborhood",
+ size_threshold = 10)
+
+plotSpatial(spe,
+ node_color_by = "spatial_community",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ theme(legend.position = "none") +
+ ggtitle("Spatial tumor communities") +
+ scale_color_manual(values = rev(colors()))
+```
+
+
+
+The example shown above might not be of interest if different tissue structures
+exist within which spatial communities should be computed. In the following
+example, we perform spatial community detection separately for tumor and stromal
+cells.
+
+The general procedure is as follows:
+
+1. create a `colData(spe)` entry that specifies if a cell is part of the tumor
+or stroma compartment.
+
+2. use the `detectCommunity` function of the `imcRtools`
+package to cluster cells within the tumor or stroma compartment solely based on
+their spatial interaction graph as constructed by the `steinbock` package.
+
+Both tumor and stromal spatial communities are stored in the `colData` of
+the `SpatialExperiment` object under the `spatial_community` identifier.
+
+**Of note:** Here, and in contrast to the function call above, we set the seed
+argument within the `SerialParam` function for reproducibility
+purposes. We need this here due to the way the `detectCommunity` function
+is implemented when setting the `group_by` parameter.
+
+
+```r
+spe$tumor_stroma <- ifelse(spe$celltype == "Tumor", "Tumor", "Stroma")
+
+library(BiocParallel)
+spe <- detectCommunity(spe,
+ colPairName = "neighborhood",
+ size_threshold = 10,
+ group_by = "tumor_stroma",
+ BPPARAM = SerialParam(RNGseed = 220819))
+```
+
+We can now separately visualize the tumor and stromal communities.
+
+
+```r
+plotSpatial(spe[,spe$celltype == "Tumor"],
+ node_color_by = "spatial_community",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ theme(legend.position = "none") +
+ ggtitle("Spatial tumor communities") +
+ scale_color_manual(values = rev(colors()))
+```
+
+
+
+```r
+plotSpatial(spe[,spe$celltype != "Tumor"],
+ node_color_by = "spatial_community",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ theme(legend.position = "none") +
+ ggtitle("Spatial non-tumor communities") +
+ scale_color_manual(values = rev(colors()))
+```
+
+
+
+The example data was acquired using a panel that mainly focuses on immune cells.
+We are therefore unable to detect many tumor sub-phenotypes and will
+focus on the stromal communities.
+
+In the next step, the fraction of cell types within each
+spatial stromal community is displayed.
+
+
+```r
+library(pheatmap)
+library(viridis)
+
+cur_spe <- spe[,spe$celltype != "Tumor"]
+
+for_plot <- prop.table(table(cur_spe$spatial_community,
+ cur_spe$celltype),
+ margin = 1)
+
+pheatmap(for_plot,
+ color = colorRampPalette(c("dark blue", "white", "dark red"))(100),
+ show_rownames = FALSE,
+ scale = "column")
+```
+
+
+
+We observe that many spatial stromal communities are made up of myeloid cells or
+"stromal" (non-immune) cells. Other communities are mainly made up of B cells
+and BnT cells indicating tertiary lymphoid structures (TLS). While plasma cells,
+CD4$^+$ or CD8$^+$ T cells tend to aggregate, only in few spatial stromal
+communities consists of mainly neutrophils.
+
+## Cellular neighborhood analysis
+
+The following section highlights the use of the `imcRtools` package to
+detect cellular neighborhoods. This approach has been proposed by
+[@Goltsev2018] and [@Schurch2020] to group cells based on information
+contained in their direct neighborhood.
+
+[@Goltsev2018] perfomed Delaunay triangulation-based graph construction,
+neighborhood aggregation and then clustered cells. [@Schurch2020] on the
+other hand constructed a 10-nearest neighbor graph before aggregating
+information across neighboring cells.
+
+In the following code chunk we will use the 20-nearest neighbor graph as
+constructed above to define the direct cellular neighborhood. The
+[aggregateNeighbors](https://bodenmillergroup.github.io/imcRtools/reference/aggregateNeighbors.html)
+function allows neighborhood aggregation in 2 different ways:
+
+1. For each cell the function computes the fraction of cells of a
+ certain type (e.g., cell type) among its neighbors.
+2. For each cell it aggregates (e.g., mean) the expression counts
+ across all neighboring cells.
+
+Based on these measures, cells can now be clustered into cellular
+neighborhoods. We will first compute the fraction of the different cell
+types among the 20-nearest neighbors and use kmeans clustering to group
+cells into 6 cellular neighborhoods.
+
+**Of note:** constructing a 20-nearest neighbor graph and clustering
+using kmeans with `k=6` is only an example. Similar to the analysis done
+in Section \@ref(snn-graph), it is recommended to perform a parameter
+sweep across different graph construction algorithms and different
+parmaters `k` for kmeans clustering. Finding the best CN detection
+settings is also subject to the question at hand. Constructing graphs
+with more neighbors usually results in larger CNs.
+
+
+```r
+# By celltypes
+spe <- aggregateNeighbors(spe,
+ colPairName = "knn_interaction_graph",
+ aggregate_by = "metadata",
+ count_by = "celltype")
+
+set.seed(220705)
+
+cn_1 <- kmeans(spe$aggregatedNeighbors, centers = 6)
+spe$cn_celltypes <- as.factor(cn_1$cluster)
+
+plotSpatial(spe,
+ node_color_by = "cn_celltypes",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_brewer(palette = "Set3")
+```
+
+
+
+The next code chunk visualizes the cell type compositions of the
+detected cellular neighborhoods (CN).
+
+
+```r
+for_plot <- prop.table(table(spe$cn_celltypes, spe$celltype),
+ margin = 1)
+
+pheatmap(for_plot,
+ color = colorRampPalette(c("dark blue", "white", "dark red"))(100),
+ scale = "column")
+```
+
+
+
+CN 1 and CN 6 are mainly composed of tumor cells with CN 6 forming the
+tumor/stroma border. CN 3 is mainly composed of B and BnT cells
+indicating TLS. CN 5 is composed of aggregated plasma cells and most T
+cells.
+
+We will now detect cellular neighborhoods by computing the mean
+expression across the 20-nearest neighbor prior to kmeans clustering
+(k=6).
+
+
+```r
+# By expression
+spe <- aggregateNeighbors(spe,
+ colPairName = "knn_interaction_graph",
+ aggregate_by = "expression",
+ assay_type = "exprs",
+ subset_row = rowData(spe)$use_channel)
+
+set.seed(220705)
+
+cn_2 <- kmeans(spe$mean_aggregatedExpression, centers = 6)
+spe$cn_expression <- as.factor(cn_2$cluster)
+
+plotSpatial(spe,
+ node_color_by = "cn_expression",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_brewer(palette = "Set3")
+```
+
+
+
+Also here, we can visualize the cell type composition of each cellular
+neighborhood.
+
+
+```r
+for_plot <- prop.table(table(spe$cn_expression, spe$celltype),
+ margin = 1)
+
+pheatmap(for_plot,
+ color = colorRampPalette(c("dark blue", "white", "dark red"))(100),
+ scale = "column")
+```
+
+
+
+When clustering cells based on the mean expression within the direct
+neighborhood, tumor cells are split across CN 6, CN 1 and CN 4 without
+forming a clear tumor/stroma interface. This result reflects
+patient-to-patient differences in the expression of tumor markers.
+
+CN 3 again contains B cells and BnT cells but also CD8 and undefined
+cells, therefore it is less representative of TLS compared to CN 3 in
+previous CN approach. CN detection based on mean marker expression is
+therefore sensitive to staining/expression differences between samples
+as well as lateral spillover due to imperfect segmentation.
+
+An alternative to the `aggregateNeighbors` function is provided by the
+[lisaClust](https://bioconductor.org/packages/release/bioc/html/lisaClust.html)
+Bioconductor package [@Patrick2023]. In contrast to `imcRtools`, the
+`lisaClust` package computes local indicators of spatial associations
+(LISA) functions and clusters cells based on those. More precise, the
+package summarizes L-functions from a Poisson point process model to
+derive numeric vectors for each cell which can then again be clustered
+using kmeans. All steps are supported by the `lisaClust` function which
+can be applied to a `SingleCellExperiment` and `SpatialExperiment` object.
+
+
+In the following example, we calculate the LISA curves within a 10µm, 20µm and
+50µm neighborhood around each cell. Increasing these radii will lead to broader
+and smoother spatial clusters. However, a number of parameter settings should be
+tested to estimate the robustness of the results.
+
+
+```r
+library(lisaClust)
+
+set.seed(220705)
+spe <- lisaClust(spe,
+ k = 6,
+ Rs = c(10, 20, 50),
+ spatialCoords = c("Pos_X", "Pos_Y"),
+ cellType = "celltype",
+ imageID = "sample_id")
+
+plotSpatial(spe,
+ node_color_by = "region",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_brewer(palette = "Set3")
+```
+
+
+
+Similar to the example above, we can now observe the cell type
+composition per spatial cluster.
+
+
+```r
+for_plot <- prop.table(table(spe$region, spe$celltype),
+ margin = 1)
+
+pheatmap(for_plot,
+ color = colorRampPalette(c("dark blue", "white", "dark red"))(100),
+ scale = "column")
+```
+
+
+
+In this case, CN 1 and 4 contain tumor cells but no CN is forming the
+tumor/stroma interface. CN 3 represents TLS. CN 2 indicates T cell
+subtypes and plasma cells are aggregated to CN 5.
+
+As an alternative way of visualizing the enrichment of cell types within the
+detected CNs, the `lisaClust` package provides the `regionMap` function.
+
+
+```r
+regionMap(spe,
+ cellType = "celltype",
+ region = "region")
+```
+
+
+
+## Spatial context analysis
+
+Downstream of CN assignments, we will analyze the spatial context (SC)
+of each cell using three functions from the `imcRtools` package.
+
+While CNs can represent sites of unique local processes, the term SC was
+coined by Bhate and colleagues [@Bhate2022] and describes tissue regions
+in which distinct CNs may be interacting. Hence, SCs may be interesting
+regions of specialized biological events.
+
+Here, we will first detect SCs using the `detectSpatialContext` function. This
+function relies on CN fractions for each cell in a spatial interaction
+graph (originally a KNN graph), which we will calculate using
+`buildSpatialGraph` and `aggregateNeighbors`. We will focus on the CNs
+derived from cell type fractions but other CN assignments are possible.
+
+**Of note**, the window size (k for KNN) for `buildSpatialGraph` should
+reflect a length scale on which biological signals can be exchanged and
+depends, among others, on cell density and tissue area. In view of their
+divergent functionality, we recommend to use a larger window size for SC
+(interaction between local processes) than for CN (local processes)
+detection. Since we used a 20-nearest neighbor graph for CN assignment,
+we will use a 40-nearest neighbor graph for SC detection. As before,
+different parameters should be tested.
+
+Subsequently, the CN fractions are sorted from high-to-low and the SC of
+each cell is assigned as the minimal combination of SCs that additively
+surpass a user-defined threshold. The default threshold of 0.9 aims to
+represent the dominant CNs, hence the most prevalent signals, in a given
+window.
+
+For more details and biological validation, please refer to
+[@Bhate2022].
+
+
+```r
+library(circlize)
+library(RColorBrewer)
+
+# Construct a 40-nearest neighbor graph
+spe <- buildSpatialGraph(spe,
+ img_id = "sample_id",
+ type = "knn",
+ name = "knn_spatialcontext_graph",
+ k = 40)
+
+# Compute the fraction of cellular neighborhoods around each cell
+spe <- aggregateNeighbors(spe,
+ colPairName = "knn_spatialcontext_graph",
+ aggregate_by = "metadata",
+ count_by = "cn_celltypes",
+ name = "aggregatedNeighborhood")
+
+# Detect spatial contexts
+spe <- detectSpatialContext(spe,
+ entry = "aggregatedNeighborhood",
+ threshold = 0.90,
+ name = "spatial_context")
+
+# Define SC color scheme
+n_SCs <- length(unique(spe$spatial_context))
+col_SC <- setNames(colorRampPalette(brewer.pal(9, "Paired"))(n_SCs),
+ sort(unique(spe$spatial_context)))
+
+# Visualize spatial contexts on images
+plotSpatial(spe,
+ node_color_by = "spatial_context",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_manual(values = col_SC)
+```
+
+
+
+We detect a total of 52 distinct
+SCs across this dataset.
+
+For ease of interpretation, we will directly compare the CN and SC
+assignments for `Patient3_001`.
+
+
+```r
+library(patchwork)
+
+# Compare CN and SC for one patient
+p1 <- plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "cn_celltypes",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_brewer(palette = "Set3")
+
+p2 <- plotSpatial(spe[,spe$sample_id == "Patient3_001"],
+ node_color_by = "spatial_context",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_manual(values = col_SC, limits = force)
+
+p1 + p2
+```
+
+
+
+As expected, we can observe that interfaces between different CNs make
+up distinct SCs. For instance, interface between CN 3 (TLS region
+consisting of B and BnT cells) and CN 5 (Plasma- and T-cell dominated)
+turns to SC 3_5. On the other hand, the core of CN 3 becomes SC 3, since
+the most abundant CN of the neighborhood for these cells is just the CN
+itself.
+
+Next, we filter the SCs based on user-defined thresholds for number of
+group entries (here at least 3 patients) and/or total number of cells
+(here minimum of 100 cells) per SC using the `filterSpatialContext` function.
+
+
+```r
+## Filter spatial contexts
+# By number of group entries
+spe <- filterSpatialContext(spe,
+ entry = "spatial_context",
+ group_by = "patient_id",
+ group_threshold = 3,
+ name = "spatial_context_filtered")
+
+plotSpatial(spe,
+ node_color_by = "spatial_context_filtered",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_manual(values = col_SC, limits = force)
+```
+
+
+
+```r
+# Filter out small and infrequent spatial contexts
+spe <- filterSpatialContext(spe,
+ entry = "spatial_context",
+ group_by = "patient_id",
+ group_threshold = 3,
+ cells_threshold = 100,
+ name = "spatial_context_filtered")
+
+plotSpatial(spe,
+ node_color_by = "spatial_context_filtered",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_manual(values = col_SC, limits = force)
+```
+
+
+
+Lastly, we can use the `plotSpatialContext` function to generate *SC
+graphs*, analogous to *CN combination maps* in [@Bhate2022]. Returned
+objects are `ggplots`, which can be easily modified further. We will
+create a SC graph for the filtered SCs here.
+
+
+```r
+## Plot spatial context graph
+
+# Colored by name, size by n_cells
+plotSpatialContext(spe,
+ entry = "spatial_context_filtered",
+ group_by = "sample_id",
+ node_color_by = "name",
+ node_size_by = "n_cells",
+ node_label_color_by = "name")
+```
+
+
+
+```r
+# Colored by n_cells, size by n_group
+plotSpatialContext(spe,
+ entry = "spatial_context_filtered",
+ group_by = "sample_id",
+ node_color_by = "n_cells",
+ node_size_by = "n_group",
+ node_label_color_by = "n_cells") +
+ scale_color_viridis()
+```
+
+
+
+SC 1 (Tumor-dominated), SC 1_6 (Tumor and Tumor-Stroma interface) and SC
+4_5 (Plasma/T cell and Myeloid/Neutrophil interface) are the most
+frequent SCs in this dataset. Moreover, we may compare the degree of the
+different nodes in the SC graph. For example, we can observe that SC 1
+has only one degree (directed to SC 1_6), while SC 5 (T cells and plasma cells) has
+a much higher degree (n = 4) and potentially more CN interactions.
+
+## Patch detection
+
+The previous section focused on detecting cellular neighborhoods in a rather
+unsupervised fashion. However, the `imcRtools` package also provides methods for
+detecting spatial compartments in a supervised fashion. The
+[patchDetection](https://bodenmillergroup.github.io/imcRtools/reference/patchDetection.html)
+function allows the detection of connected sets of similar cells as proposed by
+[@Hoch2022]. In the following example, we will use the `patchDetection` function
+to detect tumor patches in three steps:
+
+1. Find connected sets of tumor cells (using the `steinbock` graph).
+2. Components which contain less than 10 cells are excluded.
+3. Expand the components by 1µm to construct a concave hull around the patch and
+include cells within the patch.
+
+
+```r
+spe <- patchDetection(spe,
+ patch_cells = spe$celltype == "Tumor",
+ img_id = "sample_id",
+ expand_by = 1,
+ min_patch_size = 10,
+ colPairName = "neighborhood",
+ BPPARAM = MulticoreParam())
+```
+
+```
+## The returned object is ordered by the 'sample_id' entry.
+```
+
+```r
+plotSpatial(spe,
+ node_color_by = "patch_id",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ theme(legend.position = "none") +
+ scale_color_manual(values = rev(colors()))
+```
+
+
+
+We can now calculate the fraction of T cells within each tumor patch to roughly
+estimate T cell infiltration.
+
+
+```r
+library(tidyverse)
+colData(spe) %>% as_tibble() %>%
+ group_by(patch_id, sample_id) %>%
+ summarize(Tcell_count = sum(celltype == "CD8" | celltype == "CD4"),
+ patch_size = n(),
+ Tcell_freq = Tcell_count / patch_size) %>%
+ filter(!is.na(patch_id)) %>%
+ ggplot() +
+ geom_point(aes(log10(patch_size), Tcell_freq, color = sample_id)) +
+ theme_classic()
+```
+
+
+
+We can now measure the size of each patch using the
+[patchSize](https://bodenmillergroup.github.io/imcRtools/reference/patchSize.html)
+function and visualize tumor patch distribution per patient.
+
+
+```r
+patch_size <- patchSize(spe, "patch_id")
+
+patch_size <- merge(patch_size,
+ colData(spe)[match(patch_size$patch_id, spe$patch_id),],
+ by = "patch_id")
+
+ggplot(as.data.frame(patch_size)) +
+ geom_boxplot(aes(patient_id, log10(size))) +
+ geom_point(aes(patient_id, log10(size)))
+```
+
+
+
+The
+[minDistToCells](https://bodenmillergroup.github.io/imcRtools/reference/minDistToCells.html)
+function can be used to calculate the minimum distance between each cell and a
+cell set of interest. Here, we highlight its use to calculate the minimum
+distance of all cells to the detected tumor patches. Negative values indicate
+the minimum distance of each tumor patch cell to a non-tumor patch cell.
+
+
+```r
+spe <- minDistToCells(spe,
+ x_cells = !is.na(spe$patch_id),
+ img_id = "sample_id")
+```
+
+```
+## The returned object is ordered by the 'sample_id' entry.
+```
+
+```r
+plotSpatial(spe,
+ node_color_by = "distToCells",
+ img_id = "sample_id",
+ node_size_fix = 0.5) +
+ scale_color_gradient2(low = "dark blue", mid = "white", high = "dark red")
+```
+
+
+
+Finally, we can observe the minimum distances to tumor patches in a cell type
+specific manner.
+
+
+```r
+library(ggridges)
+
+ggplot(as.data.frame(colData(spe))) +
+ geom_density_ridges(aes(distToCells, celltype, fill = celltype)) +
+ geom_vline(xintercept = 0, color = "dark red", linewidth = 2) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype)
+```
+
+
+
+## Interaction analysis
+
+**Bug notice: we discovered and fixed a bug in the `testInteractions` function in version below 1.5.5 which affected `SingleCellExperiment` or `SpatialExperiment` objects in which cells were not grouped by image. Please make sure you have the newest version (>= 1.6.0) installed.**
+
+The next section focuses on statistically testing the pairwise interaction
+between all cell types of the dataset. For this, the `imcRtools` package
+provides the
+[testInteractions](https://bodenmillergroup.github.io/imcRtools/reference/testInteractions.html)
+function which implements the interaction testing strategy proposed by
+[@Shapiro2017].
+
+Per grouping level (e.g., image), the `testInteractions` function computes the
+averaged cell type/cell type interaction count and compares this count against
+an empirical null distribution which is generated by permuting all cell labels (while maintaining the tissue structure).
+
+In the following example, we use the `steinbock` generated spatial interaction
+graph and estimate the interaction or avoidance between cell types in the
+dataset.
+
+
+```r
+library(scales)
+out <- testInteractions(spe,
+ group_by = "sample_id",
+ label = "celltype",
+ colPairName = "neighborhood",
+ BPPARAM = SerialParam(RNGseed = 221029))
+
+head(out)
+```
+
+```
+## DataFrame with 6 rows and 10 columns
+## group_by from_label to_label ct p_gt p_lt
+## 
+
+In the plot above the red tiles indicate cell type pairs that were detected to
+significantly interact on a large number of images. On the other hand, blue
+tiles show cell type pairs which tend to avoid each other on a large number
+of images.
+
+Here we can observe that tumor cells are mostly compartmentalized and are in
+avoidance with other cell types. As expected, B cells interact with BnT cells;
+regulatory T cells interact with CD4+ T cells and CD8+ T cells. Most cell types
+show self interactions indicating spatial clustering.
+
+The `imcRtools` package further implements an interaction testing strategy
+proposed by [@Schulz2018] where the hypothesis is tested if at least n cells of
+a certain type are located around a target cell type (`from_cell`). This type of
+testing can be performed by selecting `method = "patch"` and specifying the
+number of patch cells via the `patch_size` parameter.
+
+
+```r
+out <- testInteractions(spe,
+ group_by = "sample_id",
+ label = "celltype",
+ colPairName = "neighborhood",
+ method = "patch",
+ patch_size = 3,
+ BPPARAM = SerialParam(RNGseed = 221029))
+
+out %>% as_tibble() %>%
+ group_by(from_label, to_label) %>%
+ summarize(sum_sigval = sum(sigval, na.rm = TRUE)) %>%
+ ggplot() +
+ geom_tile(aes(from_label, to_label, fill = sum_sigval)) +
+ scale_fill_gradient2(low = muted("blue"), mid = "white", high = muted("red")) +
+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
+```
+
+
+
+These results are comparable to the interaction testing presented above. The
+main difference comes from the lack of symmetry. We can now for example see that
+3 or more myeloid cells sit around CD4$^+$ T cells while this interaction is not
+as strong when considering CD4$^+$ T cells sitting around myeloid cells.
+
+Finally, we save the updated `SpatialExperiment` object.
+
+
+```r
+saveRDS(spe, "data/spe.rds")
+```
+
+
+
+## Session Info
+
+The page you requested cannot be found (perhaps it was moved or renamed).
+You may want to try searching to find the page's new location, or use +the table of contents to find the page you are looking for.
+In Section 7.4 we observed staining/expression differences +between the individual samples. This can arise due to technical (e.g., +differences in sample processing) as well as biological (e.g., differential +expression between patients/indications) effects. However, the combination of these effects +hinders cell phenotyping via clustering as highlighted in Section 9.2.
+To integrate cells across samples, we can use computational +strategies developed for correcting batch effects in single-cell RNA sequencing +data. In the following sections, we will use functions of the +batchelor, +harmony and +Seurat +packages to correct for such batch effects.
+Of note: the correction approaches presented here aim at removing any +differences between samples. This will also remove biological differences +between the patients/indications. Nevertheless, integrating cells across samples +can facilitate the detection of cell phenotypes via clustering.
+First, we will read in the SpatialExperiment object containing the single-cell
+data.
The batchelor package provides the mnnCorrect and fastMNN functions to
+correct for differences between samples/batches. Both functions build up on
+finding mutual nearest neighbors (MNN) among the cells of different samples and
+correct expression differences between the batches (Haghverdi et al. 2018). The mnnCorrect function
+returns corrected expression counts while the fastMNN functions performs the
+correction in reduced dimension space. As such, fastMNN returns integrated
+cells in form of a low dimensional embedding.
Paper: Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors
+Documentation: batchelor
Here, we apply the fastMNN function to integrate cells between
+patients. By setting auto.merge = TRUE the function estimates the best
+batch merging order by maximizing the number of MNN pairs at each merging step.
+This is more time consuming than merging sequentially based on how batches appear in the
+dataset (default). We again select the markers defined in Section 5.2
+for sample correction.
The function returns a SingleCellExperiment object which contains corrected
+low-dimensional coordinates for each cell in the reducedDim(out, "corrected")
+slot. This low-dimensional embedding can be further used for clustering and
+non-linear dimensionality reduction. We check that the order of cells is the
+same between the input and output object and then transfer the corrected
+coordinates to the main SpatialExperiment object.
library(batchelor)
+set.seed(220228)
+out <- fastMNN(spe, batch = spe$patient_id,
+ auto.merge = TRUE,
+ subset.row = rowData(spe)$use_channel,
+ assay.type = "exprs")
+
+# Check that order of cells is the same
+stopifnot(all.equal(colnames(spe), colnames(out)))
+
+# Transfer the correction results to the main spe object
+reducedDim(spe, "fastMNN") <- reducedDim(out, "corrected")The computational time of the fastMNN function call is
+2.33 minutes.
Of note, the warnings that the fastMNN function produces can be avoided as follows:
BSPARAM = BiocSingular::ExactParam()Warning in (function (A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth = TRUE, :
+ You're computing too large a percentage of total singular values, use a standard svd instead.d = 30In check_numbers(k = k, nu = nu, nv = nv, limit = min(dim(x)) - :
+ more singular values/vectors requested than availableThe fastMNN function further returns outputs that can be used to assess the
+quality of the batch correction. The metadata(out)$merge.info entry collects
+diagnostics for each individual merging step. Here, the batch.size and
+lost.var entries are important. The batch.size entry reports the relative
+magnitude of the batch effect and the lost.var entry represents the percentage
+of lost variance per merging step. A large batch.size and low lost.var
+indicate sufficient batch correction.
## DataFrame with 3 rows and 3 columns
+## left right batch.size
+## <List> <List> <numeric>
+## 1 Patient4 Patient2 0.381635
+## 2 Patient4,Patient2 Patient1 0.581013
+## 3 Patient4,Patient2,Patient1 Patient3 0.767376## Patient1 Patient2 Patient3 Patient4
+## [1,] 0.000000000 0.031154864 0.00000000 0.046198914
+## [2,] 0.043363546 0.009772150 0.00000000 0.011931892
+## [3,] 0.005394755 0.003023119 0.07219394 0.005366304We observe that Patient4 and Patient2 are most similar with a low batch effect. +Merging cells of Patient3 into the combined batch of Patient1, +Patient2 and Patient4 resulted in the highest percentage of lost variance and +the detection of the largest batch effect. In the next paragraph we can +visualize the correction results.
+The simplest option to check if the sample effects were corrected is by using +non-linear dimensionality reduction techniques and observe mixing of cells across +samples. We will recompute the UMAP embedding using the corrected +low-dimensional coordinates for each cell.
+library(scater)
+
+set.seed(220228)
+spe <- runUMAP(spe, dimred= "fastMNN", name = "UMAP_mnnCorrected") Next, we visualize the corrected UMAP while overlaying patient IDs.
+library(cowplot)
+library(dittoSeq)
+library(viridis)
+
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP before correction")
+p2 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP after correction")
+
+plot_grid(p1, p2)
We observe an imperfect merging of Patient3 into all other samples. This +was already seen when displaying the merging information above. +We now also visualize the expression of selected markers across all cells +before and after batch correction.
+markers <- c("Ecad", "CD45RO", "CD20", "CD3", "FOXP3", "CD206", "MPO", "SMA", "Ki67")
+
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_mnnCorrected",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
We observe that immune cells across patients are merged after batch correction
+using fastMNN. However, the tumor cells of different patients still cluster
+separately.
The harmony algorithm performs batch correction by iteratively clustering and
+correcting the positions of cells in PCA space (Korsunsky et al. 2019). We will first
+perform PCA on the asinh-transformed counts and then call the RunHarmony
+function to perform data integration.
Paper: Fast, sensitive and accurate integration of single-cell data with Harmony
+Documentation: harmony
Similar to the fastMNN function, harmony returns the corrected
+low-dimensional coordinates for each cell. These can be transfered to the
+reducedDim slot of the original SpatialExperiment object.
library(harmony)
+library(BiocSingular)
+
+spe <- runPCA(spe,
+ subset_row = rowData(spe)$use_channel,
+ exprs_values = "exprs",
+ ncomponents = 30,
+ BSPARAM = ExactParam())
+
+set.seed(230616)
+out <- RunHarmony(spe, group.by.vars = "patient_id")
+
+# Check that order of cells is the same
+stopifnot(all.equal(colnames(spe), colnames(out)))
+
+reducedDim(spe, "harmony") <- reducedDim(out, "HARMONY")The computational time of the HarmonyMatrix function call is
+1.3 minutes.
We will now again visualize the cells in low dimensions after UMAP embedding.
+ +# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP before correction")
+p2 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP_harmony", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP after correction")
+
+plot_grid(p1, p2)
And we visualize selected marker expression as defined above.
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_harmony",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
We observe a more aggressive merging of cells from different patients compared
+to the results after fastMNN correction. Importantly, immune cell and epithelial
+markers are expressed in distinct regions of the UMAP.
The Seurat package provides a number of functionalities to analyze single-cell
+data. As such it also allows the integration of cells across different samples.
+Conceptually, Seurat performs batch correction similarly to fastMNN by
+finding mutual nearest neighbors (MNN) in low dimensional space before
+correcting the expression values of cells (Stuart et al. 2019).
Paper: Comprehensive Integration of Single-Cell Data
+Documentation: Seurat
To use Seurat, we will first create a Seurat object from the SpatialExperiment
+object and add relevant metadata. The object also needs to be split by patient
+prior to integration.
library(Seurat)
+library(SeuratObject)
+seurat_obj <- as.Seurat(spe, counts = "counts", data = "exprs")
+seurat_obj <- AddMetaData(seurat_obj, as.data.frame(colData(spe)))
+
+seurat.list <- SplitObject(seurat_obj, split.by = "patient_id")To avoid long run times, we will use an approach that relies on reciprocal PCA
+instead of canonical correlation analysis for dimensionality reduction and
+initial alignment. For an extended tutorial on how to use Seurat for data
+integration, please refer to their
+vignette.
We will first define the features used for integration and perform PCA on cells
+of each patient individually. The FindIntegrationAnchors function detects MNNs between
+cells of different patients and the IntegrateData function corrects the
+expression values of cells. We slightly increase the number of neighbors to be
+considered for MNN detection (the k.anchor parameter). This increases the integration
+strength.
features <- rownames(spe)[rowData(spe)$use_channel]
+seurat.list <- lapply(X = seurat.list, FUN = function(x) {
+ x <- ScaleData(x, features = features, verbose = FALSE)
+ x <- RunPCA(x, features = features, verbose = FALSE, approx = FALSE)
+ return(x)
+})
+
+anchors <- FindIntegrationAnchors(object.list = seurat.list,
+ anchor.features = features,
+ reduction = "rpca",
+ k.anchor = 20)
+
+combined <- IntegrateData(anchorset = anchors)We now select the integrated assay and perform PCA dimensionality reduction.
+The cell coordinates in PCA reduced space can then be transferred to the
+original SpatialExperiment object. Of note: by splitting the object into
+individual batch-specific objects, the ordering of cells in the integrated
+object might not match the ordering of cells in the input object. In this case,
+columns will need to be reordered. Here, we test if the ordering of cells in the
+integrated Seurat object matches the ordering of cells in the main
+SpatialExperiment object.
DefaultAssay(combined) <- "integrated"
+
+combined <- ScaleData(combined, verbose = FALSE)
+combined <- RunPCA(combined, npcs = 30, verbose = FALSE, approx = FALSE)
+
+# Check that order of cells is the same
+stopifnot(all.equal(colnames(spe), colnames(combined)))
+
+reducedDim(spe, "seurat") <- Embeddings(combined, reduction = "pca")The computational time of the Seurat function calls is
+4.29 minutes.
As above, we recompute the UMAP embeddings based on Seurat integrated results
+and visualize the embedding.
Visualize patient IDs.
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP before correction")
+p2 <- dittoDimPlot(spe, var = "patient_id",
+ reduction.use = "UMAP_seurat", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP after correction")
+
+plot_grid(p1, p2)
Visualization of marker expression.
+# Before correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
# After correction
+plot_list <- multi_dittoDimPlot(spe, var = markers, reduction.use = "UMAP_seurat",
+ assay = "exprs", size = 0.2, list.out = TRUE)
+plot_list <- lapply(plot_list, function(x) x + scale_color_viridis())
+plot_grid(plotlist = plot_list)
Similar to the methods presented above, Seurat integrates immune cells correctly.
+When visualizing the patient IDs, slight patient-to-patient differences within tumor
+cells can be detected.
Choosing the correct integration approach is challenging without having ground truth +cell labels available. It is recommended to compare different techniques and different +parameter settings. Please refer to the documentation of the individual tools +to become familiar with the possible parameter choices. Furthermore, in the following +section, we will discuss clustering and classification approaches in light of +expression differences between samples.
+In general, it appears that MNN-based approaches are less conservative in terms
+of merging compared to harmony. On the other hand, harmony could well merge
+cells in a way that regresses out biological signals.
The modified SpatialExperiment object is saved for further downstream analysis.
## R version 4.3.1 (2023-06-16)
+## Platform: x86_64-pc-linux-gnu (64-bit)
+## Running under: Ubuntu 22.04.3 LTS
+##
+## Matrix products: default
+## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
+## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
+##
+## locale:
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [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: Etc/UTC
+## tzcode source: system (glibc)
+##
+## attached base packages:
+## [1] stats4 stats graphics grDevices utils datasets methods
+## [8] base
+##
+## other attached packages:
+## [1] testthat_3.1.10 SeuratObject_4.1.4
+## [3] Seurat_4.4.0 BiocSingular_1.16.0
+## [5] harmony_1.0.1 Rcpp_1.0.11
+## [7] viridis_0.6.4 viridisLite_0.4.2
+## [9] dittoSeq_1.12.1 cowplot_1.1.1
+## [11] scater_1.28.0 ggplot2_3.4.3
+## [13] scuttle_1.10.2 SpatialExperiment_1.10.0
+## [15] batchelor_1.16.0 SingleCellExperiment_1.22.0
+## [17] SummarizedExperiment_1.30.2 Biobase_2.60.0
+## [19] GenomicRanges_1.52.0 GenomeInfoDb_1.36.3
+## [21] IRanges_2.34.1 S4Vectors_0.38.2
+## [23] BiocGenerics_0.46.0 MatrixGenerics_1.12.3
+## [25] matrixStats_1.0.0
+##
+## loaded via a namespace (and not attached):
+## [1] RcppAnnoy_0.0.21 splines_4.3.1
+## [3] later_1.3.1 bitops_1.0-7
+## [5] tibble_3.2.1 R.oo_1.25.0
+## [7] polyclip_1.10-6 lifecycle_1.0.3
+## [9] rprojroot_2.0.3 edgeR_3.42.4
+## [11] globals_0.16.2 lattice_0.21-8
+## [13] MASS_7.3-60 magrittr_2.0.3
+## [15] limma_3.56.2 plotly_4.10.2
+## [17] sass_0.4.7 rmarkdown_2.25
+## [19] jquerylib_0.1.4 yaml_2.3.7
+## [21] httpuv_1.6.11 sctransform_0.4.0
+## [23] spatstat.sparse_3.0-2 sp_2.0-0
+## [25] reticulate_1.32.0 pbapply_1.7-2
+## [27] RColorBrewer_1.1-3 ResidualMatrix_1.10.0
+## [29] pkgload_1.3.3 abind_1.4-5
+## [31] zlibbioc_1.46.0 Rtsne_0.16
+## [33] purrr_1.0.2 R.utils_2.12.2
+## [35] RCurl_1.98-1.12 GenomeInfoDbData_1.2.10
+## [37] ggrepel_0.9.3 irlba_2.3.5.1
+## [39] spatstat.utils_3.0-3 listenv_0.9.0
+## [41] pheatmap_1.0.12 goftest_1.2-3
+## [43] spatstat.random_3.1-6 dqrng_0.3.1
+## [45] fitdistrplus_1.1-11 parallelly_1.36.0
+## [47] DelayedMatrixStats_1.22.6 leiden_0.4.3
+## [49] codetools_0.2-19 DropletUtils_1.20.0
+## [51] DelayedArray_0.26.7 tidyselect_1.2.0
+## [53] farver_2.1.1 ScaledMatrix_1.8.1
+## [55] spatstat.explore_3.2-3 jsonlite_1.8.7
+## [57] BiocNeighbors_1.18.0 ellipsis_0.3.2
+## [59] progressr_0.14.0 ggridges_0.5.4
+## [61] survival_3.5-5 tools_4.3.1
+## [63] ica_1.0-3 glue_1.6.2
+## [65] gridExtra_2.3 xfun_0.40
+## [67] dplyr_1.1.3 HDF5Array_1.28.1
+## [69] withr_2.5.1 fastmap_1.1.1
+## [71] rhdf5filters_1.12.1 fansi_1.0.4
+## [73] digest_0.6.33 rsvd_1.0.5
+## [75] R6_2.5.1 mime_0.12
+## [77] colorspace_2.1-0 scattermore_1.2
+## [79] tensor_1.5 spatstat.data_3.0-1
+## [81] R.methodsS3_1.8.2 RhpcBLASctl_0.23-42
+## [83] utf8_1.2.3 tidyr_1.3.0
+## [85] generics_0.1.3 data.table_1.14.8
+## [87] httr_1.4.7 htmlwidgets_1.6.2
+## [89] S4Arrays_1.0.6 uwot_0.1.16
+## [91] pkgconfig_2.0.3 gtable_0.3.4
+## [93] lmtest_0.9-40 XVector_0.40.0
+## [95] brio_1.1.3 htmltools_0.5.6
+## [97] bookdown_0.35 scales_1.2.1
+## [99] png_0.1-8 knitr_1.44
+## [101] rstudioapi_0.15.0 reshape2_1.4.4
+## [103] rjson_0.2.21 nlme_3.1-162
+## [105] cachem_1.0.8 zoo_1.8-12
+## [107] rhdf5_2.44.0 stringr_1.5.0
+## [109] KernSmooth_2.23-21 parallel_4.3.1
+## [111] miniUI_0.1.1.1 vipor_0.4.5
+## [113] desc_1.4.2 pillar_1.9.0
+## [115] grid_4.3.1 vctrs_0.6.3
+## [117] RANN_2.6.1 promises_1.2.1
+## [119] beachmat_2.16.0 xtable_1.8-4
+## [121] cluster_2.1.4 waldo_0.5.1
+## [123] beeswarm_0.4.0 evaluate_0.21
+## [125] magick_2.8.0 cli_3.6.1
+## [127] locfit_1.5-9.8 compiler_4.3.1
+## [129] rlang_1.1.1 crayon_1.5.2
+## [131] future.apply_1.11.0 labeling_0.4.3
+## [133] plyr_1.8.8 ggbeeswarm_0.7.2
+## [135] stringi_1.7.12 deldir_1.0-9
+## [137] BiocParallel_1.34.2 munsell_0.5.0
+## [139] lazyeval_0.2.2 spatstat.geom_3.2-5
+## [141] Matrix_1.6-1.1 patchwork_1.1.3
+## [143] sparseMatrixStats_1.12.2 future_1.33.0
+## [145] Rhdf5lib_1.22.1 shiny_1.7.5
+## [147] ROCR_1.0-11 igraph_1.5.1
+## [149] bslib_0.5.1A common step during single-cell data analysis is the annotation of cells based +on their phenotype. Defining cell phenotypes is often subjective and relies +on previous biological knowledge. The Orchestrating Single Cell Analysis with Bioconductor book +presents a number of approaches to phenotype cells detected by single-cell RNA +sequencing based on reference datasets or gene set analysis.
+In highly-multiplexed imaging, target proteins or molecules are manually +selected based on the biological question at hand. It narrows down the feature +space and facilitates the manual annotation of clusters to derive cell +phenotypes. We will therefore discuss and compare a number of clustering +approaches to group cells based on their similarity in marker expression in +Section 9.2.
+Unlike single-cell RNA sequencing or CyTOF data, single-cell data derived from +highly-multiplexed imaging data often suffers from “lateral spillover” between +neighboring cells. This spillover caused by imperfect segmentation often hinders +accurate clustering to define specific cell phenotypes in multiplexed imaging +data. Tools have been developed to correct lateral spillover between cells +(Bai et al. 2021) but the approach requires careful selection of the markers to +correct. In Section 9.3 we will train and apply a random +forest classifier to classify cell phenotypes in the dataset as alternative +approach to clustering-based cell phenotyping. This approach has been previously used to +identify major cell phenotypes in metastatic melanoma and avoids clustering of +cells (Hoch et al. 2022).
+We will first read in the previously generated SpatialExperiment object and
+sample 2000 cells to visualize cluster membership.
In the first section, we will present clustering approaches to identify cellular +phenotypes in the dataset. These methods group cells based on their similarity +in marker expression or by their proximity in low dimensional space. A number of +approaches have been developed to cluster data derived from single-cell RNA +sequencing technologies (Yu et al. 2022) or CyTOF (Weber and Robinson 2016). For demonstration +purposes, we will highlight common clustering approaches that are available in R +and have been used for clustering cells obtained from IMC. Two approaches rely +on graph-based clustering and one approach uses self organizing maps (SOM).
+The PhenoGraph clustering approach was first described to group cells of a CyTOF
+dataset (Levine et al. 2015). The algorithm first constructs a graph by detecting the
+k nearest neighbours based on euclidean distance in expression space. In the
+next step, edges between nodes (cells) are weighted by their overlap in nearest
+neighbor sets. To quantify the overlap in shared nearest neighbor sets, the
+jaccard index is used. The Louvain modularity optimization approach is used to
+detect connected communities and partition the graph into clusters of cells.
+This clustering strategy was used by Jackson, Fischer et al. and Schulz et
+al. to cluster IMC data (Jackson et al. 2020; Schulz et al. 2018).
There are several different PhenoGraph implementations available in R. Here, we +use the one available at +https://github.com/i-cyto/Rphenograph. +For large datasets, +https://github.com/stuchly/Rphenoannoy +offers a more performant implementation of the algorithm.
+In the following code chunk, we select the asinh-transformed mean pixel
+intensities per cell and channel and subset the channels to the ones containing
+biological variation. This matrix is transposed to store cells in rows. Within
+the Rphenograph function, we select the 45 nearest neighbors for graph
+building and louvain community detection (default). The function returns a list
+of length 2, the first entry being the graph and the second entry containing the
+community object. Calling membership on the community object will return
+cluster IDs for each cell. These cluster IDs are then stored within the
+colData of the SpatialExperiment object. Cluster IDs are mapped on top of
+the UMAP embedding and single-cell marker expression within each cluster are
+visualized in form of a heatmap.
It is recommended to test different inputs to k as shown in the next section.
+Selecting larger values for k results in larger clusters.
library(Rphenograph)
+library(igraph)
+library(dittoSeq)
+library(viridis)
+
+mat <- t(assay(spe, "exprs")[rowData(spe)$use_channel,])
+
+set.seed(230619)
+out <- Rphenograph(mat, k = 45)
+
+clusters <- factor(membership(out[[2]]))
+
+spe$pg_clusters <- clusters
+
+dittoDimPlot(spe, var = "pg_clusters",
+ reduction.use = "UMAP", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("Phenograph clusters on UMAP")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("pg_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$pg_clusters))],
+ metadata(spe)$color_vectors$patient_id))
The Rphenograph function call took
+2.31 minutes.
We can observe that some of the clusters only contain cells of a single patient.
+This can often be observed in the tumor compartment. In the next step, we
+use the integrated cells (see Section 8) in low dimensional
+embedding for clustering. Here, the low dimensional embedding can
+be directly accessed from the reducedDim slot.
mat <- reducedDim(spe, "fastMNN")
+
+set.seed(230619)
+out <- Rphenograph(mat, k = 45)
+
+clusters <- factor(membership(out[[2]]))
+
+spe$pg_clusters_corrected <- clusters
+
+dittoDimPlot(spe, var = "pg_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("Phenograph clusters on UMAP, integrated cells")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("pg_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$pg_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
Clustering using the integrated embedding leads to clusters that contain cells +of different patients. Cluster annotation can now be performed by manually +labeling cells based on their marker expression (see Notes in Section +9.2.5).
+The bluster +package provides a simple interface to cluster cells using a number of different +clustering approaches and different metrics to access cluster stability.
+For simplicity, we will focus on graph based clustering as this is the most
+popular and a fast method for single-cell clustering. The bluster package
+provides functionalities to build k-nearest neighbor (KNN) graphs and its weighted
+version, shared nearest neighbor (SNN) graphs where nodes represent cells.
+The user can chose the number of neighbors to consider (parameter k),
+the edge weighting method (parameter type) and the community detection
+function to use (parameter cluster.fun). As all parameters affect the clustering
+results, the bluster package provides the clusterSweep function to test
+a number of parameter settings in parallel. In the following code chunk,
+we select the asinh-transformed mean pixel intensities and subset the markers
+of interest. The resulting matrix is transposed to fit to the requirements of
+the bluster package (cells in rows).
We test two different settings for k, two for type and fix the cluster.fun
+to louvain as this is one of the most common approaches for community detection.
+This function call is parallelized by setting the BPPARAM parameter.
library(bluster)
+library(BiocParallel)
+library(ggplot2)
+
+mat <- t(assay(spe, "exprs")[rowData(spe)$use_channel,])
+
+combinations <- clusterSweep(mat,
+ BLUSPARAM=SNNGraphParam(),
+ k=c(10L, 20L),
+ type = c("rank", "jaccard"),
+ cluster.fun = "louvain",
+ BPPARAM = MulticoreParam(RNGseed = 220427))We next calculate two metrics to estimate cluster stability: the average +silhouette width and the neighborhood purity.
+We use the approxSilhouette function to compute the silhouette width for each
+cell and compute the average across all cells per parameter setting. Please see
+?silhouette for more information on how the silhouette width is computed for
+each cell. A large average silhouette width indicates a cluster parameter
+setting for which cells that are well clustered.
The neighborPurity function computes the fraction of cells around each cell
+with the same cluster ID. Per parameter setting, we compute the average
+neighborhood purity across all cells. A large average neighborhood purity
+indicates a cluster parameter setting for which cells that are well clustered.
sil <- vapply(as.list(combinations$clusters),
+ function(x) mean(approxSilhouette(mat, x)$width),
+ 0)
+
+ggplot(data.frame(method = names(sil),
+ sil = sil)) +
+ geom_point(aes(method, sil)) +
+ theme_classic(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+ xlab("Cluster parameter combination") +
+ ylab("Average silhouette width")
pur <- vapply(as.list(combinations$clusters),
+ function(x) mean(neighborPurity(mat, x)$purity),
+ 0)
+
+ggplot(data.frame(method = names(pur),
+ pur = pur)) +
+ geom_point(aes(method, pur)) +
+ theme_classic(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+ xlab("Cluster parameter combination") +
+ ylab("Average neighborhood purity")
The cluster parameter sweep took +8.81 minutes.
+Performing a cluster sweep takes some time as multiple function calls are run in parallel. +We do however recommend testing a number of different parameter settings to +assess clustering performance.
+Once parameter settings are known, we can either use the clusterRows function
+of the bluster package to cluster cells or its convenient wrapper function
+exported by the
+scran package.
+The scran::clusterCells function accepts a SpatialExperiment (or
+SingleCellExperiment) object which stores cells in columns. By default, the
+function detects the 10 nearest neighbours for each cell, performs rank-based
+weighting of edges (see ?makeSNNGraph for more information) and uses the
+cluster_walktrap function to detect communities in the graph.
As we can see above, the clustering approach in this dataset with k being 20
+and rank-based edge weighting leads to the highest silhouette width and highest
+neighborhood purity.
library(scran)
+
+set.seed(220620)
+clusters <- clusterCells(spe[rowData(spe)$use_channel,],
+ assay.type = "exprs",
+ BLUSPARAM = SNNGraphParam(k=20,
+ cluster.fun = "louvain",
+ type = "rank"))
+
+spe$nn_clusters <- clusters
+
+dittoDimPlot(spe, var = "nn_clusters",
+ reduction.use = "UMAP", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SNN clusters on UMAP")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("nn_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$nn_clusters))],
+ metadata(spe)$color_vectors$patient_id))
The shared nearest neighbor graph clustering approach took +1.31 minutes.
+This function was used by (Tietscher et al. 2022) to cluster cells obtained by IMC. Setting
+type = "jaccard" performs clustering similar to Rphenograph above and Seurat.
Similar to the results obtained by Rphenograph, some of the clusters are
+patient-specific. We can now perform clustering of the integrated cells
+by directly specifying which low-dimensional embedding to use:
set.seed(220621)
+clusters <- clusterCells(spe,
+ use.dimred = "fastMNN",
+ BLUSPARAM = SNNGraphParam(k = 20,
+ cluster.fun = "louvain",
+ type = "rank"))
+
+spe$nn_clusters_corrected <- clusters
+
+dittoDimPlot(spe, var = "nn_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SNN clusters on UMAP, integrated cells")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("nn_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$nn_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
An alternative to graph-based clustering is offered by the
+CATALYST
+package. The cluster function internally uses the
+FlowSOM
+package to group cells into 100 (default) clusters based on self organizing maps
+(SOM). In the next step, the
+ConsensusClusterPlus
+package is used to perform hierarchical consensus clustering of the previously
+detected 100 SOM nodes into 2 to maxK clusters. Cluster stability for each k
+can be assessed by plotting the delta_area(spe). The optimal number
+of clusters can be found by selecting the k at which a plateau is reached.
+In the example below, an optimal k lies somewhere around 13.
library(CATALYST)
+
+# Run FlowSOM and ConsensusClusterPlus clustering
+set.seed(220410)
+spe <- cluster(spe,
+ features = rownames(spe)[rowData(spe)$use_channel],
+ maxK = 30)
+
+# Assess cluster stability
+delta_area(spe)
spe$som_clusters <- cluster_ids(spe, "meta13")
+
+dittoDimPlot(spe, var = "som_clusters",
+ reduction.use = "UMAP", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SOM clusters on UMAP")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("som_clusters", "patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$som_clusters))],
+ metadata(spe)$color_vectors$patient_id))
Running FlowSOM clustering took 0.22 minutes.
+The CATALYST package does not provide functionality to perform FlowSOM and
+ConsensusClusterPlus clustering directly on the batch-corrected, integrated cells. As an
+alternative to the CATALYST package, the bluster package provides SOM
+clustering when specifying the SomParam() parameter. Similar to the CATALYST
+approach, we will first cluster the dataset into 100 clusters (also called
+“codes”). These codes are then further clustered into a maximum of 30 clusters
+using ConsensusClusterPlus (using hierarchical clustering and euclidean
+distance). The delta area plot can be accessed using the (not exported)
+.plot_delta_area function from CATALYST. Here, it seems that the plateau is
+reached at a k of 16 and we will store the final cluster IDs within the
+SpatialExperiment object.
library(kohonen)
+library(ConsensusClusterPlus)
+
+# Select integrated cells
+mat <- reducedDim(spe, "fastMNN")
+
+# Perform SOM clustering
+set.seed(220410)
+som.out <- clusterRows(mat, SomParam(100), full = TRUE)
+
+# Cluster the 100 SOM codes into larger clusters
+ccp <- ConsensusClusterPlus(t(som.out$objects$som$codes[[1]]),
+ maxK = 30,
+ reps = 100,
+ distance = "euclidean",
+ seed = 220410,
+ plot = NULL)
# Link ConsensusClusterPlus clusters with SOM codes and save in object
+som.cluster <- ccp[[16]][["consensusClass"]][som.out$clusters]
+spe$som_clusters_corrected <- as.factor(som.cluster)
+
+dittoDimPlot(spe, var = "som_clusters_corrected",
+ reduction.use = "UMAP_mnnCorrected", size = 0.2,
+ do.label = TRUE) +
+ ggtitle("SOM clusters on UMAP, integrated cells")
dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("som_clusters_corrected","patient_id"),
+ annot.colors = c(dittoColors(1)[1:length(unique(spe$som_clusters_corrected))],
+ metadata(spe)$color_vectors$patient_id))
The FlowSOM clustering approach has been used by (Hoch et al. 2022) to sub-cluster tumor
+cells as measured by IMC.
Finally, we can compare the results of different clustering approaches. For +this, we visualize the number of cells that are shared between different +clustering results in a pairwise fashion. In the following heatmaps a high match +between clustering results can be seen for those clusters that are uniquely +detected in both approaches.
+First, we will visualize the match between the three different approaches +applied to the asinh-transformed counts.
+library(patchwork)
+library(pheatmap)
+library(gridExtra)
+
+tab1 <- table(paste("Rphenograph", spe$pg_clusters),
+ paste("SNN", spe$nn_clusters))
+tab2 <- table(paste("Rphenograph", spe$pg_clusters),
+ paste("SOM", spe$som_clusters))
+tab3 <- table(paste("SNN", spe$nn_clusters),
+ paste("SOM", spe$som_clusters))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
We observe that Rphenograph and the shared nearest neighbor (SNN) approach by
+scran show similar results (first heatmap above). For example, Rphenograph
+cluster 20 (a tumor cluster) is perfectly captured by SNN cluster 12. On the
+other hand, the Neutrophil cluster (SNN cluster 6) is split into Rphenograph
+cluster 2 and Rphenograph cluster 6. A common approach
+is to now merge clusters that contain similar cell types and annotate them by
+hand (see below).
Below, a comparison between the clustering results of the integrated cells +is shown.
+tab1 <- table(paste("Rphenograph", spe$pg_clusters_corrected),
+ paste("SNN", spe$nn_clusters_corrected))
+tab2 <- table(paste("Rphenograph", spe$pg_clusters_corrected),
+ paste("SOM", spe$som_clusters_corrected))
+tab3 <- table(paste("SNN", spe$nn_clusters_corrected),
+ paste("SOM", spe$som_clusters_corrected))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
In comparison to clustering on the non-integrated cells, the clustering results +of the integrated cells show higher overlap. The SNN approach resulted in fewer +clusters and therefore matches better with the SOM clustering approach.
+The bluster package provides a number of metrics to assess cluster stability
+here.
+For brevity we only highlighted the use of the silhouette width and the
+neighborhood purity but different metrics should be tested to assess cluster
+stability.
To assign cell types to clusters, we manually annotate clusters based on their +marker expression. For example, SNN cluster 12 (clustering of the integrated +cells) shows high, homogeneous expression of CD20 and we might therefore label +this cluster as B cells. The next chapter 10 will +highlight single-cell visualization methods that can be helpful for manual +cluster annotations.
+An example how to label clusters can be seen below:
+library(dplyr)
+cluster_celltype <- recode(spe$nn_clusters_corrected,
+ "1" = "Tumor_proliferating",
+ "2" = "Myeloid",
+ "3" = "Tumor",
+ "4" = "Tumor",
+ "5" = "Stroma",
+ "6" = "Proliferating",
+ "7" = "Myeloid",
+ "8" = "Plasma_cell",
+ "9" = "CD8",
+ "10" = "CD4",
+ "11" = "Neutrophil",
+ "12" = "Bcell",
+ "13" = "Stroma")
+
+spe$cluster_celltype <- cluster_celltypeIn this section, we will highlight a cell type classification approach based +on ground truth labeling and random forest classification. The rational for +this supervised cell phenotyping approach is to use the information contained +in the pre-defined markers to detect cells of interest. This approach was +used by Hoch et al. to classify cell types in a metastatic melanoma IMC +dataset (Hoch et al. 2022).
+The antibody panel used in the example data set mainly focuses on immune cell +types and little on tumor cell phenotypes. Therefore we will label the following +cell types:
+The “B next to T cell” phenotype (BnTcell) is commonly observed in immune
+infiltrated regions measured by IMC. We include this phenotype to account for B
+cell/T cell interactions where precise classification into B cells or T cells is
+not possible. The exact gating scheme can be seen at
+img/Gating_scheme.pdf.
As related approaches, Astir and +Garnett use pre-defined panel +information to classify cell phenotypes based on their marker expression.
+The cytomapper
+package provides the cytomapperShiny function that allows gating of cells
+based on their marker expression and visualization of selected cells directly
+on the images.
library(cytomapper)
+if (interactive()) {
+
+ images <- readRDS("data/images.rds")
+ masks <- readRDS("data/masks.rds")
+
+ cytomapperShiny(object = spe, mask = masks, image = images,
+ cell_id = "ObjectNumber", img_id = "sample_id")
+}The labeled cells for this data set can be accessed at
+10.5281/zenodo.6554544 and were downloaded
+in Section 4. Gating is performed per image and the
+cytomapperShiny function allows the export of gated cells in form of a
+SingleCellExperiment or SpatialExperiment object. The cell label is stored
+in colData(object)$cytomapper_CellLabel and the gates are stored in
+metadata(object). In the next section, we will read in and consolidate the
+labeled data.
For consistent visualization of cell types, we will now pre-define their colors:
+celltype <- setNames(c("#3F1B03", "#F4AD31", "#894F36", "#1C750C", "#EF8ECC",
+ "#6471E2", "#4DB23B", "grey", "#F4800C", "#BF0A3D", "#066970"),
+ c("Tumor", "Stroma", "Myeloid", "CD8", "Plasma_cell",
+ "Treg", "CD4", "undefined", "BnTcell", "Bcell", "Neutrophil"))
+
+metadata(spe)$color_vectors$celltype <- celltypeHere, we will read in the individual SpatialExperiment objects containing the
+labeled cells and concatenate them. In the process of concatenating the
+SpatialExperiment objects along their columns, the sample_id entry is
+appended by .1, .2, .3, ... due to replicated entries.
library(SingleCellExperiment)
+label_files <- list.files("data/gated_cells",
+ full.names = TRUE, pattern = ".rds$")
+
+# Read in SPE objects
+spes <- lapply(label_files, readRDS)
+
+# Merge SPE objects
+concat_spe <- do.call("cbind", spes)In the following code chunk we will identify cells that were labeled multiple +times. This occurs when different cell phenotypes are gated per image and can +affect immune cells that are located inside the tumor compartment.
+We will first identify those cells that were uniquely labeled. In the next step,
+we will identify those cells that were labeled twice AND were labeled as Tumor
+cells. These cells will be assigned their immune cell label. Finally, we will
+save the unique labels within the original SpatialExperiment object.
Of note: this concatenation strategy is specific for cell phenotypes contained in this +example dataset. The gated cell labels might need to be processed in a slightly +different way when working with other samples.
+For these tasks, we will define a filter function:
+filter_labels <- function(object,
+ label = "cytomapper_CellLabel") {
+ cur_tab <- unclass(table(colnames(object), object[[label]]))
+
+ cur_labels <- colnames(cur_tab)[apply(cur_tab, 1, which.max)]
+ names(cur_labels) <- rownames(cur_tab)
+
+ cur_labels <- cur_labels[rowSums(cur_tab) == 1]
+
+ return(cur_labels)
+}This function is now applied to all cells and then only non-tumor cells.
+labels <- filter_labels(concat_spe)
+
+cur_spe <- concat_spe[,concat_spe$cytomapper_CellLabel != "Tumor"]
+
+non_tumor_labels <- filter_labels(cur_spe)
+
+additional_cells <- setdiff(names(non_tumor_labels), names(labels))
+
+final_labels <- c(labels, non_tumor_labels[additional_cells])
+
+# Transfer labels to SPE object
+spe_labels <- rep("unlabeled", ncol(spe))
+names(spe_labels) <- colnames(spe)
+spe_labels[names(final_labels)] <- final_labels
+spe$cell_labels <- spe_labels
+
+# Number of cells labeled per patient
+table(spe$cell_labels, spe$patient_id)##
+## Patient1 Patient2 Patient3 Patient4
+## Bcell 152 131 234 263
+## BnTcell 396 37 240 1029
+## CD4 45 342 167 134
+## CD8 60 497 137 128
+## Myeloid 183 378 672 517
+## Neutrophil 97 4 17 16
+## Plasma_cell 34 536 87 59
+## Stroma 84 37 85 236
+## Treg 139 149 49 24
+## Tumor 2342 906 1618 1133
+## unlabeled 7214 9780 7826 9580Based on these labels, we can now train a random forest classifier to classify +all remaining, unlabeled cells.
+In this section, we will use the +caret framework for machine +learning in R. This package provides an interface to train a number of +regression and classification models in a coherent fashion. We use a random +forest classifier due to low number of parameters, high speed and an observed +high performance for cell type classification (Hoch et al. 2022).
+In the following section, we will first split the SpatialExperiment object
+into labeled and unlabeled cells. Based on the labeled cells, we split
+the data into a train (75% of the data) and test (25% of the data) dataset.
+We currently do not provide an independently labeled validation dataset.
The caret package provides the trainControl function, which specifies model
+training parameters and the train function, which performs the actual model
+training. While training the model, we also want to estimate the best model
+parameters. In the case of the chosen random forest model (method = "rf"), we
+only need to estimate a single parameters (mtry) which corresponds to the
+number of variables randomly sampled as candidates at each split. To estimate
+the best parameter, we will perform a 5-fold cross validation (set within
+trainControl) over a tune length of 5 entries to mtry. In the following
+code chunk, the createDataPartition and the train function are not deterministic,
+meaning they return different results across different runs. We therefore set
+a seed here for both functions.
library(caret)
+
+# Split between labeled and unlabeled cells
+lab_spe <- spe[,spe$cell_labels != "unlabeled"]
+unlab_spe <- spe[,spe$cell_labels == "unlabeled"]
+
+# Randomly split into train and test data
+set.seed(221029)
+trainIndex <- createDataPartition(factor(lab_spe$cell_labels), p = 0.75)
+
+train_spe <- lab_spe[,trainIndex$Resample1]
+test_spe <- lab_spe[,-trainIndex$Resample1]
+
+# Define fit parameters for 5-fold cross validation
+fitControl <- trainControl(method = "cv",
+ number = 5)
+
+# Select the arsinh-transformed counts for training
+cur_mat <- t(assay(train_spe, "exprs")[rowData(train_spe)$use_channel,])
+
+# Train a random forest classifier
+rffit <- train(x = cur_mat,
+ y = factor(train_spe$cell_labels),
+ method = "rf", ntree = 1000,
+ tuneLength = 5,
+ trControl = fitControl)
+
+rffit## Random Forest
+##
+## 10049 samples
+## 37 predictor
+## 10 classes: 'Bcell', 'BnTcell', 'CD4', 'CD8', 'Myeloid', 'Neutrophil', 'Plasma_cell', 'Stroma', 'Treg', 'Tumor'
+##
+## No pre-processing
+## Resampling: Cross-Validated (5 fold)
+## Summary of sample sizes: 8040, 8039, 8038, 8038, 8041
+## Resampling results across tuning parameters:
+##
+## mtry Accuracy Kappa
+## 2 0.9643726 0.9524051
+## 10 0.9780071 0.9707483
+## 19 0.9801973 0.9736577
+## 28 0.9787052 0.9716635
+## 37 0.9779095 0.9705890
+##
+## Accuracy was used to select the optimal model using the largest value.
+## The final value used for the model was mtry = 19.Training the classifier took +11.77 minutes.
+We next observe the accuracy of the classifer when predicting cell phenotypes +across the cross-validation and when applying the classifier to the test +dataset.
+First, we can visualize the classification accuracy during parameter +tuning:
+ggplot(rffit) +
+ geom_errorbar(data = rffit$results,
+ aes(ymin = Accuracy - AccuracySD,
+ ymax = Accuracy + AccuracySD),
+ width = 0.4) +
+ theme_classic(base_size = 15)
The best value for mtry is 19 and is used when predicting new data.
It is often recommended to visualize the variable importance of the +classifier. The following plot specifies which variables (markers) are +most important for classifying the data.
+ +
As expected, the markers that were used for gating (Ecad, CD3, CD20, HLADR, +CD8a, CD38, FOXP3) were important for classification.
+To assess the accuracy, sensitivity, specificity, among other quality measures of +the classifier, we will now predict cell phenotypes in the test data.
+# Select the arsinh-transformed counts of the test data
+cur_mat <- t(assay(test_spe, "exprs")[rowData(test_spe)$use_channel,])
+
+# Predict the cell phenotype labels of the test data
+set.seed(231019)
+cur_pred <- predict(rffit, newdata = cur_mat)While the overall classification accuracy can appear high, we also want
+to check if each cell phenotype class is correctly predicted.
+For this, we will calculate the confusion matrix between predicted and actual
+cell labels. This measure may highlight individual cell phenotype classes that
+were not correctly predicted by the classifier. When setting mode = "everything",
+the confusionMatrix function returns all available prediction measures including
+sensitivity, specificity, precision, recall and the F1 score per cell
+phenotype class.
cm <- confusionMatrix(data = cur_pred,
+ reference = factor(test_spe$cell_labels),
+ mode = "everything")
+
+cm## Confusion Matrix and Statistics
+##
+## Reference
+## Prediction Bcell BnTcell CD4 CD8 Myeloid Neutrophil Plasma_cell Stroma
+## Bcell 186 2 0 0 0 0 6 0
+## BnTcell 4 423 1 0 0 0 0 0
+## CD4 0 0 163 0 0 2 3 2
+## CD8 0 0 0 199 0 0 8 0
+## Myeloid 0 0 2 1 437 0 0 0
+## Neutrophil 0 0 0 0 0 30 0 0
+## Plasma_cell 1 0 3 2 0 0 158 0
+## Stroma 0 0 2 0 0 0 0 108
+## Treg 0 0 0 0 0 0 3 0
+## Tumor 4 0 1 3 0 1 1 0
+## Reference
+## Prediction Treg Tumor
+## Bcell 0 1
+## BnTcell 0 1
+## CD4 0 5
+## CD8 0 3
+## Myeloid 0 0
+## Neutrophil 0 0
+## Plasma_cell 1 0
+## Stroma 0 0
+## Treg 89 2
+## Tumor 0 1487
+##
+## Overall Statistics
+##
+## Accuracy : 0.9806
+## 95% CI : (0.9753, 0.985)
+## No Information Rate : 0.4481
+## P-Value [Acc > NIR] : < 2.2e-16
+##
+## Kappa : 0.9741
+##
+## Mcnemar's Test P-Value : NA
+##
+## Statistics by Class:
+##
+## Class: Bcell Class: BnTcell Class: CD4 Class: CD8
+## Sensitivity 0.95385 0.9953 0.94767 0.97073
+## Specificity 0.99714 0.9979 0.99622 0.99650
+## Pos Pred Value 0.95385 0.9860 0.93143 0.94762
+## Neg Pred Value 0.99714 0.9993 0.99716 0.99809
+## Precision 0.95385 0.9860 0.93143 0.94762
+## Recall 0.95385 0.9953 0.94767 0.97073
+## F1 0.95385 0.9906 0.93948 0.95904
+## Prevalence 0.05830 0.1271 0.05142 0.06129
+## Detection Rate 0.05561 0.1265 0.04873 0.05949
+## Detection Prevalence 0.05830 0.1283 0.05232 0.06278
+## Balanced Accuracy 0.97549 0.9966 0.97195 0.98361
+## Class: Myeloid Class: Neutrophil Class: Plasma_cell
+## Sensitivity 1.0000 0.909091 0.88268
+## Specificity 0.9990 1.000000 0.99779
+## Pos Pred Value 0.9932 1.000000 0.95758
+## Neg Pred Value 1.0000 0.999095 0.99340
+## Precision 0.9932 1.000000 0.95758
+## Recall 1.0000 0.909091 0.88268
+## F1 0.9966 0.952381 0.91860
+## Prevalence 0.1306 0.009865 0.05351
+## Detection Rate 0.1306 0.008969 0.04723
+## Detection Prevalence 0.1315 0.008969 0.04933
+## Balanced Accuracy 0.9995 0.954545 0.94024
+## Class: Stroma Class: Treg Class: Tumor
+## Sensitivity 0.98182 0.98889 0.9920
+## Specificity 0.99938 0.99846 0.9946
+## Pos Pred Value 0.98182 0.94681 0.9933
+## Neg Pred Value 0.99938 0.99969 0.9935
+## Precision 0.98182 0.94681 0.9933
+## Recall 0.98182 0.98889 0.9920
+## F1 0.98182 0.96739 0.9927
+## Prevalence 0.03288 0.02691 0.4481
+## Detection Rate 0.03229 0.02661 0.4445
+## Detection Prevalence 0.03288 0.02810 0.4475
+## Balanced Accuracy 0.99060 0.99368 0.9933To easily visualize these results, we can now plot the true positive rate +(sensitivity) versus the false positive rate (1 - specificity). The size of the +point is determined by the number of true positives divided by the total number +of cells.
+library(tidyverse)
+
+data.frame(cm$byClass) %>%
+ mutate(class = sub("Class: ", "", rownames(cm$byClass))) %>%
+ ggplot() +
+ geom_point(aes(1 - Specificity, Sensitivity,
+ size = Detection.Rate,
+ fill = class),
+ shape = 21) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme_classic(base_size = 15) +
+ ylab("Sensitivity (TPR)") +
+ xlab("1 - Specificity (FPR)")
We observe high sensitivity and specificity for most cell types. Plasma cells +show the lowest true positive rate with 88% being sufficiently high.
+Finally, to observe which cell phenotypes were wrongly classified, we can visualize +the distribution of classification probabilities per cell phenotype class:
+set.seed(231019)
+cur_pred <- predict(rffit,
+ newdata = cur_mat,
+ type = "prob")
+cur_pred$truth <- factor(test_spe$cell_labels)
+
+cur_pred %>%
+ pivot_longer(cols = Bcell:Tumor) %>%
+ ggplot() +
+ geom_boxplot(aes(x = name, y = value, fill = name), outlier.size = 0.5) +
+ facet_wrap(. ~ truth, ncol = 1) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme(panel.background = element_blank(),
+ axis.text.x = element_text(angle = 45, hjust = 1))
The boxplots indicate the classification probabilities per class. The classifier +is well trained if classification probabilities are only high for the one +specific class.
+In the final section, we will now use the tuned and tested random forest +classifier to predict the cell phenotypes of the unlabeled data.
+First, we predict the cell phenotypes and extract their classification +probabilities.
+# Select the arsinh-transformed counts of the unlabeled data for prediction
+cur_mat <- t(assay(unlab_spe, "exprs")[rowData(unlab_spe)$use_channel,])
+
+# Predict the cell phenotype labels of the unlabeled data
+set.seed(231014)
+cell_class <- as.character(predict(rffit,
+ newdata = cur_mat,
+ type = "raw"))
+names(cell_class) <- rownames(cur_mat)
+
+table(cell_class)## cell_class
+## Bcell BnTcell CD4 CD8 Myeloid Neutrophil
+## 817 979 3620 2716 6302 559
+## Plasma_cell Stroma Treg Tumor
+## 2692 4904 1170 10641# Extract prediction probabilities for each cell
+set.seed(231014)
+cell_prob <- predict(rffit,
+ newdata = cur_mat,
+ type = "prob")Each cell is assigned to the class with highest probability. There are however
+cases, where the highest probability is low meaning the cell can not be uniquely
+assigned to a class. We next want to identify these cells and label them as
+“undefined”. Here, we select a maximum classification probability threshold
+of 40% but this threshold needs to be adjusted for other datasets. The adjusted
+cell labels are then stored in the SpatialExperiment object.
library(ggridges)
+
+# Distribution of maximum probabilities
+tibble(max_prob = rowMax(as.matrix(cell_prob)),
+ type = cell_class) %>%
+ ggplot() +
+ geom_density_ridges(aes(x = max_prob, y = cell_class, fill = cell_class)) +
+ scale_fill_manual(values = metadata(spe)$color_vectors$celltype) +
+ theme_classic(base_size = 15) +
+ xlab("Maximum probability") +
+ ylab("Cell type") +
+ xlim(c(0,1.2))## Picking joint bandwidth of 0.0238
# Label undefined cells
+cell_class[rowMax(as.matrix(cell_prob)) < 0.4] <- "undefined"
+
+# Store labels in SpatialExperiment onject
+cell_labels <- spe$cell_labels
+cell_labels[colnames(unlab_spe)] <- cell_class
+spe$celltype <- cell_labels
+
+table(spe$celltype, spe$patient_id)##
+## Patient1 Patient2 Patient3 Patient4
+## Bcell 179 527 431 458
+## BnTcell 416 586 594 1078
+## CD4 391 1370 699 1385
+## CD8 518 1365 479 1142
+## Myeloid 1369 2197 1723 2731
+## Neutrophil 348 9 148 176
+## Plasma_cell 650 2122 351 274
+## Stroma 633 676 736 3261
+## Treg 553 409 243 310
+## Tumor 5560 3334 5648 2083
+## undefined 129 202 80 221We can now compare the cell labels derived by classification to the different +clustering strategies. The first comparison is against the clustering results +using the asinh-transformed counts.
+tab1 <- table(spe$celltype,
+ paste("Rphenograph", spe$pg_clusters))
+tab2 <- table(spe$celltype,
+ paste("SNN", spe$nn_clusters))
+tab3 <- table(spe$celltype,
+ paste("SOM", spe$som_clusters))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
We can see that Tumor and Myeloid cells span multiple clusters while +Neutrophiles are detected as an individual cluster by all clustering approaches.
+We next compare the cell classification against clustering results using the +integrated cells.
+tab1 <- table(spe$celltype,
+ paste("Rphenograph", spe$pg_clusters_corrected))
+tab2 <- table(spe$celltype,
+ paste("SNN", spe$nn_clusters_corrected))
+tab3 <- table(spe$celltype,
+ paste("SOM", spe$som_clusters_corrected))
+
+pheatmap(log10(tab1 + 10), color = viridis(100))
We observe a high agreement between the shared nearest neighbor clustering +approach using the integrated cells and the cell phenotypes derived by +classification.
+In the next sections, we will highlight visualization strategies to verify the +correctness of the phenotyping approach. Specifically, Section +11.2.3 shows how to outline identified cell phenotypes on +composite images.
+Finally, we save the updated SpatialExperiment object.
## R version 4.3.1 (2023-06-16)
+## Platform: x86_64-pc-linux-gnu (64-bit)
+## Running under: Ubuntu 22.04.3 LTS
+##
+## Matrix products: default
+## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
+## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
+##
+## locale:
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [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: Etc/UTC
+## tzcode source: system (glibc)
+##
+## attached base packages:
+## [1] stats4 stats graphics grDevices utils datasets methods
+## [8] base
+##
+## other attached packages:
+## [1] testthat_3.1.10 ggridges_0.5.4
+## [3] lubridate_1.9.3 forcats_1.0.0
+## [5] stringr_1.5.0 purrr_1.0.2
+## [7] readr_2.1.4 tidyr_1.3.0
+## [9] tibble_3.2.1 tidyverse_2.0.0
+## [11] caret_6.0-94 lattice_0.21-8
+## [13] cytomapper_1.12.0 EBImage_4.42.0
+## [15] dplyr_1.1.3 gridExtra_2.3
+## [17] pheatmap_1.0.12 patchwork_1.1.3
+## [19] ConsensusClusterPlus_1.64.0 kohonen_3.0.12
+## [21] CATALYST_1.24.0 scran_1.28.2
+## [23] scuttle_1.10.2 BiocParallel_1.34.2
+## [25] bluster_1.10.0 viridis_0.6.4
+## [27] viridisLite_0.4.2 dittoSeq_1.12.1
+## [29] ggplot2_3.4.3 igraph_1.5.1
+## [31] Rphenograph_0.99.1.9003 SpatialExperiment_1.10.0
+## [33] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
+## [35] Biobase_2.60.0 GenomicRanges_1.52.0
+## [37] GenomeInfoDb_1.36.3 IRanges_2.34.1
+## [39] S4Vectors_0.38.2 BiocGenerics_0.46.0
+## [41] MatrixGenerics_1.12.3 matrixStats_1.0.0
+##
+## loaded via a namespace (and not attached):
+## [1] bitops_1.0-7 RColorBrewer_1.1-3
+## [3] doParallel_1.0.17 tools_4.3.1
+## [5] backports_1.4.1 utf8_1.2.3
+## [7] R6_2.5.1 HDF5Array_1.28.1
+## [9] rhdf5filters_1.12.1 GetoptLong_1.0.5
+## [11] withr_2.5.1 sp_2.0-0
+## [13] cli_3.6.1 sandwich_3.0-2
+## [15] labeling_0.4.3 sass_0.4.7
+## [17] nnls_1.5 mvtnorm_1.2-3
+## [19] randomForest_4.7-1.1 proxy_0.4-27
+## [21] systemfonts_1.0.4 colorRamps_2.3.1
+## [23] svglite_2.1.1 R.utils_2.12.2
+## [25] scater_1.28.0 parallelly_1.36.0
+## [27] plotrix_3.8-2 limma_3.56.2
+## [29] flowCore_2.12.2 rstudioapi_0.15.0
+## [31] generics_0.1.3 shape_1.4.6
+## [33] gtools_3.9.4 car_3.1-2
+## [35] Matrix_1.6-1.1 RProtoBufLib_2.12.1
+## [37] waldo_0.5.1 ggbeeswarm_0.7.2
+## [39] fansi_1.0.4 abind_1.4-5
+## [41] R.methodsS3_1.8.2 terra_1.7-46
+## [43] lifecycle_1.0.3 multcomp_1.4-25
+## [45] yaml_2.3.7 edgeR_3.42.4
+## [47] carData_3.0-5 rhdf5_2.44.0
+## [49] recipes_1.0.8 Rtsne_0.16
+## [51] grid_4.3.1 promises_1.2.1
+## [53] dqrng_0.3.1 crayon_1.5.2
+## [55] shinydashboard_0.7.2 beachmat_2.16.0
+## [57] cowplot_1.1.1 magick_2.8.0
+## [59] pillar_1.9.0 knitr_1.44
+## [61] ComplexHeatmap_2.16.0 metapod_1.8.0
+## [63] rjson_0.2.21 future.apply_1.11.0
+## [65] codetools_0.2-19 glue_1.6.2
+## [67] data.table_1.14.8 vctrs_0.6.3
+## [69] png_0.1-8 gtable_0.3.4
+## [71] cachem_1.0.8 gower_1.0.1
+## [73] xfun_0.40 S4Arrays_1.0.6
+## [75] mime_0.12 prodlim_2023.08.28
+## [77] DropletUtils_1.20.0 survival_3.5-5
+## [79] timeDate_4022.108 iterators_1.0.14
+## [81] cytolib_2.12.1 hardhat_1.3.0
+## [83] lava_1.7.2.1 statmod_1.5.0
+## [85] ellipsis_0.3.2 TH.data_1.1-2
+## [87] ipred_0.9-14 nlme_3.1-162
+## [89] rprojroot_2.0.3 bslib_0.5.1
+## [91] irlba_2.3.5.1 svgPanZoom_0.3.4
+## [93] vipor_0.4.5 rpart_4.1.19
+## [95] colorspace_2.1-0 raster_3.6-23
+## [97] nnet_7.3-19 tidyselect_1.2.0
+## [99] compiler_4.3.1 BiocNeighbors_1.18.0
+## [101] desc_1.4.2 DelayedArray_0.26.7
+## [103] bookdown_0.35 scales_1.2.1
+## [105] tiff_0.1-11 digest_0.6.33
+## [107] fftwtools_0.9-11 rmarkdown_2.25
+## [109] XVector_0.40.0 htmltools_0.5.6
+## [111] pkgconfig_2.0.3 jpeg_0.1-10
+## [113] sparseMatrixStats_1.12.2 fastmap_1.1.1
+## [115] rlang_1.1.1 GlobalOptions_0.1.2
+## [117] htmlwidgets_1.6.2 shiny_1.7.5
+## [119] DelayedMatrixStats_1.22.6 farver_2.1.1
+## [121] jquerylib_0.1.4 zoo_1.8-12
+## [123] jsonlite_1.8.7 ModelMetrics_1.2.2.2
+## [125] R.oo_1.25.0 BiocSingular_1.16.0
+## [127] RCurl_1.98-1.12 magrittr_2.0.3
+## [129] GenomeInfoDbData_1.2.10 Rhdf5lib_1.22.1
+## [131] munsell_0.5.0 Rcpp_1.0.11
+## [133] ggnewscale_0.4.9 pROC_1.18.4
+## [135] stringi_1.7.12 brio_1.1.3
+## [137] zlibbioc_1.46.0 MASS_7.3-60
+## [139] plyr_1.8.8 listenv_0.9.0
+## [141] parallel_4.3.1 ggrepel_0.9.3
+## [143] splines_4.3.1 hms_1.1.3
+## [145] circlize_0.4.15 locfit_1.5-9.8
+## [147] ggpubr_0.6.0 ggsignif_0.6.4
+## [149] pkgload_1.3.3 reshape2_1.4.4
+## [151] ScaledMatrix_1.8.1 XML_3.99-0.14
+## [153] drc_3.0-1 evaluate_0.21
+## [155] tzdb_0.4.0 foreach_1.5.2
+## [157] tweenr_2.0.2 httpuv_1.6.11
+## [159] RANN_2.6.1 polyclip_1.10-6
+## [161] future_1.33.0 clue_0.3-65
+## [163] ggforce_0.4.1 rsvd_1.0.5
+## [165] broom_1.0.5 xtable_1.8-4
+## [167] e1071_1.7-13 rstatix_0.7.2
+## [169] later_1.3.1 class_7.3-22
+## [171] FlowSOM_2.8.0 beeswarm_0.4.0
+## [173] cluster_2.1.4 timechange_0.2.0
+## [175] globals_0.16.2The following section discusses possible quality indicators for data obtained +by IMC and other highly multiplexed imaging technologies. Here, we will focus +on describing quality metrics on the single-cell as well as image level.
+ +The first step after image segmentation is to observe its accuracy.
+Without having ground-truth data readily available, a common approach to
+segmentation quality control is to overlay segmentation masks on composite images
+displaying channels that were used for segmentation.
+The cytomapper
+package supports exactly this tasks by using the plotPixels function.
Here, we select 3 random images and perform image- and channel-wise +normalization (channels are first min-max normalized and scaled to a range of +0-1 before clipping the maximum intensity to 0.2).
+library(cytomapper)
+set.seed(20220118)
+img_ids <- sample(seq_along(images), 3)
+
+# Normalize and clip images
+cur_images <- images[img_ids]
+cur_images <- cytomapper::normalize(cur_images, separateImages = TRUE)
+cur_images <- cytomapper::normalize(cur_images, inputRange = c(0, 0.2))
+
+plotPixels(cur_images,
+ mask = masks[img_ids],
+ img_id = "sample_id",
+ missing_colour = "white",
+ colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
+ colour = list(CD163 = c("black", "yellow"),
+ CD20 = c("black", "red"),
+ CD3 = c("black", "green"),
+ Ecad = c("black", "cyan"),
+ DNA1 = c("black", "blue")),
+ image_title = NULL,
+ legend = list(colour_by.title.cex = 0.7,
+ colour_by.labels.cex = 0.7))
We can see that nuclei are centered within the segmentation masks and all cell
+types are correctly segmented (note: to zoom into the image you can right click
+and select Open Image in New Tab). A common challenge here is to segment large (e.g.,
+epithelial cells - in cyan) versus small (e.g., B cells - in red). However, the
+segmentation approach here appears to correctly segment cells across different
+sizes.
An easier and interactive way of observing segmentation quality is to use the +interactive image viewer provided by the +cytoviewer R/Bioconductor +package (Meyer, Eling, and Bodenmiller 2023). Under “Image-level” > “Basic controls”, up to six markers +can be selected for visualization. The contrast of each marker can be adjusted. +Under “Image-level” > “Advanced controls”, click the “Show cell outlines” box +to outline segmented cells on the images.
+library(cytoviewer)
+
+app <- cytoviewer(image = images,
+ mask = masks,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id")
+
+if (interactive()) {
+ shiny::runApp(app, launch.browser = TRUE)
+}An additional approach to observe cell segmentation quality and potentially also +antibody specificity issues is to visualize single-cell expression in form of a +heatmap. Here, we sub-sample the dataset to 2000 cells for visualization +purposes and overlay the cancer type from which the cells were extracted.
+library(dittoSeq)
+library(viridis)
+cur_cells <- sample(seq_len(ncol(spe)), 2000)
+
+dittoHeatmap(spe[,cur_cells],
+ genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs",
+ cluster_cols = TRUE,
+ scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = "indication",
+ annotation_colors = list(indication = metadata(spe)$color_vectors$indication))
We can differentiate between epithelial cells (Ecad+) and immune cells +(CD45RO+). Some of the markers are detected in specific cells (e.g., Ki67, CD20, +Ecad) while others are more broadly expressed across cells (e.g., HLADR, B2M, +CD4).
+Image-level quality control is often performed using tools that offer a +graphical user interface such as QuPath, +FIJI and the previously mentioned +cytoviewer package. Viewers +that were specifically developed for IMC data can be seen +here. In this +section, we will specifically focus on quantitative metrics to assess image +quality.
+It is often of interest to calculate the signal-to-noise ratio (SNR) for +individual channels and markers. Here, we define the SNR as:
+\[SNR = I_s/I_n\]
+where \(I_s\) is the intensity of the signal (mean intensity of pixels with true
+signal) and \(I_n\) is the intensity of the noise (mean intensity of pixels
+containing noise). This definition of the SNR is just one of many and other
+measures can be applied. Finding a threshold that separates pixels containing
+signal and pixels containing noise is not trivial and different approaches can
+be chosen. Here, we use the otsu thresholding approach to find pixels of the
+“foreground” (i.e., signal) and “background” (i.e., noise). The SNR is then
+defined as the mean intensity of foreground pixels divided by the mean intensity
+of background pixels. We compute this measure as well as the mean signal
+intensity per image. The plot below shows the average SNR versus the average
+signal intensity across all images.
library(tidyverse)
+library(ggrepel)
+library(EBImage)
+
+cur_snr <- lapply(names(images), function(x){
+ img <- images[[x]]
+ mat <- apply(img, 3, function(ch){
+ # Otsu threshold
+ thres <- otsu(ch, range = c(min(ch), max(ch)), levels = 65536)
+ # Signal-to-noise ratio
+ snr <- mean(ch[ch > thres]) / mean(ch[ch <= thres])
+ # Signal intensity
+ ps <- mean(ch[ch > thres])
+
+ return(c(snr = snr, ps = ps))
+ })
+ t(mat) %>% as.data.frame() %>%
+ mutate(image = x,
+ marker = colnames(mat)) %>%
+ pivot_longer(cols = c(snr, ps))
+})
+
+cur_snr <- do.call(rbind, cur_snr)
+
+cur_snr %>%
+ group_by(marker, name) %>%
+ summarize(log_mean = log2(mean(value))) %>%
+ pivot_wider(names_from = name, values_from = log_mean) %>%
+ ggplot() +
+ geom_point(aes(ps, snr)) +
+ geom_label_repel(aes(ps, snr, label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
We observe PD1, LAG3 and cleaved PARP to have high SNR but low signal intensity +meaning that in general these markers are not abundantly expressed. The Iridium +intercalator (here marked as DNA1 and DNA2) has the highest signal intensity +but low SNR. This might be due to staining differences between individual nuclei +where some nuclei are considered as background. We do however observe high +SNR and sufficient signal intensity for the majority of markers.
+Otsu thesholding and SNR calculation does not perform well if the markers are +lowly abundant. In the next code chunk, we will remove markers that have +a positive signal of below 2 per image.
+cur_snr <- cur_snr %>%
+ pivot_wider(names_from = name, values_from = value) %>%
+ filter(ps > 2) %>%
+ pivot_longer(cols = c(snr, ps))
+
+cur_snr %>%
+ group_by(marker, name) %>%
+ summarize(log_mean = log2(mean(value))) %>%
+ pivot_wider(names_from = name, values_from = log_mean) %>%
+ ggplot() +
+ geom_point(aes(ps, snr)) +
+ geom_label_repel(aes(ps, snr, label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
This visualization shows a reduces SNR for PD1, LAG3 and cleaved PARP which was +previously inflated due to low signal.
+Another quality indicator is the image area covered by cells (or biological
+tissue). This metric identifies ROIs where little cells are present, possibly
+hinting at incorrect selection of the ROI. We can compute the percentage of
+covered image area using the metadata contained in the SpatialExperiment
+object:
cell_density <- colData(spe) %>%
+ as.data.frame() %>%
+ group_by(sample_id) %>%
+ # Compute the number of pixels covered by cells and
+ # the total number of pixels
+ summarize(cell_area = sum(area),
+ no_pixels = mean(width_px) * mean(height_px)) %>%
+ # Divide the total number of pixels
+ # by the number of pixels covered by cells
+ mutate(covered_area = cell_area / no_pixels)
+
+# Visualize the image area covered by cells per image
+ggplot(cell_density) +
+ geom_point(aes(reorder(sample_id,covered_area), covered_area)) +
+ theme_minimal(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
+ ylim(c(0, 1)) +
+ ylab("% covered area") + xlab("")
We observe that two of the 14 images show unusually low cell coverage. These
+two images can now be visualized using cytomapper.
# Normalize and clip images
+cur_images <- images[c("Patient4_005", "Patient4_007")]
+cur_images <- cytomapper::normalize(cur_images, separateImages = TRUE)
+cur_images <- cytomapper::normalize(cur_images, inputRange = c(0, 0.2))
+
+plotPixels(cur_images,
+ mask = masks[c("Patient4_005", "Patient4_007")],
+ img_id = "sample_id",
+ missing_colour = "white",
+ colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
+ colour = list(CD163 = c("black", "yellow"),
+ CD20 = c("black", "red"),
+ CD3 = c("black", "green"),
+ Ecad = c("black", "cyan"),
+ DNA1 = c("black", "blue")),
+ legend = list(colour_by.title.cex = 0.7,
+ colour_by.labels.cex = 0.7))
These two images display less dense tissue structure but overall the images are +intact and appear to be segmented correctly.
+Finally, it can be beneficial to visualize the mean marker expression per image
+to identify images with outlying marker expression. This check does not
+indicate image quality per se but can highlight biological differences. Here,
+we will use the aggregateAcrossCells function of the
+scuttle package to compute the mean expression per
+image. For visualization purposes, we again asinh transform the mean expression
+values.
library(scuttle)
+
+image_mean <- aggregateAcrossCells(spe,
+ ids = spe$sample_id,
+ statistics="mean",
+ use.assay.type = "counts")
+assay(image_mean, "exprs") <- asinh(counts(image_mean))
+
+dittoHeatmap(image_mean, genes = rownames(spe)[rowData(spe)$use_channel],
+ assay = "exprs", cluster_cols = TRUE, scale = "none",
+ heatmap.colors = viridis(100),
+ annot.by = c("indication", "patient_id", "ROI"),
+ annotation_colors = list(indication = metadata(spe)$color_vectors$indication,
+ patient_id = metadata(spe)$color_vectors$patient_id,
+ ROI = metadata(spe)$color_vectors$ROI),
+ show_colnames = TRUE)
We observe extensive biological variation across the 14 images specifically for +some of the cell phenotype markers including the macrophage marker CD206, the B +cell marker CD20, the neutrophil marker CD15, and the proliferation marker Ki67. +These differences will be further studied in the following chapters.
+In the following paragraphs we will look at different metrics and visualization +approaches to assess data quality (as well as biological differences) on the +single-cell level.
+Related to the signal-to-noise ratio (SNR) calculated above on the pixel-level, +a similar measure can be derived on the single-cell level. Here, we will use +a two component Gaussian mixture model for each marker to find cells +with positive and negative expression. The SNR is defined as:
+\[SNR = I_s/I_n\]
+where \(I_s\) is the intensity of the signal (mean intensity of cells with
+positive signal) and \(I_n\) is the intensity of the noise (mean intensity of
+cells lacking expression). To define cells with positive and negative marker
+expression, we fit the mixture model across the transformed counts of all cells
+contained in the SpatialExperiment object. Next, for each marker we calculate
+the mean of the non-transformed counts for the positive and the negative cells.
+The SNR is then the ratio between the mean of the positive signal and the mean
+of the negative signal.
library(mclust)
+
+set.seed(220224)
+mat <- sapply(seq_len(nrow(spe)), function(x){
+ cur_exprs <- assay(spe, "exprs")[x,]
+ cur_counts <- assay(spe, "counts")[x,]
+
+ cur_model <- Mclust(cur_exprs, G = 2)
+ mean1 <- mean(cur_counts[cur_model$classification == 1])
+ mean2 <- mean(cur_counts[cur_model$classification == 2])
+
+ signal <- ifelse(mean1 > mean2, mean1, mean2)
+ noise <- ifelse(mean1 > mean2, mean2, mean1)
+
+ return(c(snr = signal/noise, ps = signal))
+})
+
+cur_snr <- t(mat) %>% as.data.frame() %>%
+ mutate(marker = rownames(spe))
+
+cur_snr %>% ggplot() +
+ geom_point(aes(log2(ps), log2(snr))) +
+ geom_label_repel(aes(log2(ps), log2(snr), label = marker)) +
+ theme_minimal(base_size = 15) + ylab("Signal-to-noise ratio [log2]") +
+ xlab("Signal intensity [log2]")
Next, we observe the distributions of cell size across the individual images. +Differences in cell size distributions can indicate segmentation biases due to +differences in cell density or can indicate biological differences due to cell +type compositions (tumor cells tend to be larger than immune cells).
+dittoPlot(spe, var = "area",
+ group.by = "sample_id",
+ plots = "boxplot") +
+ ylab("Cell area") + xlab("")
## Min. 1st Qu. Median Mean 3rd Qu. Max.
+## 3.00 47.00 70.00 76.38 98.00 466.00The median cell size is 70 pixels with a median major axis
+length of 11.3. The largest cell
+has an area of 466 pixels which relates to a diameter of
+21.6 pixels assuming a circular shape.
+Overall, the distribution of cell sizes is similar across images with images from
+Patient4_005 and Patient4_007 showing a reduced average cell size. These
+images contain fewer tumor cells which can explain the smaller average cell size.
We detect very small cells in the dataset and will remove them. +The chosen threshold is arbitrary and needs to be adjusted per dataset.
+ +## [1] 65Another quality indicator can be an absolute measure of cell density often +reported in cells per mm\(^2\).
+cell_density <- colData(spe) %>%
+ as.data.frame() %>%
+ group_by(sample_id) %>%
+ summarize(cell_count = n(),
+ no_pixels = mean(width_px) * mean(height_px)) %>%
+ mutate(cells_per_mm2 = cell_count/(no_pixels/1000000))
+
+ggplot(cell_density) +
+ geom_point(aes(reorder(sample_id,cells_per_mm2), cells_per_mm2)) +
+ theme_minimal(base_size = 15) +
+ theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
+ ylab("Cells per mm2") + xlab("")
The number of cells per mm\(^2\) varies across images which also depends on the +number of tumor/non-tumor cells. As we can see in the following sections, some +immune cells appear in cell dense regions while other stromal regions are less +dense.
+The data presented here originate from samples from different locations with +potential differences in pre-processing and each sample was stained individually. +These (and other) technical aspects can induce staining differences between +samples or batches of samples. Observing potential staining differences can be +crucial to assess data quality. We will use ridgeline visualizations to check +differences in staining patterns:
+multi_dittoPlot(spe, vars = rownames(spe)[rowData(spe)$use_channel],
+ group.by = "patient_id", plots = "ridgeplot",
+ assay = "exprs",
+ color.panel = metadata(spe)$color_vectors$patient_id)
We observe variations in the distributions of marker expression across patients. +These variations may arise partly from different abundances of cells in +different images (e.g., Patient3 may have higher numbers of CD11c+ and PD1+ +cells) as well as staining differences between samples. While most of the +selected markers are specifically expressed in immune cell subtypes, we can see +that E-Cadherin (a marker for epithelial (tumor) cells) shows a similar +expression range across all patients.
+Finally, we will use non-linear dimensionality reduction methods to project
+cells from a high-dimensional (40) down to a low-dimensional (2) space. For this
+the scater package provides the runUMAP and
+runTSNE function. To ensure reproducibility, we will need to set a seed;
+however different seeds and different parameter settings (e.g., the perplexity
+parameter in the runTSNE function) need to be tested to avoid
+over-interpretation of visualization artefacts. For dimensionality reduction, we
+will use all channels that show biological variation across the dataset.
+However, marker selection can be performed with different biological questions
+in mind. Here, both the runUMAP and runTSNE function are not deterministic,
+meaning they produce different results across different runs. We therefore
+set a seed in this chunk for reproducibility purposes.
library(scater)
+
+set.seed(220225)
+spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs")
+spe <- runTSNE(spe, subset_row = rowData(spe)$use_channel, exprs_values = "exprs") After dimensionality reduction, the low-dimensional embeddings are stored in the
+reducedDim slot.
## List of length 2
+## names(2): UMAP TSNE## UMAP1 UMAP2
+## Patient1_001_1 -4.810167 -3.777362
+## Patient1_001_2 -4.397347 -3.456036
+## Patient1_001_3 -4.369883 -3.445561
+## Patient1_001_4 -4.081614 -3.162119
+## Patient1_001_5 -6.234012 -2.433976
+## Patient1_001_6 -5.666597 -3.428058Visualization of the low-dimensional embedding facilitates assessment of
+potential “batch effects”. The dittoDimPlot
+function allows flexible visualization. It returns ggplot objects which
+can be further modified.
library(patchwork)
+
+# visualize patient id
+p1 <- dittoDimPlot(spe, var = "patient_id", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on UMAP")
+p2 <- dittoDimPlot(spe, var = "patient_id", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$patient_id) +
+ ggtitle("Patient ID on TSNE")
+
+# visualize region of interest id
+p3 <- dittoDimPlot(spe, var = "ROI", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$ROI) +
+ ggtitle("ROI ID on UMAP")
+p4 <- dittoDimPlot(spe, var = "ROI", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$ROI) +
+ ggtitle("ROI ID on TSNE")
+
+# visualize indication
+p5 <- dittoDimPlot(spe, var = "indication", reduction.use = "UMAP", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$indication) +
+ ggtitle("Indication on UMAP")
+p6 <- dittoDimPlot(spe, var = "indication", reduction.use = "TSNE", size = 0.2) +
+ scale_color_manual(values = metadata(spe)$color_vectors$indication) +
+ ggtitle("Indication on TSNE")
+
+(p1 + p2) / (p3 + p4) / (p5 + p6)
# visualize marker expression
+p1 <- dittoDimPlot(spe, var = "Ecad", reduction.use = "UMAP",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "Ecad") +
+ ggtitle("E-Cadherin expression on UMAP")
+p2 <- dittoDimPlot(spe, var = "CD45RO", reduction.use = "UMAP",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "CD45RO") +
+ ggtitle("CD45RO expression on UMAP")
+p3 <- dittoDimPlot(spe, var = "Ecad", reduction.use = "TSNE",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "Ecad") +
+ ggtitle("Ecad expression on TSNE")
+p4 <- dittoDimPlot(spe, var = "CD45RO", reduction.use = "TSNE",
+ assay = "exprs", size = 0.2) +
+ scale_color_viridis(name = "CD45RO") +
+ ggtitle("CD45RO expression on TSNE")
+
+(p1 + p2) / (p3 + p4)
We observe a strong separation of tumor cells (Ecad+ cells) between the +patients. Here, each patient was diagnosed with a different tumor type. The +separation of tumor cells could be of biological origin since tumor cells tend +to display differences in expression between patients and cancer types and/or of +technical origin: the panel only contains a single tumor marker (E-Cadherin) and +therefore slight technical differences in staining causes visible separation +between cells of different patients. Nevertheless, the immune compartment +(CD45RO+ cells) mix between patients and we can rule out systematic staining +differences between patients.
+The modified SpatialExperiment object is saved for further downstream analysis.
## R version 4.3.1 (2023-06-16)
+## Platform: x86_64-pc-linux-gnu (64-bit)
+## Running under: Ubuntu 22.04.3 LTS
+##
+## Matrix products: default
+## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
+## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
+##
+## locale:
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [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: Etc/UTC
+## tzcode source: system (glibc)
+##
+## attached base packages:
+## [1] stats4 stats graphics grDevices utils datasets methods
+## [8] base
+##
+## other attached packages:
+## [1] testthat_3.1.10 patchwork_1.1.3
+## [3] scater_1.28.0 mclust_6.0.0
+## [5] scuttle_1.10.2 ggrepel_0.9.3
+## [7] lubridate_1.9.3 forcats_1.0.0
+## [9] stringr_1.5.0 dplyr_1.1.3
+## [11] purrr_1.0.2 readr_2.1.4
+## [13] tidyr_1.3.0 tibble_3.2.1
+## [15] tidyverse_2.0.0 viridis_0.6.4
+## [17] viridisLite_0.4.2 dittoSeq_1.12.1
+## [19] ggplot2_3.4.3 cytoviewer_1.0.1
+## [21] cytomapper_1.12.0 SingleCellExperiment_1.22.0
+## [23] SummarizedExperiment_1.30.2 Biobase_2.60.0
+## [25] GenomicRanges_1.52.0 GenomeInfoDb_1.36.3
+## [27] IRanges_2.34.1 S4Vectors_0.38.2
+## [29] BiocGenerics_0.46.0 MatrixGenerics_1.12.3
+## [31] matrixStats_1.0.0 EBImage_4.42.0
+##
+## loaded via a namespace (and not attached):
+## [1] RColorBrewer_1.1-3 rstudioapi_0.15.0
+## [3] jsonlite_1.8.7 magrittr_2.0.3
+## [5] ggbeeswarm_0.7.2 magick_2.8.0
+## [7] farver_2.1.1 rmarkdown_2.25
+## [9] zlibbioc_1.46.0 vctrs_0.6.3
+## [11] memoise_2.0.1 DelayedMatrixStats_1.22.6
+## [13] RCurl_1.98-1.12 terra_1.7-46
+## [15] svgPanZoom_0.3.4 htmltools_0.5.6
+## [17] S4Arrays_1.0.6 BiocNeighbors_1.18.0
+## [19] raster_3.6-23 Rhdf5lib_1.22.1
+## [21] rhdf5_2.44.0 sass_0.4.7
+## [23] bslib_0.5.1 desc_1.4.2
+## [25] htmlwidgets_1.6.2 fontawesome_0.5.2
+## [27] cachem_1.0.8 mime_0.12
+## [29] lifecycle_1.0.3 pkgconfig_2.0.3
+## [31] rsvd_1.0.5 colourpicker_1.3.0
+## [33] Matrix_1.6-1.1 R6_2.5.1
+## [35] fastmap_1.1.1 GenomeInfoDbData_1.2.10
+## [37] shiny_1.7.5 digest_0.6.33
+## [39] colorspace_2.1-0 shinycssloaders_1.0.0
+## [41] rprojroot_2.0.3 irlba_2.3.5.1
+## [43] dqrng_0.3.1 pkgload_1.3.3
+## [45] beachmat_2.16.0 labeling_0.4.3
+## [47] timechange_0.2.0 fansi_1.0.4
+## [49] nnls_1.5 abind_1.4-5
+## [51] compiler_4.3.1 withr_2.5.1
+## [53] tiff_0.1-11 BiocParallel_1.34.2
+## [55] HDF5Array_1.28.1 R.utils_2.12.2
+## [57] DelayedArray_0.26.7 rjson_0.2.21
+## [59] tools_4.3.1 vipor_0.4.5
+## [61] beeswarm_0.4.0 httpuv_1.6.11
+## [63] R.oo_1.25.0 glue_1.6.2
+## [65] rhdf5filters_1.12.1 promises_1.2.1
+## [67] grid_4.3.1 Rtsne_0.16
+## [69] generics_0.1.3 gtable_0.3.4
+## [71] tzdb_0.4.0 R.methodsS3_1.8.2
+## [73] hms_1.1.3 ScaledMatrix_1.8.1
+## [75] BiocSingular_1.16.0 sp_2.0-0
+## [77] utf8_1.2.3 XVector_0.40.0
+## [79] RcppAnnoy_0.0.21 pillar_1.9.0
+## [81] limma_3.56.2 later_1.3.1
+## [83] lattice_0.21-8 tidyselect_1.2.0
+## [85] locfit_1.5-9.8 miniUI_0.1.1.1
+## [87] knitr_1.44 gridExtra_2.3
+## [89] bookdown_0.35 edgeR_3.42.4
+## [91] svglite_2.1.1 xfun_0.40
+## [93] shinydashboard_0.7.2 brio_1.1.3
+## [95] DropletUtils_1.20.0 pheatmap_1.0.12
+## [97] stringi_1.7.12 fftwtools_0.9-11
+## [99] yaml_2.3.7 evaluate_0.21
+## [101] codetools_0.2-19 archive_1.1.6
+## [103] BiocManager_1.30.22 cli_3.6.1
+## [105] uwot_0.1.16 xtable_1.8-4
+## [107] systemfonts_1.0.4 munsell_0.5.0
+## [109] jquerylib_0.1.4 Rcpp_1.0.11
+## [111] png_0.1-8 parallel_4.3.1
+## [113] ellipsis_0.3.2 jpeg_0.1-10
+## [115] sparseMatrixStats_1.12.2 bitops_1.0-7
+## [117] SpatialExperiment_1.10.0 scales_1.2.1
+## [119] ggridges_0.5.4 crayon_1.5.2
+## [121] BiocStyle_2.28.1 rlang_1.1.1
+## [123] cowplot_1.1.1The following section describes how to visualize the abundance of biomolecules +(e.g., protein or RNA) as well as cell-specific metadata on images. Section +11.1 focuses on visualizing pixel-level information +including the generation of pseudo-color composite images. Section +11.2 highlights the visualization of cell metadata (e.g., +cell phenotype) as well as summarized pixel intensities on cell segmentation +masks.
+The +cytomapper +R/Bioconductor package was developed to support the handling and visualization +of multiple multi-channel images and segmentation masks (Eling et al. 2020). The main +data object for image handling is the +CytoImageList +container which we used in Section 5 to store multi-channel +images and segmentation masks.
+We will first read in the previously processed data and randomly select 3 images +for visualization purposes.
+library(SpatialExperiment)
+library(cytomapper)
+spe <- readRDS("data/spe.rds")
+images <- readRDS("data/images.rds")
+masks <- readRDS("data/masks.rds")
+
+# Sample images
+set.seed(220517)
+cur_id <- sample(unique(spe$sample_id), 3)
+cur_images <- images[names(images) %in% cur_id]
+cur_masks <- masks[names(masks) %in% cur_id]The following section gives examples for visualizing individual channels or
+multiple channels as pseudo-color composite images. For this the cytomapper
+package exports the plotPixels function which expects a CytoImageList object
+storing one or multiple multi-channel images. In the simplest use case, a
+single channel can be visualized as follows:
The plot above shows the tissue expression of the epithelial tumor marker
+E-cadherin on the 3 selected images. The bcg parameter (default c(0, 1, 1))
+stands for “background”, “contrast”, “gamma” and controls these attributes of
+the image. This parameter takes a named list where each entry specifies these
+attributes per channel. The first value of the numeric vector will be added to
+the pixel intensities (background); pixel intensities will be multiplied by the
+second entry of the vector (contrast); pixel intensities will be exponentiated
+by the third entry of the vector (gamma). In most cases, it is sufficient to
+adjust the second (contrast) entry of the vector.
The following example highlights the visualization of 6 markers (maximum allowed +number of markers) at once per image. The markers indicate the spatial +distribution of tumor cells (E-cadherin), T cells (CD3), B cells (CD20), CD8+ T +cells (CD8a), plasma cells (CD38) and proliferating cells (Ki67).
+plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1),
+ CD8a = c(0, 5, 1),
+ CD38 = c(0, 8, 1),
+ Ki67 = c(0, 5, 1)))
The default colors for visualization are chosen by the additive RGB (red, green,
+blue) color model. For six markers the default colors are: red, green, blue,
+cyan (green + blue), magenta (red + blue), yellow (green + red). These colors
+are the easiest to distinguish by eye. However, you can select other colors for
+each channel by setting the colour parameter:
plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1)),
+ colour = list(Ecad = c("black", "burlywood1"),
+ CD3 = c("black", "cyan2"),
+ CD20 = c("black", "firebrick1")))
The colour parameter takes a named list in which each entry specifies the
+colors from which a color gradient is constructed via colorRampPalette. These
+are usually vectors of length 2 in which the first entry is "black" and the
+second entry specifies the color of choice. Although not recommended, you can
+also specify more than two colors to generate a more complex color gradient.
As an alternative to setting the bcg parameter, images can first be
+normalized. Normalization here means to scale the pixel intensities per channel
+between 0 and 1 (or a range specified by the ft parameter in the normalize
+function). By default, the normalize function scales pixel intensities across
+all images contained in the CytoImageList object (separateImages = FALSE).
+Each individual channel is scaled independently (separateChannels = TRUE).
After 0-1 normalization, maximum pixel intensities can be clipped to enhance the
+contrast of the image (setting the inputRange parameter). In the following
+example, the clipping to 0 and 0.2 is the same as multiplying the pixel
+intensities by a factor of 5.
# 0 - 1 channel scaling across all images
+norm_images <- cytomapper::normalize(cur_images)
+
+# Clip channel at 0.2
+norm_images <- cytomapper::normalize(norm_images, inputRange = c(0, 0.2))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
The default setting of scaling pixel intensities across all images ensures +comparable intensity levels across images. Pixel intensities can also be +scaled per image therefore correcting for staining/expression differences +between images:
+# 0 - 1 channel scaling per image
+norm_images <- cytomapper::normalize(cur_images, separateImages = TRUE)
+
+# Clip channel at 0.2
+norm_images <- cytomapper::normalize(norm_images, inputRange = c(0, 0.2))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
As we can see, the marker Ki67 appears brighter on image 2 and 3 in comparison +to scaling the channel across all images.
+Finally, the normalize function also accepts a named list input for the
+inputRange argument. In this list, the clipping range per channel can be set
+individually:
# 0 - 1 channel scaling per image
+norm_images <- cytomapper::normalize(cur_images,
+ separateImages = TRUE,
+ inputRange = list(Ecad = c(0, 50),
+ CD3 = c(0, 30),
+ CD20 = c(0, 40),
+ CD8a = c(0, 50),
+ CD38 = c(0, 10),
+ Ki67 = c(0, 70)))
+
+plotPixels(norm_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"))
In the following section, we will show examples on how to visualize single +cells either as segmentation masks or outlined on composite images. This type +of visualization allows to observe the spatial distribution of cell phenotypes, +the visual assessment of morphological features and quality control in terms +of cell segmentation and phenotyping.
+The cytomapper package provides the plotCells function that accepts a
+CytoImageList object containing segmentation masks. These are defined as
+single channel images where sets of pixels with the same integer ID identify
+individual cells. This integer ID can be found as an entry in the colData(spe)
+slot and as pixel information in the segmentation masks. The entry in
+colData(spe) needs to be specified via the cell_id argument to the
+plotCells function. In that way, data contained in the SpatialExperiment
+object can be mapped to segmentation masks. For the current dataset, the cell
+IDs are stored in colData(spe)$ObjectNumber.
As cell IDs are only unique within a single image, plotCells also requires
+the img_id argument. This argument specifies the colData(spe) as well as the
+mcols(masks) entry that stores the unique image name from which each cell was
+extracted. In the current dataset the unique image names are stored in
+colData(spe)$sample_id and mcols(masks)$sample_id.
Providing these two entries that allow mapping between the SpatialExperiment
+object and segmentation masks, we can now color individual cells based on their
+cell type:
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype")
For consistent visualization, the plotCells function takes a named list as
+color argument. The entry name must match the colour_by argument.
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype))
If only individual cell types should be visualized, the SpatialExperiment
+object can be subsetted (e.g., to only contain CD8+ T cells). In the following
+example CD8+ T cells are colored in red and all other cells that are not
+contained in the dataset are colored in white (as set by the missing_color
+argument).
CD8 <- spe[,spe$celltype == "CD8"]
+
+plotCells(cur_masks,
+ object = CD8,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = c(CD8 = "red")),
+ missing_colour = "white")
In terms of visualizing metadata, any entry in the colData(spe) slot can be
+visualized. The plotCells function automatically detects if the entry
+is continuous or discrete. In this fashion, we can now visualize the area of each
+cell:
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "area")
Similar to visualizing single-cell metadata on segmentation masks, we can
+use the plotCells function to visualize the aggregated pixel intensities
+per cell. In the current dataset pixel intensities were aggregated by computing
+the mean pixel intensity per cell and per channel. The plotCells function
+accepts the exprs_values argument (default counts) that allows selecting
+the assay which stores the expression values that should be visualized.
In the following example, we visualize the asinh-transformed mean pixel +intensities of the epithelial marker E-cadherin on segmentation masks.
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "Ecad",
+ exprs_values = "exprs")
We will now visualize the maximum number of +allowed markers as composites on the segmentation masks. As above the markers +indicate the spatial distribution of tumor cells (E-cadherin), T cells (CD3), B +cells (CD20), CD8+ T cells (CD8a), plasma cells (CD38) and proliferating cells +(Ki67).
+plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ exprs_values = "exprs")
While visualizing 6 markers on the pixel-level may still allow the distinction +of different tissue structures, observing single-cell expression levels is +difficult when visualizing many markers simultaneously due to often overlapping +expression.
+Similarly to adjusting marker colors when visualizing pixel intensities, we
+can change the color gradients per marker by setting the color argument:
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ exprs_values = "exprs",
+ colour = list(Ecad = c("black", "burlywood1"),
+ CD3 = c("black", "cyan2"),
+ CD20 = c("black", "firebrick1")))
The following section highlights the combined visualization of pixel- and
+cell-level information at once. For this, besides the SpatialExperiment object,
+the plotPixels function accepts two CytoImageList objects. One for the
+multi-channel images and one for the segmentation masks. By specifying the
+outline_by parameter, the outlines of cells can now be colored based on their
+metadata.
The following example first generates a 3-channel composite images displaying +the expression of E-cadherin, CD3 and CD20 before coloring the cells’ outlines +by their cell phenotype.
+plotPixels(image = cur_images,
+ mask = cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ outline_by = "celltype",
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1)),
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ thick = TRUE)
Distinguishing individual cell phenotypes is nearly impossible in the images +above.
+However, the SpatialExperiment object can be subsetted to only contain cells
+of a single or few phenotypes. This allows the selective visualization of cell
+outlines on composite images.
Here, we select all CD8+ T cells from the dataset and outline them on a 2-channel +composite image displaying the expression of CD3 and CD8a.
+CD8 <- spe[,spe$celltype == "CD8"]
+
+plotPixels(image = cur_images,
+ mask = cur_masks,
+ object = CD8,
+ cell_id = "ObjectNumber", img_id = "sample_id",
+ colour_by = c("CD3", "CD8a"),
+ outline_by = "celltype",
+ bcg = list(CD3 = c(0, 5, 1),
+ CD8a = c(0, 5, 1)),
+ colour = list(celltype = c("CD8" = "white")),
+ thick = TRUE)
This type of visualization allows the quality control of two things: 1. +segmentation quality of individual cell types can be checked and 2. cell +phenotyping accuracy can be visually assessed against expected marker expression.
+The cytomapper package provides a number of function arguments to adjust the
+visual appearance of figures that are shared between the plotPixels and
+plotCells function.
For a full overview of the arguments please refer to ?plotting-param.
We use the following example to highlight how to adjust the scale bar, the image +title, the legend appearance and the margin between images.
+plotPixels(cur_images,
+ colour_by = c("Ecad", "CD3", "CD20", "CD8a", "CD38", "Ki67"),
+ bcg = list(Ecad = c(0, 5, 1),
+ CD3 = c(0, 5, 1),
+ CD20 = c(0, 5, 1),
+ CD8a = c(0, 5, 1),
+ CD38 = c(0, 8, 1),
+ Ki67 = c(0, 5, 1)),
+ scale_bar = list(length = 100,
+ label = expression("100 " ~ mu * "m"),
+ cex = 0.7,
+ lwidth = 10,
+ colour = "grey",
+ position = "bottomleft",
+ margin = c(5,5),
+ frame = 3),
+ image_title = list(text = mcols(cur_images)$indication,
+ position = "topright",
+ colour = "grey",
+ margin = c(5,5),
+ font = 2,
+ cex = 2),
+ legend = list(colour_by.title.cex = 0.7,
+ margin = 10),
+ margin = 40)
By default, all images are displayed on the same graphics device. This can be
+useful when saving all images at once (see next section) to zoom into the
+individual images instead of opening each image individually. However, when
+displaying images in a markdown document these are more accessible when
+visualized individually. For this, the plotPixels and plotCells function
+accepts the display parameter that when set to "single" displays each
+resulting image in its own graphics device:
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ display = "single",
+ legend = NULL)
The final section addresses how to save composite images and how to return them +for integration with other plots.
+The plotPixels and plotCells functions accept the save_plot argument which
+takes a named list of the following entries: filename indicates the location
+and file type of the image saved to disk; scale adjusts the resolution of the
+saved image (this only needs to be adjusted for small images).
plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ save_plot = list(filename = "data/celltype_image.png"))The composite images (together with their annotation) can also be returned. In
+the following code chunk we save two example plots to variables (out1 and
+out2).
out1 <- plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = "celltype",
+ colour = list(celltype = metadata(spe)$color_vectors$celltype),
+ return_plot = TRUE)out2 <- plotCells(cur_masks,
+ object = spe,
+ cell_id = "ObjectNumber",
+ img_id = "sample_id",
+ colour_by = c("Ecad", "CD3", "CD20"),
+ exprs_values = "exprs",
+ return_plot = TRUE)The composite images are stored in out1$plot and out2$plot and can be
+converted into a graph object recognized by the
+cowplot
+package.
The final function call of the following chunk plots both object next to each +other.
+library(cowplot)
+library(gridGraphics)
+p1 <- ggdraw(out1$plot, clip = "on")
+p2 <- ggdraw(out2$plot, clip = "on")
+
+plot_grid(p1, p2)
The +cytoviewer +package allows the interactive visualization of multi-channel images and +segmentation masks. It also allows to map cellular metadata onto segmentation +masks and outlining of cells on composite images. For a full introduction to the +package, please refer to the vignette.
+ +## R version 4.3.1 (2023-06-16)
+## Platform: x86_64-pc-linux-gnu (64-bit)
+## Running under: Ubuntu 22.04.3 LTS
+##
+## Matrix products: default
+## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
+## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
+##
+## locale:
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [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: Etc/UTC
+## tzcode source: system (glibc)
+##
+## attached base packages:
+## [1] grid stats4 stats graphics grDevices utils datasets
+## [8] methods base
+##
+## other attached packages:
+## [1] cytoviewer_1.0.1 gridGraphics_0.5-1
+## [3] cowplot_1.1.1 cytomapper_1.12.0
+## [5] EBImage_4.42.0 SpatialExperiment_1.10.0
+## [7] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
+## [9] Biobase_2.60.0 GenomicRanges_1.52.0
+## [11] GenomeInfoDb_1.36.3 IRanges_2.34.1
+## [13] S4Vectors_0.38.2 BiocGenerics_0.46.0
+## [15] MatrixGenerics_1.12.3 matrixStats_1.0.0
+##
+## loaded via a namespace (and not attached):
+## [1] splines_4.3.1 later_1.3.1
+## [3] bitops_1.0-7 tibble_3.2.1
+## [5] R.oo_1.25.0 svgPanZoom_0.3.4
+## [7] polyclip_1.10-6 XML_3.99-0.14
+## [9] lifecycle_1.0.3 rstatix_0.7.2
+## [11] edgeR_3.42.4 doParallel_1.0.17
+## [13] lattice_0.21-8 MASS_7.3-60
+## [15] backports_1.4.1 magrittr_2.0.3
+## [17] limma_3.56.2 sass_0.4.7
+## [19] rmarkdown_2.25 plotrix_3.8-2
+## [21] jquerylib_0.1.4 yaml_2.3.7
+## [23] httpuv_1.6.11 sp_2.0-0
+## [25] RColorBrewer_1.1-3 ConsensusClusterPlus_1.64.0
+## [27] multcomp_1.4-25 abind_1.4-5
+## [29] zlibbioc_1.46.0 Rtsne_0.16
+## [31] purrr_1.0.2 R.utils_2.12.2
+## [33] RCurl_1.98-1.12 TH.data_1.1-2
+## [35] tweenr_2.0.2 sandwich_3.0-2
+## [37] circlize_0.4.15 GenomeInfoDbData_1.2.10
+## [39] ggrepel_0.9.3 irlba_2.3.5.1
+## [41] CATALYST_1.24.0 terra_1.7-46
+## [43] dqrng_0.3.1 svglite_2.1.1
+## [45] DelayedMatrixStats_1.22.6 codetools_0.2-19
+## [47] DropletUtils_1.20.0 DelayedArray_0.26.7
+## [49] scuttle_1.10.2 ggforce_0.4.1
+## [51] tidyselect_1.2.0 shape_1.4.6
+## [53] raster_3.6-23 farver_2.1.1
+## [55] ScaledMatrix_1.8.1 viridis_0.6.4
+## [57] jsonlite_1.8.7 BiocNeighbors_1.18.0
+## [59] GetoptLong_1.0.5 ellipsis_0.3.2
+## [61] scater_1.28.0 ggridges_0.5.4
+## [63] survival_3.5-5 iterators_1.0.14
+## [65] systemfonts_1.0.4 foreach_1.5.2
+## [67] tools_4.3.1 ggnewscale_0.4.9
+## [69] Rcpp_1.0.11 glue_1.6.2
+## [71] gridExtra_2.3 xfun_0.40
+## [73] dplyr_1.1.3 HDF5Array_1.28.1
+## [75] shinydashboard_0.7.2 withr_2.5.1
+## [77] fastmap_1.1.1 rhdf5filters_1.12.1
+## [79] fansi_1.0.4 rsvd_1.0.5
+## [81] digest_0.6.33 R6_2.5.1
+## [83] mime_0.12 colorspace_2.1-0
+## [85] gtools_3.9.4 jpeg_0.1-10
+## [87] R.methodsS3_1.8.2 utf8_1.2.3
+## [89] tidyr_1.3.0 generics_0.1.3
+## [91] data.table_1.14.8 htmlwidgets_1.6.2
+## [93] S4Arrays_1.0.6 pkgconfig_2.0.3
+## [95] gtable_0.3.4 ComplexHeatmap_2.16.0
+## [97] RProtoBufLib_2.12.1 XVector_0.40.0
+## [99] htmltools_0.5.6 carData_3.0-5
+## [101] bookdown_0.35 fftwtools_0.9-11
+## [103] clue_0.3-65 scales_1.2.1
+## [105] png_0.1-8 colorRamps_2.3.1
+## [107] knitr_1.44 rstudioapi_0.15.0
+## [109] reshape2_1.4.4 rjson_0.2.21
+## [111] cachem_1.0.8 zoo_1.8-12
+## [113] rhdf5_2.44.0 GlobalOptions_0.1.2
+## [115] stringr_1.5.0 shinycssloaders_1.0.0
+## [117] miniUI_0.1.1.1 parallel_4.3.1
+## [119] vipor_0.4.5 pillar_1.9.0
+## [121] vctrs_0.6.3 promises_1.2.1
+## [123] ggpubr_0.6.0 BiocSingular_1.16.0
+## [125] car_3.1-2 cytolib_2.12.1
+## [127] beachmat_2.16.0 xtable_1.8-4
+## [129] cluster_2.1.4 archive_1.1.6
+## [131] beeswarm_0.4.0 evaluate_0.21
+## [133] magick_2.8.0 mvtnorm_1.2-3
+## [135] cli_3.6.1 locfit_1.5-9.8
+## [137] compiler_4.3.1 rlang_1.1.1
+## [139] crayon_1.5.2 ggsignif_0.6.4
+## [141] FlowSOM_2.8.0 plyr_1.8.8
+## [143] flowCore_2.12.2 ggbeeswarm_0.7.2
+## [145] stringi_1.7.12 viridisLite_0.4.2
+## [147] BiocParallel_1.34.2 nnls_1.5
+## [149] munsell_0.5.0 tiff_0.1-11
+## [151] colourpicker_1.3.0 Matrix_1.6-1.1
+## [153] sparseMatrixStats_1.12.2 ggplot2_3.4.3
+## [155] Rhdf5lib_1.22.1 shiny_1.7.5
+## [157] fontawesome_0.5.2 drc_3.0-1
+## [159] memoise_2.0.1 igraph_1.5.1
+## [161] broom_1.0.5 bslib_0.5.1Compiled: 2023-10-19
+This workflow highlights the use of common R/Bioconductor packages +to analyze single-cell data obtained from segmented multi-channel images. We will not perform multi-channel image processing and segmentation in R +but rather link to available approaches in Section 3. While we +use imaging mass cytometry (IMC) data as an example, the concepts presented here can be applied to images +obtained by other highly-multiplexed imaging technologies (e.g. CODEX, MIBI, +mIF, etc.).
+We will give an introduction to IMC in Section 2 and highlight +strategies to extract single-cell data from multi-channel images in Section +3.
+Reproducible code written in R is available from Section 4 +onwards and the workflow can be largely divided into the following parts:
+Multi-channel image and spatial, single-cell analysis is complex and we +highlight an example workflow here. However, this workflow is not complete and +does not cover all possible aspects of exploratory data analysis. Instead, we +demonstrate this workflow as a solid basis that supports other aspects of data +analysis. It offers interoperability with other packages for single-cell and +spatial analysis and the user will need to become familiar with the general +framework to efficiently analyse data obtained from multiplexed imaging +technologies.
+We provide the workflow as an open-source resource. It does not mean that +this workflow is tested on all possible datasets or biological questions.
+If you notice an issue or missing information, please report an issue +here. We also +welcome contributions in form of pull requests or feature requests in form of +issues. Have a look at the source code at:
+ +The workflow has been published in
+https://www.nature.com/articles/s41596-023-00881-0
+which you can cite as follows:
Windhager, J., Zanotelli, V.R.T., Schulz, D. et al. An end-to-end workflow for multiplexed image processing and analysis.
+ Nat Protoc (2023). Version 1.0.0 [2023-06-30]
+Version 1.0.1 [2023-10-19]
+predict call after training a classifier* nils.eling@uzh.ch
+1: Department for Quantitative Biomedicine, University of Zurich
+2: Institute for Molecular Health Sciences, ETH Zurich
Version 1.0.0 [2023-06-30]
+Version 1.0.1 [2023-10-19]
+predict call after training a classifierHighly multiplexed imaging (HMI) enables the simultaneous detection of dozens of +biological molecules (e.g., proteins, transcripts; also referred to as +“markers”) in tissues. Recently established multiplexed tissue imaging +technologies rely on cyclic staining with fluorescently-tagged antibodies +(Lin et al. 2018; Gut, Herrmann, and Pelkmans 2018), or the use of oligonucleotide-tagged (Goltsev et al. 2018; Saka et al. 2019) or metal-tagged (Giesen et al. 2014; Angelo et al. 2014) antibodies, among others. +The key strength of these technologies is that they allow in-depth analysis of +single cells within their spatial tissue context. As a result, these methods +have enabled analysis of the spatial architecture of the tumor microenvironment +(Lin et al. 2018; Jackson et al. 2020; Ali et al. 2020; Schürch et al. 2020), determination of nucleic acid +and protein abundances for assessment of spatial co-localization of cell types +and chemokines (Hoch et al. 2022) and spatial niches of virus infected cells (Jiang et al. 2022), +and characterization of pathological features during COVID-19 infection +(Rendeiro et al. 2021; Mitamura et al. 2021), Type 1 diabetes progression (Damond et al. 2019) and +autoimmune disease (Ferrian et al. 2021).
+Imaging mass cytometry (IMC) utilizes metal-tagged antibodies to detect over 40 +proteins and other metal-tagged molecules in biological samples. IMC can be used +to perform highly multiplexed imaging and is particularly suited to profiling +selected areas of tissues across many samples.
+
+Overview of imaging mass cytometry data acquisition. Taken from (Giesen et al. 2014)
IMC has first been published in 2014 (Giesen et al. 2014) and has been commercialized by +Standard BioToolsTM to be distributed as the Hyperion Imaging +SystemTM (documentation is available +here). +Similar to other HMI technologies such as MIBI (Angelo et al. 2014), CyCIF (Lin et al. 2018), +4i (Gut, Herrmann, and Pelkmans 2018), CODEX (Goltsev et al. 2018) and SABER (Saka et al. 2019), IMC captures the spatial +expression of multiple proteins in parallel. With a nominal 1 μm resolution, +IMC is able to detect cytoplasmic and nuclear localization of proteins. The +current ablation frequency of IMC is 200Hz, meaning that a 1 mm\(^2\) area +can be imaged within about 2 hours.
+Technical aspects of how data acquisition works can be found in the original +publication (Giesen et al. 2014). Briefly, antibodies to detect targets in biological +material are labeled with heavy metals (e.g., lanthanides) that do not occur in +biological systems and thus can be used upon binding to their target as a +readout similar to fluorophores in fluorescence microscopy. Thin sections of the +biological sample on a glass slide are stained with an antibody cocktail. +Stained microscopy slides are mounted on a precise motor-driven stage inside the +ablation chamber of the IMC instrument. A high-energy UV laser is focused on the +tissue, and each individual laser shot ablates tissue from an area of roughly 1 +μm\(^2\). The energy of the laser is absorbed by the tissue resulting +in vaporization followed by condensation of the ablated material. The ablated +material from each laser shot is transported in the gas phase into the plasma of +the mass cytometer, where first atomization of the particles and then ionization +of the atoms occurs. The ion cloud is then transferred into a vacuum, and all +ions below a mass of 80 m/z are filtered using a quadrupole mass filter. The +remaining ions (mostly those used to tag antibodies) are analyzed in a +time-of-flight mass spectrometer to ultimately obtain an accumulated mass +spectrum from all ions that correspond to a single laser shot. One can regard +this spectrum as the information underlying a 1 μm\(^2\) pixel. With +repetitive laser shots (e.g., at 200 Hz) and a simultaneous lateral sample +movement, a tissue can be ablated pixel by pixel. Ultimately an image is +reconstructed from each pixel mass spectrum.
+In principle, IMC can be applied to the same type of samples as conventional +fluorescence microscopy. The largest distinction from fluorescence microscopy is +that for IMC, primary-labeled antibodies are commonly used, whereas in +fluorescence microscopy secondary antibodies carrying fluorophores are widely +applied. Additionally, for IMC, samples are dried before acquisition and can be +stored for years. Formalin-fixed and paraffin-embedded (FFPE) samples are widely +used for IMC. The FFPE blocks are cut to 2-5 μm thick sections and are +stained, dried, and analyzed with IMC.
+Metal-labeled antibodies are used to stain molecules in tissues enabling to +delineate tissue structures, cells, and subcellular structures. Metal-conjugated +antibodies can either be purchased directly from Standard BioToolsTM (MaxPar IMC Antibodies), +or antibodies can be purchased and labeled individually (MaxPar Antibody +Labeling). +Antibody labeling using the MaxPar kits is performed via TCEP antibody reduction +followed by crosslinking with sulfhydryl-reactive maleimide-bearing metal +polymers. For each antibody it is essential to validate its functionality, +specificity and optimize its usage to provide optimal signal to noise. To +facilitate antibody handling, a database is highly useful. +Airlab is such a platform; it +allows antibody lot tracking, validation data uploads, and panel generation for +subsequent upload to the IMC acquisition software from Standard BioToolsTM
+Depending on the sample type, different staining protocols can be used. +Generally, once antibodies of choice have been conjugated to a metal tag, +titration experiments are performed to identify the optimal staining +concentration. For FFPE samples, different staining protocols have been +described, and different antibodies show variable staining with different +protocols. Protocols such as the one provided by Standard BioToolsTM or the one describe by +(Ijsselsteijn et al. 2019) are recommended. Briefly, for FFPE tissues, a dewaxing +step is performed to remove the paraffin used to embed the material, followed by +a graded re-hydration of the samples. Thereafter, heat-induced epitope retrieval +(HIER), a step aiming at the reversal of formalin-based fixation, is used to +unmask epitopes within tissues and make them accessible to antibodies. Epitope +unmasking is generally performed in either basic, EDTA-based buffers (pH 9.2) or +acidic, citrate-based buffers (pH 6). Next, a buffer containing bovine serum +albumin (BSA) is used to block non-specific binding. This buffer is also used to +dilute antibody stocks for the actual antibody staining. Staining time and +temperature may vary and optimization must be performed to ensure that each +single antibody performs well. However, overnight staining at 4°C or 3-5 +hours at room temperature seem to be suitable in many cases.
+Following antibody incubation, unbound antibodies are washed away and a +counterstain comparable to DAPI is applied to enable the identification of +nuclei. The Iridium intercalator +from Standard BioToolsTM is a reagent of choice and applied in a brief 5 minute staining. +Finally, the samples are washed again and then dried under an airflow. Once +dried, the samples are ready for analysis using IMC and are +usually stable for a long period of time (at least one year).
+Data is acquired using the CyTOF software from Standard BioToolsTM (see manuals +here).
+The regions of interest are selected by providing coordinates for ablation. To +determine the region to be imaged, so called “panoramas” can be generated. These +are stitched images of single fields of views of about 200 μm in diameter. +Panoramas provide an optical overview of the tissue with a resolution similar to +10x in microscopy and are intended to help with the selection of regions of +interest for ablation. The tissue should be centered on the glass side, since +the imaging mass cytometer cannot access roughly 5 mm from each of the slide +edges. Currently, the instruments can process one slide at a time and usually one MCD +file per sample slide is generated.
+Many regions of interest can be defined on a single slide and acquisition +parameters such as channels to acquire, acquisition speed (100 Hz or 200 Hz), +ablation energy, and other parameters are user-defined. It is recommended that +all isotope channels are recorded. This will result in larger raw data files but valuable information such as +potential contamination of the argon gas (e.g., Xenon) or of the samples (e.g., +lead, barium) is stored.
+To process a large number of slides or to select regions on whole-slide samples, +panoramas may not provide sufficient information. If this is the case, +multi-color immunofluorescence of the same slide prior to staining with +metal-labeled antibodies may be performed. To allow for region selection based +on immunofluorescence images and to align those images with a panorama of the +same or consecutive sections of the sample, we developed +napping.
+Acquisition time is directly proportional to the total size of ablation, and run +times for samples of large area or for large sample numbers can roughly be calculated by +dividing the ablation area in square micrometer by the ablation speed (e.g., +200Hz). In addition to the proprietary MCD file format, TXT files can also +be generated for each region of interest. This is recommended as a back-up +option in case of errors that may corrupt MCD files but not TXT files.
+Upon completion of the acquisition an MCD file of variable size is generated. A +single MCD file can hold raw acquisition data for multiple regions of interest, +optical images providing a slide level overview of the sample (“panoramas”), and +detailed metadata about the experiment. Additionally, for each acquisition a +TXT file is generated which holds the same pixel information as the matched +acquisition in the MCD file.
+The Hyperion Imaging SystemTM produces files in the following folder structure:
+.
++-- {XYZ}_ROI_001_1.txt
++-- {XYZ}_ROI_002_2.txt
++-- {XYZ}_ROI_003_3.txt
++-- {XYZ}.mcdHere, {XYZ} defines the filename, ROI_001, ROI_002, ROI_003 are
+user-defined names (descriptions) for the selected regions of interest (ROI),
+and 1, 2, 3 indicate the unique acquisition identifiers. The ROI
+description entry can be specified in the Standard BioTools software when
+selecting ROIs. The MCD file contains the raw imaging data and the full metadata
+of all acquired ROIs, while each TXT file contains data of a single ROI without
+metadata. To follow a consistent naming scheme and to bundle all metadata, we
+recommend to zip the folder. Each ZIP file should only contain data from a
+single MCD file, and the name of the ZIP file should match the name of the MCD
+file.
We refer to this data as raw data and the further +processing of this data is described in Section 3.
+ +