--- title: "Auditable gene mapping with SpaMTPdb and SpaMTPData" output: BiocStyle::html_document: toc: true vignette: > %\VignetteIndexEntry{Auditable gene mapping with SpaMTPdb and SpaMTPData} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "##", error = FALSE, cache = FALSE) ``` ## Reference ownership and reproducibility SpaMTPdb provides RaMP annotation tables and a separate versioned HGNC gene reference. SpaMTPData supplies native experiment data and declared species. SpaMTP resolves identifiers, annotates row metadata and runs pathway analyses. Use SpaMTP >= 0.99.6, SpaMTPdb >= 0.99.4 and SpaMTPData >= 0.99.5 together. ```{r registry} library(SpaMTP) SpaMTPdb::spaMTPdbGeneReference(metadata = TRUE) SpaMTPdb::spaMTPdbVersion() ``` The HGNC reference uses a fixed quarterly archive, with its own version and checksum. It does not change the existing RaMP 3.0.7 or experiment 1.1.0 files. HGNC mapping is for **Homo sapiens**; it does not convert mouse genes to human orthologues. Native experiments loaded by SpaMTPData carry their registry species in `metadata(x)$SpaMTPData`, including their nested experiments. Species checks therefore still apply after extracting RNA with `altExp()`. ## Build once, inspect before testing The following examples run during vignette builds, using a bundled 240-location subset of the published human-brain dataset. The locations nearest the coordinate-wise median were selected without consulting expression or mapping results. All 339 measured RNA features and their paired MS features are retained; see `help("vignetteData")` for the recipe and source checksums. Pinned annotation references are downloaded if not already available locally or in a verified cache. Local files and verified caches can also be used with `offline = TRUE`. ```{r reference, eval=TRUE} database <- loadSpaMTPDatabase( c("source_df", "analyte", "analytehaspathway", "pathway"), version = "3.0.7" ) reference <- SpaMTPdb::spaMTPdbGeneReference( version = "2026.7.7" ) index <- buildGeneMappingIndex(database$source_df, gene_reference = reference) mapping <- mapGeneIdentifiers( c("gene_symbol:AATK", "entrez:9625", "gene_symbol:BANK1"), index = index ) mapping ``` The index records a fingerprint of the source table's gene identifiers and common names. Pathway workflows reject a cached index built from different identity assignments, even if the RaMP ID set is unchanged. Rebuild older saved indices with `buildGeneMappingIndex()`; changing row order or duplicating identical source rows does not invalidate the fingerprint. Symbols, previous symbols and aliases are resolved to an HGNC gene. An approved symbol takes priority over other genes' aliases. Entrez, HGNC, UniProt and versioned Ensembl gene IDs are supported. A shared protein ID can be resolved by a compatible gene-specific ID in the same RaMP record. Conflicting specific IDs quarantine the record rather than joining different genes. `gene_id` is the stable HGNC identity. `ramp_ids` lists its non-conflicting RaMP records, and `canonical_ramp_id` is one representative used to count the gene once. For example, AATK's UniProt and Entrez records in RaMP 3.0.7 belong to one HGNC gene, so their pathway memberships are united before enrichment. Inspect `status`, `method`, `candidate_gene_ids` and `candidate_ramp_ids`. The latter exposes raw RaMP candidates for direct IDs and crossreferences. A shared crossreference with unresolved or conflicting competitors remains ambiguous even when only one candidate has a known HGNC identity. `not_in_ramp` means HGNC recognizes the gene but no usable RaMP record was found; `ambiguous` and `conflicting_ramp_records` require review. `mapped_ramp_only` denotes a raw RaMP record whose HGNC identity was not resolved. A warning is not evidence that the gene has no biological function or pathway membership. ## Annotate experiment features without altering expression ```{r experiment, eval=TRUE} spe <- readRDS(system.file("extdata", "human_demo.rds", package = "SpaMTP")) SingleCellExperiment::altExpNames(spe) # Select the actual transcriptome altExp shown above. If identifiers live in # a rowData column, set id_column; otherwise row names are used. spe <- annotateGeneIdentifiers( spe, assay = "SPT", index = index, organism = "Homo sapiens" ) rna <- SingleCellExperiment::altExp(spe, "SPT") head(SummarizedExperiment::rowData(rna)) table(SummarizedExperiment::rowData(rna)$spamtp_gene_status) unresolved <- SummarizedExperiment::rowData(rna)$spamtp_gene_status != "mapped" SummarizedExperiment::rowData(rna)[unresolved, c("spamtp_gene_symbol", "spamtp_gene_status", "spamtp_gene_ramp_ids")] stopifnot(ncol(rna) == ncol(spe), nrow(rna) == 339L) ``` Mapping columns remain aligned with original features. Counts, feature names, coordinates and paired modalities are preserved. If multiple expression rows map to the same gene, annotation records that fact without summing expression. Resolve those features before differential analysis. Regional pathway and network workflows reject different differential statistics for the same gene within a cluster rather than selecting an arbitrary row. For these pinned inputs, 329 genes map and 10 produce an intentional unresolved identifier warning: two have conflicting RaMP records and eight are absent from RaMP. Their original expression remains in the object, but they are excluded from pathway testing. The table above exposes those exclusions rather than silently treating the entire 339-gene panel as mapped. ## Use the whole measured panel as the universe For a runnable foreground, we declare the measured genes AATK and BANK1 as genes of interest. This is an illustrative, predeclared query, **not** a list of statistically significant genes inferred from this single section. ```{r enrichment, eval=TRUE} genes_of_interest <- c("AATK", "BANK1") measured_panel_genes <- rownames(rna) stopifnot(all(genes_of_interest %in% measured_panel_genes)) result <- fishersPathwayAnalysis( Analyte = list(genes = genes_of_interest), universe = list(genes = measured_panel_genes), database = database, gene_index = index, min_path_size = 5, pval_cutoff = NULL ) head(attr(result, "enrichment")$gene_mapping$inputs) head(result) ``` Foreground and universe use the same complete identity index. All RaMP records of a resolved gene contribute to its union of pathway memberships, but the gene counts once in K, M, U and overlap. The universe is fixed before size filtering; zero-overlap eligible pathways remain in the BH family. `pval_cutoff` filters FDR after correction. Missing identifiers are excluded with an explicit audit; do not set U to the raw panel length to hide mapping failures. Official SpaMTPdb workflows harmonize human genes automatically. A custom bundle can provide `gene_reference`, or an index can be supplied explicitly. Custom bundles without a reference retain RaMP-only behavior without a network request. `gene_mapping = "ramp"` explicitly reproduces historical mapping or uses a separately curated species-specific database; it is not an orthology conversion switch. Previously saved gene-ID mappings and p-values may change after identity reconciliation even though Fisher arithmetic is unchanged. Use the same `gene_index` in `findRegionalPathways()` and `pathwayNetworkPlots()` to keep annotation identity and provenance consistent. ## One membership index for enrichment, expression and plots `buildPathwayIndex()` joins the gene identity index to RaMP memberships once. A pathway's `database_size` counts unique identities after reconciliation and conflict exclusion; `database_raw_size` counts original RaMP records. `measured_size` and `coverage_fraction` describe the selected assay or measured universe. `used_size` and `used_members` record what actually contributes to a score or test. In regional analysis these can be smaller than the measured set because of the requested ranking filters. A pathway score is an expression summary, not evidence that the whole pathway is activated. This small example runs offline. Its fictional gene has two RaMP records with complementary pathway memberships; both contribute while the gene counts once. ```{r shared-membership-example} reference_example <- data.frame( hgnc_id = c("HGNC:1", "HGNC:2", "HGNC:3"), symbol = c("GENEA", "GENEB", "GENEC"), entrez_id = c("101", "102", "103"), uniprot_ids = c("U1", "U2", "U3") ) database_example <- list( source_df = data.frame(rampId = paste0("RAMP_G_", c(1, 2, 3, 4)), sourceId = c("entrez:101", "uniprot:U1", "entrez:102", "entrez:103")), analytehaspathway = data.frame(rampId = paste0("RAMP_G_", c(1, 3, 2, 4)), pathwayRampId = c("P1", "P1", "P2", "P2")), pathway = data.frame(pathwayRampId = c("P1", "P2"), pathwayName = c("First", "Second")), gene_reference = reference_example ) pathway_index_example <- buildPathwayIndex(database_example) example <- SpaMTPData::spaMTPExampleData() E_example <- matrix(seq_len(36), 3, dimnames = list(c("GENEA", "GENEB", "GENEC"), colnames(example))) SingleCellExperiment::altExp(example, "transcriptome") <- SingleCellExperiment::SingleCellExperiment(assays = list(logcounts = log1p(E_example))) example <- createPathwayAssay(example, analyte_type = "genes", assay = "transcriptome", slot = "logcounts", new_assay = "gene_ids", pathway_index = pathway_index_example) example <- createPathwayObject(example, assay = "gene_ids", new.assay = "gene_scores", pathway_index = pathway_index_example) SummarizedExperiment::rowData(SingleCellExperiment::altExp(example, "gene_scores"))[ , c("pathwayRampId", "database_size", "measured_size", "coverage_fraction", "used_members")] ``` For the real human subset loaded above, use the full measured RNA panel and the same index for every entry point. This chunk runs during package builds on all 339 RNA features and 240 retained observations. The [end-to-end guide](End_to_End_Workflow.html) also uses these data for descriptive regional rankings. Neither example treats locations as biological replicates. ```{r human-pathway-closure, eval=TRUE} pathway_index <- buildPathwayIndex(database = database, gene_index = index) spe <- createPathwayAssay(spe, analyte_type = "genes", assay = "SPT", slot = "counts", new_assay = "gene_ids", pathway_index = pathway_index) spe <- createPathwayObject(spe, assay = "gene_ids", new.assay = "gene_pathways", pathway_index = pathway_index, min_path_size = 5, max_path_size = 500) pathways <- SingleCellExperiment::altExp(spe, "gene_pathways") coverage <- S4Vectors::metadata(pathways)$pathway_mapping$coverage head(coverage[, c("pathwayRampId", "database_size", "database_raw_size", "measured_size", "coverage_fraction", "excluded_conflict_count")]) head(coverage$used_members) head(coverage$excluded_conflict_records) selected <- head(coverage$pathwayRampId[coverage$eligible], 2) plots <- plotPathwaysSpatially(selected, spe, assay = "SPT", slot = "counts", pathway_index = pathway_index) attr(plots, "pathway_coverage") for (plot in plots) print(plot) stopifnot(ncol(pathways) == ncol(spe), length(plots) == 2L) ``` The original assays, feature names, coordinates and paired modalities are preserved. A derived gene assay keeps its selected layer name: logcounts do not become counts. Different expression rows mapping to one gene error by default; choose `duplicate_genes = "mean"` or `"sum"` only after considering the measurement scale (summing log-expression is generally inappropriate). Identical repeated rows are deduplicated. Named pathway plots accept exact `pathwayRampId` values; ambiguous display names error instead of choosing the first match. Gene and pathway reference fingerprints are retained in outputs. Fisher's complete coverage table is `attr(result, "pathway_coverage")`; its per-row size columns also appear in the result. Pathway score objects put coverage in rowData and the complete audit in metadata. Named plots return coverage attributes and the exact plotted score vector. Unmeasured and conflict-only pathways stay in the complete audit even when no score is drawn. ```{r session-info} sessionInfo() ```