---
title: "SegStability: correspondence from shared molecules"
author: "Rama Krishna Pudota"
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{SegStability: correspondence from shared molecules}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# Introduction

Different cell boundaries can assign the same detected molecule to different
cells. SegStability describes that correspondence and its effect on categorical
cell-state conclusions. It compares alternative segmentations of **one
immutable molecule universe**. Each molecule needs an exact, stable ID, and
each segmentation needs a complete cell census and one explicit status per
molecule: `assigned`, `unassigned`, or `ambiguous`. An ambiguous molecule is
never assigned to a candidate cell by the package.

The examples below use 24 synthetic molecules. Coordinates and state labels are
illustrative and have no biological interpretation. The example includes empty
cells, an unassigned molecule on each side, and one ambiguous molecule with two
recorded candidates. See `?segmentationExample` for its tables.

# Installation

Before Bioconductor acceptance, install the development package from its
public GitHub repository. The commands below are for an interactive R session;
this vignette does not run them.

```{r installation, eval=FALSE}
install.packages(c("BiocManager", "remotes"))
BiocManager::install(c("S4Vectors", "BiocStyle"))
remotes::install_github("ramakrishna-p21/SegStability", build_vignettes = TRUE)
```

```{r data}
library(SegStability)
data("segmentationExample", package = "SegStability")
d <- segmentationExample
names(d)
head(d$molecules[, c("molecule_id", "feature_id")])
```

# Build and inspect a segmentation set

The constructor matches assignment rows by molecule ID. Their row order has no
identity meaning. `cells` lists every cell, including those with no assigned
molecules. Candidate evidence is separate from the assignment's unresolved
status.

```{r construction}
make_partition <- function(side) {
    a <- d$assignments[[side]]
    ambiguity <- d$ambiguity[[side]]
    if (!is.null(ambiguity)) {
        ambiguity <- S4Vectors::SimpleList(
            evidence = S4Vectors::DataFrame(ambiguity$evidence),
            candidates = S4Vectors::DataFrame(ambiguity$candidates))
    }
    SegmentationPartition(d$molecules, side, d$cells[[side]],
                          a$molecule_id, a$status, a$cell_id,
                          ambiguity = ambiguity)
}
ss <- SegmentationSet(d$molecules,
    list(first = make_partition("first"),
         second = make_partition("second")),
    specimen_id = "synthetic-example")
segmentationNames(ss)
segmentationInfo(ss)[, c("name", "n_cells", "n_molecules")]
```

`[[` retrieves a named, universe-bound partition. `assignmentData()` returns
the original external molecule IDs, semantic statuses, and typed cell IDs in
shared-universe order. Missing cell IDs for unresolved assignments are retained.
Coordinate declarations can be read from either set or comparison summaries;
reading a declaration does not certify spatial compatibility.

```{r access}
cellCensus(ss[["first"]])
tail(assignmentData(ss, "first"), 4)
summary(ss)$coordinate_info
```

# Correspondence across segmentations

```{r comparison}
cmp <- compareSegmentations(ss, "first", "second")
summary(cmp)$coordinate_info
edge <- correspondenceEdges(cmp)
edge_table <- data.frame(
    first = cellSummary(cmp, "first")$cell_id[S4Vectors::from(edge)],
    second = cellSummary(cmp, "second")$cell_id[S4Vectors::to(edge)],
    weight = S4Vectors::mcols(edge)$weight)
edge_table
componentSummary(cmp)[, c("component_id", "component_class")]
```

An edge weight counts molecules uniquely assigned to its two endpoint cells.
For a first-side cell, $N_X$ is its complete uniquely assigned molecule count;
for a second-side cell, the corresponding count is $N_A$. Endpoint fractions
divide the edge weight by that endpoint's own denominator. The denominator is
independent of how many molecules the other segmentation assigned. For example,
the first-side `a1` has five assigned molecules, including `m21`, which is
unassigned on the second side;
its edge to `b1` therefore has fraction `4/5`.

Dominance uses the complete graph. Exact maximum ties remain ties; a reciprocal
anchor requires one unique maximum at each endpoint. Connected components
include empty and unmatched census cells. Component labels such as `split` and
`merge` describe graph shape, rather than biological division or a preferred
segmentation.

```{r accounting}
cellSummary(cmp, "first")[, c("cell_id", "assigned_mass", "matched_mass",
                               "unmatched_mass", "dominance_count")]
assignmentSummary(cmp)[assignmentSummary(cmp)$n_molecules > 0, ]
```

# Categorical cell-state stability

Cell states are supplied separately and keyed by census cell ID. Comparable
state names should use one shared vocabulary. A missing label is unresolved,
and never counted as disagreement. The example has 20 shared molecules, of
which 18 have labels on both endpoint cells: 15 share a label and three differ.
The original first and second assigned denominators are 21 and 22.

```{r states}
state <- stateStability(cmp, d$states)
state$global[, c("matched_mass", "resolved_mass", "same_state_mass",
                 "different_state_mass", "unresolved_label_mass")]
```

# Reading support sensitivity

`supportSensitivity()` selects edges in each requested direction when their
weight divided by the **original** assigned mass of that direction's endpoint
meets the threshold. The directions are separate calculations. Filtering does
not change the primary graph, its dominance, components, census, or
denominators.

```{r support}
support <- supportSensitivity(cmp, thresholds = c(0, 0.6),
                              states = d$states)
as.data.frame(support[, c("direction", "threshold", "original_matched_mass",
                         "retained_mass", "excluded_matched_mass",
                         "unmatched_mass", "assigned_mass")])
```

`retained_mass` counts shared molecules on selected edges;
`excluded_matched_mass` counts originally shared molecules on filtered edges.
`unmatched_mass` is the original assigned mass with no partner in the other
segmentation. These are different losses. `assigned_mass` stays at the original
denominator. The `cells` and `components` columns hold detailed tables for
each direction and threshold; they retain original census and component rows.
State columns describe the selected edges when labels were supplied. With
`states = NULL`, state-dependent fields are `NA`, rather than zero. Missing
labels remain unresolved. `endpoint_unavailable_cells` counts cells with a
missing supplied state. `total_nonestimable_cells` counts zero-assigned cells;
`conditional_nonestimable_cells` counts cells with no selected resolved mass.
These fields distinguish an unavailable ratio from a measured value of zero.

```{r support-detail}
support$cells[[2]][["first"]][, c("cell_id", "assigned_mass",
                                 "selected_retained_mass",
                                 "selected_removed_mass")]
support$components[[2]][, c("component_id", "component_class",
                            "original_matched_mass", "selected_retained_mass")]
```

# Optional spatial evidence

`anchorTopology()` requires an observation window, polygon geometry and a
qualified common coordinate declaration. Eligible reciprocal anchors and their
fan classes are derived from the complete correspondence graph. Missing or
ineligible geometry remains explicit. `compareGeometry()` measures polygon
overlap as separate evidence; its maximum-IoU matching is a derived alignment
and does not replace shared-molecule correspondence. Neither function repairs,
transforms or clips source geometry.

This toy spatial example is run only when the optional `sf` package is present.
The coordinate declaration here is supplied for the constructed toy polygons;
it should never be inferred from an unqualified real dataset.

```{r geometry}
if (requireNamespace("sf", quietly = TRUE)) {
    ring <- rbind(c(1, 1), c(2, 1), c(2, 2), c(1, 2), c(1, 1))
    polygon <- sf::st_sfc(sf::st_polygon(list(ring)))
    geometry <- list(
        first = sf::st_sf(cell_id = "a1", geometry = polygon),
        second = sf::st_sf(cell_id = "b1", geometry = polygon))
    info <- S4Vectors::SimpleList(
        units = "micrometres", origin = c(0, 0),
        axis_directions = c("right", "up"),
        transforms = S4Vectors::SimpleList(),
        registration = "synthetic common XY",
        specimen_id = "synthetic-example")
    window <- c(xmin = 0, xmax = 10, ymin = 0, ymax = 10)
    topology <- anchorTopology(cmp, geometry, window, info)
    overlap <- compareGeometry(cmp, geometry, info)
    head(topology$anchors)
    head(overlap$pairs)
}
```

# Agreement among segmentation runs

When the same molecule universe has multiple stochastic segmentation runs,
`replicateAgreement()` reports all nine assignment-status combinations and an
adjusted Rand index (ARI) on **jointly assigned molecules only**. Run-local cell
names are never treated as corresponding identities. This illustration relabels
the synthetic first partition to show the label-invariant calculation; it is
not a claim about variability of a real stochastic algorithm.

```{r replicates}
a <- d$assignments$first
a$cell_id <- ifelse(is.na(a$cell_id), NA_character_,
                    paste0("runA_", a$cell_id))
b <- d$assignments$first
b$cell_id <- ifelse(is.na(b$cell_id), NA_character_,
                    paste0("runB_", b$cell_id))
runA_cells <- data.frame(cell_id = paste0("runA_", d$cells$first$cell_id))
runB_cells <- data.frame(cell_id = paste0("runB_", d$cells$first$cell_id))
runA <- SegmentationPartition(d$molecules, "runA", runA_cells,
    a$molecule_id, a$status, a$cell_id)
runB <- SegmentationPartition(d$molecules, "runB", runB_cells,
    b$molecule_id, b$status, b$cell_id)
runs <- SegmentationSet(d$molecules, list(runA = runA, runB = runB),
                        "synthetic-example")
replicateAgreement(runs, c("runA", "runB"))$pairwise_summary
```

# Interoperability and interpretation

`fromMoleculeExperiment()` extracts retained molecule IDs from an explicit
assay/sample and requires separately supplied assignments and a complete cell
census. It returns arguments for `SegmentationSet()`; source row order and
coordinates cannot reconstruct a lost molecule ID. The small example in
`?fromMoleculeExperiment` runs when that optional package is installed.
`geometryFromSpatialFeatureExperiment()` extracts a keyed polygon sidecar with
an explicit geometry-to-cell map and census; see its help page for a runnable
guarded example. Geometry does not establish molecule correspondence or
registration on its own.

## A qualified Xenium application

The [ovrlpy study](https://doi.org/10.1038/s41587-026-03004-8) released saved
segmentations of Xenium mouse-brain data. A development-time replay used
`Xenium_V1_FF_Mouse_Brain_MultiSection_1`, replicate 1, and source-defined
500-µm tile (9,9): x and y in [4500,5000) µm. The immutable ROI universe has
199,448 original molecule IDs. The comparison is between the original vendor
assignment and saved Baysor v0.7.1 output. It exercised
`compareSegmentations()` and `supportSensitivity()`; it did not evaluate cell
states, paired geometry, topology, or stochastic replicates.

The small machine-readable record is installed at
`extdata/validation/ovrlpy_validation_summary.csv` and documented in
`?ovrlpy_validation_summary`. The optional development-time scripts
`scripts/prepare_ovrlpy_validation.py` and
`scripts/reproduce_ovrlpy_validation.R` reconstruct the inputs from the
specified public source files and check the reported results. They are not run
while building or checking this package. No raw molecules, large segmentation
outputs, or benchmark working files are bundled. This evaluated chunk reads
only the bundled summary:

```{r real-data-summary}
path <- system.file("extdata", "validation",
                    "ovrlpy_validation_summary.csv", package = "SegStability")
validation <- read.csv(path, colClasses = "character")
validation[validation$key %in% c("roi_universe_size", "correspondence_edges",
                                 "shared_mass"), c("key", "value")]
```

Each method's census includes native-center cells and cells assigned an ROI
molecule; a cell crossing the tile has a denominator restricted to its own ROI
molecules. The saved outputs differed in upstream preparation, so this is not
an isolated boundary-method comparison. Correspondence and directional
support describe this qualified input pair; they neither rank segmentation
accuracy nor establish biological ground truth.

SegStability does not run a segmentation method, infer missing provenance,
choose a preferred segmentation, or establish biological ground truth. State
agreement is conditional on supplied labels; geometry overlap is different
evidence from molecule sharing. Interpret results within those inputs and
their validation reports.

```{r session-info}
sessionInfo()
```