---
title: "Designing and evaluating a spatial multi-omics workflow"
author: "Tianyao Lu"
output:
  BiocStyle::html_document:
    toc: true
vignette: >
  %\VignetteIndexEntry{Designing and evaluating a spatial multi-omics workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "##", fig.width = 9,
                     fig.height = 5, error = FALSE, cache = FALSE)
library(SpaMTP)
```

# Analysis questions determine the workflow

A comprehensive report should expose the consequences of analytical choices.
Its purpose is not to produce an attractive embedding or maximize the number of
significant features. The workflow follows the analysis chain of the historical
mouse-brain vignette: inspect measurements and correspondence; compare individual,
spatial and integrated representations; characterize regions; examine paired
features; then connect measured identities to pathways and networks.

| Question | Native interface | Evidence to inspect |
|:--|:--|:--|
| What was measured and retained? | `normalizeSMData()`, native containers and image access | Selected layer, QC exclusions, scale and acquisition provenance |
| What changes when spatial information is added? | `runMetabolicPCA()`, `scaleSMData()`, `runSpatialGraphPCA()` | Same selected features and all retained observations in both representations |
| Does integration improve the representation? | `multiOmicIntegration()`, `getKmeanClusters()` | Individual versus joint structure; external-reference ARI, approximate silhouettes and sensitivity across K |
| Which features characterize a region? | `findAllDEMs(method="markers")`, `demsHeatmap()` | Pairwise AUC, Cohen effects, detection, mean expression and spatial-block sensitivity |
| Which differences reproduce across biological samples? | `findAllDEMs(method="replicate")` | Explicit design, replicate means, effect intervals and the complete BH family |
| Is an association explained by regional composition? | `findCorrelatedFeatures()`, `plot3DFeature()` | Raw versus covariate-residual correlations and matched spatial layers |
| Which members support a pathway? | `buildPathwayIndex()`, `createPathwayAssay()`, `findRegionalPathways()` | Measured coverage, ambiguity exclusions, complete ranks and leading edges |
| Are pathway members coordinated and where are they? | `runRAMPGeseca()`, `plotRegionalPathways()`, `pathwayNetworkPlots()` | The same identity assay and membership index in every downstream view |

The main example uses a bundled **real DHB/Visium mouse-brain subset** from
SpaMTPData 1.1.0. It retains 160 paired locations, all 2,075 measured mass features
and all 32,285 RNA features. On each lesion side, the 80 locations nearest the
archived ACB coordinate centroid are retained. Selection uses coordinates and
archived labels, not expression or statistical results. The ACB, CP and nearby
regions remain available for descriptive comparisons.

The [gene-mapping guide](Gene_Identifier_Mapping.html) and
[annotation guide](Metabolite_Annotation_Pipeline.html) use the same bundled
real-data demonstrators for their respective tasks. The installed full-data
recipes remain available for analyses of the complete sections. Only the
replicate-design example below uses simulated observations.

# Supply modalities and study design explicitly

```{r inputs}
x <- readRDS(system.file("extdata", "dhb_demo.rds", package = "SpaMTP"))
spec <- list(
    main = list(type = "metabolomics", layer = "counts", species = "Mus musculus"),
    SPT = list(type = "transcriptomics", layer = "counts", species = "Mus musculus")
)
stopifnot(identical(dim(x), c(2075L, 160L)),
          nrow(SingleCellExperiment::altExp(x, "SPT")) == 32285L)
table(x$RegionLoupe, x$lesion)
S4Vectors::metadata(x)$vignette_subset[c("rule", "features", "source")]
```

ACB denotes nucleus accumbens, CP caudoputamen, OT olfactory tubercle and aco
anterior commissure, olfactory limb, following the archived region labels.
These two sides belong to **one section**, not independent biological replicates.
The very small aco compartment is retained in the input but does not meet the
workflow's minimum of three observations for regional marker scoring. Cropping
also changes boundary neighborhoods, so these results describe this subset.
The selected sides contain different proportions of ACB and CP. A difference
associated with lesion status can therefore reflect regional composition; the
raw and region-adjusted associations below illustrate this distinction.

```{r input-spatial-layout}
coordinates <- as.data.frame(SpatialExperiment::spatialCoords(x))
colnames(coordinates)[1:2] <- c("x", "y")
coordinates$region <- x$RegionLoupe
coordinates$lesion <- x$lesion
ggplot2::ggplot(coordinates, ggplot2::aes(x, y, colour = region, shape = lesion)) +
    ggplot2::geom_point(size = 2) + ggplot2::scale_y_reverse() + ggplot2::coord_equal() +
    ggplot2::theme_classic() +
    ggplot2::labs(x = "Original image x", y = "Original image y")
```

The workflow preserves the supplied layers. Its declared baseline uses native
library-size normalization followed by conversion to log2 units for raw MS/RNA
counts, and log2(1 + value) for protein counts. For already transformed data,
choose `normalization="none"`. This baseline does not replace vendor-specific
spectral processing, RNA count QC or a platform-appropriate normalization design.
In particular, logged values must not be logged again or labelled raw counts.

Biological replication is an explicit metadata field. A pixel, technical pool,
spatial block, image channel or integration cluster does not create an independent
specimen. Transcriptomics, metabolomics and proteomics can be declared as three
modalities; H&E, DAPI and RGB channels remain contextual images unless their
measurements support a separate stated assay model.

# Compare representations using the same observations and features

```{r workflow}
fit <- runSpaMTPWorkflow(x, spec, group = "lesion", regions = "RegionLoupe",
    clusters = 3, npcs = 8, max_features = 1000, association_features = 4,
    coordinate_units = "original image pixels",
    structure = list(primary = "spatial", k_grid = c(3, 4),
                     reference = "RegionLoupe", umap = FALSE,
                     stability_blocks = 3))
fit$structure$metrics
```

For this short demonstration, PCA uses at most 1,000 variable features per
modality, eight components, two candidate values of K and three spatial-block
omissions. All measured features remain in the object and regional marker
analysis. Displays use the first two components (`umap = FALSE`); set `umap = TRUE`
to add UMAP displays. These are explicit teaching settings, not newly recommended
analysis defaults.

`runMetabolicPCA(features=...)` records all selected features explicitly.
Ordinary and graph PCA use the same variable features, centred and scaled in the
same way. Graph PCA uses a sparse Laplacian for every retained observation and
keeps different samples in separate spatial graphs. It is not fitted on a small
report preview and then discarded: graph embeddings feed the spatial joint
representation and can be selected as the downstream primary representation.

The joint representation uses `multiOmicIntegration()` to concatenate standardized
modality embeddings. This native implementation is **not Seurat WNN**. Recovering
and benchmarking WNN requires a separate, validated method implementation.
UMAP is a visualization of a representation; clusters and metrics use its complete
retained dimensions. A nicer-looking UMAP is not a performance test.

The number of clusters and primary method are chosen before examining reference
ARI. Approximate silhouettes assess separation, while ARI compares to supplied
external annotation. Spatial-neighbor agreement is a smoothness diagnostic; spatial
regularization increases it by construction. Spatial-block omission refits
clustering **conditional on full-data embeddings**. It does not validate feature
selection, normalization or generalization to new specimens.

# Region characterization is richer than a mean difference

```{r markers}
markers <- fit$region_analysis$modalities$main$markers
head(markers$DEMs[, c("gene", "cluster", "logFC", "mean_auc",
                     "min_auc", "mean_cohen", "direction_stability")])
demsHeatmap(markers, n = 2, only.pos = TRUE, logfc.threshold = 0,
    order.by = "mean_auc", fontsize_row = 9, fontsize_col = 10)
```

The native marker mode delegates pairwise scoring to `scran::scoreMarkers()`.
AUC measures pairwise concordance, including half credit for ties. Mean and minimum
AUC distinguish broadly specific markers from features that separate only some
regions. Cohen effects account for within-region variability. Summaries equally
weight other eligible regions so that the largest anatomical compartment does not
automatically dominate the reference. Detection fractions refer to non-zero values
on the explicitly selected layer; residualized layers need different interpretation.

Coordinate-only blocks support a leave-one-block-out sensitivity analysis. The
result records the range of Cohen effects, marker-rank changes, direction agreement
and omitted fits that were ineligible because a region became too small. These
ranges are **not confidence intervals**. Region labels are held fixed, including
when learned from these data, and no pixel-based P values are manufactured.

`demsHeatmap()` uses the analysis expression matrix directly. Its FDR ordering now
sorts numerical FDR; it no longer sorts a regulation label or applies CPM to an
already transformed marker matrix. The heatmap, detection plot, spatial feature
maps and interactive explorer use the same ranked feature identifiers.

# Keep biological-replicate inference separate

The published brain sections do not provide independent biological replicates.
This design example therefore explicitly simulates four paired patients, each
with three observations per condition; it is not an analysis of the real sections.

```{r replicate-design, eval=TRUE}
design <- expand.grid(pixel = seq_len(3), condition = c("control", "treated"),
                      patient_id = paste0("simulated_patient", seq_len(4)))
set.seed(1234)
expression <- matrix(rnorm(6 * nrow(design), mean = 5, sd = 0.3), nrow = 6,
    dimnames = list(paste0("feature", seq_len(6)),
                    paste0("observation", seq_len(nrow(design)))))
expression[1:2, design$condition == "treated"] <-
    expression[1:2, design$condition == "treated"] + 1.2
replicate_example <- SingleCellExperiment::SingleCellExperiment(
    assays = list(logcounts = expression), colData = S4Vectors::DataFrame(design))
replicate_de <- findAllDEMs(replicate_example, ident = "condition", assay = "main",
    slot = "logcounts", method = "replicate", replicate = "patient_id",
    contrasts = list(treatment_control = list(
        numerator = "treated", denominator = "control", paired = TRUE)))
stopifnot(identical(replicate_de$tests$treatment_control$status, "completed"))
head(replicate_de$tests$treatment_control$table)
```

This mode fits limma to equally weighted biological-replicate means of the selected
values. Paired contrasts retain complete pairs; at least three complete pairs (or
three independent replicates per arm) are required. The result retains the design,
coefficient, residual degrees of freedom, replicate expression matrix, pointwise
moderated intervals, and every tested feature. BH correction is within modality
and contrast. A `logFC` column is a difference on the selected workflow scale;
it is not automatically a fold change in raw abundance. RNA count models may
instead require independently designed count pseudobulking and edgeR/DESeq2.

The historical random-pool mode remains available for compatibility and emits a
warning about technical partitions. It is never used by the automatic workflow.

# Ask whether colocalization survives regional composition

```{r associations}
head(fit$associations[[1]])
```

`findCorrelatedFeatures(covariates=...)` retains raw Pearson correlation and the
correlation between both features after regression on the declared sample,
replicate and region fields. The workflow queries regional marker anchors against
all PCA-selected target features. It records those selection sets and the matched
observations. Strong raw association can disappear after removing a common regional
shift. A residual association can still reflect spatial autocorrelation or other
unmeasured structure, and selected pairs require independent validation.

Exact feature identifiers work for RNA and protein as well as metabolites; no
artificial m/z metadata is needed. The report includes paired scatter plots and
native `plot3DFeature()` views of the same selected pairs.

# Map identities before pathway statistics

Here the bundled 240-location human-brain subset supplies the complete measured
339-gene panel and its
archived `cout_labels` cell labels. These are descriptive groups from a single
section, not independent replicates or newly validated anatomical boundaries.

```{r pathway-chain, eval=TRUE}
database <- loadSpaMTPDatabase(
    c("source_df", "analyte", "analytehaspathway", "pathway"), version = "3.0.7")
gene_reference <- SpaMTPdb::spaMTPdbGeneReference(version = "2026.7.7")
gene_index <- buildGeneMappingIndex(database$source_df,
                                  gene_reference = gene_reference)
index <- buildPathwayIndex(database, gene_index = gene_index)
real_data <- readRDS(system.file("extdata", "human_demo.rds", package = "SpaMTP"))
real_data <- normalizeSMData(real_data, "LogNormalize", assay = "SPT", verbose = FALSE)
identity <- createPathwayAssay(real_data, analyte_type = "genes",
    assay = "SPT", slot = "logcounts", new_assay = "gene_ids",
    pathway_index = index, organism = "Homo sapiens")
identity_markers <- findAllDEMs(identity, ident = "cout_labels", assay = "gene_ids",
    slot = "logcounts", method = "markers", spatial_blocks = 0)
# Build a feature-by-region matrix of the complete identity-level effects.
ranks <- identity_markers$expression
for (region in colnames(ranks)) {
    tab <- identity_markers$DEMs[identity_markers$DEMs$cluster == region, ]
    ranks[, region] <- tab$logFC[match(rownames(ranks), tab$gene)]
}
set.seed(1234)
regional <- findRegionalPathways(identity, ident = "cout_labels",
    ranks = list(genes = ranks), pathway_index = index,
    min_path_size = 5, max_path_size = 500)
stopifnot(nrow(regional) > 0, ncol(identity) == ncol(real_data))
plotRegionalPathways(regional, database = database)
```

The ranks interface uses the complete measured feature ranking, with no FDR or
foreground prefilter. It retains every used rank and the complete eligible pathway
family. Competitive feature-set P values do not establish reproducibility across
specimens. Joint gene/metabolite ranks can be RMS-scaled per modality, but pathways
with many measured genes can still be dominated by transcriptomic coverage.

For metabolomics, `createPathwayAssay(metabolite_ambiguity="exclude")` excludes
mass features with multiple candidate compound identities. Original candidates,
exclusions, aggregation weights and measured members remain inspectable. The
workflow uses this policy by default; the public native interface retains the
historical `expand` option. Sparse matrices stay sparse during identity aggregation.
The network viewer can consume this identity assay directly, preventing re-expansion
of excluded mass candidates. GESECA, scores and networks therefore use the same
identity definition as regional enrichment.

The pathway dot plot uses **adjusted** P values for significance outlines and
`1 - Jaccard similarity` for pathway ordering. It handles a single pathway. Mouse
genes are not silently uppercased into a human HGNC index; the registered mouse
recipes currently show compound pathways until an explicit compatible gene identity
or orthology reference is supplied.

# Apply the workflow to complete reference sections

The installed scripts retain the full DHB/Visium and FMP10/Visium analyses,
including annotation, spatial integration, pathways and report generation. The
vignette builds the small demonstration above; running the complete recipes is
an explicit analysis step with larger resource and runtime requirements.

```{r reference-recipes}
prepare <- system.file("workflows", "prepare_mouse_brain.R", package = "SpaMTP")
run <- system.file("workflows", "run_mouse_brain.R", package = "SpaMTP")
source(prepare)
source(run)
formals(prepareMouseBrain)
formals(runMouseBrain)
```

Call `prepareMouseBrain(outputDir = a_new_input_directory)` to resolve the
pinned native SpaMTPData 1.1.0 resources, then
`runMouseBrain(a_new_input_directory, outputDir = a_new_report_directory,
case = "dhb")`. An explicit `nativeDir` selects verified local files and defaults
to offline mode. Existing output directories are not overwritten. Source
`examples/pathway_network_mouse_brain_demo.R` from the installed package for
the corresponding full-data network recipe.

For FMP10/Visium, preparation additionally requires a measured target spot radius
in the original image coordinate system through `spatialConfig`. For the exact
registered V11L12-109_B1 pair, `extdata/mouse_brain_scalefactors.json` retains the
Space Ranger value `spot_diameter_fullres = 376.942`, from
[Vicari et al.](https://doi.org/10.17632/w7nw4km7xd.1); use half that value as
`target_radius`. Do not reuse it for another sample. H&E is optional and omitted
from the native resources. Independently supplied images require explicit scale
and orientation checks.

Pairing is documented for these registered specimens; species alone does not
establish pairing. Both full-data runners re-annotate with the pinned chemical
database and an explicit positive-mode DHB or FMP-10 index. Ambiguous compound
assignments are excluded; mouse genes are not converted to human genes. The
[annotation guide](Metabolite_Annotation_Pipeline.html) executes annotation,
compound pathways and network generation on the small real DHB example.

To regenerate the bundled subsets, source `workflows/prepare_vignette_data.R`
and call `prepareVignetteData()` with a new output directory. See
`help("vignetteData")` for the deterministic selection rules and source checksums.

# Render and retain a reviewable analysis

```{r report, eval=TRUE}
report <- renderSpaMTPReport(fit, tempfile("dhb-demo-report-"))
stopifnot(file.exists(report))
```

These panels come from the report just computed on all 160 retained locations.
Compare the archived region labels in the embeddings with the spatial partitions;
smoothness alone does not establish better recovery of anatomy. Cluster numbers
are local to each panel and do not identify matching cell types across methods.

```{r report-panels, out.width="100%"}
figures <- utils::read.csv(file.path(dirname(report), "figure_index.csv"))
selected <- match(c("structure_embedding", "structure_spatial"), figures$name)
stopifnot(!anyNA(selected))
knitr::include_graphics(file.path(dirname(report), figures$png[selected]))
```

Reports follow the analytical sequence above. They include a linked region/DE
explorer with feature search, directional AUC/effect/stability/FDR filters, spatial
zoom and pan, coordinated feature profiles and complete CSV exports. Native 3D and
network views are available when their input requirements are satisfied. Figures
are embedded as PNGs and also exported as PDFs. Open the HTML in a browser or an IDE
HTML preview that permits JavaScript; Markdown previews may disable interaction.
The output directory also contains the complete workflow RDS, settings, inputs used
by each stage, source and database provenance, and file checksums.

For methods, see the [Bioconductor marker-scoring documentation](https://bioconductor.org/books/release/OSCA.basic/marker-detection.html)
and the [GraphPCA paper](https://doi.org/10.1186/s13059-024-03429-x).

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