Paururo · Oct 17, 2026 · Oct 17, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,23 @@ and Recore uses semantic versioning.
 
 ### Fixed
 
+- Memory no longer grows by megabytes with every genome. Each genome's
+  annotation was held as a String per name in up to six maps, and its
+  gene-to-cluster lookup as a map of every name, copied again for the corrected
+  matrix: 3.4 MB per genome with minimal GFFs and 6.2 MB with PGAP-style ones
+  (six names per CDS), so 5,000 genomes would have needed about 17 and 32 GB.
+  Each name is now stored once, in one buffer per genome with a sorted table
+  of hashes; the lookup keeps one cluster per name and is corrected in place;
+  matrix cells are one string per cluster; and `cluster_synteny_blocks.tsv` and
+  `cluster_locus_context.tsv` are built in chunks bounded by genome-cluster
+  pairs instead of by clusters, the locus contexts now in parallel too. On 800
+  synthetic genomes peak memory went from 2.8 to 0.56 GB (5.5 to 1.0 GB with
+  PGAP-style GFFs) and run time from 33 to 22 seconds. Every output is byte for
+  byte the same.
+- When two corrections in one genome named the same gene ID, the cluster its
+  neighbours took in `cluster_locus_context.tsv` and
+  `cluster_synteny_blocks.tsv` depended on hash-map order and could change from
+  run to run. The first cluster in matrix order now wins.
 - Large runs are faster and use less memory, with identical outputs.
   `cluster_synteny_blocks.tsv` held every cluster's blocks before writing the
   first and was built in one thread: on 200 genomes of 4,000 genes, half the run

diff --git a/benchmarks/README.md b/benchmarks/README.md
@@ -149,6 +149,12 @@ equal the false absences planted). The datasets are generated afresh and
 removed afterwards unless `--keep-workdir` is given. The numbers belong to the
 machine they ran on, which the table records; the shape is what carries over.
 
+The GFFs carry an ID and a locus tag. `--rich-gff` writes PGAP-style records
+instead (a name, a gene symbol, an old locus tag, a protein ID,
+cross-references, a product and a long note), six names per CDS, which is what
+an annotation index has to hold for most real assemblies; that run goes to
+`benchmarks/search_scaling_rich_metrics.tsv`.
+
 ## AMAP performance benchmark
 
 `amap_performance_benchmark.py` profiles Recore on the full local AMAP

diff --git a/benchmarks/search_scaling.py b/benchmarks/search_scaling.py
@@ -33,6 +33,7 @@
 
 REPO_ROOT = Path(__file__).resolve().parents[1]
 DEFAULT_METRICS = REPO_ROOT / "benchmarks" / "search_scaling_metrics.tsv"
+DEFAULT_RICH_METRICS = REPO_ROOT / "benchmarks" / "search_scaling_rich_metrics.tsv"
 BASES = np.frombuffer(b"ACGT", dtype=np.uint8)
 STOPS = {b"TAA", b"TAG", b"TGA"}
 
@@ -44,7 +45,33 @@ def random_orf(rng: np.random.Generator, codons: int) -> bytes:
     return b"ATG" + b"".join(sense[:codons]) + b"TAA"
 
 
-def build_dataset(out: Path, n_genomes: int, n_genes: int, gaps: int, seed: int) -> dict:
+PRODUCTS = [
+    "DNA-directed RNA polymerase subunit beta",
+    "ATP-dependent Clp protease proteolytic subunit",
+    "PE family protein",
+    "hypothetical protein",
+    "ABC transporter ATP-binding protein",
+    "transcriptional regulator, TetR family",
+]
+
+
+def gff_attributes(genome: str, i: int, rich: bool) -> str:
+    """The CDS attributes: an ID and a locus tag, or a PGAP-style record with a
+    name, a product, a long note, cross-references and a protein ID, which is
+    what most real annotations carry and what an annotation index holds."""
+    tag = f"{genome}_{i:05d}"
+    if not rich:
+        return f"ID={tag};locus_tag={tag}"
+    return (
+        f"ID=cds-{tag};Parent=gene-{tag};Name=WP_{i:09d}.1;gene=gen{i:05d};locus_tag={tag};"
+        f"old_locus_tag=Rv{i:04d};protein_id=WP_{i:09d}.1;Dbxref=GenBank:WP_{i:09d}.1,GeneID:{1000000 + i};"
+        f"product={PRODUCTS[i % len(PRODUCTS)]};"
+        "Note=Derived by automated computational analysis using gene prediction method%3A Protein Homology. "
+        "GO_function%3A GO%3A0003677 - DNA binding [Evidence IEA]"
+    )
+
+
+def build_dataset(out: Path, n_genomes: int, n_genes: int, gaps: int, seed: int, rich_gff: bool = False) -> dict:
     """Write pa.csv, pan_ref.fa, genomes/ and annotations/ under `out`."""
     rng = np.random.default_rng(seed)
     genes = [random_orf(rng, int(rng.integers(250, 500))) for _ in range(n_genes)]
@@ -82,7 +109,7 @@ def build_dataset(out: Path, n_genomes: int, n_genes: int, gaps: int, seed: int)
             if i not in absent[g]:
                 gff_lines.append(
                     f"{genome}_1\tsynthetic\tCDS\t{pos + 1}\t{pos + len(gene)}\t.\t+\t0\t"
-                    f"ID={genome}_{i:05d};locus_tag={genome}_{i:05d}"
+                    + gff_attributes(genome, i, rich_gff)
                 )
             pos += len(gene)
         contig = bytearray(b"".join(seq_parts))
@@ -161,11 +188,19 @@ def main() -> int:
     parser.add_argument("--gaps", type=int, default=80, help="matrix gaps per genome, half of them false")
     parser.add_argument("--threads", type=int, default=8)
     parser.add_argument("--seed", type=int, default=1)
+    parser.add_argument(
+        "--rich-gff",
+        action="store_true",
+        help="PGAP-style annotations (names, products, notes, cross-references) instead of an ID and a locus tag",
+    )
     parser.add_argument("--minimap2")
     parser.add_argument("--recore", help="recore binary (default: target/release/recore, built if missing)")
     parser.add_argument("--workdir", help="where datasets and outputs go (default: a temporary directory)")
     parser.add_argument("--keep-workdir", action="store_true")
-    parser.add_argument("--metrics", default=str(DEFAULT_METRICS))
+    parser.add_argument(
+        "--metrics",
+        help="output TSV (default: search_scaling_metrics.tsv, or search_scaling_rich_metrics.tsv with --rich-gff)",
+    )
     args = parser.parse_args()
 
     minimap2 = find_minimap2(args.minimap2)
@@ -187,7 +222,7 @@ def main() -> int:
             if data.exists():
                 shutil.rmtree(data)
             t0 = time.monotonic()
-            info = build_dataset(data, n, args.genes, args.gaps, args.seed)
+            info = build_dataset(data, n, args.genes, args.gaps, args.seed, args.rich_gff)
             print(f"{n} genomes: dataset in {time.monotonic() - t0:.0f}s ({info['genome_mb']} Mb each)", flush=True)
             cache = workdir / f"cache_{n}"
             if cache.exists():
@@ -199,6 +234,7 @@ def main() -> int:
                 "genomes": n,
                 "genome_mb": info["genome_mb"],
                 "genes_per_genome": args.genes,
+                "annotation": "rich" if args.rich_gff else "minimal",
                 "gaps_searched": plain["searched"],
                 "false_absences": info["false_absences"],
                 "rescued": plain["rescued"],
@@ -223,7 +259,7 @@ def main() -> int:
         if not args.keep_workdir and not args.workdir:
             shutil.rmtree(workdir, ignore_errors=True)
 
-    metrics = Path(args.metrics)
+    metrics = Path(args.metrics) if args.metrics else DEFAULT_RICH_METRICS if args.rich_gff else DEFAULT_METRICS
     metrics.parent.mkdir(parents=True, exist_ok=True)
     with open(metrics, "w", newline="") as fh:
         writer = csv.DictWriter(fh, fieldnames=list(rows[0].keys()), delimiter="\t")

diff --git a/benchmarks/search_scaling_metrics.tsv b/benchmarks/search_scaling_metrics.tsv
@@ -1,5 +1,5 @@
-genomes	genome_mb	genes_per_genome	gaps_searched	false_absences	rescued	threads	seconds_no_cache	seconds_cold_cache	seconds_warm_cache	seconds_per_genome	peak_rss_mb	recore	minimap2	machine
-100	4.91	4000	8000	4000	3999	8	3.4	3.4	2.4	0.034	461	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
-200	4.91	4000	16000	8000	8000	8	7.1	7.9	6.1	0.035	762	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
-400	4.91	4000	32000	16000	16000	8	18.0	18.5	13.0	0.045	1456	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
-800	4.91	4000	64000	32000	32000	8	48.1	50.1	38.1	0.06	2570	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+genomes	genome_mb	genes_per_genome	annotation	gaps_searched	false_absences	rescued	threads	seconds_no_cache	seconds_cold_cache	seconds_warm_cache	seconds_per_genome	peak_rss_mb	recore	minimap2	machine
+100	4.91	4000	minimal	8000	4000	3999	8	2.4	2.5	1.3	0.024	236	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+200	4.91	4000	minimal	16000	8000	8000	8	4.7	5.1	2.8	0.024	300	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+400	4.91	4000	minimal	32000	16000	16000	8	10.7	10.5	6.3	0.027	377	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+800	4.91	4000	minimal	64000	32000	32000	8	21.4	22.1	13.5	0.027	556	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
diff --git a/benchmarks/search_scaling_rich_metrics.tsv b/benchmarks/search_scaling_rich_metrics.tsv
@@ -0,0 +1,5 @@
+genomes	genome_mb	genes_per_genome	annotation	gaps_searched	false_absences	rescued	threads	seconds_no_cache	seconds_cold_cache	seconds_warm_cache	seconds_per_genome	peak_rss_mb	recore	minimap2	machine
+100	4.91	4000	rich	8000	4000	3999	8	2.2	2.2	1.1	0.022	322	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+200	4.91	4000	rich	16000	8000	8000	8	4.2	4.3	2.3	0.021	430	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+400	4.91	4000	rich	32000	16000	16000	8	8.8	9.1	5.0	0.022	637	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
+800	4.91	4000	rich	64000	32000	32000	8	20.7	22.0	12.7	0.026	1053	recore 0.2.0	2.31-r1302	arm64 Darwin 14 cpus
diff --git a/docs/how-it-works/algorithm.md b/docs/how-it-works/algorithm.md
@@ -83,7 +83,11 @@ For each genome with a GFF file, builds a `GffIndex`:
 
 1. Genes sorted by (contig, start position)
 2. Per-contig range index for O(log n) overlap queries via binary search (`partition_point`)
-3. Gene ID index for O(1) lookup by name
+3. Every name of every CDS (its gene ID, then its aliases) stored once, back to back in one buffer, with a sorted table of name hashes for lookup by name
+4. Products and notes read once to recognise mobile elements, then dropped
+
+A CDS costs 32 bytes plus its names, so the index grows with a genome's gene
+count and the length of its names, not with how many maps refer to them.
 
 The index supports:
 - `overlapping(contig, start, end)` → all CDS features overlapping a query region
@@ -94,14 +98,19 @@ The index supports:
 
 ## Phase 4: Gene-to-Cluster Maps
 
-For each genome involved in the search, builds a HashMap from gene ID → ClusterIdx. This maps every gene ID found in the CSV cells to its Panaroo cluster assignment.
+For each genome with an annotation, builds a map from gene ID → ClusterIdx. This maps every gene ID found in the CSV cells to its Panaroo cluster assignment.
 
 Panaroo appends suffixes to gene IDs that don't exist in the original GFF:
 - `_pseudo` for pseudogene annotations
 - `_dup_N` for duplicate gene copies
 
 The map indexes both the original ID and stripped versions to maximize GFF matching success.
 
+The map is then read into one cluster for each name of the genome's
+`GffIndex` (4 bytes each, a `GeneClusters`) and dropped. Every later lookup is
+for a name the annotation uses, so this answers each one as the map did. The
+genomes are done in parallel.
+
 ---
 
 ## Phase 5: Search (minimap2)

diff --git a/docs/how-it-works/limits.md b/docs/how-it-works/limits.md
@@ -108,28 +108,36 @@ Measured with [`benchmarks/search_scaling.py`](https://github.com/PathoGenOmics-
 on an Apple M4 Pro (14 cores, 24 GB) with 8 threads and minimap2 2.31:
 synthetic MTBC-sized genomes (one 4.9 Mb contig and 4,000 genes each), 80
 matrix gaps per genome of which 40 are false absences, GFFs and `--synteny`.
+The GFFs give each CDS an ID and a locus tag; the last two columns are the
+same runs with PGAP-style GFFs (`--rich-gff`), six names per CDS plus a
+product and a note.
 
-| Genomes | Gaps searched | Run time | Warm `--search-cache` | Peak memory |
-|---------|---------------|----------|-----------------------|-------------|
-| 100 | 8,000 | 3.4 s | 2.4 s | 0.5 GB |
-| 200 | 16,000 | 7.1 s | 6.1 s | 0.8 GB |
-| 400 | 32,000 | 18.0 s | 13.0 s | 1.5 GB |
-| 800 | 64,000 | 48.1 s | 38.1 s | 2.6 GB |
+| Genomes | Gaps searched | Run time | Warm `--search-cache` | Peak memory | PGAP-style GFFs: run time | PGAP-style GFFs: peak memory |
+|---------|---------------|----------|-----------------------|-------------|---------------------------|------------------------------|
+| 100 | 8,000 | 2.4 s | 1.3 s | 0.24 GB | 2.2 s | 0.32 GB |
+| 200 | 16,000 | 4.7 s | 2.8 s | 0.30 GB | 4.2 s | 0.43 GB |
+| 400 | 32,000 | 10.7 s | 6.3 s | 0.38 GB | 8.8 s | 0.64 GB |
+| 800 | 64,000 | 21.4 s | 13.5 s | 0.56 GB | 20.7 s | 1.05 GB |
 
 Every planted false absence was recovered, except one at 100 genomes that
 synteny rejected because the random deletions had changed its neighbourhood.
 
-Memory grows by about 3 MB per genome, mostly the annotations synteny keeps for
-every genome. Time grows a little faster than the genome count, from 0.03 s
-per genome at 100 genomes to 0.06 s at 800, and the minimap2 search is the
-smaller part of it (10 of the 48 seconds at 800 genomes). The rest is synteny,
-and the audit tables in `cluster_synteny_blocks.tsv`, whose size follows the
-number of distinct gene neighbourhoods. These synthetic genomes are a worst
-case for that: every gap reshuffles the context of ten neighbours, and at 800
-genomes the table reached 2.3 GB. Assemblies with conserved synteny write far
-less.
-
-`--search-cache` removes the search from a rerun (10 seconds to 1 at 800
+Memory grows by about 0.5 MB per genome, or 1 MB with PGAP-style GFFs, nearly
+all of it the annotation synteny reads for every genome: 5,000 genomes like
+these would need about 2.5 GB, or 5.5 GB with PGAP-style GFFs. Each name is
+stored once per genome and the gene-to-cluster lookup keeps one cluster per
+name; when every name was a String in several maps, memory grew by 3.4 MB per
+genome (6.2 MB with PGAP-style GFFs). Time grows with the genome count, about
+0.027 s per genome, and the minimap2 search is about 9 of the 21 seconds at 800
+genomes. The rest is synteny, and the audit tables in
+`cluster_synteny_blocks.tsv`, whose size follows the number of distinct gene
+neighbourhoods. These synthetic genomes are a worst case for that: every gap
+reshuffles the context of ten neighbours, and at 800 genomes the table reached
+2.3 GB. Assemblies with conserved synteny write far less, and the table is
+written a bounded chunk of clusters at a time, so its size does not show in the
+memory.
+
+`--search-cache` removes the search from a rerun (21 seconds to 13.5 at 800
 genomes), so it pays off most where the search is the expensive part: many gaps
 per genome (`--min-prevalence 0`), `--multi-allele`, or larger genomes. The
 cache directory can be deleted at any time; the next run simply searches again.

diff --git a/src/parser.rs b/src/parser.rs
@@ -8,6 +8,39 @@ use crate::types::{ClusterIdx, GenomeIdx};
 
 // ── Presence/absence matrix (encapsulated) ──
 
+/// One cluster's cells, back to back in one string. A String per cell cost 24
+/// bytes plus an allocation, which with gene IDs a dozen characters long was
+/// most of what the matrix held for each genome.
+struct CellRow {
+    text: Box<str>,
+    /// Where each genome's cell ends in `text`; each starts where the one
+    /// before ends.
+    ends: Box<[u32]>,
+}
+
+impl CellRow {
+    fn new<'a>(cells: impl IntoIterator<Item = &'a str>) -> Self {
+        let mut text = String::new();
+        let mut ends = Vec::new();
+        for cell in cells {
+            text.push_str(cell);
+            ends.push(u32::try_from(text.len()).expect("a matrix row under 4 GB"));
+        }
+        CellRow {
+            text: text.into_boxed_str(),
+            ends: ends.into_boxed_slice(),
+        }
+    }
+
+    fn get(&self, gi: usize) -> &str {
+        let start = match gi {
+            0 => 0,
+            _ => self.ends[gi - 1] as usize,
+        };
+        &self.text[start..self.ends[gi] as usize]
+    }
+}
+
 /// Presence/absence matrix from Panaroo/Roary output.
 ///
 /// Fields are private to enforce consistent access and prevent
@@ -17,7 +50,7 @@ pub struct PresenceAbsence {
     genomes: Vec<String>,
     present_genomes: Vec<Vec<GenomeIdx>>,
     present_bits: Vec<Vec<bool>>,
-    cell_values: Vec<Vec<String>>,
+    cell_values: Vec<CellRow>,
     header: Vec<String>,
     metadata: Vec<Vec<String>>,
     genome_col_indices: Vec<usize>,
@@ -80,7 +113,7 @@ impl PresenceAbsence {
 
     #[inline]
     pub fn cell(&self, ci: ClusterIdx, gi: GenomeIdx) -> &str {
-        &self.cell_values[ci.val()][gi.val()]
+        self.cell_values[ci.val()].get(gi.val())
     }
 
     // ── CSV structure (for output reconstruction) ──
@@ -162,6 +195,10 @@ impl PresenceAbsence {
         let genome_col_indices: Vec<usize> = (0..n_ge).collect();
         let mut header = vec!["Gene".to_string()];
         header.extend(genomes.iter().cloned());
+        let cell_values = cell_values
+            .iter()
+            .map(|row| CellRow::new(row.iter().map(String::as_str)))
+            .collect();
         PresenceAbsence {
             clusters,
             genomes,
@@ -271,16 +308,14 @@ pub fn parse_presence_absence(path: &Path) -> Result<PresenceAbsence> {
 
         let mut present_row = Vec::new();
         let mut bit_row = vec![false; genomes.len()];
-        let mut row_values = Vec::with_capacity(genomes.len());
-
-        for (gi, &col_idx) in genome_indices.iter().enumerate() {
-            let val = record.get(col_idx).unwrap_or("").trim().to_string();
+        let row_values = CellRow::new(genome_indices.iter().enumerate().map(|(gi, &col_idx)| {
+            let val = record.get(col_idx).unwrap_or("").trim();
             if !val.is_empty() {
                 present_row.push(GenomeIdx(gi));
                 bit_row[gi] = true;
             }
-            row_values.push(val);
-        }
+            val
+        }));
 
         let meta_values: Vec<String> = meta_col_indices
             .iter()
@@ -394,15 +429,14 @@ fn parse_pirate_gene_families(path: &Path) -> Result<PresenceAbsence> {
         let field = |col: Option<usize>| col.and_then(|i| record.get(i)).unwrap_or("").to_string();
         let mut present_row = Vec::new();
         let mut bit_row = vec![false; genomes.len()];
-        let mut row_values = Vec::with_capacity(genomes.len());
-        for (gi, &col) in genome_cols.iter().enumerate() {
-            let value = record.get(col).unwrap_or("").trim().to_string();
+        let row_values = CellRow::new(genome_cols.iter().enumerate().map(|(gi, &col)| {
+            let value = record.get(col).unwrap_or("").trim();
             if !value.is_empty() {
                 present_row.push(GenomeIdx(gi));
                 bit_row[gi] = true;
             }
-            row_values.push(value);
-        }
+            value
+        }));
         metadata.push(vec![
             field(Some(family_col)),
             field(gene_name_col),
@@ -1373,13 +1407,30 @@ mod tests {
             Accessory Order,QC,Min group size nuc,Max group size nuc,\
             Avg group size nuc,g1,g2,g3\n\
             gA,gA,prod,3,3,1.0,0,0,0,,100,100,100,id1,id2,id3\n\
-            gB,gB,prod,2,2,1.0,0,0,0,,200,200,200,id4,,id5\n",
+            gB,gB,prod,2,2,1.0,0,0,0,,200,200,200,id4,,id5\n\
+            gC,gC,prod,1,2,2.0,0,0,0,,300,300,300,, id6;id7 ,\n",
         );
         let pa = parse_presence_absence(csv.path()).unwrap();
-        assert_eq!(pa.n_clusters(), 2);
+        assert_eq!(pa.n_clusters(), 3);
         assert_eq!(pa.n_genomes(), 3);
         assert_eq!(pa.core_size(), 1);
-        assert_eq!(pa.gaps().len(), 1);
+        assert_eq!(pa.gaps().len(), 3);
+        // Every cell comes back as written, trimmed, empty ones included.
+        let cells: Vec<Vec<&str>> = (0..3)
+            .map(|ci| {
+                (0..3)
+                    .map(|gi| pa.cell(ClusterIdx(ci), GenomeIdx(gi)))
+                    .collect()
+            })
+            .collect();
+        assert_eq!(
+            cells,
+            vec![
+                vec!["id1", "id2", "id3"],
+                vec!["id4", "", "id5"],
+                vec!["", "id6;id7", ""],
+            ]
+        );
     }
 
     #[test]