diff --git a/CHANGELOG.md b/CHANGELOG.md index e6bdfda..1375bf3 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 93ef65f..055eac7 100644 --- a/README.md +++ b/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 | diff --git a/benchmarks/README.md b/benchmarks/README.md index b79e17f..6f75d30 100644 --- a/benchmarks/README.md +++ b/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 diff --git a/benchmarks/search_scaling.py b/benchmarks/search_scaling.py new file mode 100644 index 0000000..edc52d7 --- /dev/null +++ b/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()) diff --git a/benchmarks/search_scaling_metrics.tsv b/benchmarks/search_scaling_metrics.tsv new file mode 100644 index 0000000..64941c3 --- /dev/null +++ b/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 diff --git a/docs/assets/example-report.html b/docs/assets/example-report.html index 6600d2d..c9622b5 100644 --- a/docs/assets/example-report.html +++ b/docs/assets/example-report.html @@ -9564,7 +9564,7 @@

${escapeHtml(name)}

window.recoreExplorerRunning = true; })(); -}catch(error){window.recoreExplorerFailed(error);}