Read PPanGGOLiN, PanTA and PIRATE; write an Rtab and rescued sequences; cache searches; faster large runs

CHANGELOG.md

@@ -9,6 +9,14 @@ and Recore uses semantic versioning.
 
 ### Fixed
 
+- 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
+  and 1.2 GB. It is now built in parallel a chunk of clusters at a time and
+  written in the same order, byte for byte. Synteny validation runs per cluster
+  in parallel (a rescue only reads its own cluster's reference contexts), and
+  GFF loading and the internal-stop scan per genome. On that dataset a run went
+  from 16.0 to 7.4 seconds and its peak memory from 3.3 to 0.8 GB.
 - `evolutionary_events.tsv` no longer presents a gain or loss built on gaps
   Recore never searched as an audited call. New columns `absent_not_searched`
   and `not_searched_reason` count those absences and say why; a call none of
@@ -223,6 +231,24 @@ and Recore uses semantic versioning.
 
 ### Added
 
+- `gene_presence_absence_recore.Rtab`, the corrected matrix as 1/0 in Panaroo's
+  Rtab layout, written in the same pass as the CSV so the two cannot disagree.
+- `--write-rescued-fasta` writes `rescued_genes.fna` and `rescued_genes.faa`:
+  every cell the binary matrix gained, read in the gene's direction and
+  translated with the bacterial code from the phase its alignment fixes, under
+  the corrected GFF's feature IDs.
+- Input from PPanGGOLiN (`matrix.csv` and its family FASTA), PanTA and PIRATE
+  (`PIRATE.gene_families.tsv` and `representative_sequences.ffn`), each checked
+  on real output. PIRATE's table is read into Panaroo's layout, its fissions
+  `(a:b)` are understood, and PanTA's `sample-contig-gene` IDs now match the
+  GFF; before, synteny validated nothing on PanTA input and no rescue reached
+  the matrix.
+- `--search-cache DIR` keeps minimap2's results between runs, keyed by genome,
+  query sequence, the search settings and the minimap2 version. A rerun with
+  another policy, coverage, synteny or prevalence setting aligns only what it
+  has not seen, and writes exactly what an uncached run writes.
+- `benchmarks/search_scaling.py` and its table: run time and memory from 100 to
+  800 MTBC-sized genomes.
 - The report can pin two genomes: their rows stay in a band under the column
   header whatever the scroll, and the genome card counts the clusters where
   they differ (noting differences that rest on a gap never searched), lists

README.md

@@ -26,9 +26,15 @@ Recore searches every absent cell against the genome assemblies with minimap2,
 classifies the evidence, validates orthology against annotation and synteny when
 they are available, and writes both:
 
-- a **Panaroo-compatible corrected binary matrix** for downstream analyses
+- a **Panaroo-compatible corrected binary matrix** for downstream analyses,
+  as CSV and as the 1/0 `.Rtab` that pyseer and most R code read
 - a **multistate evidence matrix** for audit, plotting, and manuscript review
 
+It reads the output of Panaroo, Roary, PPanGGOLiN, PanTA and PIRATE, can write
+the rescued genes as nucleotide and protein FASTA for functional annotation,
+and with `--search-cache` reruns with another policy or threshold without
+aligning again.
+
 The central design choice is conservative. Recore separates *sequence evidence*
 from *clean binary presence*: full-length evidence can be corrected in the binary
 matrix, while partial homology, pseudogenes, paralog traps and
@@ -72,6 +78,7 @@ Four outputs answer most questions:
 | File | For |
 |---|---|
 | `gene_presence_absence_recore.csv` | the corrected binary matrix, for downstream tools |
+| `gene_presence_absence_recore.Rtab` | the same matrix as 1/0, for pyseer and Rtab readers |
 | `gene_presence_absence_recore_states.csv` | audit and figures, one biological state per cell |
 | `orthology_evidence.tsv` | publication wording and manual review |
 | `recore_report.html` | the whole run in one self-contained page |

benchmarks/README.md

@@ -129,6 +129,26 @@ For new runs, prefer the integrated report flag:
 recore ... --write-html-report -o recore_output
 ```
 
+## Search scaling benchmark
+
+`search_scaling.py` measures how a run grows with the number of genomes, on
+synthetic MTBC-sized datasets: one 4.9 Mb contig and 4,000 genes per genome,
+and per genome 40 false absences (gene in the assembly, not in the matrix or
+the GFF) and 40 true ones (gene gone). Each size runs three times with the
+release binary and `--synteny`: without a search cache, with a cold one and
+with the warm one.
+
+```bash
+python3 benchmarks/search_scaling.py --genomes 100 200 400 800 --threads 8 \
+  --minimap2 /path/to/minimap2
+```
+
+It writes `benchmarks/search_scaling_metrics.tsv`: wall seconds for each of the
+three runs, seconds per genome, peak memory and the rescue count (which should
+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.
+
 ## AMAP performance benchmark
 
 `amap_performance_benchmark.py` profiles Recore on the full local AMAP

benchmarks/search_scaling.py

@@ -0,0 +1,237 @@
+#!/usr/bin/env python3
+"""How long the search takes as the number of genomes grows.
+
+Generates MTBC-sized synthetic datasets (one 4.4 Mb contig and about 4,000
+genes per genome by default), with a Panaroo-style matrix in which every genome
+has some false absences (the gene is in the assembly but not in the matrix and
+not annotated) and some true ones (the gene is gone from the assembly). Each
+size is run three times with the release binary: without a search cache, with
+a cold cache and with the warm one, and wall time, peak memory and the rescue
+count are written to a TSV.
+
+The numbers are for the machine they ran on; the shape (time per genome, and
+what the cache saves) is what carries over.
+
+    python3 benchmarks/search_scaling.py --genomes 100 200 400 --threads 8
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import os
+import platform
+import re
+import shutil
+import subprocess
+import sys
+import tempfile
+import time
+from pathlib import Path
+
+import numpy as np
+
+REPO_ROOT = Path(__file__).resolve().parents[1]
+DEFAULT_METRICS = REPO_ROOT / "benchmarks" / "search_scaling_metrics.tsv"
+BASES = np.frombuffer(b"ACGT", dtype=np.uint8)
+STOPS = {b"TAA", b"TAG", b"TGA"}
+
+
+def random_orf(rng: np.random.Generator, codons: int) -> bytes:
+    """ATG, `codons` sense codons, TAA."""
+    body = BASES[rng.integers(0, 4, size=(codons + 16) * 3)].tobytes()
+    sense = [body[i:i + 3] for i in range(0, len(body), 3) if body[i:i + 3] not in STOPS]
+    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:
+    """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)]
+    names = [f"cl{i:05d}" for i in range(n_genes)]
+    spacer_len = 120
+    (out / "genomes").mkdir(parents=True, exist_ok=True)
+    (out / "annotations").mkdir(parents=True, exist_ok=True)
+    with open(out / "pan_ref.fa", "w") as fa:
+        for name, seq in zip(names, genes):
+            fa.write(f">{name}\n{seq.decode()}\n")
+
+    genomes = [f"G{g:05d}" for g in range(n_genomes)]
+    # absent[g] = genes missing from the matrix in genome g; deleted[g] = the
+    # subset also missing from the assembly. Half of each genome's gaps are
+    # false absences, the other half true ones.
+    absent: list[set[int]] = []
+    deleted: list[set[int]] = []
+    for _ in genomes:
+        chosen = rng.choice(n_genes, size=gaps, replace=False)
+        absent.append(set(int(i) for i in chosen))
+        deleted.append(set(int(i) for i in chosen[: gaps // 2]))
+
+    genome_bp = 0
+    for g, genome in enumerate(genomes):
+        seq_parts = []
+        gff_lines = ["##gff-version 3"]
+        pos = 0
+        for i, gene in enumerate(genes):
+            spacer = BASES[rng.integers(0, 4, size=spacer_len)].tobytes()
+            seq_parts.append(spacer)
+            pos += spacer_len
+            if i in deleted[g]:
+                continue
+            seq_parts.append(gene)
+            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}"
+                )
+            pos += len(gene)
+        contig = bytearray(b"".join(seq_parts))
+        # About one difference in a thousand, so alignments are not trivial.
+        sites = rng.integers(0, len(contig), size=len(contig) // 1000)
+        for s in sites:
+            contig[s] = BASES[(BASES.tolist().index(contig[s]) + 1) % 4]
+        genome_bp = len(contig)
+        with open(out / "genomes" / f"{genome}.fasta", "w") as fh:
+            fh.write(f">{genome}_1\n")
+            text = contig.decode()
+            for k in range(0, len(text), 80):
+                fh.write(text[k:k + 80] + "\n")
+        (out / "annotations" / f"{genome}.gff").write_text("\n".join(gff_lines) + "\n")
+
+    header = (
+        "Gene,Non-unique Gene name,Annotation,No. isolates,No. sequences,"
+        "Avg sequences per isolate,Genome Fragment,Order within Fragment,"
+        "Accessory Order,QC,Min group size nuc,Max group size nuc,Avg group size nuc,"
+        + ",".join(genomes)
+    )
+    with open(out / "pa.csv", "w") as fh:
+        fh.write(header + "\n")
+        for i, name in enumerate(names):
+            cells = ["" if i in absent[g] else f"{genomes[g]}_{i:05d}" for g in range(n_genomes)]
+            n = sum(1 for c in cells if c)
+            fh.write(f"{name},{name},synthetic,{n},{n},1,0,0,0,,0,0,0," + ",".join(cells) + "\n")
+    return {
+        "genome_mb": round(genome_bp / 1e6, 2),
+        "false_absences": sum(len(a) - len(d) for a, d in zip(absent, deleted)),
+        "true_absences": sum(len(d) for d in deleted),
+    }
+
+
+def find_minimap2(explicit: str | None) -> Path:
+    for candidate in [explicit, shutil.which("minimap2")]:
+        if candidate and Path(candidate).exists():
+            return Path(candidate)
+    raise SystemExit("ERROR: minimap2 not found; pass --minimap2")
+
+
+def run_recore(binary: Path, data: Path, out: Path, threads: int, cache: Path | None, env: dict) -> dict:
+    """One run under /usr/bin/time; wall seconds, peak RSS and rescues."""
+    cmd = [
+        str(binary), "-p", str(data / "pa.csv"), "-r", str(data / "pan_ref.fa"),
+        "-g", str(data / "genomes"), "-a", str(data / "annotations"), "--synteny",
+        "-t", str(threads), "-o", str(out),
+    ]
+    if cache is not None:
+        cmd += ["--search-cache", str(cache)]
+    timer = ["/usr/bin/time", "-l"] if platform.system() == "Darwin" else ["/usr/bin/time", "-v"]
+    start = time.monotonic()
+    proc = subprocess.run(timer + cmd, env=env, capture_output=True, text=True)
+    seconds = time.monotonic() - start
+    if proc.returncode != 0:
+        sys.stderr.write(proc.stderr[-3000:])
+        raise SystemExit(f"ERROR: recore failed on {data}")
+    rss = 0
+    m = re.search(r"(\d+)\s+maximum resident set size", proc.stderr)
+    if m:
+        rss = int(m.group(1)) // (1024 * 1024)
+    else:
+        m = re.search(r"Maximum resident set size \(kbytes\): (\d+)", proc.stderr)
+        if m:
+            rss = int(m.group(1)) // 1024
+    summary = (out / "summary.txt").read_text()
+    rescued = int(re.search(r"Rescued:\s+(\d+)", summary).group(1))
+    searched = int(re.search(r"Gaps searched:\s+(\d+)", summary).group(1))
+    return {"seconds": round(seconds, 1), "peak_rss_mb": rss, "rescued": rescued, "searched": searched}
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument("--genomes", type=int, nargs="+", default=[50, 100, 200])
+    parser.add_argument("--genes", type=int, default=4000, help="genes per genome")
+    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("--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))
+    args = parser.parse_args()
+
+    minimap2 = find_minimap2(args.minimap2)
+    env = os.environ.copy()
+    env["PATH"] = f"{minimap2.parent}{os.pathsep}{env.get('PATH', '')}"
+    binary = Path(args.recore) if args.recore else REPO_ROOT / "target" / "release" / "recore"
+    if not binary.exists():
+        subprocess.run(["cargo", "build", "--release"], cwd=REPO_ROOT, check=True)
+    version = subprocess.run([str(binary), "--version"], capture_output=True, text=True).stdout.strip()
+    mm2_version = subprocess.run([str(minimap2), "--version"], capture_output=True, text=True).stdout.strip()
+
+    workdir = Path(args.workdir) if args.workdir else Path(tempfile.mkdtemp(prefix="recore_scaling_"))
+    workdir.mkdir(parents=True, exist_ok=True)
+    machine = f"{platform.machine()} {platform.system()} {os.cpu_count()} cpus"
+    rows = []
+    try:
+        for n in args.genomes:
+            data = workdir / f"data_{n}"
+            if data.exists():
+                shutil.rmtree(data)
+            t0 = time.monotonic()
+            info = build_dataset(data, n, args.genes, args.gaps, args.seed)
+            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():
+                shutil.rmtree(cache)
+            plain = run_recore(binary, data, workdir / f"out_{n}_plain", args.threads, None, env)
+            cold = run_recore(binary, data, workdir / f"out_{n}_cold", args.threads, cache, env)
+            warm = run_recore(binary, data, workdir / f"out_{n}_warm", args.threads, cache, env)
+            row = {
+                "genomes": n,
+                "genome_mb": info["genome_mb"],
+                "genes_per_genome": args.genes,
+                "gaps_searched": plain["searched"],
+                "false_absences": info["false_absences"],
+                "rescued": plain["rescued"],
+                "threads": args.threads,
+                "seconds_no_cache": plain["seconds"],
+                "seconds_cold_cache": cold["seconds"],
+                "seconds_warm_cache": warm["seconds"],
+                "seconds_per_genome": round(plain["seconds"] / n, 3),
+                "peak_rss_mb": max(plain["peak_rss_mb"], cold["peak_rss_mb"], warm["peak_rss_mb"]),
+                "recore": version,
+                "minimap2": mm2_version,
+                "machine": machine,
+            }
+            rows.append(row)
+            print("  " + ", ".join(f"{k}={v}" for k, v in row.items() if k not in {"recore", "minimap2", "machine"}), flush=True)
+            if not args.keep_workdir:
+                shutil.rmtree(data, ignore_errors=True)
+                for suffix in ("plain", "cold", "warm"):
+                    shutil.rmtree(workdir / f"out_{n}_{suffix}", ignore_errors=True)
+                shutil.rmtree(cache, ignore_errors=True)
+    finally:
+        if not args.keep_workdir and not args.workdir:
+            shutil.rmtree(workdir, ignore_errors=True)
+
+    metrics = Path(args.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")
+        writer.writeheader()
+        writer.writerows(rows)
+    print(f"Wrote {metrics}")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

benchmarks/search_scaling_metrics.tsv

@@ -0,0 +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

docs/getting-started/inputs.md

@@ -2,14 +2,35 @@
 
 | Input | Required | Description |
 |-------|----------|-------------|
-| `gene_presence_absence.csv` | yes | Panaroo-style matrix: metadata columns followed by one column per genome |
+| `gene_presence_absence.csv` | yes | Panaroo-style matrix: metadata columns followed by one column per genome. Roary, PPanGGOLiN and PanTA write this layout, and PIRATE's own table is read too (see below) |
 | `pan_genome_reference.fa` | yes | Cluster representative sequences; FASTA IDs must match matrix cluster names |
 | genome FASTAs | yes | One assembly per matrix genome column |
 | GFF annotations | optional, recommended | Enables annotation-aware `SPLIT`, `START_VARIANT`, `PARALOG`, pseudogene and synteny classification |
 | reference genome name | optional | Used by `--reference-genome` to audit whether clusters are absent from a chosen reference |
 | phylogenetic tree | optional | Newick, for `--tree`: parsimony on the corrected state, and the phylogeny pane of the HTML report |
 | mobile-element FASTA | optional | For `--mobile-element-reference`: sequence-based IS support instead of product names |
 
+## Which files each pangenome tool gives you
+
+Pass `-p` the matrix and `-r` the representative sequences, and `-a` the same
+GFF files the tool was run on, so the gene IDs in the matrix cells are the ones
+in the annotations. Each pair was checked on real output of the version shown.
+
+| Tool | `-p` | `-r` |
+|------|------|------|
+| Panaroo | `gene_presence_absence.csv` | `pan_genome_reference.fa` |
+| Roary | `gene_presence_absence.csv` | `pan_genome_reference.fa` |
+| PPanGGOLiN (2.3.2) | `matrix.csv`, from `ppanggolin all` or `ppanggolin write_pangenome --csv` | `all_nucleotide_families.fasta`, from `ppanggolin fasta --gene_families all` |
+| PanTA (1.0.1) | `gene_presence_absence.csv` | `representative_clusters_nucl.fasta` |
+| PIRATE (1.0.5) | `PIRATE.gene_families.tsv` | `representative_sequences.ffn` |
+
+For PIRATE, give recore its own table, not the output of `PIRATE_to_roary.pl`:
+that script names every row `allele__family`, which matches no sequence in
+PIRATE's representative FASTA. The table is read into Panaroo's layout, so the
+corrected files are the same as for a Panaroo run, and a fission written
+`(locus1:locus2)` counts its pieces as one gene. PanTA writes its gene IDs as
+`sample-contig-gene`; recore matches them to the GFF by the gene part.
+
 ## Name matching
 
 This is where most failed runs come from.

docs/getting-started/quickstart.md

@@ -56,6 +56,7 @@ change what a run means, rather than how fast it is:
 | Read this | For |
 |-----------|-----|
 | `gene_presence_absence_recore.csv` | the corrected binary matrix, for downstream tools |
+| `gene_presence_absence_recore.Rtab` | the same matrix as 1/0, for pyseer and Rtab readers |
 | `gene_presence_absence_recore_states.csv` | audit and figures, one biological state per cell |
 | `orthology_evidence.tsv` | publication wording and manual review |
 | `recore_report.html` | the whole run in one page, offline |