--- title: "Confidence-weighted spatial analysis with fuzzyCell" author: - name: Shreya Rao affiliation: The University of Sydney format: html: toc: true toc-depth: 3 number-sections: true embed-resources: true code-tools: true knitr: opts_chunk: collapse: true comment: "#>" message: false warning: false vignette: > %\VignetteIndexEntry{Confidence-weighted spatial analysis with fuzzyCell} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- # Overview Cell type annotation in spatial transcriptomics data is subject to considerable uncertainty. Transcript diffusion and imperfect segmentation assign transcripts to neighbouring cells, and limited gene panels yield sparse expression profiles, so closely related subtypes such as CD4+ and CD8+ T cells are frequently misclassified. Downstream analyses nonetheless typically treat annotations as fixed, and this uncertainty propagates, unaccounted for, into estimates of cell type proportions, co-localisation and niche composition. *fuzzyCell* quantifies this uncertainty by bootstrapping. Within each group of closely related subtypes, a classifier is trained on the original annotations, and each cell is repeatedly reclassified after multinomial resampling of its transcript counts. The proportion of resamples in which a cell is assigned to a subtype defines its **confidence** in that subtype. These confidences serve as **soft cell type weights**: each cell contributes to every subtype of its group in proportion to its confidence, rather than to a single label. Each analysis can be performed with the original annotations (*unweighted*) and with the confidence weights (*weighted*), and comparing the two identifies conclusions that are sensitive to annotation uncertainty. | If you want to know... | Use | Built on | |:-------------------------------------|:----------------|:----------------| | whether a cell type's proportion differs between outcomes | `runProportionsPipeline()`, `runProportionsRelativeToParentPipeline()` | linear models of image-level proportions | | whether two cell types co-localise more or less between outcomes | `runColocalisationModels()` | [spicyR](https://bioconductor.org/packages/spicyR)'s L-function, and Statial's Kontextual test for relative co-localisation | | whether tissue neighbourhoods (niches) differ between outcomes | `runNicheIdentification()`, `testNicheOutcome()` | [lisaClust](https://bioconductor.org/packages/lisaClust) | The rest of this vignette follows the workflow on a real dataset. # Installation ```{r} #| label: install #| eval: false if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("fuzzyCell") ``` # Example data: lung biopsies We follow the whole workflow on Xenium lung biopsies from patients with pulmonary fibrosis and from controls (Vannan et al., 2025). The full data are available as `SpatialDatasets::spe_Vannan_2025()`. The package includes a small subset, `lungSubset`: five control and five disease biopsies, each cropped to a window of about 3,000 cells at its centre. Every image comes from a different patient, and `disease_status` records whether the patient is a control or has pulmonary fibrosis. ```{r} #| label: data library(fuzzyCell) library(SpatialExperiment) library(ggplot2) data(lungSubset) spe = lungSubset spe # one image (and one outcome) per patient imageOutcome = unique(colData(spe)[, c("imageID", "disease_status")]) table(imageOutcome$disease_status) ``` ::: callout-note Any `SummarizedExperiment` (such as a `SpatialExperiment`) works as long as it has a `"counts"` assay and the cell metadata columns `cellID`, `imageID`, `cellType`, `x` and `y` in its `colData`. If your columns are named differently, the analysis functions take arguments (such as `cellType` and `imageID`) giving the column names. ::: # The cell type tree Confidence scores are estimated only between closely related cell types, so *fuzzyCell* requires a hierarchy specifying which subtypes may be confused with one another. The hierarchy is supplied as a nested list. A group of subtypes to be refined together is specified as `list(children = ...)`; every other entry is either an organisational level (such as the lineage `"Immune"`) or a leaf cell type that is not refined (such as `"SMCs/Pericytes"`). Each top-level entry must be named. The cell type tree for the lung data is defined as follows: ```{r} #| label: tree tree = list( "Immune" = list( "Myeloid" = list( "Macrophage lineage" = list(children = c("Alveolar Macrophages", "Interstitial Macrophages", "SPP1+ Macrophages", "Macrophages - IFN-activated", "Monocytes/MDMs")), "Dendritic cells" = list(children = c("cDCs", "Migratory DCs", "Langerhans cells", "pDCs")), "Granulocytes" = list(children = c("Mast", "Basophils", "Neutrophils")), "Proliferating Myeloid"), "Lymphoid" = list( "T cells" = list(children = c("CD4+ T-cells", "CD8+ T-cells", "Tregs", "Proliferating T-cells")), "NK lineage" = list(children = c("NK/NKT", "Proliferating NK/NKT")), "B lineage" = list(children = c("B cells", "Plasma", "Proliferating B cells")))), "Mesenchymal" = list( "Fibroblast lineage" = list(children = c("Alveolar FBs", "Adventitial FBs", "Activated Fibrotic FBs", "Inflammatory FBs", "Subpleural FBs", "Proliferating FBs", "Myofibroblasts")), "SMCs/Pericytes", "Mesothelial"), "Endothelial lineage" = list(children = c("Arteriole", "Venous", "Capillary", "Lymphatic")), "Epithelial lineage" = list( "Alveolar" = list(children = c("AT1", "AT2", "Transitional AT2", "Proliferating AT2")), "Airway" = list(children = c("Basal", "Secretory", "Goblet", "Multiciliated", "RASC", "PNEC", "KRT5-/KRT17+", "Proliferating Airway")))) buildCellTypeMap(tree)$groupsWithChildren ``` The same tree is included in the package and can be loaded with `data(lungCellTypeTree)`. Each of the groups listed above can receive confidence scores. Cell types outside a group, such as `SMCs/Pericytes` and `Mesothelial`, are not refined and retain their original labels. # Bootstrap confidence `bootstrapLineageConfidence()` does the work, group by group, as described in the Overview. Within each group it uses ridge multinomial logistic regression on the first 50 principal components of the group's cells (`classifier = "pcaLogistic"`, the default, and what we use here). `"logistic"` fits it on all genes. Other available classifiers include random forest (`"rf"`), linear discriminant analysis (`"lda"`) and support vector machines (`"svm"`). The classifier is trained with the same number of cells of each subtype, which gives every subtype the same prior. `priors = "observed"` instead weights each subtype by its share of the group's cells. We use 20 resamples to keep the vignette quick, but we encourage at least 100 for your own analysis. ```{r} #| label: bootstrap #| results: hide confList = bootstrapLineageConfidence(spe, tree, nBoot = 20, classifier = "pcaLogistic", seed = 1) ``` The result, `confList`, is a list with one data frame per group. Each data frame has one row per cell in the group, a `cellID` column, and one column per subtype holding that cell's confidence in it (between 0 and 1, adding up to 1 across the subtypes): ```{r} #| label: confhead names(confList) head(confList[["Endothelial lineage"]]) ``` How spread out are the confidences? Here are the scores for the subtypes of two groups, the B lineage and the endothelial cells. ```{r} #| label: confplot #| fig-height: 6 #| fig-width: 6 #| fig-alt: "Histograms of each subtype's confidence scores, for the B lineage and the endothelial cells." groups = c("B lineage", "Endothelial lineage") scores = do.call(rbind, lapply(groups, function(g) { conf = confList[[g]] do.call(rbind, lapply(setdiff(names(conf), "cellID"), function(subtype) data.frame(panel = paste0(g, ": ", subtype), confidence = conf[[subtype]]))) })) ggplot(scores, aes(confidence)) + geom_histogram(bins = 20, fill = "#4c78a8", colour = "white") + facet_wrap(~panel, ncol = 2, scales = "free_y") + labs(x = "confidence in the subtype", y = "number of cells") + theme_classic() ``` In the B lineage every cell scores close to 0 or 1 for each subtype, so the classifier is sure which one it is. The endothelial cells are different: many capillary cells, and some venous cells, have scores spread between 0 and 1, so they are split between subtypes. Few cells score close to 1 for arterioles, but many get a small arteriole score. Cells with scores in the middle count partly towards each subtype in the weighted analyses, which is where weighted and unweighted results can differ. # Cell type proportions Does the proportion of each cell type differ between controls and patients with pulmonary fibrosis? `runProportionsPipeline()` fits, for each cell type, a linear model of its per-image proportion on the outcome, once with the confidences as soft counts (weighted) and once with the original labels (unweighted). It takes as input the `SpatialExperiment` (`spe`), the cell type tree and the confidence list (`confList`). The outcome of interest is given by `outcomeCol`, the name of a column of `colData(spe)` with exactly two levels and a single value per image. `subject` is the name of the `colData(spe)` column identifying the patient each image comes from: every image must belong to a single subject, and every subject must contribute only one image. In the results, `estimate` is the difference in proportion for the second level of the outcome relative to the first (alphabetical order for text, factor order for factors); here that is `Disease` minus `Control`. `p.adj` is the Benjamini-Hochberg adjusted p-value across the cell types tested. ```{r} #| label: props #| results: hide props = runProportionsPipeline(spe, tree, confList, outcomeCol = "disease_status", subject = "subject") ``` ```{r} #| label: propstables head(props$modelWeighted) head(props$modelUnweighted) ``` For most cell types the weighted and unweighted analyses agree. Arterioles are one where they do not: with the original labels their proportion does not differ clearly between the outcomes, but with the confidences it is lower in disease. ```{r} #| label: propscompare cols = c("cellType", "estimate", "p.value", "p.adj") rbind( unweighted = subset(props$modelUnweighted, cellType == "Arteriole")[, cols], weighted = subset(props$modelWeighted, cellType == "Arteriole")[, cols]) ``` We can use `plotProportionBoxplot()` to visualise the per-patient proportions behind a cell type's estimate. ```{r} #| label: propsplot #| fig-height: 4 #| fig-width: 5 plotProportionBoxplot(props, "Arteriole", spe, "disease_status") ``` The plot is interactive: hovering over a point tells you which sample it is, so you can find the patients behind an estimate. As an example, let's look more closely at one image, the control biopsy `VUHD113`. `plotOriginalLabelSpatial()` shows the endothelial cells of the image with their hard labels, and `plotWeightedSpatial()` shows the same cells coloured by their confidence of being arterioles. ```{r} #| label: imgplots #| fig-height: 4 #| fig-width: 6 # the endothelial lineage with its original labels plotOriginalLabelSpatial(spe, tree, "VUHD113", "Endothelial lineage") # the same cells, weighted by their confidence of being arterioles plotWeightedSpatial(spe, tree, confList, "VUHD113", "Endothelial lineage", "Arteriole") ``` Marker expression helps to judge whether the confidences are sensible. Here is the arterial marker `HEY1` in the endothelial cells of this image: ```{r} #| label: imgmarkers #| fig-height: 4 #| fig-width: 6 plotOneMarkerSpatial(spe, "VUHD113", "Endothelial lineage", "HEY1", cellTypeTree = tree) ``` ```{r} #| label: hey1 #| echo: false endo = confList[["Endothelial lineage"]] endo = endo[endo$cellID %in% spe$cellID[spe$imageID == "VUHD113"], ] hey1 = counts(spe)["HEY1", endo$cellID] > 0 label = spe$cellType[match(endo$cellID, spe$cellID)] pct = function(x) paste0(round(100 * mean(x)), "%") ``` HEY1 expression is concentrated in the same vessels as the cells with high arteriole confidence, and is largely absent from the rest of the endothelium. In this image HEY1 is detected in `r pct(hey1[label == "Arteriole"])` of the cells labelled as arterioles and in `r pct(hey1[label != "Arteriole"])` of the other endothelial cells. It is detected in `r pct(hey1[endo$Arteriole >= 0.5])` of the cells with an arteriole confidence of at least 0.5, compared with `r pct(hey1[endo$Arteriole == 0])` of the cells with none, so the confidences are consistent with the marker. ## Proportions relative to the parent group The proportion of a subtype among all cells can change simply because its whole group changes. `runProportionsRelativeToParentPipeline()` instead models each subtype as a proportion of its group (for example, AT2 cells among all alveolar epithelial cells). It takes the same inputs as `runProportionsPipeline()` and returns weighted and unweighted models in the same format: ```{r} #| label: relprops #| eval: false relProps = runProportionsRelativeToParentPipeline(spe, tree, confList, outcomeCol = "disease_status", subject = "subject") ``` # Co-localisation `runColocalisationModels()` tests, for every pair of cell types and groups, whether their spatial co-localisation differs between the two outcomes. For an ordered pair `from` → `to`, co-localisation in each image is quantified by the cross-type L-function, computed by spicyR at a radius `r` (50 by default, in the units of the spatial coordinates). Positive values indicate that more `to` cells lie within `r` of `from` cells than expected if the cells were placed at random (attraction), negative values indicate fewer than expected (avoidance), and values near 0 indicate neither. For each pair, the per-image values are then modelled against the outcome with a linear model, once with the original annotations (unweighted) and once with the confidence weights (weighted). In the weighted analysis, each cell contributes to a label in proportion to its confidence. *Relative* co-localisation, which compares a pair with the parent population of the `to` cell type, is assessed with Statial's Kontextual test using the same weights. The confidence weights themselves are returned by `confidenceWeights()` as a matrix with one column per cell type and group, and can be passed directly to spicyR and lisaClust through their `cellWeights` argument. ```{r} #| label: coloc #| results: hide coloc = runColocalisationModels(spe, tree, confList, outcomeCol = "disease_status", subject = "subject", r = 50) ``` The weighted and unweighted results are returned together. Mast cells and proliferating fibroblasts are a pair where they disagree: with the original labels their co-localisation does not differ significantly between the outcomes, while with the confidences it is lower in disease. ```{r} #| label: colocresults subset(coloc$modelUnweighted, from == "Mast" & to == "Proliferating FBs") subset(coloc$modelWeighted, from == "Mast" & to == "Proliferating FBs") ``` `plotMetricBoxplot()` shows the per-image values behind the pair: ```{r} #| label: colocplot #| fig-height: 4 #| fig-width: 5 plotMetricBoxplot(coloc, spe, "Mast", "Proliferating FBs", "disease_status", metric = "lfunction") ``` The difference is largely driven by the `Control` samples showing increased attraction after weighting. We can take a closer look at one such sample, `VUHD113`: ```{r} #| label: vuhd113plots #| fig-height: 4.5 #| fig-width: 6.5 # the two cell types with their original labels plotHardSpatialPair(spe, "VUHD113", "Mast", "Proliferating FBs", tree) # the same two cell types, weighted by their confidences plotWeightedSpatialPair(spe, tree, confList, "VUHD113", "Granulocytes", "Mast", "Fibroblast lineage", "Proliferating FBs") ``` ```{r} #| label: vuhd113 #| echo: false img = spe$imageID == "VUHD113" xy = SpatialExperiment::spatialCoords(spe)[img, ] label = spe$cellType[img] mast = label == "Mast" fb = confList[["Fibroblast lineage"]] weight = setNames(fb[["Proliferating FBs"]], fb$cellID)[spe$cellID[img]] weight[is.na(weight)] = 0 nearMast = apply(xy, 1, function(p) any(sqrt(colSums((t(xy[mast, ]) - p)^2)) <= 50)) ``` With the original labels, the image contains only `r sum(label == "Proliferating FBs")` proliferating fibroblasts, none of which lie within 50 µm of any of its `r sum(mast)` mast cells, so the L-function for the pair is at its minimum. With the confidences, proliferating-fibroblast weight (blue) is spread across `r sum(weight > 0)` fibroblast-lineage cells. Of the `r sum(weight >= 0.5)` cells with a confidence of at least 0.5, `r sum(weight >= 0.5 & label == "Alveolar FBs")` are labelled as alveolar fibroblasts, and `r sum(weight > 0 & nearMast)` of the weighted cells lie within 50 µm of a mast cell, which produces the apparent attraction. The mast cells (red) keep high confidences after weighting, so the change comes almost entirely from the fibroblasts. # Niches Niches are recurring local neighbourhoods: groups of cells that share a similar mix of surrounding cell types. `runNicheIdentification()` describes every cell by its local indicators of spatial association from lisaClust, at several radii, and clusters those profiles into `k` niches. It does this twice, with the hard labels and with the weights, in which each neighbour counts towards a cell type in proportion to its confidence. ```{r} #| label: niches #| results: hide niches = runNicheIdentification(spe, tree, confList, r = c(20, 50), k = 8) ``` The niche plots take the niche assignments returned by `runNicheIdentification()`. `nicheCompositionMap()` shows which cell types are enriched in each niche, here for the unweighted (left) and weighted (right) niches: ```{r} #| label: nichemap #| fig-height: 7 #| fig-width: 5.5 #| layout-ncol: 2 nicheCompositionMap(spe, niches$nicheUnweighted) nicheCompositionMap(spe, niches$nicheWeighted) ``` Most niches have similar compositions with and without weights. The main exception is the alveolar niche (`niche8`): without weights it combines AT1 and AT2 cells with capillaries and neutrophils, whereas with weights it is split into an epithelial niche of AT2, transitional AT2 and AT1 cells (`niche1`) and a vascular niche of venous, arteriole and capillary cells (`niche8`). An alluvial plot, drawn here with the ggalluvial package, shows how cells move between the unweighted and weighted niches: ```{r} #| label: nichealluvial #| fig-height: 5 #| fig-width: 6.5 library(ggalluvial) flows = as.data.frame(table(unweighted = niches$nicheUnweighted, weighted = niches$nicheWeighted), responseName = "cells") ggplot(flows, aes(axis1 = unweighted, axis2 = weighted, y = cells)) + geom_alluvium(aes(fill = unweighted), width = 1/4) + geom_stratum(width = 1/4) + geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 3) + scale_x_discrete(limits = c("Unweighted", "Weighted"), expand = c(0.1, 0.1)) + labs(y = "Cells", fill = "Unweighted niche") + theme_minimal() ``` ```{r} #| label: nichemoves #| echo: false moves = table(niches$nicheUnweighted, niches$nicheWeighted) share = function(from, to) paste0(round(100 * moves[from, to] / sum(moves[from, ])), "%") ``` The split of the alveolar niche is visible as flows of cells: the unweighted `niche8` divides between the weighted `niche8` (`r share("niche8", "niche8")`), `niche1` (`r share("niche8", "niche1")`) and `niche4` (`r share("niche8", "niche4")`). `plotNicheSpatial()` shows where a niche's cells are in an image: ```{r} #| label: nichespatial #| fig-height: 5 #| fig-width: 9 plotNicheSpatial(spe, "VUHD038", "niche1", nicheCol = niches$nicheWeighted) ``` The proportion of each niche can be tested against the outcome in the same way as the cell types: ```{r} #| label: nichetest nicheTests = testNicheOutcome(spe, niches$nicheUnweighted, niches$nicheWeighted, outcomeCol = "disease_status", subject = "subject") nicheTests$testWeighted nicheTests$testUnweighted ``` A positive estimate means that the niche is more common in disease. Without weights, the alveolar niche (`niche8`) is less common in disease. With weights, both parts of the divided alveolar niche, the epithelial `niche1` and the vascular `niche8`, tend to be less common in disease. Neither difference remains significant after adjusting for the number of niches tested. # Practical notes - **Use at least 100 bootstraps** (`nBoot`) for real analyses, and set `seed` so the confidences can be reproduced. - **Large datasets:** `cores` runs the bootstrap on several cores, and `subprocess = TRUE` runs each group in a fresh R session, which keeps the workers small. The analysis functions take `cores` too. - **An alternative to the bootstrap:** `crossfitLineageProbabilities()` returns cross-fitted class probabilities in the same format as `bootstrapLineageConfidence()`, so they can be used in exactly the same way. - **Your own pipeline:** `confidenceWeights()` gives the cells-by-labels weight matrix for any function that takes `cellWeights`, such as `spicyR::spicy()` and `lisaClust::lisa()`. # References - Vannan et al. (2025). Xenium spatial transcriptomics of pulmonary fibrosis and control lung biopsies, available as `SpatialDatasets::spe_Vannan_2025()`; the subset used here is `lungSubset`. - [spicyR](https://bioconductor.org/packages/spicyR): spatial linear and mixed-effects modelling of cell co-localisation. - [lisaClust](https://bioconductor.org/packages/lisaClust): clustering of local indicators of spatial association. - [Statial](https://bioconductor.org/packages/Statial): the Kontextual measure of relative co-localisation. # Session info ```{r} #| label: sessioninfo sessionInfo() ```