bio-expression-matrix-gene-id-mapping

What this skill does

## Version Compatibility

Reference examples tested with: biomaRt 2.58+, AnnotationDbi 1.66+, org.Hs.eg.db 3.18+, org.Mm.eg.db 3.18+, GenomicFeatures 1.54+, mygene 1.38+ (Python), pyensembl 2.3+, pandas 2.2+, rtracklayer 1.62+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Gene ID Mapping

**"Convert gene IDs from X to Y"** -> Query the appropriate annotation source (local org.db for speed, biomaRt for Ensembl-specific attributes, mygene for cross-database aliases, Ensembl REST for low-level access), with version pinning for reproducibility and explicit handling of one-to-many mappings, withdrawn symbols, and species-specific naming.

## The Single Most Important Modern Insight -- Excel autocorrect renamed real genes

Ziemann, Eren, El-Osta 2016 *Genome Biol* 17:177 scanned 18 leading genomics journals and found ~20% of papers with Excel-attached supplementary gene lists had silently mangled symbols (`SEPT2` -> `2-Sep`, `MARCH1` -> `1-Mar`, ...). Five years later the problem persisted. HGNC's response (Bruford, Braschi, Denny, Jones, Seal, Tweedie 2020 *Nat Genet* 52:754) was to rename the affected genes:

| Old | New | Affected |
|-----|-----|----------|
| `SEPT#` | `SEPTIN#` | SEPT1 - SEPT14 -> SEPTIN1 - SEPTIN14 |
| `MARCH#` | `MARCHF#` | MARCH1 - MARCH11 -> MARCHF1 - MARCHF11 |
| `MARC#` | `MTARC#` | MARC1, MARC2 -> MTARC1, MTARC2 |
| `DEC1` | `DELEC1` | DEC1 -> DELEC1 |

Code that hard-codes old symbols silently drops these genes when joined against post-2020 annotations. Detection on import: if a gene column contains `^\d{1,2}-(Jan|Feb|Mar|...|Dec)$` patterns, the file was Excel-corrupted. Always `read.csv(colClasses=c(gene='character'))` (R) or `pd.read_csv(dtype={'gene': str})` (Python) -- but the damage is at Excel-save time, not import time.

Two related insights that determine half the practical work:

1. **Ensembl version suffixes matter sometimes and not others.** `ENSG00000123456.7` is release-specific; the unversioned `ENSG00000123456` is the stable cross-release ID. STRIP for cross-release joins, MSigDB lookups, gene-set databases. KEEP for intra-release reproducibility and clinical reports. CRITICAL: the naive `sub('\\..*', '', x)` regex ALSO strips the GENCODE `_PAR_Y` suffix in releases 25-43, collapsing chrY PAR duplicates onto their chrX counterparts. Use `sub('\\.[0-9]+(_PAR_Y)?$', '\\1', x)`.

2. **Never use HGNC symbols as the primary computational key.** Symbols change. Use Ensembl or Entrez as keys; carry symbols only as display labels in the final results table.

## Algorithmic Taxonomy

| Tool | Source | Speed | Strength | Use for |
|------|--------|-------|----------|---------|
| AnnotationDbi + `org.Hs.eg.db` / `org.Mm.eg.db` | NCBI Gene snapshot, pinned at Bioc install | Fast, local | Stable, version-pinned | Default for Ensembl <-> Entrez <-> Symbol within Bioconductor |
| biomaRt | Ensembl BioMart over HTTP | Slow for >5k queries; timeouts | Ensembl-specific attributes (biotype, transcript versions, paralogs, orthologs) | Need Ensembl-specific fields; archive endpoints for release pinning |
| mygene.info / mygene (Python) | REST API to a curated meta-database | Server-side batching of 1000 IDs | Best for symbol/alias/prev_symbol resolution | Cross-database; HGNC withdrawn symbol resolution; non-R environments |
| Ensembl REST | Direct REST API to Ensembl | Rate-limited (15 req/sec) | Low-level access to variant consequence, sequence, etc. | Specialized queries not covered by biomaRt |
| pyensembl | Local Ensembl database (Python) | Fast, local, version-pinned | Reproducible offline; gene objects with transcript and exon access | Python pipelines needing rich annotation |
| HGNC API direct | https://rest.genenames.org | REST | Authoritative source for HGNC | Symbol provenance, prev/alias detection |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| R Bioconductor pipeline, Ensembl <-> Entrez <-> Symbol | `AnnotationDbi::mapIds(org.Hs.eg.db, ...)` | Fastest, version-pinned, stable |
| Need Ensembl-only attributes (biotype, paralog, ortholog) | `biomaRt::useEnsembl(version=N)` | Only biomaRt exposes these |
| Cross-database with alias and withdrawn-symbol fallback | mygene `querymany(scopes='symbol,alias,prev_symbol')` | Designed for this case |
| Python pipeline, reproducible | `pyensembl` with pinned release | Offline, version-locked |
| Clinical report needing canonical transcript per gene | MANE Select (Morales 2022 *Nature* 604:310) | Cross-database consensus (RefSeq + Ensembl) |
| Cross-species mouse <-> human | Ensembl Compara `getLDS` filtered to `one2one` | Compara has best coverage; one2one most defensible |
| Building tx2gene for tximport | `GenomicFeatures::makeTxDbFromGFF` on the SAME GTF used in quantification | Annotation pinning matters |
| Need to reproduce a 2023 analysis exactly | `useEnsembl(version=109)` (or whichever release was used) | Without `version=`, biomaRt floats to current release |
| GRCh37 (legacy clinical) | `useEnsembl(GRCh=37)` dedicated permanent endpoint | GRCh37 -> GRCh38 mappings are not 1:1 |

## AnnotationDbi + org.db

**Goal:** Map Ensembl gene IDs to symbols, Entrez IDs, or descriptions using a local Bioconductor annotation package.

**Approach:** `mapIds()` with the source keytype and target column; handle one-to-many via `multiVals`.

```r
library(org.Hs.eg.db)
library(AnnotationDbi)

ensembl_ids <- sub('\\.[0-9]+(_PAR_Y)?$', '\\1', rownames(counts))

symbols  <- mapIds(org.Hs.eg.db, keys = ensembl_ids,
                    keytype = 'ENSEMBL', column = 'SYMBOL',
                    multiVals = 'first')
entrez   <- mapIds(org.Hs.eg.db, keys = ensembl_ids,
                    keytype = 'ENSEMBL', column = 'ENTREZID',
                    multiVals = 'first')
descrips <- mapIds(org.Hs.eg.db, keys = ensembl_ids,
                    keytype = 'ENSEMBL', column = 'GENENAME',
                    multiVals = 'first')

keytypes(org.Hs.eg.db)
```

`multiVals` options: `'first'` (silent), `'asNA'` (NA for ambiguous), `'list'` (preserve all). For DE results tables, `'first'` is typical but the mapping rate should be reported.

For mouse: `org.Mm.eg.db`. For other organisms: check Bioconductor `AnnotationData -> OrgDb` list.

## biomaRt with Version Pinning

**Goal:** Query Ensembl BioMart with the EXACT release version, for reproducibility.

**Approach:** `useEnsembl(version=N)` pins; `listEnsemblArchives()` lists available archives.

```r
library(biomaRt)

ensembl <- useEnsembl(biomart = 'genes',
                       dataset = 'hsapiens_gene_ensembl',
                       version = 110)

ensembl_grch37 <- useEnsembl(biomart = 'genes',
                              dataset = 'hsapiens_gene_ensembl',
                              GRCh = 37)

mapping <- getBM(
    attributes = c('ensembl_gene_id', 'hgnc_symbol', 'entrezgene_id',
                   'gene_biotype', 'description'),
    filters    = 'ensembl_gene_id',
    values     = ensembl_ids,
    mart       = ensembl
)
```

The `filters=` argument is PLURAL. The singular `filter=` may work via R's partial matching but breaks unpredictably if another argument starts with `f`. Always spell `filters=` and `values=` fully.

Multiple filters:

```r
genes_in_region <- getBM(
    attributes = c('ensembl_gene_id', 'hgnc_symbol'),
    filters    = c('chromosome_name', 'start', 'end'),
    values     = list('16', 1100000, 1250000),
    mart       = ensembl
)
```

Without `version=`, biomaRt floats to the current release -- a script written in 2023 against Ensembl 109 produces different mappings in 2026 against Ensembl 113. ALWAYS pi