Workflow: CosMX data analysis
This workflow is part of a series of vignettes illustrating analysis workflows for spatial transcriptomics data using the OSTA ecosystem. OSTA is a textbook aimed at providing practical guidance for the analysis of spatial transcriptomics data using Bioconductor packages. The book is available at this link, a preprint of the book is available at this link, and the code repository is at this link.
This workflow extends the existing workflow present in OSTA book, with additional details specifically suited for the BBCC25 conference.
Instructors:
- Dario Righelli1
- Francesca Calanca2,3
- Emanuela Scafuri2,3
1 Department of Electrical Engineer and Information Technology,
University of Naples “Federico II”, Italy
2 Department of Chemical, Materials and Production Engineering,
University of Naples “Federico II”, Italy
3 Institute for Applied Mathematics “M. Picone”, IAC-CNR, Naples,
Italy
Introduction
In this demo, we will analyze a 1k-plex CosMx data (Bruker) of a mouse brain sample, see other datasets at this link for additional examples and exercising yourself. Main shown steps are quality control and preprocessing, cell type annotation, identification of marker genes, and cell-cell interaction analysis.
Dependencies
We now load all the required libraries.
library(DESpace)
library(ggplot2)
library(OSTA.data)
library(patchwork)
library(pheatmap)
library(scater)
library(scrapper)
library(scRNAseq)
library(sf)
library(SingleR)
library(SpaceTrooper)
library(SpatialExperiment)
library(SpatialExperimentIO)
library(spatialFDA)
library(SpatialFeatureExperiment)
library(spdep)
library(Voyager)
# set seed for random number generation
# in order to make results reproducible
set.seed(112358)Load data
CosMx data
We load the dataset which was pre-stored in the
OSTA.data package. The dataset contains about 50,000 cells
profiled for about 1,000 RNA targets. We then use the
readCosmxSXE function from the
SpatialExperimentIO package to read the data into a
SpatialExperiment object Righelli et
al. (2022).
Then, it is possible to adapt the dataset to work with
SpaceTrooper by updating the object and reading the cell
polygons (this process is optional if we don’t need the polygons).
# retrieve CosMx dataset from OSF repo
id <- "CosMx1k_MouseBrain2"
pa <- OSTA.data_load(id, mol=FALSE)
dir.create(td <- tempfile())
unzip(pa, exdir=td)
cos <- readCosmxSXE(td, addTx=FALSE)
# prepare data for 'SpaceTrooper'
cos <- updateCosmxSPE(cos, td, sampleName="CosMx")
cos <- readAndAddPolygonsToSPE(cos)
cos$in_tissue <- TRUE # artifact to work with ggspavis
cos## class: SpatialExperiment
## dim: 950 48556
## metadata(4): fov_positions fov_dim polygons technology
## assays(1): counts
## rownames(950): Chrna4 Slc6a1 ... Cck Aqp4
## rowData names(0):
## colnames(48556): f1_c1 f1_c10 ... f99_c98 f99_c99
## colData names(22): fov cellID ... polygons in_tissue
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(1): NegPrb
## spatialCoords names(2) : CenterX_global_px CenterY_global_px
## imgData names(1): sample_idIn alternative, it is possible to directly load imaging-based
technologies dataset using the SpaceTrooper package with one of the
following commands. Note that for adding the poligons to the
SpatialExperiment object, you can use the
readAndAddPolygonsToSPE function as shown above, or just
use the keepPolygons=TRUE argument (for xenium and merfish) in the
reading functions.
NB: these are only for example purpose. If you want to download xenium or merfish data, you need to find a dataset of one of these technologies (i.e. in OSTA.data) and then put the path of the dataset to the functions.
cos <- readCosmxSPE(td, sampleName = "CosMx")
cospro <- readCosmxProteinSPE(td, sampleName = "CosMxProtein")
xen <- readXeniumSPE(dirName = "path/to/xenium/data", sampleName = "Xenium", keepPolygons = TRUE) # do not run as it is
mer <- readMerfishSPE(dirName = "path/to/merscope/data", sampleName = "MERSCOPE", keepPolygons = TRUE) # do not run as it isOnce loaded, you can access additional metadata of the
SpatialExperiment object cos as follows: - colData(cos):
cell-level metadata, including FOV ID, cell area, and QC metrics. -
alternatively, you can directly access a specific metadata column with
cos$<column_name>. (e.g., cos$fov gives
the FOV IDs for each cell). - rowData(cos): gene-level metadata,
including gene names and biotypes. - assays(cos): names of the assays
available in the object (e.g., “counts” for raw counts, “logcounts” for
log-normalized counts if available). - spatialCoords(cos): spatial
coordinates of the cells.
Basic Visualization
We can have an overview of the data, by plotting the spatial distribution of the cells, represented by their centroids. This is useful to have a global bird’s-eye-view of the tissue.
plotCentroids(cos)
Unlike, Xenium and MERSCOPE, CosMx does not stitch fields of view (FOVs) in a single image, but rather treats each FOV independently. It is therefore useful to visualize a map of the FOVs to understand their spatial arrangement.
plotCellsFovs(cos)
Excercises:
- Explore the
cosobject usingcolData(cos),rowData(cos), andassays(cos). - How many cells and genes are in the dataset? (hint: use
dim(cos)) - Explore the cell metadata. What information is available for each cell? Do you recognize something specific to CosMx data?
- Check the distribution of counts across cells. Are there any cells with very low counts?
- Extract the spatial coordinates of the cells using
spatialCoords(cos). What is the range of the x and y coordinates? - plot the centroids of the cell using
plotCentroids(cos). Colouring them for a specific cell metadata. -
BONUS: Load a different Xenium dataset from the
OSTA.datapackage for further analysis (use theJanesick_breastCancer_Xenium_rep1). This could be heavy to load for your laptop if you think to do it together with the cosmx dataset.
Quality Control
Quality control (QC) is a crucial step in imaging-based spatial transcriptomics, just as it is in scRNA-seq, where it helps remove low-quality cells and reveal technical or biological biases at the sample, field-of-view (FOV), or slide-area level.
A key difference, however, is that imaging-based ST platforms typically do not quantify the whole transcriptome—or even an unbiased subset of it—but rely instead on predefined gene panels. These panels, depending on platform and experimental goals, may be pan-tissue or focused on specific tissues or cell types. Their size varies widely, from a few hundred genes (e.g., 313 genes in the first Xenium panel) to several thousand (e.g., the Xenium 5k panel or the CosMx 6k and 18k panels). As a result, unlike sequencing-based assays, imaging platforms may not provide enough ribosomal or mitochondrial transcripts to compute standard QC metrics. They do, however, include sets of negative probes that allow detection of non-specific RNA capture.
Finally, imaging-based ST technologies also offer rich morphological information (e.g., cell area, eccentricity) and signals from antibody stains—such as nuclear or membrane markers—which can further support QC and downstream analyses.
Here we compute the cell-level QC metrics using the
SpaceTrooper package.
old <- names(colData(cos))
lys <- c("NegPrb", "Negative", "SystemControl")
cos <- spatialPerCellQC(cos, rmZeros=TRUE, negProbList=lys)This produces the following new QC metrics in the object’s
colData:
## [1] "sum" "detected"
## [3] "altexps_NegPrb_sum" "altexps_NegPrb_detected"
## [5] "altexps_NegPrb_percent" "total"
## [7] "control_sum" "control_detected"
## [9] "target_sum" "target_detected"
## [11] "CenterX_global_px" "CenterY_global_px"
## [13] "ctrl_total_ratio" "log2Ctrl_total_ratio"
## [15] "CenterX_global_um" "CenterY_global_um"
## [17] "Area_um" "dist_border_x"
## [19] "dist_border_y" "dist_border"
## [21] "log2AspectRatio" "CountArea"
## [23] "log2CountArea"If we want to explore the distribution of a specific QC metric, e.g.,
the log2 of counts per area, we can use the plotMetricHist
function. When computing the spatial outliers, we can use either the
scuttle (sc) or medcaple (mc) methods, or both. Then we can visualize
the two outlier thresholds on the histogram.
cos <- computeSpatialOutlier(cos, computeBy = "log2CountArea", method="both")
plotMetricHist(cos, "log2CountArea", useFences = "log2CountArea_outlier_sc")
plotMetricHist(cos, "log2CountArea", useFences = "log2CountArea_outlier_mc")
Or visualize the metric in the spatial context:
plotCentroids(cos, colourBy="log2CountArea")## Coordinate system already present.
## ℹ Adding new coordinate system, which will replace the existing one.
Now, we will simply compute SpaceTrooper’s aggregated QC score and keep the cells that pass the threshold. More details will be provided once the paper will be published (Banzi , Righelli et al., in preparation).
cos <- computeQCScore(cos)## Metric log2Ctrl_total_ratio produced NA skewness values.It will be removed from the QC score formula.
plotCentroids(cos, colourBy="QC_score")## Coordinate system already present.
## ℹ Adding new coordinate system, which will replace the existing one.
After that we can compute the QC flags, visualize them and filter out low-quality cells.
cos <- computeQCScoreFlags(cos, qsThreshold=0.5)
plotCentroids(cos, colourBy="low_qcscore") +
theme(legend.key.size=rel(0)) +
guides(col=guide_legend(override.aes=list(size=2)))
SpaceTrooper::plotMetricHist(cos, "QC_score")
cos <- cos[, !cos$low_qcscore]Excercises:
- Explore the newly added QC metrics in
colData(cos). What do they represent? - Visualize the distribution of another QC metric, e.g., `log2NegativeArea
- Are you able to identify outlier cells based on this metric? how you visualize them?
- Visualize the spatial distribution of another QC metric.
- Compute the outliers for the QC_score (before filtering low qc cells) using both methods (sc and mc). How many cells are flagged as low quality by each method?
Processing
For the sake of runtime and laptop phisical limitations, we will perform downstream analyses only on a small portion of the tissue, defined by the following FOVs:
fs <- c(72:74, 90:92, 97:99, 114:116)
plotZoomFovsMap(cos, fovs=fs)
sub <- cos[, cos$fov %in% fs]Next, we’ll log-normalize counts by area, and perform principal component analysis (PCA) on all 950 RNA targets:
# cell area-based normalization
logcounts(sub) <- sweep(assay(sub), 2, sub$Area, `/`)
# principal component analysis
sub <- runPCA(sub)Let’s visualize the expression of the first two principal components (PCs) in the spatial context.
# add PCs as cell metadata
pcs <- reducedDim(sub, "PCA")
colData(sub) <- cbind(colData(sub), pcs)
# visualize PCs 1 & 2 in space
plotCentroids(sub, colourBy="PC1") +
plotCentroids(sub, colourBy="PC2")
Cell type Annotation
Cell type annotation is a key step in the analysis of spatial transcriptomics data, as it allows us to assign biological meaning to the observed spatial patterns of gene expression.
Here, we use the Zeisel et al. (2015) dataset as a reference to annotate the cell types with SingleR. This dataset is composed by single-cell RNA-seq data from the mouse brain, and contains 3000 cells from a variety of annotated cell types at two different levels of granularity (level 1 and level 2). We will use the coarse level 1 annotations for this demo.
sceZ <- ZeiselBrainData()
sceZ <- logNormCounts(sceZ)
table(sceZ$level1class)##
## astrocytes_ependymal endothelial-mural interneurons
## 224 235 290
## microglia oligodendrocytes pyramidal CA1
## 98 820 939
## pyramidal SS
## 399
table(sceZ$level2class)##
## (none) Astro1 Astro2 CA1Pyr1 CA1Pyr2 CA1PyrInt CA2Pyr2 Choroid
## 189 68 61 380 447 49 41 10
## ClauPyr Epend Int1 Int10 Int11 Int12 Int13 Int14
## 5 20 12 21 10 21 15 22
## Int15 Int16 Int2 Int3 Int4 Int5 Int6 Int7
## 18 20 24 10 15 20 22 23
## Int8 Int9 Mgl1 Mgl2 Oligo1 Oligo2 Oligo3 Oligo4
## 26 11 17 16 45 98 87 106
## Oligo5 Oligo6 Peric Pvm1 Pvm2 S1PyrDL S1PyrL23 S1PyrL4
## 125 359 21 32 33 81 74 26
## S1PyrL5 S1PyrL5a S1PyrL6 S1PyrL6b SubPyr Vend1 Vend2 Vsmc
## 16 28 39 21 22 32 105 62
pred <- SingleR(test=sub, ref=sceZ, labels=sceZ$level1class, de.method="wilcox")
sub$SingleR_label <- factor(pred$pruned.labels)
nk <- nlevels(sub$SingleR_label)
pal <- hcl.colors(nk, "Spectral")
plotPolygons(sub,
colourBy="SingleR_label") +
scale_fill_manual(values=pal)
Note that the annotation performed here is a transfer label from a reference dataset that only partially matches with the tissue under study and purposely annotated at a coarse level (level 1). Ideally, one can use a more appropriate reference dataset, or better yet, build a custom reference from single-cell RNA-seq data of the same tissue.
A different approach consists in using, in addition to the gene expression profiles, the image data and the spatial context of the cells. InSituType (Danaher et al. 2022) provides a framework to perform such analysis in unsupervised, semi-supervised, or supervised manners. Please refer to its documentation for more details.
Marker genes
Once that we have identified the cell types, we can identify cluster-related spatially-variable genes using the DESpace package (Cai et al. 2024). In this workflow, we use this as a way to identify marker genes for the cell types annotated above.
## gene_id LR logCPM PValue FDR
## Itm2a Itm2a 2596.761 11.03859 0 0
## Flt1 Flt1 2601.543 11.05386 0 0
## Cldn5 Cldn5 2825.269 11.11885 0 0
## Gja1 Gja1 4982.193 11.29054 0 0
## Gpr37l1 Gpr37l1 4989.928 11.30344 0 0
## Bcan Bcan 3567.078 11.61163 0 0We focus on the genes with the smallest p-values, and compute the average expression of each gene in each cell type. We can visualize this with a heatmap.
df <- res$gene_results
gs <- df$gene_id[rank(df$PValue, ties.method="min") == 1]
mu <- t(apply(counts(sub)[gs, ], 1, tapply, sub$SingleR_label, mean))
pheatmap(t(mu),
scale="column", show_colnames=FALSE,
cluster_rows=TRUE, cluster_cols=TRUE,
main="Average expression of top SVGs per cell type")
We can also visualize the spatial distribution of some of these genes, e.g., Mbp that marks oligodendrocytes, and Calm1, Calm2, and Snap25 that mark different subsets of neurons. Note also how some genes show a difference of expression within the pyramidal neuron population, suggesting the presence of subtypes.
gs <- c("Mbp", "Calm1", "Calm2", "Snap25")
ps <- lapply(gs, \(g) {
sub[[g]] <- logcounts(sub)[g, ]
plotPolygons(sub, colourBy=g) + ggtitle(g)
})
dy <- range(logcounts(sub)[gs, ])
wrap_plots(ps, guides="collect") &
scale_fill_viridis_c(
"expression",
limits=dy, breaks=dy,
labels=c("low", "high")) &
theme(plot.title=element_text(hjust=0.5))
Neighborhood-based gene expression analyses
We start by creating a cell neighborhood graph, based on a distance threshold of 50 microns. To do so, we use the SpatialFeatureExperiment package (Moses et al. 2023).
Note that an alternative approach uses the nearest neighbors (kNN) method to define the neighborhood graph. Here, we prefer to use a distance-based approach, as it better captures the morphological structure of this tissue area.
# convert to SFE & remove cells with NA labels
sfe <- as(sub, "SpatialFeatureExperiment")
sfe <- sfe[, !is.na(sfe$SingleR_label)]
# this is needed for Voyager's plotting functions
# sfe$polygons<-.renameGeometry(sfe$polygons, from="global", to="geomery")
colnames(sfe$polygons)[6] <- "geometry"
st_geometry(sfe$polygons) <- "geometry"
colGeometry(sfe, "cellSeg") <- sfe$polygons
colGraph(sfe, "poly2nb") <-
findSpatialNeighbors(sfe,
type="cellSeg",
method="poly2nb",
style="W", snap=50)
plotColGraph(sfe,
colGraphName="poly2nb",
colGeometryName="cellSeg") +
theme_void()
Given the neighborhood graph, we can calculate local measures of spatial autocorrelation for specific genes in the dataset.
Moran’s I is a measure of spatial autocorrelation that quantifies the degree to which observed values (e.g., gene expression) are similar across neighboring spatial locations. The core idea is to compare the value of each cell with the average value of its spatial neighbors, using a spatial-weights matrix that defines what “neighboring” means.
When Moran’s I is:
positive → similar values cluster together (High–High or Low–Low regions)
negative → dissimilar values tend to be adjacent Low-High (checkerboard-likeLow patterns)
≈ 0 → no spatial structure (random pattern)
Here, we focus on Snap25 and Calm1, neuronal marker genes identified above that show interesting spatial expression distributions (see above). First, we compute the local Moran’s I coefficient.
gs <- c("Snap25", "Calm1")
sfe <- runUnivariate(sfe,
features=gs,
type="localmoran",
zero.policy=TRUE,
colGraphName="poly2nb",
colGeometryName="cellSeg")
plotLocalResult(sfe,
name="localmoran",
features=gs, ncol=2,
colGeometryName="cellSeg",
divergent=TRUE, diverge_center=0)
From the plots, we can see that Snap25 shows positive spatial autocorrelation in a very localized part of the tissue, while Calm1 shows a more diffuse autocorrelation structure.
We focus on Snap25 and visualize its autocorrelation structure using a Moran scatter plot, which shows the relationship between the expression of the gene in each cell and the average expression in its neighbors.
sfe <- runUnivariate(sfe,
feature=gs[1],
zero.policy=TRUE,
type="moran.plot",
colGraphName="poly2nb",
colGeometryName="cellSeg")
moranPlot(sfe,
feature=gs[1],
graphName="poly2nb",
swap_rownames="symbol")
We can then classify cells based on the significance of their local Moran’s I values, and visualize the resulting clusters in the spatial context.
res <- localResults(sfe)[["localmoran"]][[1]]
res$locClust <- with(res, ifelse(
`-log10p_adj` > -log10(0.05),
as.character(mean), "non-siginificant"))
localResults(sfe)[["localmoran"]][[1]] <- res
plotLocalResult(sfe,
features=gs[1],
name="localmoran",
attribute="locClust",
colGeometryName="cellSeg")
In this case, we can see a significant positive autocorrelation in the lower left part of the structure, while the rest of the tissue does not show significant autocorrelation.
Note that other local spatial autocorrelation measures, such as Geary’s C coefficient and Getis-Ord statistic, as well as their global versions can be used for similar analyses Moses et al. (2023).
Cell-cell interactions
Joint counts
We can also use the neighborhood graph to investigate cell-cell interactions between different cell types annotated above. Here, we use the join count statistic implemented in the spdep package to quantify the tendency of cells of different types to be neighbors.
In practice, join counts evaluate, for each pair of labels, how frequently two adjacent cells belong to those categories, and compare this to an expectation derived from a permutation-based null model. The resulting z-scores indicate the strength and direction of the interaction:
Positive z-values → the two cell types co-locate more than expected (attraction / spatial clustering).
Negative z-values → the two cell types co-locate less than expected (avoidance / spatial segregation).
Values near zero → no significant deviation from spatial randomness.
jc <- joincount.multi(sfe$SingleR_label, colGraph(sfe, "poly2nb"), zero.policy=TRUE)
head(jc[order(abs(jc[, "z-value"]), decreasing=TRUE), ])## Joincount Expected Variance z-value
## Jtot 1132.15707 2494.52543 265.87369 -83.55200
## pyramidal CA1:pyramidal CA1 735.66944 351.21205 62.92504 48.46592
## pyramidal CA1:oligodendrocytes 22.39901 303.28546 75.62290 -32.30012
## oligodendrocytes:oligodendrocytes 196.71667 65.35907 19.34973 29.86192
## pyramidal CA1:astrocytes_ependymal 97.55301 379.29404 92.67327 -29.26664
## pyramidal SS:pyramidal CA1 195.38757 445.19310 106.83889 -24.16781We can see that there is a strong tendency for cells of the same type to be neighbors (high positive z-values). There are also some interesting negative interactions between different cell types: e.g., pyramidal CA1 tend to be away from oligodendrocytes and astrocytes. This highlights the compartimentalized structure of the brain tissue.
Point-pattern analysis
We complement the graph-based interaction analysis with point-pattern statistics, which treat each cell as a spatial point and quantify deviations from complete spatial randomness (CSR). Using the spatialFDA package, we compute Besag’s L function, a variance-stabilized form of Ripley’s K function, to characterize both within-type and between-type spatial interactions.
Besag’s L(r) summarizes how many neighboring cells occur
within distance r of a reference cell, compared to the
expectation under CSR. When the empirical curve lies above the CSR line,
the point pattern exhibits clustering, whereas values below the line
indicate dispersion or spatial inhibition.
Keep in mind that Ripley’s K function, and the pair correlation function (Emons et al. 2025) makes the quite strong assuption of homogeneity of the point pattern, which is irrealistic in tissue biology. In particular, CSR is a naive null model that does not take into account tissue structure and cell density variations across the tissue. Hence, these plots should be considered as exploratory analyses rather than formal hypothesis tests.
res <- calcMetricPerFov(
spe=sub, fun="Lest",
subsetby="sample_id",
marks="SingleR_label",
selection="pyramidal CA1",
rSeq=seq(0, 500, l=100),
by="sample_id", ncores=1)
plotMetricPerFov(res,
theo=TRUE, correction="iso",
x="r", imageId="sample_id")
In this type of plots, the dashed line represents the theoretical value of the metric under complete spatial randomness (CSR). Values above this line indicate clustering, while values below indicate dispersion. Here, we can see that pyramidal CA1 neurons tend to cluster together, as expected.
Note that the L-function lays below the CSR line for the first 50 microns, which is a consequence of the cell segmentation that does not allow cells to overlap, indicating the impossibility of finding two cells at a distance smaller than their size.
While in the plot above we focused on a single cell type, we can also investigate the interaction between two different cell types, e.g., pyramidal CA1 neurons and oligodendrocytes, using the cross-L function.
res <- calcMetricPerFov(
spe=sub, fun="Lcross",
subsetby="sample_id",
marks="SingleR_label",
selection=c("pyramidal CA1", "oligodendrocytes"),
rSeq=seq(0, 500, l=100),
by="sample_id", ncores=1)
plotMetricPerFov(res,
theo=TRUE, correction="iso",
x="r", imageId="sample_id")
In this case, we can see that pyramidal CA1 neurons and oligodendrocytes tend to be spatially segregated, as indicated by the cross-L function being below the CSR line.
We only briefly discussed cell-cell interactions and spatial exploratory analyses here: we refer interested readers to the comprehensive documentation and tutorials available at the pasta and Voyager websites.
Excercises
- Preprocessing: what percentage of variance is explained by PC1 and PC2? Explore the explained percentage of variance with a plot (i.e. elbow plot, histogram).
- Annotation: create a frequency table of the annotation column. Which is the most abundant cell type in this tissue subset? Which is the rarest?
- Do you know any specific markers for microglia? Try to explore their expression. (i.e. extract the latest column from the heatmap)
- The neighborhood graph was built using a distance snap of 50 microns. Try re-running, changing the distance snap e.g., 80 and compare the two settings.