SpaMTP: Spatial Metabolite Annotation

Overview

SpaMTP provides analysis and visualisation methods for spatial metabolomics and paired spatial multi-omics experiments. Its primary containers are Cardinal’s MSImagingArrays/MSImagingExperiment and Bioconductor’s SpatialExperiment. Seurat is not an analysis dependency. Only the optional boundary-conversion example uses SeuratObject; all analysis runs on native containers. Explicit conversion functions connect existing Seurat workflows. Versioned annotation resources are supplied by the companion SpaMTPdb package and experiment data by SpaMTPData.

This vignette demonstrates a small, fully reproducible accurate-mass annotation workflow. It deliberately uses an in-memory database so that package checking does not require a network connection or a large data download.

See modern data access for replacing historical container-internal reads and writes, and the indexed annotation pipeline for the complete native annotation tutorial. The companion workflow demonstrates how SpaMTPData, SpaMTPdb and SpaMTP work together on published native resources; only its explicitly labelled historical conversion uses SeuratObject.

library(SpaMTP)

Bioconductor container model

The container depends on the data’s processing stage:

Processing stage Container SpaMTP mapping
Raw or continuous spectra Cardinal MSImagingArrays / MSImagingExperiment File-backed spectra and unequal mass axes remain in Cardinal/matter
Aligned or binned MSI SpatialExperiment Intensities in assays(), m/z and annotations in rowData(), pixel metadata in colData()
Pixel coordinates SpatialExperiment spatialCoords()
Optical images SpatialExperiment imgData()
Paired transcriptome SpatialExperiment altExp()

Raw imzML should be read and processed in Cardinal without materialising the spectra. Once features share a common axis, conversion is explicit:

For a self-contained file-I/O example, first write two synthetic unequal-axis spectra to a temporary imzML/ibd pair. These are not a biological acquisition. The same reader accepts a real imzML path; do not discard its backing ibd file until file-backed processing has finished.

binned_raw <- local({
  directory <- tempfile("imzml-example-")
  dir.create(directory)
  on.exit(unlink(directory, recursive = TRUE))
  input <- Cardinal::MSImagingArrays(
    centroided = TRUE, continuous = FALSE,
    spectraData = list(intensity = list(c(1, 2, 3), c(4, 5, 6, 7)),
                      mz = list(c(100, 101, 102), c(100, 100.25, 101, 104))),
    pixelData = Cardinal::PositionDataFrame(
      coord = data.frame(x = c(1, 2), y = c(1, 1)),
      run = factor(c("section", "section"))))
  path <- file.path(directory, "section.imzML")
  Cardinal::writeImzML(input, path, bundle = FALSE, verbose = FALSE)
  raw <- Cardinal::readImzML(path, memory = FALSE, as = "MSImagingArrays")
  raw <- normalizeSMData(raw) # Cardinal queues TIC processing.
  binSpaMTP(raw, resolution = 0.01, units = "mz")
})
## queued: intensity normalization
stopifnot(methods::is(binned_raw, "SpatialExperiment"), ncol(binned_raw) == 2L,
  max(SummarizedExperiment::rowData(binned_raw)$mz) >= 104)

The S4 binning method executes Cardinal’s queued processing and obtains a common axis before creating an SPE. The raw-array mass range is scanned across all spectra, not inferred from the first pixel alone. Calling asSpatialExperiment() on unaligned MSImagingArrays gives an explicit request to bin first; no resolution is guessed during conversion.

The reverse coercion is also registered, so either the SpaMTP generic or the standard S4 coercion can be used:

cardinal <- asCardinal(binned_raw)
cardinal <- methods::as(binned_raw, "MSImagingExperiment")
stopifnot(methods::is(cardinal, "MSImagingExperiment"),
          identical(dim(cardinal), dim(binned_raw)))

A binned MSI workflow

The following small object illustrates the native layout without external files.

intensity <- matrix(
  seq_len(24),
  nrow = 4,
  dimnames = list(
    paste0("mz-", c(100.01, 100.02, 101.01, 101.02)),
    paste0("pixel", seq_len(6))
  )
)
coordinates <- cbind(
  x = rep(seq_len(3), each = 2),
  y = rep(seq_len(2), times = 3)
)
spe <- SpatialExperiment::SpatialExperiment(
  assays = list(counts = intensity),
  rowData = S4Vectors::DataFrame(
    mz = c(100.01, 100.02, 101.01, 101.02)
  ),
  colData = S4Vectors::DataFrame(
    region = rep(c("edge", "core"), 3)
  ),
  spatialCoords = coordinates,
  sample_id = rep("sample01", ncol(intensity))
)

spe <- normalizeSMData(spe, normalisation.type = "LogNormalize")
## Stored LogNormalize values in assay `logcounts`.
plotSpatialFeature(spe, features = "mz-100.01", assayName = "logcounts")

For paired measurements, transcriptomic data are stored as an alternative experiment and matched by pixel name. Integration uses Bioconductor reduced dimensions rather than Seurat assay internals.

transcriptome_sce <- SingleCellExperiment::SingleCellExperiment(
  assays = list(counts = matrix(
    c(1, 2, 3, 4, 4, 2, 1, 3, 2, 5, 4, 1,
      3, 1, 5, 2, 1, 4, 2, 5, 5, 3, 1, 4),
    nrow = 4,
    dimnames = list(paste0("gene", seq_len(4)), colnames(spe))
  ))
)
spe <- addTranscriptome(spe, transcriptome_sce, name = "transcriptome")
spe <- normalizeSMData(spe, "LogNormalize", assay = "transcriptome")
## Stored LogNormalize values in assay `logcounts`.
set.seed(42)
spe <- multiOmicIntegration(spe, dims.list = list(1, 1))
SingleCellExperiment::altExpNames(spe)
## [1] "transcriptome"
SingleCellExperiment::reducedDim(spe, "integrated")
##           main_PC1 transcriptome_PC1
## pixel1 -1.98227963        -0.5992166
## pixel2 -0.03539622         1.1278169
## pixel3  0.32889099        -1.0438482
## pixel4  0.48584137        -0.6564457
## pixel5  0.57348117        -0.1763852
## pixel6  0.62946232         1.3480789

# Reuse scater directly for standard expression visualisation.
scater::plotExpression(
  spe,
  features = "mz-100.01",
  x = "region",
  exprs_values = "logcounts"
)

The joint embedding concatenates standardized PCA scores with equal weight per modality. It is not Seurat WNN, and does not claim equivalent results. Only main and transcriptome are included by default; supply modalities explicitly to include another measured modality. Derived pathway scores are not included automatically.

Alignment and pixel mapping

Spatial mapping uses matching sample_id values and a common coordinate system. Independent slices are never connected simply because their x/y coordinates overlap. The example below maps MSI to paired transcriptomic pixel centres; lower-resolution spots require an explicit radius or polygon.

target <- SpatialExperiment::SpatialExperiment(
  assays = SummarizedExperiment::assays(transcriptome_sce),
  colData = SummarizedExperiment::colData(spe),
  spatialCoords = SpatialExperiment::spatialCoords(spe)
)
mapped <- mapSpatialOmics(
  spe, target, ST.hires = TRUE, SM.pixel.width = 1
)
stopifnot(all(mapped$mapped), all(mapped$n_source_pixels == 1))
SingleCellExperiment::altExpNames(mapped)
## [1] "transcriptome"

affine <- diag(3)
affine[1:2, 3] <- c(0.2, -0.1)
aligned <- alignSpatialOmics(spe, target, transformation = affine)
checkAlignment(aligned, target)

MSI values are averaged over matched source pixels. Unmapped targets are retained with mapped = FALSE and zero-filled intensities unless dropUnmapped = TRUE. These missing measurements must not be interpreted as observed biological zeros. Optical images are stored with addSpatialImage() and imgData(); their scale factors change display dimensions, not the stored spatial coordinates.

Binning and statistical reuse

binned <- binSpaMTP(spe, resolution = 0.05, units = "mz")
stopifnot(isTRUE(all.equal(
  as.numeric(Matrix::colSums(SummarizedExperiment::assay(binned, "counts"))),
  as.numeric(Matrix::colSums(SummarizedExperiment::assay(spe, "counts")))
)))

Aligned-matrix binning assigns each feature once and conserves total intensity under sum aggregation. Raw Cardinal input retains Cardinal’s reference-window binning semantics. Re-annotate the newly binned features before using annotation-dependent analyses.

TMM factors are computed by edgeR::normLibSizes(). Differential testing uses limma::lmFit() and limma::treat(), while runPooling() uses edgeR::sumTechReps() for technical partitions. Random pixel pools are not independent biological replicates and cannot support population-level inference as if they were. Use independent biological samples for such inference.

Native feature plots and annotation ranking

MSI plots read masses and annotations from rowData(), including when feature IDs are not formatted as mz-.... A plusminus window includes every peak within the tolerance of the nearest measured mass. Annotation plots sum all matching peaks; overlapping windows count each peak once.

SummarizedExperiment::rowData(spe)$all_IsomerNames <-
  c("Candidate A", "Candidate A", "Candidate B", "Candidate C")
imageMZAnnotationPlot(spe, "Candidate A", slot = "counts",
                      dark.background = FALSE)

massIntensityPlot(spe, group.by = "region", mz.labels = 100.01)

Use spatialMZPlot(spe, 100.01, images="optical") after attaching an image with addSpatialImage(). Its imgData() scale factor relates raster pixels to spatialCoords() units. Samples are displayed separately. Former Seurat segmentation/molecule/blending arguments are not supported. Interactive mass windows are available without Seurat; creating the app does not launch it:

app <- interactiveSpatialPlot(spe)
class(app)
## [1] "shiny.appobj"
# shiny::runApp(app)

The following synthetic example ranks three candidates for one feature using precomputed paired pathway expression. Real workflows can create these scores with createPathwayObject(), or let calculateAnnotationStatistics() reuse it on an alternative experiment indexed by RaMP analyte IDs.

ranking_demo <- spe
SummarizedExperiment::rowData(ranking_demo)$all_Ramp_IDs <-
  rep("RAMP_C_1; RAMP_C_2; RAMP_C_3", nrow(ranking_demo))
pathway_values <- rbind(pA = 1:6, pB = c(2, 1, 5, 3, 6, 4), pC = 6:1)
colnames(pathway_values) <- colnames(ranking_demo)
SingleCellExperiment::altExp(ranking_demo, "pathway") <-
  SingleCellExperiment::SingleCellExperiment(
    assays = list(pathwayScores = pathway_values))
rank_db <- list(
  source_df = data.frame(rampId = paste0("RAMP_C_", 1:3),
                         commonName = c("A", "B", "C")),
  analytehaspathway = data.frame(
    rampId = c("RAMP_C_1", "RAMP_C_1", "RAMP_C_2", "RAMP_C_3"),
    pathwayRampId = c("pA", "pB", "pB", "pC")))
calculateSingleAnnotationStatistics(100.01, ranking_demo, database = rank_db)
## # A tibble: 3 × 7
##   metabolite ramp_id  n_sig_path max_cor z_score   pval pval_adj
##   <chr>      <chr>         <int>   <dbl>   <dbl>  <dbl>    <dbl>
## 1 A          RAMP_C_1          2   1       1.58  0.0574    0.172
## 2 C          RAMP_C_3          0  -1      -0.423 0.664     0.876
## 3 B          RAMP_C_2          1   0.657  -1.15  0.876     0.876

This retains the historical weighted z-score formula, now calculated directly from paired Pearson correlations. Its standard-normal tail and BH-adjusted values are heuristic ranking measures, not calibrated identification p-values. Correlated pathways and spatially dependent pixels violate the simple interpretation of that tail probability. Constant expression has undefined correlation; a constant candidate-score component contributes zero.

Optional Seurat conversion boundary

SpaMTP analysis does not load Seurat or SeuratObject. Only these explicit converters need SeuratObject. Convert existing data once, then continue with native normalization, scaling, reductions and plotting.

if (requireNamespace("SeuratObject", quietly = TRUE)) {
    seurat_object <- spatialExperimentToSeurat(spe, includeAltExps = FALSE)
    spe_again <- seuratToSpatialExperiment(seurat_object)
    sce_again <- seuratToSingleCellExperiment(seurat_object)
    stopifnot(identical(colnames(spe_again), colnames(sce_again)))
}

Use seuratToSingleCellExperiment() for non-spatial transcriptomics; it never invents spatial coordinates. For spatial input, coordinates must have x/y columns and named pixels, or be supplied by a suitable Seurat image. Multiple images require an explicit selection or coordinate table. A single image/FOV is read through its public spatial accessor and takes precedence over cached metadata coordinates, which may predate alignment. The selected source is recorded in metadata(spe_again)$seurat_interoperability. Only one exact layer is transferred; split layers are not silently partially imported. Optical rasters, reductions and graphs are not transferred. Paired alternative assays are matched by pixel IDs, and unpaired assays are skipped with a warning.

For an existing Seurat MSI workflow, the native sequence is:

scaled_spe <- normalizeSMData(spe, "LogNormalize", verbose = FALSE)
scaled_spe <- scaleSMData(scaled_spe)
SummarizedExperiment::assayNames(scaled_spe)
## [1] "counts"    "logcounts" "scaled"

scaleSMData() centres/scales each feature across pixels and stores a separate scaled assay. This creates a dense matrix. Unlike the historical WNN workflow, multiOmicIntegration() combines equal-weight PCA embeddings; the algorithms must not be described as numerically equivalent.

Inspect available database resources

spaMTPDatabaseInfo() reads the lightweight SpaMTPdb registry without downloading any annotation table.

database_registry <- spaMTPDatabaseInfo(version = "latest")
database_registry[
  database_registry$resource %in% c("chem_props", "smiles_features"),
  c("resource", "version", "category", "rows", "columns")
]
##           resource version  category   rows columns
## 1       chem_props   3.0.7      core 289754      11
## 16 smiles_features   3.0.7 structure 284818      36

Create a small structure-aware database

The example contains D-2-hydroxyglutaric acid and glucose. SMILES connectivity is decomposed into interpretable functional groups and ionisation sites before candidate generation.

demo_database <- data.frame(
  ramp_id = c("RAMP_C_2HG", "RAMP_C_GLUCOSE"),
  chem_source_id = c("HMDB:HMDB0000606", "HMDB:HMDB0000122"),
  common_name = c("D-2-hydroxyglutaric acid", "D-glucose"),
  monoisotop_mass = c(148.037173, 180.063388),
  mol_formula = c("C5H8O5", "C6H12O6"),
  iso_smiles = c("O=C(O)CC(O)C(=O)O", "OCC1OC(O)C(O)C(O)C1O"),
  stringsAsFactors = FALSE
)

demo_database <- annotateSMILESStructure(demo_database)
demo_database[, c(
  "common_name", "carboxyl_sites", "hydroxyl_sites",
  "max_exchangeable_protons", "positive_mode_score",
  "negative_mode_score", "alkali_affinity_score"
)]
##                common_name carboxyl_sites hydroxyl_sites max_exchangeable_protons
## 1 D-2-hydroxyglutaric acid              2              3                        2
## 2                D-glucose              0              5                        0
##   positive_mode_score negative_mode_score alkali_affinity_score
## 1                0.52                0.94                  0.99
## 2                0.47                0.45                  0.87

Build and query an indexed search space

Candidate masses are generated once and sorted. Repeated observed peaks can then be queried using a binary range search rather than scanning the full database for every peak.

annotation_index <- buildMZAnnotationIndex(
  db = demo_database,
  polarity = "negative",
  adducts = c("[M-H]-", "[2M-H]-"),
  infer_structure = "always"
)

observed_mz <- 147.02935
candidates <- queryMZAnnotationIndex(
  observed_mz = observed_mz,
  index = annotation_index,
  ppm = 5
)

candidates[, c(
  "observed_mz", "expected_mz", "metabolite_names", "adduct",
  "ppm_error", "score"
)]
##   observed_mz expected_mz         metabolite_names adduct ppm_error      score
## 1    147.0293    147.0299 D-2-hydroxyglutaric acid    M-H  3.717158 0.07816111

Reproducibility

For production analyses, record the database version returned by spaMTPDatabaseInfo() together with ion mode, matrix profile, adduct rules, mass tolerance, and annotation-score threshold.

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] rmarkdown_2.32   SpaMTP_0.99.11   BiocStyle_2.41.0
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.6.1               later_1.4.8                 filelock_1.0.3             
##   [4] tibble_3.3.1                lifecycle_1.0.5             httr2_1.3.0                
##   [7] sf_1.1-3                    edgeR_4.99.7                globals_0.20.0             
##  [10] lattice_0.23-1              crosstalk_1.2.2             magrittr_2.0.5             
##  [13] SpaMTPdb_0.99.4             limma_3.99.2                plotly_4.12.1              
##  [16] sass_0.4.10                 jquerylib_0.1.4             yaml_2.3.12                
##  [19] metapod_1.21.0              httpuv_1.6.17               otel_0.2.0                 
##  [22] spam_2.11-4                 sp_2.2-4                    cowplot_1.2.0              
##  [25] DBI_1.3.0                   buildtools_1.0.0            RColorBrewer_1.1-3         
##  [28] abind_1.4-8                 GenomicRanges_1.65.4        purrr_1.2.2                
##  [31] BiocGenerics_0.59.12        rappdirs_0.3.4              SpaMTPData_0.99.5          
##  [34] IRanges_2.47.5              S4Vectors_0.51.10           ggrepel_0.9.8              
##  [37] irlba_2.4.1                 listenv_1.1.0               maketools_1.3.2            
##  [40] pheatmap_1.0.13             units_1.0-1                 RSpectra_0.16-2            
##  [43] matter_2.15.0               dqrng_0.4.1                 parallelly_1.48.0          
##  [46] codetools_0.2-20            DelayedArray_0.39.8         scuttle_1.23.2             
##  [49] tidyselect_1.2.1            farver_2.1.2                ScaledMatrix_1.21.0        
##  [52] viridis_0.6.5               matrixStats_1.5.0           stats4_4.6.1               
##  [55] BiocFileCache_3.3.0         Seqinfo_1.3.2               jsonlite_2.0.0             
##  [58] CardinalIO_1.11.0           BiocNeighbors_2.7.3         e1071_1.7-17               
##  [61] progressr_1.0.0             scater_1.41.2               tools_4.6.1                
##  [64] Rcpp_1.1.2                  glue_1.8.1                  gridExtra_2.3.1            
##  [67] SparseArray_1.13.4          BiocBaseUtils_1.15.1        xfun_0.61                  
##  [70] MatrixGenerics_1.25.0       dplyr_1.2.1                 withr_3.0.3                
##  [73] BiocManager_1.30.27         fastmap_1.2.0               bluster_1.23.1             
##  [76] digest_0.6.39               rsvd_1.0.5                  R6_2.6.1                   
##  [79] mime_0.13                   RSQLite_3.53.3              utf8_1.2.6                 
##  [82] tidyr_1.3.2                 generics_0.1.4              data.table_1.18.6.1        
##  [85] FNN_1.1.4.1                 class_7.3-24                httr_1.4.9                 
##  [88] htmlwidgets_1.6.4           S4Arrays_1.13.2             ontologyIndex_2.12         
##  [91] pkgconfig_2.0.3             gtable_0.3.6                blob_1.3.0                 
##  [94] S7_0.2.2                    SingleCellExperiment_1.35.2 XVector_0.53.0             
##  [97] sys_3.4.3                   htmltools_0.5.9             dotCall64_1.2              
## [100] fgsea_1.39.2                ProtGenerics_1.45.0         SeuratObject_5.4.0         
## [103] scales_1.4.0                Biobase_2.73.2              png_0.1-9                  
## [106] SpatialExperiment_1.23.0    scran_1.41.1                knitr_1.52                 
## [109] rjson_0.2.23                nlme_3.1-171                curl_8.0.0                 
## [112] proxy_0.4-29                cachem_1.1.0                stringr_1.6.0              
## [115] Cardinal_3.15.0             BiocVersion_3.24.0          KernSmooth_2.23-27         
## [118] parallel_4.6.1              vipor_0.4.7                 AnnotationDbi_1.75.2       
## [121] pillar_1.11.1               grid_4.6.1                  vctrs_0.7.3                
## [124] promises_1.5.0              BiocSingular_1.29.2         dbplyr_2.6.0               
## [127] beachmat_2.29.3             xtable_1.8-8                cluster_2.1.8.3            
## [130] beeswarm_0.4.0              evaluate_1.0.5              magick_2.9.1               
## [133] cli_3.6.6                   locfit_1.5-9.12             compiler_4.6.1             
## [136] rlang_1.3.0                 crayon_1.5.3                future.apply_1.20.2        
## [139] labeling_0.4.3              mclust_6.1.3                classInt_0.4-11            
## [142] naturalsort_0.1.3           ggbeeswarm_0.7.3            stringi_1.8.9              
## [145] viridisLite_0.4.3           BiocParallel_1.47.0         Biostrings_2.81.9          
## [148] Matrix_1.7-6                ExperimentHub_3.3.2         bit64_4.8.6                
## [151] future_1.76.0               ggplot2_4.0.3               KEGGREST_1.53.6            
## [154] statmod_1.5.2               shiny_1.14.0                SummarizedExperiment_1.43.0
## [157] AnnotationHub_4.3.2         igraph_2.3.4                memoise_2.0.1              
## [160] bslib_0.12.0                fastmatch_1.1-8             bit_4.6.0