---
title: "Imaging-based Spatial Transcriptomics Data Quality Control with SpaceTrooper"
author: 
    - name: "Benedetta Banzi"
    - name: "Dario Righelli"
date: "`r BiocStyle::doc_date()`"
output:
    BiocStyle::html_document:
        toc: true
        toc_float: true
vignette: >
    %\VignetteIndexEntry{Imaging-based Spatial Transcriptomics Data Quality Control with SpaceTrooper}
    %\VignetteEncoding{UTF-8}
    %\VignetteEngine{knitr::rmarkdown}
editor_options:
    chunk_output_type: console
---

```{r setup, include=FALSE}
# Set chunk options: suppress echo, messages, and warnings in code output
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```
## Introduction

`SpaceTrooper` is an `R/Bioconductor` package for Quality Control (QC)
of imaging-based spatial transcriptomics and proteomics data. It provides 
multi-platform data harmonization, cell-level QC, and visualization utilities. 
The package leverages [SpatialExperiment](https://www.bioconductor.org/packages/release/bioc/html/SpatialExperiment.html) objects to support data from 
**CosMx**, **Xenium**, and **MERFISH** technologies.


## Quality Control pipeline

In this section, we show how to run the standard `SpaceTrooper` workflow on 
Spatial Transcriptomics data, to perform QC by computing the Quality Score (QS). 

For demonstration purposes, we showcase the `CosMx_DBKero_Tiny` dataset as the 
main tutorial example, a small subset of the CosMx 1k breast cancer dataset 
analyzed also in the [paper](https://www.biorxiv.org/content/10.64898/2025.12.24.696336v1). 

### Load example data

`SpaceTrooper` provides three different reading functions, one per each 
technology, that load data into a `SpatialExperiment` object. For more details 
on how to correctly specify the inputs, please refer to the [SpaceTrooper utilities](https://bioconductor.org/packages/devel/bioc/vignettes/SpaceTrooper/inst/doc/SpaceTrooper_utilities.html) vignette.

```{r load-cosmx}
library(SpaceTrooper)
library(ggplot2)

cosmxFolder <- system.file(file.path("extdata", "CosMx_DBKero_Tiny"), 
                        package="SpaceTrooper")
spe <- readCosmxSPE(cosmxFolder, sampleName="DBKero_CosMx")

spe
```

In addition, we provide the code for loading **Xenium**, and **MERFISH** example
data. Please, note that for these technologies the parameters are configured to 
automatically load polygons and compute missing metrics, so these steps don't 
need to be run again. To proceed with these technologies, skip directly to
[Add QC metrics](#add-qc-metrics).

```{r load-xenmer}
xeniumFolder <- system.file( "extdata", "Xenium_small", package="SpaceTrooper")
xen <- readXeniumSPE(xeniumFolder, computeMissingMetrics=TRUE,
                        keepPolygons=TRUE)

merfishFolder <- system.file("extdata", "Merfish_Tiny", package="SpaceTrooper")
mer <- readMerfishSPE(merfishFolder, boundariesType="parquet",
                        computeMissingMetrics=TRUE, keepPolygons=TRUE)
```

### Field of View (FOVs) visualization

In this section, we show how to plot each cell position onto the map of the 
Fields of View (FOVs), whose coordinates are provided only for CosMx technology.

**IMPORTANT**: we have noticed that according to CosMx version, there could be a 
misalignment of FOVs and cell centroids, that can be easily corrected with a 
single line of code. Therefore, it is crucial to generate this plot to check for 
any spatial shift, as such misalignment directly affects the QS computation.

```{r plot-fovs}
# to check misalignment
plotCellsFovs(spe, size = 3, alpha = 0.7)
```

The dataset is a subset with just a single FoV, whose number is displayed at the
center. Cell centroids, shown in dark red, are all contained by FoV boundaries,
thus no shift correction is needed here.

When an experiment has multiple FoVs, you can see the map and the topological
organization of the FoVs, together with their numbers.

If any misalignment is observed, it can be corrected by adjusting the FoV 
coordinates directly. For instance, if the FoVs are shifted upward by one FoV 
height (which, for CosMx technology, corresponds to 4,256 pixels), you can 
correct this by subtracting that value from the FoVs y coordinates:

```{r fov-correction}
# code line for shift correction
# metadata(spe)$fov_positions$y_global_px <- metadata(spe)$fov_positions$y_global_px - 4256
```

### Load polygons

In this section we show how to load polygons after `SpatialExperiment` creation, 
only for visualization purposes. They are stored as an `sf` object within 
`colData`.

This step is not mandatory for CosMx, because the pipeline can be executed even
without them. 

```{r load-poly,  message = TRUE}
# polygon loading
spe <- readAndAddPolygonsToSPE(spe, boundariesType="csv")
```

Please pay attention to any warnings. Cell polygons with fewer than four 
vertices cannot be handled by the geometry packages used by `SpaceTrooper`. 
Therefore, the corresponding cells are discarded from the `SpatialExperiment` 
object.

### Add QC metrics

The following sections will work the same way for all types of data.

The `spatialPerCellQC` function computes additional metrics per each cell, that
are saved inside the `SpatialExperiment`and accessible through `colData(spe)`. 
It is mandatory to run this function before computing QS.

The `negProbList` parameter is, by default, a vector containing all control probe 
patterns encountered so far across the supported technologies. Because these 
patterns continue to evolve, some may not yet be included. If you find that your
control probe patterns are missing from the default list, you can define a 
custom vector, as shown below.

By default, the function automatically removes 0 count cells, but this can be 
handled with the `rmZeros` parameter.

```{r cosmx-analysis-qc,  message = TRUE}
spe <- spatialPerCellQC(spe, rmZeros=TRUE,
        negProbList=c("NegPrb", "Negative", "SystemControl"))

colnames(colData(spe))
```

Several metrics are added, both derived from expression assay and cell 
morphology. Some of them are directly used to compute QS:

- **log2SignalDensity**: log2-transformed ratio between total probe counts per 
    cell and cell area in µm² (or volume in µm³ for MERFISH technology);
- **Area_um**: cell area in µm² (or volume in µm³ in the case of MERFISH 
    technology);
- **log2Ctrl_total_ratio**: log2-transformed ratio between total negative 
    control probe counts and total probe count per cell;
- **log2AspectRatio**: log2-transformed ratio between cell maximum length along 
    the x dimension and cell maximum length along the y dimension (pixels). This
    metric is taken in absolute value and combined with **dist_border** to 
    consider only cells within 50 pixels from the nearest FoV border (only for 
    CosMx technology).

For a better detailed explanation of the other metrics, please refer to 
the [SpaceTrooper utilities](https://bioconductor.org/packages/devel/bioc/vignettes/SpaceTrooper/inst/doc/SpaceTrooper_utilities.html) vignette.

### Compute Quality Score

`computeQCScore` function calculates QS per each cell. QS combines several
metrics into a formula: **log2SignalDensity**, corresponding to signal density, 
**Area_um** or **volume**, i.e. cell size, **log2Ctrl_total_ratio**, namely 
background signal and **log2AspectRatio** combined with **dist_border**, which 
jointly correspond to the border effect (only for CosMx technology). 

`glmnet` package is used to estimate the coefficients of the formula, resulting 
in a robust score that captures low-quality cells. During model training, the
selected cells are used as good or bad examples to learn how each term in the 
formula contributes to cell quality.

**IMPORTANT**: please, pay attention to any warning. If a term has too few bad 
examples (fewer than 0.1% of the total number of cells), it is excluded from the 
formula and therefore not used in the QS computation.

The QS (stored as `QC_score` in `coldata`) ranges from 0 to 1, with 0 meaning 
low-quality and 1 high-quality.

We are setting a seed to ensure reproducibility in this tutorial, because there 
are stochastic processes underlying QS computation.

```{r cosmx-analysis-score,  message = TRUE}
set.seed(1713)

spe <- computeQCScore(spe)

summary(spe$QC_score)
```

In this case, all the terms were used as no warning appeared.

Subsequently, it is possible also to assess which cells have a QS lower than a 
certain threshold (default is 0.5) with the following function. It creates a
new column called `low_qcscore` inside of `colData`.

```{r cosmx-analysis-score2, message = TRUE}
spe <- computeQCScoreFlags(spe, qsThreshold=0.5)

table(spe$low_qcscore)
```

Using this threshold, 97 cells are flagged as low-quality. We do not suggest a 
fixed default threshold, but it is advisable to check the QS distribution before
setting any threshold.

### Data visualization

SpaceTrooper comes with several functions to view cells and metrics.

To view the distribution of whatever quantitative metric, `plotMetricHist` comes 
in handy.

```{r plot-hist}
# view quantitative metric distribution
plotMetricHist(spe, metric = "QC_score")
```

The QS distribution exhibits a left tail starting around 0.75. 

Cell visualization can be obtained by using either centroids (recommended when 
the dataset has a large number of cells) or polygons. 

`plotCentroids` plots cell centroids that can be colored by a certain metric
contained in the `colData` slot, by using the `colour_by` parameter.

Additionally, if you have a palette column in `colData`, containing colors for
each cell, it can be given to `palette` parameter, so that it automatically 
matches the column passed in `colour_by`. As an example, we are using the 
cell types obtained as described in the [paper](https://www.biorxiv.org/content/10.64898/2025.12.24.696336v1) with their own color palette.

```{r plot-centroids-labels}
labf <- system.file(file.path("extdata", "CosMx_DBKero_Tiny",
                            "labels_tiny.tsv"), package="SpaceTrooper")
labs <- read.table(file=labf, sep="\t", header=TRUE)

spe$labels <- as.factor(labs[match(spe$cell_id, labs$cell_id),]$label)
spe$labels_colors <- as.factor(labs[match(spe$cell_id, labs$cell_id),]$lab_color)

plotCentroids(spe, colourBy="labels", size=3, palette="labels_colors")
```

Cell centroids colored by cell types allow to view the spatial distribution of
cell populations in the sample.

When possible, polygon visualization gives a better overview of the 
cells' spatial organization and morphological characteristics.
Polygons are stored as an `sf` object within `colData`, so they can be
viewed using standard `ggplot2` functions. For CosMx technology, the polygons
must be explicitly loaded as described in [Load polygons](#Load polygons).
SpaceTrooper provides `plotPolygons` function that works just like 
`plotCentroids` but takes polygons instead of centroids.

```{r plot-polygons-fov-1}
plotPolygons(spe, colourBy="log2SignalDensity")
plotPolygons(spe, colourBy="Area_um")
plotPolygons(spe, colourBy="log2Ctrl_total_ratio")
plotPolygons(spe, colourBy="log2AspectRatio")
```

We can see in `yellow` and `darkviolet` that there are few cells with extreme 
values of either `log2SignalDensity`, `Area_um`, `log2Ctrl_total_ratio` and 
`log2AspectRatio`.

Since all plotting functions are based on `ggplot2`, you can easily customize 
the graphical outputs by adding standard ggplot2 components.

```{r plot-polygons-fov-2}
plotPolygons(spe, colourBy="QC_score") + 
    scale_fill_viridis_c(option = "plasma")
```

We can see that the QS is able to detect both the aspects highlighted by 
`log2SignalDensity`,`Area_um`, `log2Ctrl_total_ratio` or 
`log2AspectRatio`. Cells that showed either lower signal 
density, bigger size, higher background signal or border effect also display 
low QS (darker color).

It's up to the user to choose an appropriate threshold to flag cells according 
to the observed QS distribution. 

```{r plot-polygons-fov-3}
plotPolygons(spe, colourBy="low_qcscore") + 
    scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf"))
```

You can rerun`computeQCScoreFlags` to check how many cells would be flagged 
using another threshold.

```{r cosmx-analysis-score3, message = TRUE}
spe <- computeQCScoreFlags(spe, qsThreshold=0.75)

table(spe$low_qcscore)
```

Using this threshold, 193 cells are flagged as low-quality.

```{r plot-polygons-fov-4}
plotPolygons(spe, colourBy="low_qcscore") + 
    scale_fill_manual(values=c("TRUE"="red", "FALSE" = "#c0c8cf"))
```

The threshold is more stringent and more cells are flagged compared to the
previous one. It's up to you to choose a threshold that better suits your
analysis needs.

## Conclusion

In this vignette, we explored the main functionalities of the `SpaceTrooper`
package for imaging-based spatial omics data QC. 

For further insights on SpaceTrooper package usage, please refer to the [SpaceTrooper utilities](https://bioconductor.org/packages/devel/bioc/vignettes/SpaceTrooper/inst/doc/SpaceTrooper_utilities.html) vignette.

Main steps shown are:

- data loading: CosMx RNA, Xenium, MERFISH
- Quality Control: 
    - per-cell QC metrics
    - Quality Score: a score combining **signal density**, **cell size**, 
    **background noise** and  **border effect**
- visualization: 
    - centroids: with ggplot2
    - polygons: sf + ggplot2

## Session Information

```{r sessionInfo}
sessionInfo()
```
