From a3610a51713709aff3d8164d739cf7e2a41999c9 Mon Sep 17 00:00:00 2001 From: Paula Ruiz-Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Sat, 17 Oct 2026 19:06:41 +0200 Subject: [PATCH 1/2] feat: read PPanGGOLiN, PanTA and PIRATE, write an Rtab and rescued sequences, and cache searches Recore now reads PPanGGOLiN output, PanTA output (whose sample-contig-gene IDs never matched a GFF, so synteny validated nothing) and PIRATE's own gene family table, writes the corrected matrix as an Rtab, and with --write-rescued-fasta the rescued genes as nucleotide and protein FASTA. --search-cache keeps minimap2's results between runs, with byte-identical outputs. Measuring at scale showed the search is the smaller part of a large run, so the synteny block table, synteny validation, GFF loading and the stop scan now run in parallel: 200 genomes went from 16 to 7 seconds and from 3.3 to 0.8 GB, with identical outputs. --- CHANGELOG.md | 26 ++ README.md | 9 +- benchmarks/README.md | 20 + benchmarks/search_scaling.py | 237 ++++++++++ benchmarks/search_scaling_metrics.tsv | 5 + docs/getting-started/inputs.md | 23 +- docs/getting-started/quickstart.md | 1 + docs/guide/outputs.md | 34 ++ docs/guide/parameters.md | 2 + docs/how-it-works/limits.md | 32 ++ docs/index.md | 1 + src/cli.rs | 14 + src/parser.rs | 221 ++++++++- src/pipeline.rs | 455 +++++++++++------- src/pipeline/evolution.rs | 269 ++++++----- src/pipeline/output.rs | 165 +++++++ src/pipeline/scale_harness.rs | 2 + src/pipeline/tests.rs | 90 ++++ src/relationships.rs | 4 +- src/rescue.rs | 6 +- src/rescue/genes.rs | 33 +- src/rescue/orf.rs | 45 +- src/rescue/synteny.rs | 6 +- src/rescue/tests.rs | 32 ++ src/search.rs | 650 +++++++++++++++++++++++--- tests/integration.rs | 155 ++++++ 26 files changed, 2177 insertions(+), 360 deletions(-) create mode 100644 benchmarks/search_scaling.py create mode 100644 benchmarks/search_scaling_metrics.tsv 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/getting-started/inputs.md b/docs/getting-started/inputs.md index f236ab3..2c516aa 100644 --- a/docs/getting-started/inputs.md +++ b/docs/getting-started/inputs.md @@ -2,7 +2,7 @@ | 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 | @@ -10,6 +10,27 @@ | 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. diff --git a/docs/getting-started/quickstart.md b/docs/getting-started/quickstart.md index 832de15..4ac14b7 100644 --- a/docs/getting-started/quickstart.md +++ b/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 | diff --git a/docs/guide/outputs.md b/docs/guide/outputs.md index 1f24ef8..885f216 100644 --- a/docs/guide/outputs.md +++ b/docs/guide/outputs.md @@ -332,6 +332,18 @@ This file can be used directly as input to downstream tools that accept Panaroo' --- +## gene_presence_absence_recore.Rtab + +Written for every run, in the same pass as `gene_presence_absence_recore.csv`, +so the two agree on every cell. It is the corrected matrix as 1 and 0 in +Panaroo's `gene_presence_absence.Rtab` layout: a `Gene` column, then one column +per genome, tab-separated. A cell is 1 when the corrected CSV has anything in +it (input presence, `RESCUED:...` or `NSC_MERGED`), and clusters merged away by +NSC are left out, as in the CSV. This is the file pyseer's `--pres`, and most R +code written for Panaroo, reads. + +--- + ## gene_presence_absence_recore_states.csv Written for every run. Same row/column layout as the corrected matrix, but each @@ -602,6 +614,28 @@ feature records the rescue evidence explicitly. --- +## rescued_genes.fna / rescued_genes.faa + +Written only when `--write-rescued-fasta` is enabled: the sequence of every +cell the binary matrix gained, for functional annotation (eggNOG-mapper, +InterProScan, BLAST). Record IDs are the feature IDs of +`corrected_gff/*.recore.gff3`, and the description carries the genome, the +cluster, the locus, the rescue type and the phase: + +``` +>recore_Rv0003_L1_AS00101_contig1_3333_4517 genome=L1_AS00101 gene_cluster=Rv0003 locus=L1_AS00101_contig1:3333-4517(+) rescue_type=DROPPED phase=0 +``` + +`rescued_genes.fna` has the nucleotide span read in the gene's direction +(reverse-complemented on the minus strand). `rescued_genes.faa` translates it +from the phase the alignment fixes, with the bacterial code (NCBI table 11); a +first codon that is the gene's own start reads `M` whichever start it is, and +the final stop is left off. A rescue built from several alignments has no +single reading frame, so it is in the `.fna` with `phase=unknown` and not in +the `.faa`. + +--- + ## recore_report.html Written only when `--write-html-report` is enabled. diff --git a/docs/guide/parameters.md b/docs/guide/parameters.md index 7897249..7e05647 100644 --- a/docs/guide/parameters.md +++ b/docs/guide/parameters.md @@ -47,6 +47,8 @@ NSC k-mer size must be at least 1, and `--nsc-engine` must be `minhash` or |------|---------|-------------| | `--synteny` | `false` | Validate each rescue by adaptive 1-to-5-neighbour flanking context against reference genomes. Requires `--annotations`. Boosts synteny-consistent rescues (×1.2), reports quantitative synteny metrics, and rejects explicit synteny failures from the corrected matrix while writing them to `synteny_rejected_report.tsv`. Independently, PARALOG candidates without explicit synteny support are rejected to `paralog_rejected_report.tsv`. | | `--write-corrected-gff` | `false` | Write a `corrected_gff/*.recore.gff3` file for each genome with a matrix-corrected rescue, with those CDS features appended. Original GFF records are preserved when `--annotations` is provided. | +| `--write-rescued-fasta` | `false` | Write `rescued_genes.fna` and `rescued_genes.faa`: the sequence of every cell the binary matrix gained, read in the gene's direction, and its protein (bacterial code, table 11) where the reading frame is known. IDs match the features of `--write-corrected-gff`. | +| `--search-cache` | none | Directory to keep minimap2 results in and reuse them from. A rerun that changes the matrix policy, coverage, synteny or the prevalence floor aligns only the sequences it has not seen before; a changed genome, query sequence, preset, `--min-identity` or minimap2 version is searched afresh. Results are identical with and without it. | | `--write-html-report` | `false` | Write a self-contained `recore_report.html` with headline QC, correction counts, evolutionary mechanisms, top candidate clusters, and links to the underlying audit files. | | `--matrix-policy` | `balanced` | Controls which accepted rescues can enter the binary corrected matrix. `strict` corrects only `STRICT_PRESENT` and `ALIAS_CLUSTER`; `balanced` also corrects full-length `SEQUENCE_ONLY`; `sensitive` also corrects partial/disrupted accepted evidence for exploratory analyses. In all modes, the adaptive orthology-quality guard can keep weak accepted evidence as `audit_only_quality`. | | `--reference-genome` | none | Audit one genome column/reference background by exact matrix name or embedded `ASxxxxx` identifier. Runs a reference-only shadow search and writes `reference_genome_audit.tsv` with `new_to_reference`, `reference_like_or_present`, and `reference_like_ambiguous` classifications, plus `publication_confidence` labels for clean, complex/paralogous, event-derived, fragmented, pseudogene, and unresolved reference evidence. Also writes `reference_cluster_events.tsv` for reference-absent clusters whose representative sequence is best explained as an allele, fragment, or partial homologue of a reference gene. These rows are audit-only and do not change the binary matrix. | diff --git a/docs/how-it-works/limits.md b/docs/how-it-works/limits.md index 6fd58c0..f2dc93e 100644 --- a/docs/how-it-works/limits.md +++ b/docs/how-it-works/limits.md @@ -102,6 +102,38 @@ The synthetic benchmark validates policy behavior on matrix ground truth; the biological benchmark validates the locus-context guardrails. Neither is a species-specific biological truth benchmark. +## Run time and memory + +Measured with [`benchmarks/search_scaling.py`](https://github.com/PathoGenOmics-Lab/recore/blob/main/benchmarks/search_scaling.py) +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`. + +| 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 | + +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 +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. + ## Reproducible review Three Python utilities ship in `benchmarks/` for checking a run rather than diff --git a/docs/index.md b/docs/index.md index 2adb64d..fd7c0af 100644 --- a/docs/index.md +++ b/docs/index.md @@ -86,6 +86,7 @@ recore \ | Read this first | 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` | a self-contained viewer for the whole run | diff --git a/src/cli.rs b/src/cli.rs index 8ea5049..3ffaa17 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -72,6 +72,20 @@ pub struct Args { #[arg(long, default_value_t = false)] pub write_corrected_gff: bool, + /// Write rescued_genes.fna and rescued_genes.faa: the sequence of every cell + /// the binary matrix gained, for functional annotation. The IDs are those of + /// the features --write-corrected-gff writes. + #[arg(long, default_value_t = false)] + pub write_rescued_fasta: bool, + + /// Keep minimap2 results in this directory and reuse them on later runs. A + /// rerun that changes the matrix policy, coverage, synteny or prevalence + /// settings searches only what was not searched before; a changed genome, + /// query sequence, preset, --min-identity or minimap2 version is searched + /// afresh. + #[arg(long, value_name = "DIR")] + pub search_cache: Option, + /// Write a self-contained HTML run report. #[arg(long, default_value_t = false)] pub write_html_report: bool, diff --git a/src/parser.rs b/src/parser.rs index 1387a25..ff58710 100644 --- a/src/parser.rs +++ b/src/parser.rs @@ -135,7 +135,7 @@ impl PresenceAbsence { .count() } - /// Constructor for tests only — minimal PA with all genes absent. + /// Constructor for tests only: a minimal matrix with all genes absent. #[cfg(test)] pub fn new_for_test( clusters: Vec, @@ -227,7 +227,13 @@ const METADATA_COLS: &[&str] = &[ "Avg group size nuc", ]; +/// The first header field of PIRATE's own gene family table. +const PIRATE_FIRST_COLUMNS: &str = "allele_name\tgene_family\t"; + pub fn parse_presence_absence(path: &Path) -> Result { + if first_line_starts_with(path, PIRATE_FIRST_COLUMNS)? { + return parse_pirate_gene_families(path); + } let mut reader = csv::Reader::from_path(path).with_context(|| format!("Cannot open {}", path.display()))?; @@ -314,6 +320,119 @@ pub fn parse_presence_absence(path: &Path) -> Result { }) } +fn first_line_starts_with(path: &Path, prefix: &str) -> Result { + use std::io::{BufRead, BufReader}; + let file = fs::File::open(path).with_context(|| format!("Cannot open {}", path.display()))?; + let mut first = Vec::new(); + BufReader::new(file) + .read_until(b'\n', &mut first) + .with_context(|| format!("Cannot read {}", path.display()))?; + Ok(first.starts_with(prefix.as_bytes())) +} + +/// PIRATE's gene family table, `PIRATE.gene_families.tsv`. +/// +/// Tab-separated, one row per gene family: its metadata, then one column per +/// genome after `average_length(bp)`. A cell holds the family's loci in that +/// genome, `;` between copies and `(a:b)` around the pieces of a fission. The +/// table is read into Panaroo's layout, each family under its `gene_family` +/// name, so a PIRATE run writes the same files a Panaroo run does. Converting +/// it with PIRATE_to_roary.pl instead names every row `allele__family`, which +/// matches no sequence in PIRATE's representative FASTA. +fn parse_pirate_gene_families(path: &Path) -> Result { + let mut reader = csv::ReaderBuilder::new() + .delimiter(b'\t') + .quoting(false) + .flexible(true) + .from_path(path) + .with_context(|| format!("Cannot open {}", path.display()))?; + let headers: Vec = reader.headers()?.iter().map(str::to_string).collect(); + let column = |name: &str| headers.iter().position(|h| h == name); + let family_col = column("gene_family").with_context(|| { + format!( + "{}: PIRATE table without a gene_family column", + path.display() + ) + })?; + let last_metadata = column("average_length(bp)").with_context(|| { + format!( + "{}: PIRATE table without an average_length(bp) column, so where the genome columns start is unknown", + path.display() + ) + })?; + let gene_name_col = column("consensus_gene_name"); + let product_col = column("consensus_product"); + let loci_col = column("no_loci"); + let dose_col = column("average_dose"); + let genome_cols: Vec = (last_metadata + 1..headers.len()).collect(); + let genomes: Vec = genome_cols.iter().map(|&i| headers[i].clone()).collect(); + if genomes.is_empty() { + anyhow::bail!("{}: PIRATE table has no genome columns", path.display()); + } + const LAYOUT: [&str; 6] = [ + "Gene", + "Non-unique Gene name", + "Annotation", + "No. isolates", + "No. sequences", + "Avg sequences per isolate", + ]; + let header: Vec = LAYOUT + .iter() + .map(|name| name.to_string()) + .chain(genomes.iter().cloned()) + .collect(); + let genome_col_indices: Vec = (LAYOUT.len()..header.len()).collect(); + + let mut clusters = Vec::new(); + let mut present_genomes = Vec::new(); + let mut present_bits = Vec::new(); + let mut cell_values = Vec::new(); + let mut metadata = Vec::new(); + for result in reader.records() { + let record = result.with_context(|| format!("Cannot read {}", path.display()))?; + let field = |col: Option| 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(); + if !value.is_empty() { + present_row.push(GenomeIdx(gi)); + bit_row[gi] = true; + } + row_values.push(value); + } + metadata.push(vec![ + field(Some(family_col)), + field(gene_name_col), + field(product_col), + present_row.len().to_string(), + field(loci_col), + field(dose_col), + ]); + clusters.push(field(Some(family_col))); + present_genomes.push(present_row); + present_bits.push(bit_row); + cell_values.push(row_values); + } + info!( + "Parsed PIRATE gene families: {} families × {} genomes, read into Panaroo's layout", + clusters.len(), + genomes.len() + ); + Ok(PresenceAbsence { + clusters, + genomes, + present_genomes, + present_bits, + cell_values, + header, + metadata, + genome_col_indices, + }) +} + // ── FASTA / GFF parsing ── /// Parse a reference FASTA keyed by the first token of each header. @@ -331,13 +450,15 @@ pub fn parse_pan_reference( /// Parse the pan-genome reference, naming each record after the matrix cluster /// it belongs to. /// -/// Three header layouts are in use, and each record is resolved against the +/// Four header layouts are in use, and each record is resolved against the /// matrix's own cluster names in this order: /// /// 1. `>cluster ...`, Panaroo's layout: the first token is the cluster. /// 2. `>seqid cluster`, Roary's layout: the second token is the cluster. This /// used to match nothing, and a Roary run exited 0 with no rescues at all. -/// 3. `>cluster~~~ALLELE~n`: allele naming, where the cluster is the part before +/// 3. `>family;key=value;...`, PIRATE's representative sequences: the family is +/// the part of the first token before `;`. +/// 4. `>cluster~~~ALLELE~n`: allele naming, where the cluster is the part before /// `~~~`. Only as a fallback: Panaroo writes `~~~` inside real cluster names /// when it merges genes (`pstB~~~phoT`), and cutting every name at `~~~` /// either lost those clusters or, worse, filed them under another cluster @@ -367,6 +488,12 @@ fn resolve_reference_name(header: &str, clusters: Option<&HashSet<&str>>) -> Str return second.to_string(); } } + // PIRATE: `>family;representative_genome=...;locus_tag=...`. + if let Some((family, _)) = first.split_once(';') { + if clusters.contains(family) { + return family.to_string(); + } + } if let Some(pos) = first.find("~~~") { let prefix = &first[..pos]; if clusters.contains(prefix) { @@ -1685,4 +1812,92 @@ mod tests { assert_eq!(contigs["ctg1"], b"ATCGATCG"); assert_eq!(contigs["ctg2"], b"TTTTCCCC"); } + + /// PIRATE's own gene family table, a verbatim excerpt of PIRATE 1.0.5 output + /// on the tutorial data, with one cell rewritten in its fission notation + /// `(a:b)` (link_clusters.pl) and one emptied. It is read into Panaroo's + /// layout, each row under its gene_family name. + #[test] + fn pirate_gene_family_table_is_read_into_panaroo_layout() { + let table = [ + "allele_name gene_family consensus_gene_name consensus_product threshold alleles_at_maximum_threshold number_genomes average_dose min_dose max_dose genomes_containing_fissions genomes_containing_duplications number_fission_loci number_duplicated_loci no_loci products gene_names min_length(bp) max_length(bp) average_length(bp) L1_AS00101 L1_AS00102 L2_AS00201 L2_AS00202 L2_AS00203 L3_AS00301 L4_AS00401 L4_AS00402 L4_AS00403 L4_AS00404 L5_AS00501 L8_AS00801", + "g000000011_000000005 g000000011 Rv0005 DNA gyrase subunit B 95 3 12 1 1 1 0 0 0 0 12 DNA gyrase subunit B(12) Rv0005(12) 603 603 603.00 L1_AS00101_00005 L1_AS00102_00005 L2_AS00201_00005 L2_AS00202_00005 L2_AS00203_00005 L3_AS00301_00005 L4_AS00401_00005 L4_AS00402_00005 L4_AS00403_00005 L4_AS00404_00005 L5_AS00501_00005 L8_AS00801_00005", + "g000000003_000000006 g000000003 Rv0013 transcription elongation factor GreA 98 1 12 1 1 1 0 0 0 0 12 transcription elongation factor GreA(12) Rv0013(12) 1164 1164 1164.00 L1_AS00101_00012 L1_AS00102_00013 L2_AS00201_00013 L2_AS00202_00013 L2_AS00203_00013 L3_AS00301_00012 L4_AS00401_00013 L4_AS00402_00013 L4_AS00403_00013 L4_AS00404_00013 (L8_AS00801_00013a:L8_AS00801_00013b)", + ] + .join("\n"); + let f = write_temp(&format!("{table}\n")); + let pa = parse_presence_absence(f.path()).unwrap(); + assert_eq!( + pa.cluster_names(), + &["g000000011".to_string(), "g000000003".to_string()] + ); + assert_eq!(pa.n_genomes(), 12); + assert_eq!(pa.genome_name(GenomeIdx(0)), "L1_AS00101"); + assert_eq!(pa.genome_name(GenomeIdx(11)), "L8_AS00801"); + assert_eq!( + &pa.header()[..7], + &[ + "Gene", + "Non-unique Gene name", + "Annotation", + "No. isolates", + "No. sequences", + "Avg sequences per isolate", + "L1_AS00101" + ] + ); + assert_eq!( + pa.metadata(ClusterIdx(0))[..4], + ["g000000011", "Rv0005", "DNA gyrase subunit B", "12"] + ); + // Cells keep PIRATE's loci, the fission included; the emptied one is absent. + assert_eq!( + pa.cell(ClusterIdx(1), GenomeIdx(11)), + "(L8_AS00801_00013a:L8_AS00801_00013b)" + ); + assert!(!pa.is_present(ClusterIdx(1), GenomeIdx(10))); + assert_eq!(pa.present_count(ClusterIdx(1)), 11); + } + + /// PIRATE's representative_sequences.ffn headers, verbatim: the family is the + /// first token up to its first `;`. + #[test] + fn pirate_representative_headers_name_their_family() { + let f = write_temp( + ">g000000001;representative_genome=L1_AS00101;locus_tag=L1_AS00101_00001;gene_name=Rv0001;gene_product=chromosomal replication initiator DnaA;number_genomes=12\nATGAAACCC\n\ + >g000000002;representative_genome=L1_AS00101;locus_tag=L1_AS00101_00004;gene_name=Rv0004;gene_product=hypothetical protein;number_genomes=12\nATGGGG\n", + ); + let clusters: HashSet<&str> = ["g000000001", "g000000002"].into_iter().collect(); + let refs = parse_pan_reference_for(f.path(), false, &clusters).unwrap(); + assert_eq!(refs["g000000001"], vec![b"ATGAAACCC".to_vec()]); + assert_eq!(refs["g000000002"], vec![b"ATGGGG".to_vec()]); + } + + /// PanTA's gene_presence_absence.csv (PanTA 1.0.1) and PPanGGOLiN's + /// matrix.csv (2.3.2), verbatim: eight and fourteen metadata columns, every + /// field quoted. Both already read as Panaroo-style matrices. + #[test] + fn panta_and_ppanggolin_matrices_read_as_panaroo_style() { + for (lines, first_genome, cluster, cell) in [ + ( + ["\"Gene\",\"Annotation\",\"No. isolates\",\"No. sequences\",\"Avg sequences per isolate\",\"Min group size nuc\",\"Max group size nuc\",\"Avg group size nuc\",\"L1_AS00101\",\"L1_AS00102\",\"L2_AS00201\",\"L2_AS00202\",\"L2_AS00203\",\"L3_AS00301\",\"L4_AS00401\",\"L4_AS00402\",\"L4_AS00403\",\"L4_AS00404\",\"L5_AS00501\",\"L8_AS00801\"", "\"Rv0019c\",\"FHA-domain protein\",\"12\",\"12\",\"1.0\",\"1281\",\"1281\",\"1281.0\",\"L1_AS00101-L1_AS00101_contig1-L1_AS00101_00019\",\"L1_AS00102-L1_AS00102_contig1-L1_AS00102_00019\",\"L2_AS00201-L2_AS00201_contig1-L2_AS00201_00019\",\"L2_AS00202-L2_AS00202_contig1-L2_AS00202_00019\",\"L2_AS00203-L2_AS00203_contig1-L2_AS00203_00019\",\"L3_AS00301-L3_AS00301_contig1-L3_AS00301_00019\",\"L4_AS00401-L4_AS00401_contig1-L4_AS00401_00019\",\"L4_AS00402-L4_AS00402_contig1-L4_AS00402_00019\",\"L4_AS00403-L4_AS00403_contig1-L4_AS00403_00019\",\"L4_AS00404-L4_AS00404_contig1-L4_AS00404_00019\",\"L5_AS00501-L5_AS00501_contig1-L5_AS00501_00018\",\"L8_AS00801-L8_AS00801_contig1-L8_AS00801_00019\""], + "L1_AS00101", + "Rv0019c", + "L1_AS00101-L1_AS00101_contig1-L1_AS00101_00019", + ), + ( + ["\"Gene\",\"Non-unique Gene name\",\"Annotation\",\"No. isolates\",\"No. sequences\",\"Avg sequences per isolate\",\"Accessory Fragment\",\"Genome Fragment\",\"Order within Fragment\",\"Accessory Order with Fragment\",\"QC\",\"Min group size nuc\",\"Max group size nuc\",\"Avg group size nuc\",\"L1_AS00101\",\"L1_AS00102\",\"L2_AS00201\",\"L2_AS00202\",\"L2_AS00203\",\"L3_AS00301\",\"L4_AS00401\",\"L4_AS00402\",\"L4_AS00403\",\"L4_AS00404\",\"L5_AS00501\",\"L8_AS00801\"", "\"L1_AS00102_00001\",\"persistent\",\"chromosomal replication initiator DnaA\",\"12\",\"12\",\"1.0\",\"NA\",\"NA\",\"\",\"\",\"\",\"1238\",\"1238\",\"1238.0\",\"L1_AS00101_00001\",\"L1_AS00102_00001\",\"L2_AS00201_00001\",\"L2_AS00202_00001\",\"L2_AS00203_00001\",\"L3_AS00301_00001\",\"L4_AS00401_00001\",\"L4_AS00402_00001\",\"L4_AS00403_00001\",\"L4_AS00404_00001\",\"L5_AS00501_00001\",\"L8_AS00801_00001\""], + "L1_AS00101", + "L1_AS00102_00001", + "L1_AS00101_00001", + ), + ] { + let f = write_temp(&format!("{}\n{}\n", lines[0], lines[1])); + let pa = parse_presence_absence(f.path()).unwrap(); + assert_eq!(pa.n_genomes(), 12, "{cluster}"); + assert_eq!(pa.genome_name(GenomeIdx(0)), first_genome); + assert_eq!(pa.cluster_name(ClusterIdx(0)), cluster); + assert_eq!(pa.cell(ClusterIdx(0), GenomeIdx(0)), cell); + } + } } diff --git a/src/pipeline.rs b/src/pipeline.rs index 212dfc7..f1074df 100644 --- a/src/pipeline.rs +++ b/src/pipeline.rs @@ -15,7 +15,7 @@ use crate::parser::{ }; use crate::relationships::{write_cluster_components, write_cluster_relationships}; use crate::rescue::{ - build_gene_to_cluster_with_gff_aliases, check_internal_stops_in_cds, + build_gene_to_cluster_with_gff_aliases, cell_gene_id_variants, check_internal_stops_in_cds, check_internal_stops_in_hit, check_internal_stops_in_span, classify_hits, panaroo_gene_id_suffix_variants, resolve_conflicts, validate_synteny_with_cache, GffGene, GffIndex, PseudogeneReason, RescueEvent, RescueType, SyntenyEvidence, SyntenyReferenceCache, @@ -37,8 +37,8 @@ use self::evolution::{ use self::output::{ clear_optional_outputs, count_core_after, count_soft_core_after, write_corrected_csv, write_corrected_gffs, write_html_report, write_multistate_matrix, write_nsc_report, - write_run_manifest, write_summary, MultistateMatrixInputs, NotSearched, SearchCoverage, - SummaryStats, + write_rescued_fasta, write_run_manifest, write_summary, MultistateMatrixInputs, NotSearched, + SearchCoverage, SummaryStats, }; use self::quality::{ cluster_prevalences, cluster_sequence_lengths, orthology_quality, @@ -282,81 +282,52 @@ pub fn run(args: &Args) -> Result<()> { let mut parse_failures = 0usize; let mut unusable = 0usize; let annotation_files = resolve_annotation_files_for_pa(&pa, ann_dir); - for gi_val in 0..n_genomes { - let gi = GenomeIdx(gi_val); - let genome_name = pa.genome_name(gi); - let Some(path) = annotation_files.get(genome_name).cloned() else { - continue; - }; - match parse_gff_genes(&path) { - Ok(mut raw_genes) => { - if let Some(fasta_path) = genome_files.get(genome_name) { - let contig_names = - load_genome_contig_names(fasta_path).with_context(|| { - format!("Cannot read FASTA contig names for {}", genome_name) - })?; - if let Some((from, to, n_records)) = - normalize_gff_records_to_fasta(&mut raw_genes, &contig_names) - { - gff_contig_remaps.push((genome_name.to_string(), from, to, n_records)); - } - // A GFF that annotates none of the assembly's contigs, or - // has no CDS at all, used to count as loaded: every hit - // then sat on "unannotated" sequence, became MISSING_CDS, - // and nothing was corrected, with no warning that said so. - let fasta_set: ahash::HashSet<&str> = - contig_names.iter().map(String::as_str).collect(); - let on_assembly = raw_genes - .iter() - .filter(|record| fasta_set.contains(record.0.as_str())) - .count(); - if raw_genes.is_empty() || on_assembly == 0 { - unusable += 1; - if raw_genes.is_empty() { - warn!( - "GFF for {} ({}) has no CDS records; it is not used", - genome_name, - path.display() - ); - } else { - warn!( - "GFF for {} annotates none of its assembly's contigs \ - (GFF seqid '{}', FASTA contig '{}'); it is not used", - genome_name, - raw_genes[0].0, - contig_names.first().map(String::as_str).unwrap_or("") - ); - } - continue; + // Each genome's GFF is read and indexed on its own, in parallel, and the + // outcomes are folded in genome order, so the warnings, the remaps and + // the index are what the one-at-a-time loop produced. + let gff_pool = rayon::ThreadPoolBuilder::new() + .num_threads(args.threads.max(1)) + .build() + .context("Failed to build the GFF thread pool")?; + let outcomes: Vec> = gff_pool.install(|| { + use rayon::prelude::*; + (0..n_genomes) + .into_par_iter() + .map(|gi_val| { + let genome_name = pa.genome_name(GenomeIdx(gi_val)); + match annotation_files.get(genome_name) { + Some(path) => { + load_genome_gff(genome_name, path, genome_files.get(genome_name)) } + None => Ok(GffOutcome::NoFile), + } + }) + .collect() + }); + for (gi_val, outcome) in outcomes.into_iter().enumerate() { + let gi = GenomeIdx(gi_val); + match outcome? { + GffOutcome::Loaded { index, path, remap } => { + if let Some((from, to, n_records)) = remap { + gff_contig_remaps.push(( + pa.genome_name(gi).to_string(), + from, + to, + n_records, + )); } - let genes: Vec<(GffGene, Vec, Vec)> = raw_genes - .into_iter() - .map(|(contig, start, end, strand, gene_id, aliases, labels)| { - ( - GffGene { - contig, - start, - end, - strand, - gene_id, - }, - aliases, - labels, - ) - }) - .collect(); - data.insert(gi, GffIndex::build_with_aliases_and_labels(genes)); - gff_paths.insert(gi, path.clone()); + data.insert(gi, *index); + gff_paths.insert(gi, path); + } + GffOutcome::Unusable(message) => { + unusable += 1; + warn!("{}", message); } - // `{:#}` for the same reason as the FASTA warning above, and it - // matters more here: this is the one line telling a reader that - // a genome just lost its whole annotation, and the run goes on - // to a different answer without it. - Err(e) => { + GffOutcome::ParseFailure(message) => { parse_failures += 1; - warn!("Cannot parse GFF for {}: {:#}", genome_name, e) + warn!("{}", message); } + GffOutcome::NoFile => {} } } let loaded = data.len(); @@ -464,6 +435,7 @@ pub fn run(args: &Args) -> Result<()> { &args.preset, true, 500, + args.search_cache.as_deref(), )?; let reference_hit_candidates = if let Some(reference_genome) = args.reference_genome.as_deref() @@ -495,6 +467,7 @@ pub fn run(args: &Args) -> Result<()> { &args.preset, true, 5000, + args.search_cache.as_deref(), )?; reference_hit_shadow_candidates( &reference_hits, @@ -578,84 +551,36 @@ pub fn run(args: &Args) -> Result<()> { rescues_by_genome.len() ); let mut upgraded = 0usize; - - for (gi, rescue_indices) in rescues_by_genome { - let genome_name = pa.genome_name(gi); - let Some(fasta_path) = genome_files.get(genome_name) else { - continue; - }; - // Propagate, never skip: a genome that cannot be scanned used to - // be passed over in silence, and its pseudogenes went into the - // matrix as intact genes. - let contigs = - crate::parser::load_genome_contigs(fasta_path).with_context(|| { - format!( - "Cannot scan genome {} for internal stop codons", - genome_name + // One genome at a time was the slowest part of a large run after + // the search: each genome is read and scanned on its own, in + // parallel, and the findings applied afterwards. Genome order keeps + // the first error the same from run to run. + let mut genomes: Vec<(GenomeIdx, Vec)> = rescues_by_genome.into_iter().collect(); + genomes.sort_by_key(|(gi, _)| *gi); + let scans: Vec>> = pool.install(|| { + use rayon::prelude::*; + genomes + .par_iter() + .map(|(gi, rescue_indices)| { + scan_genome_for_internal_stops( + *gi, + rescue_indices, + &rescues, + &pa, + &genome_files, + &gff_indices, ) - })?; - - for rescue_idx in rescue_indices { + }) + .collect() + }); + for found in scans { + for (rescue_idx, n_stops) in found? { let rescue = &mut rescues[rescue_idx]; - // A rescue on an annotated CDS is scanned in that CDS, in its - // frame and to its stop; see check_internal_stops_in_cds. - let annotated_cds = rescue - .annotation_match - .as_deref() - .and_then(|id| gff_indices.get(&gi)?.get_by_id(id)) - .filter(|gene| { - gene.contig == rescue.contig - && gene.strand == rescue.strand - && gene.start < rescue.end - && rescue.start < gene.end - }); - let stops = if let Some(gene) = annotated_cds { - check_internal_stops_in_cds( - &gene.contig, - gene.start, - gene.end, - gene.strand, - &contigs, - ) - } else if let Some(query_start) = rescue.query_start { - // No annotation, but one alignment: its start on the - // representative fixes the reading frame. - check_internal_stops_in_span( - &rescue.contig, - rescue.start, - rescue.end, - rescue.strand, - query_start, - &contigs, - ) - } else { - check_internal_stops_in_hit( - // Build a minimal hit from rescue coords for the check - &crate::search::Hit { - query: String::new(), - cluster: String::new(), - query_length: 0, - query_start: 0, - query_end: 0, - strand: rescue.strand, - target: rescue.contig.clone(), - target_start: rescue.start, - target_end: rescue.end, - target_length: 0, - matches: 0, - alignment_length: 0, - mapq: 0, - }, - Some(&contigs), - ) + rescue.rescue_type = RescueType::Pseudogene { + reason: PseudogeneReason::InternalStop { n_stops }, }; - if let Some(n_stops) = stops { - rescue.rescue_type = RescueType::Pseudogene { - reason: PseudogeneReason::InternalStop { n_stops }, - }; - rescue.confidence *= 0.5; // Pseudogenes get lower confidence - upgraded += 1; - } + rescue.confidence *= 0.5; // Pseudogenes get lower confidence + upgraded += 1; } } @@ -671,17 +596,51 @@ pub fn run(args: &Args) -> Result<()> { // ── 8. Synteny validation ── if args.synteny && !gff_indices.is_empty() { info!("Validating synteny..."); - let mut synteny_cache = SyntenyReferenceCache::default(); - for rescue in &mut rescues { - validate_synteny_with_cache( - rescue, - &pa, - &gff_indices, - &gene_to_cluster_maps, - &mut synteny_cache, - args.verbose, - ); + // A rescue reads only its own cluster's reference contexts, so the + // rescues are grouped by cluster and the clusters validated in + // parallel, each with its own cache, then put back in their order. One + // rescue at a time it was the slowest stage of a large run. + let mut slots: Vec> = + std::mem::take(&mut rescues).into_iter().map(Some).collect(); + let mut order: Vec = (0..slots.len()).collect(); + order.sort_by_key(|&i| (slots[i].as_ref().map(|r| r.cluster_idx), i)); + let mut groups: Vec> = Vec::new(); + for i in order { + let rescue = slots[i].take().expect("each rescue is moved once"); + let same_cluster = groups + .last() + .is_some_and(|group| group[0].1.cluster_idx == rescue.cluster_idx); + if !same_cluster { + groups.push(Vec::new()); + } + groups + .last_mut() + .expect("a group was just pushed") + .push((i, rescue)); } + pool.install(|| { + use rayon::prelude::*; + groups.par_iter_mut().for_each(|group| { + let mut synteny_cache = SyntenyReferenceCache::default(); + for (_, rescue) in group.iter_mut() { + validate_synteny_with_cache( + rescue, + &pa, + &gff_indices, + &gene_to_cluster_maps, + &mut synteny_cache, + args.verbose, + ); + } + }); + }); + for (i, rescue) in groups.into_iter().flatten() { + slots[i] = Some(rescue); + } + rescues = slots + .into_iter() + .map(|rescue| rescue.expect("every rescue comes back")) + .collect(); let validated = rescues .iter() .filter(|r| r.synteny_ok == Some(true)) @@ -702,7 +661,7 @@ pub fn run(args: &Args) -> Result<()> { warn!(" Rerun with -v to see detailed synteny debug output"); } } else if args.synteny && gff_indices.is_empty() { - warn!("Synteny requested (--synteny) but no GFF annotations loaded — skipping"); + warn!("Synteny requested (--synteny) but no GFF annotations loaded, so it is skipped"); } let reference_audit_candidates = if args.reference_genome.is_some() { @@ -889,6 +848,7 @@ pub fn run(args: &Args) -> Result<()> { min_identity: args.mobile_element_min_identity, min_coverage: args.mobile_element_min_coverage, preset: &args.preset, + search_cache: args.search_cache.as_deref(), })? } else { Vec::new() @@ -993,6 +953,7 @@ pub fn run(args: &Args) -> Result<()> { &gff_indices, &corrected_gene_to_cluster_maps, synteny_reference_gi, + args.threads, &args.output, )?; info!( @@ -1061,6 +1022,11 @@ pub fn run(args: &Args) -> Result<()> { coverage: &coverage, output_dir: &args.output, })?; + if args.write_rescued_fasta { + let mut matrix_rescues: Vec<&RescueEvent> = matrix_rescue_map.values().copied().collect(); + matrix_rescues.sort_by_key(|r| (r.genome_idx, r.cluster_idx)); + write_rescued_fasta(&matrix_rescues, &pa, &genome_files, &args.output)?; + } if args.write_corrected_gff { let matrix_rescues: Vec = matrix_rescue_map.values().map(|r| (*r).clone()).collect(); @@ -1421,9 +1387,176 @@ fn stripped_gene_id(gene_id: &str) -> Option<&str> { panaroo_gene_id_suffix_variants(gene_id).into_iter().nth(1) } +/// The rescues of one genome that read through an internal stop, with how +/// many stops each has. +fn scan_genome_for_internal_stops( + gi: GenomeIdx, + rescue_indices: &[usize], + rescues: &[RescueEvent], + pa: &crate::parser::PresenceAbsence, + genome_files: &HashMap, + gff_indices: &HashMap, +) -> Result> { + let genome_name = pa.genome_name(gi); + let Some(fasta_path) = genome_files.get(genome_name) else { + return Ok(Vec::new()); + }; + // Propagate, never skip: a genome that cannot be scanned used to be passed + // over in silence, and its pseudogenes went into the matrix as intact genes. + let contigs = crate::parser::load_genome_contigs(fasta_path).with_context(|| { + format!( + "Cannot scan genome {} for internal stop codons", + genome_name + ) + })?; + let mut found = Vec::new(); + for &rescue_idx in rescue_indices { + let rescue = &rescues[rescue_idx]; + // A rescue on an annotated CDS is scanned in that CDS, in its frame and + // to its stop; see check_internal_stops_in_cds. + let annotated_cds = rescue + .annotation_match + .as_deref() + .and_then(|id| gff_indices.get(&gi)?.get_by_id(id)) + .filter(|gene| { + gene.contig == rescue.contig + && gene.strand == rescue.strand + && gene.start < rescue.end + && rescue.start < gene.end + }); + let stops = if let Some(gene) = annotated_cds { + check_internal_stops_in_cds(&gene.contig, gene.start, gene.end, gene.strand, &contigs) + } else if let Some(query_start) = rescue.query_start { + // No annotation, but one alignment: its start on the representative + // fixes the reading frame. + check_internal_stops_in_span( + &rescue.contig, + rescue.start, + rescue.end, + rescue.strand, + query_start, + &contigs, + ) + } else { + check_internal_stops_in_hit( + // Build a minimal hit from rescue coords for the check + &crate::search::Hit { + query: String::new(), + cluster: String::new(), + query_length: 0, + query_start: 0, + query_end: 0, + strand: rescue.strand, + target: rescue.contig.clone(), + target_start: rescue.start, + target_end: rescue.end, + target_length: 0, + matches: 0, + alignment_length: 0, + mapq: 0, + }, + Some(&contigs), + ) + }; + if let Some(n_stops) = stops { + found.push((rescue_idx, n_stops)); + } + } + Ok(found) +} + +/// What became of one genome's GFF. +enum GffOutcome { + Loaded { + index: Box, + path: PathBuf, + /// A single GFF seqid mapped onto a differently named FASTA contig. + remap: Option<(String, String, usize)>, + }, + /// Read, but annotating nothing of its assembly: the warning to give. + Unusable(String), + /// Not readable as GFF: the warning to give. + ParseFailure(String), + NoFile, +} + +/// Read and index one genome's GFF against its assembly's contig names. +fn load_genome_gff( + genome_name: &str, + path: &Path, + fasta_path: Option<&PathBuf>, +) -> Result { + let mut raw_genes = match parse_gff_genes(path) { + Ok(raw_genes) => raw_genes, + // `{:#}` for the same reason as the FASTA warning, and it matters more + // here: this is the one line telling a reader that a genome just lost + // its whole annotation, and the run goes on to a different answer + // without it. + Err(e) => { + return Ok(GffOutcome::ParseFailure(format!( + "Cannot parse GFF for {}: {:#}", + genome_name, e + ))) + } + }; + let mut remap = None; + if let Some(fasta_path) = fasta_path { + let contig_names = load_genome_contig_names(fasta_path) + .with_context(|| format!("Cannot read FASTA contig names for {}", genome_name))?; + remap = normalize_gff_records_to_fasta(&mut raw_genes, &contig_names); + // A GFF that annotates none of the assembly's contigs, or has no CDS at + // all, used to count as loaded: every hit then sat on "unannotated" + // sequence, became MISSING_CDS, and nothing was corrected, with no + // warning that said so. + let fasta_set: ahash::HashSet<&str> = contig_names.iter().map(String::as_str).collect(); + let on_assembly = raw_genes + .iter() + .filter(|record| fasta_set.contains(record.0.as_str())) + .count(); + if raw_genes.is_empty() { + return Ok(GffOutcome::Unusable(format!( + "GFF for {} ({}) has no CDS records; it is not used", + genome_name, + path.display() + ))); + } + if on_assembly == 0 { + return Ok(GffOutcome::Unusable(format!( + "GFF for {} annotates none of its assembly's contigs \ + (GFF seqid '{}', FASTA contig '{}'); it is not used", + genome_name, + raw_genes[0].0, + contig_names.first().map(String::as_str).unwrap_or("") + ))); + } + } + let genes: Vec<(GffGene, Vec, Vec)> = raw_genes + .into_iter() + .map(|(contig, start, end, strand, gene_id, aliases, labels)| { + ( + GffGene { + contig, + start, + end, + strand, + gene_id, + }, + aliases, + labels, + ) + }) + .collect(); + Ok(GffOutcome::Loaded { + index: Box::new(GffIndex::build_with_aliases_and_labels(genes)), + path: path.to_path_buf(), + remap, + }) +} + fn gff_gene_for_id<'a>(gff: &'a GffIndex, gene_id: &str) -> Option<&'a GffGene> { - gff.get_by_id(gene_id) - .or_else(|| stripped_gene_id(gene_id).and_then(|id| gff.get_by_id(id))) + cell_gene_id_variants(gene_id) + .into_iter() + .find_map(|id| gff.get_by_id(id)) } fn gff_gene_for_cell<'a>(gff: &'a GffIndex, cell: &str) -> Option<&'a GffGene> { diff --git a/src/pipeline/evolution.rs b/src/pipeline/evolution.rs index 10d613c..3b80fd6 100644 --- a/src/pipeline/evolution.rs +++ b/src/pipeline/evolution.rs @@ -768,6 +768,7 @@ pub(super) struct MobileElementSearchInput<'a> { pub(super) min_identity: f64, pub(super) min_coverage: f64, pub(super) preset: &'a str, + pub(super) search_cache: Option<&'a Path>, } pub(super) fn search_mobile_element_reference( @@ -817,6 +818,7 @@ pub(super) fn search_mobile_element_reference( input.preset, true, 5000, + input.search_cache, )?; let mut mobile_hits = Vec::new(); @@ -2456,74 +2458,27 @@ pub(super) fn write_cluster_locus_contexts( Ok((written_context_rows, summaries)) } +/// Unique flanking-gene blocks per cluster, in `cluster_synteny_blocks.tsv`. +/// +/// Every cluster's blocks used to be gathered before the first row was written, +/// in one thread: on 200 genomes of 4,000 genes that held 1.2 GB and took half +/// the run. Clusters do not depend on each other, so they are now built in +/// parallel a chunk at a time and written in cluster order, the same rows in +/// the same order, holding one chunk's blocks at a time. pub(super) fn write_cluster_synteny_blocks( pa: &crate::parser::PresenceAbsence, gff_indices: &HashMap, gene_to_cluster_maps: &HashMap>, reference_gi: Option, + threads: usize, output_dir: &std::path::Path, ) -> Result { - let mut block_rows: BTreeMap<(usize, String), SyntenyBlockAccumulator> = BTreeMap::new(); let reference_boundary_index = reference_gi .and_then(|gi| gff_indices.get(&gi).zip(gene_to_cluster_maps.get(&gi))) .map(|(gff, g2c)| build_reference_boundary_index(gff, g2c, pa, 5)) .unwrap_or_default(); let cluster_index = pa.cluster_index(); - for ci_val in 0..pa.n_clusters() { - let ci = ClusterIdx(ci_val); - let mut genomes = pa.present_genomes(ci).to_vec(); - genomes.sort_by_key(|gi| gi.val()); - - for gi in genomes { - let Some(gff) = gff_indices.get(&gi) else { - continue; - }; - let Some(g2c) = gene_to_cluster_maps.get(&gi) else { - continue; - }; - let Some(gene) = gff_gene_for_cell(gff, pa.cell(ci, gi)) else { - continue; - }; - - let (upstream, downstream) = gff.context_window( - &gene.contig, - gene.start, - gene.end, - Some(gene.gene_id.as_str()), - 5, - ); - let upstream_sig = context_signature(gff, &upstream, g2c, pa); - let downstream_sig = context_signature(gff, &downstream, g2c, pa); - let representative_genes = - ordered_synteny_gene_records(gff, &upstream, Some(gene), &downstream, g2c, pa); - let block_signature = synteny_gene_identity_signature(&representative_genes); - let row = block_rows.entry((ci.val(), block_signature)).or_default(); - row.count += 1; - if row.upstream_signature.is_empty() { - row.upstream_signature = upstream_sig; - row.downstream_signature = downstream_sig; - } - row.carrier_genomes.push(pa.genome_name(gi).to_string()); - if row.examples.len() < 5 { - row.examples - .push(format!("{}:{}", pa.genome_name(gi), gene.gene_id)); - } - if row.representative.is_none() { - row.representative = Some(SyntenyBlockRepresentative { - genome: pa.genome_name(gi).to_string(), - gene_id: gene.gene_id.clone(), - contig: gene.contig.clone(), - start: gene.start + 1, - end: gene.end, - strand: gene.strand, - length: gene.end.saturating_sub(gene.start), - genes: representative_genes, - }); - } - } - } - let path = output_dir.join("cluster_synteny_blocks.tsv"); let file = fs::File::create(&path)?; let mut w = BufWriter::new(file); @@ -2532,21 +2487,108 @@ pub(super) fn write_cluster_synteny_blocks( "gene_cluster\tblock_rank\tn_genomes\tfraction_of_present\tupstream_block_5\tdownstream_block_5\tblock_interpretation\texample_genomes\tcarrier_genomes\trepresentative_genome\trepresentative_gene\trepresentative_contig\trepresentative_start\trepresentative_end\trepresentative_strand\trepresentative_length_bp\trepresentative_block_genes\treference_block_genes" )?; - let mut rows_by_cluster: BTreeMap> = - BTreeMap::new(); - for ((ci, signature), acc) in block_rows { - rows_by_cluster - .entry(ci) - .or_default() - .push((signature, acc)); + const CLUSTERS_PER_CHUNK: usize = 256; + let pool = rayon::ThreadPoolBuilder::new() + .num_threads(threads.max(1)) + .build() + .context("Failed to build the synteny block thread pool")?; + let clusters: Vec = (0..pa.n_clusters()).collect(); + let mut written = 0usize; + for chunk in clusters.chunks(CLUSTERS_PER_CHUNK) { + let rows: Vec>> = pool.install(|| { + use rayon::prelude::*; + chunk + .par_iter() + .map(|&ci_val| { + cluster_synteny_block_rows( + ClusterIdx(ci_val), + pa, + gff_indices, + gene_to_cluster_maps, + &reference_boundary_index, + &cluster_index, + ) + }) + .collect() + }); + for cluster_rows in rows { + for fields in cluster_rows { + write_tsv_fields(&mut w, &fields)?; + written += 1; + } + } } + w.flush()?; + Ok(written) +} - let mut written = 0usize; - for (ci_val, mut rows) in rows_by_cluster { - rows.sort_by(|a, b| b.1.count.cmp(&a.1.count).then(a.0.cmp(&b.0))); - let present = pa.present_count(ClusterIdx(ci_val)).max(1); - let n_blocks = rows.len(); - for (rank, (_signature, acc)) in rows.into_iter().enumerate() { +/// One cluster's synteny block rows, ranked by carriers. +fn cluster_synteny_block_rows( + ci: ClusterIdx, + pa: &crate::parser::PresenceAbsence, + gff_indices: &HashMap, + gene_to_cluster_maps: &HashMap>, + reference_boundary_index: &ReferenceBoundaryIndex, + cluster_index: &HashMap, +) -> Vec> { + let mut blocks: BTreeMap = BTreeMap::new(); + let mut genomes = pa.present_genomes(ci).to_vec(); + genomes.sort_by_key(|gi| gi.val()); + for gi in genomes { + let Some(gff) = gff_indices.get(&gi) else { + continue; + }; + let Some(g2c) = gene_to_cluster_maps.get(&gi) else { + continue; + }; + let Some(gene) = gff_gene_for_cell(gff, pa.cell(ci, gi)) else { + continue; + }; + + let (upstream, downstream) = gff.context_window( + &gene.contig, + gene.start, + gene.end, + Some(gene.gene_id.as_str()), + 5, + ); + let upstream_sig = context_signature(gff, &upstream, g2c, pa); + let downstream_sig = context_signature(gff, &downstream, g2c, pa); + let representative_genes = + ordered_synteny_gene_records(gff, &upstream, Some(gene), &downstream, g2c, pa); + let block_signature = synteny_gene_identity_signature(&representative_genes); + let row = blocks.entry(block_signature).or_default(); + row.count += 1; + if row.upstream_signature.is_empty() { + row.upstream_signature = upstream_sig; + row.downstream_signature = downstream_sig; + } + row.carrier_genomes.push(pa.genome_name(gi).to_string()); + if row.examples.len() < 5 { + row.examples + .push(format!("{}:{}", pa.genome_name(gi), gene.gene_id)); + } + if row.representative.is_none() { + row.representative = Some(SyntenyBlockRepresentative { + genome: pa.genome_name(gi).to_string(), + gene_id: gene.gene_id.clone(), + contig: gene.contig.clone(), + start: gene.start + 1, + end: gene.end, + strand: gene.strand, + length: gene.end.saturating_sub(gene.start), + genes: representative_genes, + }); + } + } + + let mut rows: Vec<(String, SyntenyBlockAccumulator)> = blocks.into_iter().collect(); + rows.sort_by(|a, b| b.1.count.cmp(&a.1.count).then(a.0.cmp(&b.0))); + let present = pa.present_count(ci).max(1); + let n_blocks = rows.len(); + rows.into_iter() + .enumerate() + .map(|(rank, (_signature, acc))| { let fraction = acc.count as f64 / present as f64; let interpretation = if n_blocks <= 1 { "single_synteny_block" @@ -2559,55 +2601,46 @@ pub(super) fn write_cluster_synteny_blocks( let target_clusters = parse_context_cluster_indices( &acc.upstream_signature, &acc.downstream_signature, - &cluster_index, + cluster_index, ); - let reference_genes = best_reference_context( - &reference_boundary_index, - ClusterIdx(ci_val), - &target_clusters, - ); - write_tsv_fields( - &mut w, - &[ - pa.cluster_name(ClusterIdx(ci_val)).to_string(), - (rank + 1).to_string(), - acc.count.to_string(), - format!("{:.3}", fraction), - acc.upstream_signature, - acc.downstream_signature, - interpretation.to_string(), - acc.examples.join(";"), - acc.carrier_genomes.join(";"), - representative - .map(|rep| rep.genome.clone()) - .unwrap_or_default(), - representative - .map(|rep| rep.gene_id.clone()) - .unwrap_or_default(), - representative - .map(|rep| rep.contig.clone()) - .unwrap_or_default(), - representative - .map(|rep| rep.start.to_string()) - .unwrap_or_default(), - representative - .map(|rep| rep.end.to_string()) - .unwrap_or_default(), - representative - .map(|rep| rep.strand.to_string()) - .unwrap_or_default(), - representative - .map(|rep| rep.length.to_string()) - .unwrap_or_default(), - representative - .map(|rep| encode_synteny_gene_records(&rep.genes)) - .unwrap_or_default(), - encode_synteny_gene_records(&reference_genes), - ], - )?; - written += 1; - } - } - w.flush()?; - Ok(written) + let reference_genes = + best_reference_context(reference_boundary_index, ci, &target_clusters); + vec![ + pa.cluster_name(ci).to_string(), + (rank + 1).to_string(), + acc.count.to_string(), + format!("{:.3}", fraction), + acc.upstream_signature.clone(), + acc.downstream_signature.clone(), + interpretation.to_string(), + acc.examples.join(";"), + acc.carrier_genomes.join(";"), + representative + .map(|rep| rep.genome.clone()) + .unwrap_or_default(), + representative + .map(|rep| rep.gene_id.clone()) + .unwrap_or_default(), + representative + .map(|rep| rep.contig.clone()) + .unwrap_or_default(), + representative + .map(|rep| rep.start.to_string()) + .unwrap_or_default(), + representative + .map(|rep| rep.end.to_string()) + .unwrap_or_default(), + representative + .map(|rep| rep.strand.to_string()) + .unwrap_or_default(), + representative + .map(|rep| rep.length.to_string()) + .unwrap_or_default(), + representative + .map(|rep| encode_synteny_gene_records(&rep.genes)) + .unwrap_or_default(), + encode_synteny_gene_records(&reference_genes), + ] + }) + .collect() } diff --git a/src/pipeline/output.rs b/src/pipeline/output.rs index 2d56686..cd8d9c8 100644 --- a/src/pipeline/output.rs +++ b/src/pipeline/output.rs @@ -38,6 +38,18 @@ pub(super) fn write_corrected_csv( csv_writer.write_record(pa.header())?; + // The same matrix as 1/0, Panaroo's gene_presence_absence.Rtab layout, which + // is what pyseer, Scoary-style scripts and most R code read. Written in the + // same pass, so it cannot disagree with the CSV on any cell. + let rtab_path = output_dir.join("gene_presence_absence_recore.Rtab"); + let rtab_file = fs::File::create(&rtab_path) + .with_context(|| format!("Cannot create {}", rtab_path.display()))?; + let mut rtab = csv::WriterBuilder::new() + .delimiter(b'\t') + .from_writer(BufWriter::with_capacity(256 * 1024, rtab_file)); + rtab.write_record(std::iter::once("Gene").chain(pa.genome_names().iter().map(String::as_str)))?; + let mut rtab_row: Vec<&str> = Vec::with_capacity(pa.n_genomes() + 1); + let genome_idx_set: ahash::HashSet = pa.genome_col_indices().iter().copied().collect(); let meta_col_indices: Vec = (0..pa.header().len()) .filter(|i| !genome_idx_set.contains(i)) @@ -73,6 +85,8 @@ pub(super) fn write_corrected_csv( let mut isolate_count = 0usize; let mut sequence_count = 0usize; + rtab_row.clear(); + rtab_row.push(pa.cluster_name(ci)); for (gi_val, &col_idx) in pa.genome_col_indices().iter().enumerate() { let gi = GenomeIdx(gi_val); if let Some(r) = rescue_map.get(&(ci, gi)) { @@ -85,6 +99,7 @@ pub(super) fn write_corrected_csv( isolate_count += 1; sequence_count += 1; counts.rescued_cells += 1; + rtab_row.push("1"); } else if nsc_presence.contains(&(ci_val, gi_val)) { // Gene absent in this cluster but present in a merged NSC member row[col_idx].clear(); @@ -92,6 +107,7 @@ pub(super) fn write_corrected_csv( isolate_count += 1; sequence_count += 1; counts.nsc_cells += 1; + rtab_row.push("1"); } else { row[col_idx].clear(); let cell = pa.cell(ci, gi); @@ -99,11 +115,14 @@ pub(super) fn write_corrected_csv( if !cell.trim().is_empty() { isolate_count += 1; sequence_count += cell_sequence_count(cell); + rtab_row.push("1"); } else { counts.absent_cells += 1; + rtab_row.push("0"); } } } + rtab.write_record(&rtab_row)?; use std::fmt::Write as FmtWrite; if let Some(col) = isolates_col { @@ -127,6 +146,8 @@ pub(super) fn write_corrected_csv( csv_writer.write_record(row.iter().map(|s| s.as_str()))?; } csv_writer.flush()?; + rtab.flush() + .with_context(|| format!("Cannot write {}", rtab_path.display()))?; Ok(counts) } @@ -142,6 +163,8 @@ pub(super) fn clear_optional_outputs(output_dir: &std::path::Path) -> Result<()> "cluster_synteny_blocks.tsv", "recore_report.html", "recore_report.md", + "rescued_genes.fna", + "rescued_genes.faa", ] { let path = output_dir.join(filename); match fs::remove_file(&path) { @@ -650,6 +673,133 @@ pub(super) fn write_corrected_gffs( Ok(()) } +/// The sequence of every cell the binary matrix gained, as nucleotides read in +/// the gene's direction and, where the reading frame is known, as protein. +/// +/// IDs are the corrected GFF's feature IDs, so a sequence can be traced to its +/// feature and back; coordinates, strand and phase are the GFF's too. The +/// protein starts at the phase the alignment fixes and uses the bacterial code +/// (table 11). A rescue built from several alignments has no single frame and +/// gets no protein; its nucleotide record says `phase=unknown`. +pub(super) fn write_rescued_fasta( + rescues: &[&RescueEvent], + pa: &crate::parser::PresenceAbsence, + genome_files: &HashMap, + output_dir: &std::path::Path, +) -> Result<(usize, usize)> { + let mut by_genome: std::collections::BTreeMap> = + std::collections::BTreeMap::new(); + for rescue in rescues { + by_genome.entry(rescue.genome_idx).or_default().push(rescue); + } + let fna_path = output_dir.join("rescued_genes.fna"); + let faa_path = output_dir.join("rescued_genes.faa"); + let mut fna = BufWriter::new( + fs::File::create(&fna_path) + .with_context(|| format!("Cannot create {}", fna_path.display()))?, + ); + let mut faa = BufWriter::new( + fs::File::create(&faa_path) + .with_context(|| format!("Cannot create {}", faa_path.display()))?, + ); + let (mut n_nucleotide, mut n_protein) = (0usize, 0usize); + for (gi, mut rows) in by_genome { + let genome_name = pa.genome_name(gi); + let Some(fasta_path) = genome_files.get(genome_name) else { + continue; + }; + let contigs = crate::parser::load_genome_contigs(fasta_path) + .with_context(|| format!("Cannot read genome {} for its rescued genes", genome_name))?; + rows.sort_by(|a, b| { + a.contig + .cmp(&b.contig) + .then(a.start.cmp(&b.start)) + .then(a.end.cmp(&b.end)) + .then(a.cluster_idx.cmp(&b.cluster_idx)) + }); + for rescue in rows { + let Some(contig) = contigs.get(&rescue.contig) else { + anyhow::bail!( + "Rescue on contig {} is not in genome {}", + rescue.contig, + fasta_path.display() + ); + }; + let (start, end) = (rescue.start as usize, rescue.end as usize); + if start >= end || end > contig.len() { + anyhow::bail!( + "Rescue {}:{}-{} lies outside contig {} of genome {}", + rescue.contig, + start + 1, + end, + rescue.contig, + genome_name + ); + } + let span = &contig[start..end]; + let oriented = if rescue.strand == '-' { + crate::rescue::reverse_complement(span) + } else { + span.to_ascii_uppercase() + }; + let cluster = pa.cluster_name(rescue.cluster_idx); + let id = format!( + "recore_{}_{}_{}_{}", + safe_id_component(cluster), + safe_id_component(&rescue.contig), + rescue.start + 1, + rescue.end + ); + let phase = rescue.query_start.map(|qs| ((3 - qs % 3) % 3) as usize); + let description = format!( + "genome={} gene_cluster={} locus={}:{}-{}({}) rescue_type={} phase={}", + genome_name, + cluster, + rescue.contig, + rescue.start + 1, + rescue.end, + rescue.strand, + rescue.rescue_type, + phase.map_or_else(|| "unknown".to_string(), |p| p.to_string()) + ); + writeln!(fna, ">{id} {description}")?; + write_wrapped(&mut fna, &oriented)?; + n_nucleotide += 1; + if let Some(phase) = phase.filter(|&p| p < oriented.len()) { + // The representative begins at its start codon, so an alignment + // from its first base starts on the gene's own start codon. + let protein = crate::rescue::translate_bacterial( + &oriented[phase..], + rescue.query_start == Some(0), + ); + if !protein.is_empty() { + writeln!(faa, ">{id} {description}")?; + write_wrapped(&mut faa, protein.as_bytes())?; + n_protein += 1; + } + } + } + } + fna.flush()?; + faa.flush()?; + info!( + "Rescued gene sequences: {} nucleotide and {} protein records in {}", + n_nucleotide, + n_protein, + output_dir.display() + ); + Ok((n_nucleotide, n_protein)) +} + +/// FASTA sequence lines of 60 characters. +fn write_wrapped(w: &mut impl Write, seq: &[u8]) -> Result<()> { + for chunk in seq.chunks(60) { + w.write_all(chunk)?; + w.write_all(b"\n")?; + } + Ok(()) +} + pub(super) fn count_core_after( pa: &crate::parser::PresenceAbsence, rescue_map: &HashMap<(ClusterIdx, GenomeIdx), &RescueEvent>, @@ -1037,6 +1187,19 @@ pub(super) fn write_run_manifest( " \"write_corrected_gff\": {},", args.write_corrected_gff )?; + writeln!( + w, + " \"write_rescued_fasta\": {},", + args.write_rescued_fasta + )?; + writeln!( + w, + " \"search_cache\": {},", + args.search_cache + .as_deref() + .map(|dir| json_escape(&dir.display().to_string())) + .unwrap_or_else(|| "null".to_string()) + )?; writeln!(w, " \"write_html_report\": {},", args.write_html_report)?; if let Some(reference_genome) = args.reference_genome.as_deref() { writeln!( @@ -1763,6 +1926,7 @@ fn count_items_html( fn output_file_description(name: &str) -> &'static str { match name { "gene_presence_absence_recore.csv" => "Corrected binary matrix", + "gene_presence_absence_recore.Rtab" => "Corrected binary matrix as 1/0, for Rtab readers", "gene_presence_absence_recore_states.csv" => { "Multistate matrix: input, corrected, NSC, kept-out and not-searched codes" } @@ -1804,6 +1968,7 @@ pub(super) fn output_file_entries(output_dir: &std::path::Path) -> Vec { let rank = |name: &str| { [ "gene_presence_absence_recore.csv", + "gene_presence_absence_recore.Rtab", "orthology_evidence.tsv", "gene_presence_absence_recore_states.csv", "rescue_report.tsv", diff --git a/src/pipeline/scale_harness.rs b/src/pipeline/scale_harness.rs index 09ad2c4..d55a576 100644 --- a/src/pipeline/scale_harness.rs +++ b/src/pipeline/scale_harness.rs @@ -944,6 +944,8 @@ fn scale_report() { multi_allele: false, synteny: synteny_mode != "none", write_corrected_gff: false, + write_rescued_fasta: false, + search_cache: None, write_html_report: true, matrix_policy: MatrixPolicy::Balanced, reference_genome: reference_name.clone(), diff --git a/src/pipeline/tests.rs b/src/pipeline/tests.rs index 34e6b89..5e63d1e 100644 --- a/src/pipeline/tests.rs +++ b/src/pipeline/tests.rs @@ -16,6 +16,8 @@ fn default_args_for_validation() -> Args { multi_allele: false, synteny: false, write_corrected_gff: false, + write_rescued_fasta: false, + search_cache: None, write_html_report: false, matrix_policy: MatrixPolicy::Balanced, reference_genome: None, @@ -2160,6 +2162,94 @@ fn corrected_matrix_counts_only_what_is_written() { ); assert_eq!(counts.absent_cells, empty); assert_eq!(counts.absent_cells, 2); + + // The Rtab is the same matrix as 1/0: same rows, same cells, the merged + // member left out, the rescue and the NSC cell both 1. + let rtab = + std::fs::read_to_string(dir.path().join("gene_presence_absence_recore.Rtab")).unwrap(); + assert_eq!(rtab, "Gene\tg1\tg2\tg3\nrep\t1\t1\t0\nother\t1\t1\t0\n"); + for (csv_line, rtab_line) in written.lines().skip(1).zip(rtab.lines().skip(1)) { + let from_csv: Vec<&str> = csv_line + .split(',') + .skip(13) + .map(|cell| if cell.is_empty() { "0" } else { "1" }) + .collect(); + let from_rtab: Vec<&str> = rtab_line.split('\t').skip(1).collect(); + assert_eq!(from_csv, from_rtab, "{csv_line}"); + } +} + +/// Rescued genes come out in the gene's own direction, from the phase the +/// alignment fixes, with the IDs of the corrected GFF; a rescue without a +/// single frame gets a nucleotide record and no protein. +#[test] +fn rescued_fasta_reads_each_gene_in_its_direction_and_frame() { + let dir = tempfile::tempdir().unwrap(); + // GTG AAA CCC GGG TAA on the plus strand; ATG TTT TGG TAG on the minus. + let plus = "GTGAAACCCGGGTAA"; + let minus_rc = "CTACCAAAACAT"; + let contig = format!("AAAAA{plus}CCCCC{minus_rc}GGGGG"); + let fasta = dir.path().join("g1.fasta"); + std::fs::write(&fasta, format!(">chr description\n{contig}\n")).unwrap(); + let mut genome_files = HashMap::default(); + genome_files.insert("g1".to_string(), fasta); + let pa = PresenceAbsence::new_for_test( + vec!["c1".into(), "c2".into(), "c3".into(), "c4".into()], + vec!["g1".into()], + vec![vec!["".into()]; 4], + ); + let rescue = + |ci: usize, start: u64, end: u64, strand: char, query_start: Option| RescueEvent { + cluster_idx: ClusterIdx(ci), + start, + end, + strand, + query_start, + ..rescue_with_type(RescueType::MissingCds, None) + }; + let on_plus = rescue(0, 5, 20, '+', Some(0)); + let on_minus = rescue(1, 25, 37, '-', Some(0)); + let mid_codon = rescue(2, 6, 20, '+', Some(1)); + let no_frame = rescue(3, 5, 20, '+', None); + + let (n_nucleotide, n_protein) = write_rescued_fasta( + &[&on_plus, &on_minus, &mid_codon, &no_frame], + &pa, + &genome_files, + dir.path(), + ) + .unwrap(); + assert_eq!((n_nucleotide, n_protein), (4, 3)); + + let fna = std::fs::read_to_string(dir.path().join("rescued_genes.fna")).unwrap(); + let faa = std::fs::read_to_string(dir.path().join("rescued_genes.faa")).unwrap(); + assert!(fna.contains( + ">recore_c1_chr_6_20 genome=g1 gene_cluster=c1 locus=chr:6-20(+) rescue_type=MISSING_CDS phase=0\nGTGAAACCCGGGTAA\n" + )); + // The minus-strand gene is read as the gene, not as the contig. + assert!(fna.contains("locus=chr:26-37(-) rescue_type=MISSING_CDS phase=0\nATGTTTTGGTAG\n")); + assert!(fna.contains("locus=chr:6-20(+) rescue_type=MISSING_CDS phase=unknown\n")); + // GTG opens the gene, so it reads M; the final stop is left off. + assert!( + faa.contains("gene_cluster=c1 locus=chr:6-20(+) rescue_type=MISSING_CDS phase=0\nMKPG\n") + ); + assert!( + faa.contains("gene_cluster=c2 locus=chr:26-37(-) rescue_type=MISSING_CDS phase=0\nMFW\n") + ); + // Starting one base into a codon, the frame is two bases further in. + assert!( + faa.contains("gene_cluster=c3 locus=chr:7-20(+) rescue_type=MISSING_CDS phase=2\nKPG\n") + ); + assert!(!faa.contains("gene_cluster=c4")); +} + +#[test] +fn bacterial_translation_reads_starts_stops_and_ambiguous_codons() { + use crate::rescue::translate_bacterial; + assert_eq!(translate_bacterial(b"TTGAAATGAGGG", true), "MK*G"); + assert_eq!(translate_bacterial(b"TTGAAATAA", false), "LK"); + assert_eq!(translate_bacterial(b"atgNAAgggTAG", true), "MXG"); + assert_eq!(translate_bacterial(b"ATGAA", true), "M"); } #[test] diff --git a/src/relationships.rs b/src/relationships.rs index 9acac29..9d745c5 100644 --- a/src/relationships.rs +++ b/src/relationships.rs @@ -14,7 +14,7 @@ use std::path::Path; use crate::nscluster::NsClusterResult; use crate::parser::PresenceAbsence; -use crate::rescue::{panaroo_gene_id_suffix_variants, RescueEvent, RescueType}; +use crate::rescue::{cell_gene_id_variants, RescueEvent, RescueType}; use crate::types::ClusterIdx; #[derive(Debug, Clone)] @@ -578,7 +578,7 @@ fn build_gene_cluster_map( continue; } map.insert(gene_id.to_string(), ci_val); - for stripped in panaroo_gene_id_suffix_variants(gene_id).into_iter().skip(1) { + for stripped in cell_gene_id_variants(gene_id).into_iter().skip(1) { map.entry(stripped.to_string()).or_insert(ci_val); } } diff --git a/src/rescue.rs b/src/rescue.rs index 8cf5041..58d8b9b 100644 --- a/src/rescue.rs +++ b/src/rescue.rs @@ -12,10 +12,12 @@ mod types; pub use self::conflicts::resolve_conflicts; pub use self::genes::build_gene_to_cluster_with_gff_aliases; -pub(crate) use self::genes::panaroo_gene_id_suffix_variants; +pub(crate) use self::genes::{cell_gene_id_variants, panaroo_gene_id_suffix_variants}; pub use self::gff::{GffGene, GffIndex}; +pub(crate) use self::orf::reverse_complement; pub use self::orf::{ check_internal_stops_in_cds, check_internal_stops_in_hit, check_internal_stops_in_span, + translate_bacterial, }; pub(crate) use self::synteny::{validate_synteny_with_cache, SyntenyReferenceCache}; pub use self::types::{ @@ -27,8 +29,6 @@ use self::genes::build_gene_to_cluster; #[cfg(test)] use self::orf::count_internal_stops; #[cfg(test)] -use self::orf::reverse_complement; -#[cfg(test)] use self::synteny::{ assess_context_support, compatible_reference_support, has_recurrent_context_support, validate_synteny, weighted_side_score, ReferenceSyntenyContext, SyntenyContext, diff --git a/src/rescue/genes.rs b/src/rescue/genes.rs index 675671b..64cf8c7 100644 --- a/src/rescue/genes.rs +++ b/src/rescue/genes.rs @@ -16,8 +16,8 @@ pub fn build_gene_to_cluster(pa: &PresenceAbsence, gi: GenomeIdx) -> HashMap Vec<&str> { variants } +/// Every name a gene ID from a matrix cell may go by in the GFF, the ID itself +/// first. Then Panaroo's `_pseudo` and `_dup_N` suffixes stripped; the loci of +/// a PIRATE fission, which it writes `(locus1:locus2)`; and the parts after +/// each `-` of a PanTA ID, which is `sample-contig-gene`. The cells keep the +/// tool's own IDs, and a later variant only answers when the ones before it do +/// not, so an exact match always wins. Without the PanTA parts no cell gene +/// matched a GFF, and synteny validated nothing on PanTA input. +pub(crate) fn cell_gene_id_variants(gene_id: &str) -> Vec<&str> { + let mut variants = panaroo_gene_id_suffix_variants(gene_id); + if let Some(inner) = gene_id + .strip_prefix('(') + .and_then(|rest| rest.strip_suffix(')')) + { + for locus in inner.split(':').map(str::trim) { + if !locus.is_empty() && !variants.contains(&locus) { + variants.push(locus); + } + } + } + let mut rest = gene_id; + while let Some((_, tail)) = rest.split_once('-') { + if !tail.is_empty() && !variants.contains(&tail) { + variants.push(tail); + } + rest = tail; + } + variants +} + fn gene_id_variants(gene_id: &str) -> impl Iterator { panaroo_gene_id_suffix_variants(gene_id).into_iter() } diff --git a/src/rescue/orf.rs b/src/rescue/orf.rs index b3567f0..275a82e 100644 --- a/src/rescue/orf.rs +++ b/src/rescue/orf.rs @@ -128,8 +128,51 @@ pub fn check_internal_stops_in_span( (n_stops > 0).then_some(n_stops) } +/// Translate a coding sequence with the bacterial code (NCBI table 11), from +/// its first base, one amino acid per whole codon, with `*` for stops and `X` +/// for a codon with an ambiguous base. A final stop codon is left off. When +/// `from_start_codon`, the first codon is the gene's start and reads as `M` +/// whatever table 11 start it is (GTG and TTG included), as a CDS does. +pub fn translate_bacterial(seq: &[u8], from_start_codon: bool) -> String { + const STARTS: [&[u8; 3]; 7] = [b"ATG", b"GTG", b"TTG", b"CTG", b"ATT", b"ATC", b"ATA"]; + let mut protein = String::with_capacity(seq.len() / 3); + for (i, codon) in seq.chunks_exact(3).enumerate() { + let upper = [ + codon[0].to_ascii_uppercase(), + codon[1].to_ascii_uppercase(), + codon[2].to_ascii_uppercase(), + ]; + if i == 0 && from_start_codon && STARTS.contains(&&upper) { + protein.push('M'); + continue; + } + protein.push(codon_amino_acid(upper)); + } + if protein.ends_with('*') { + protein.pop(); + } + protein +} + +/// The standard genetic code, which table 11 shares outside start codons. +fn codon_amino_acid(codon: [u8; 3]) -> char { + let index = |base: u8| match base { + b'T' | b'U' => Some(0), + b'C' => Some(1), + b'A' => Some(2), + b'G' => Some(3), + _ => None, + }; + // Codons in TCAG order: first base slowest, third fastest. + const TABLE: &[u8; 64] = b"FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"; + match (index(codon[0]), index(codon[1]), index(codon[2])) { + (Some(a), Some(b), Some(c)) => TABLE[a * 16 + b * 4 + c] as char, + _ => 'X', + } +} + /// Reverse complement a DNA sequence. -pub(super) fn reverse_complement(seq: &[u8]) -> Vec { +pub(crate) fn reverse_complement(seq: &[u8]) -> Vec { seq.iter() .rev() .map(|&b| match b { diff --git a/src/rescue/synteny.rs b/src/rescue/synteny.rs index 6730153..91a9bb8 100644 --- a/src/rescue/synteny.rs +++ b/src/rescue/synteny.rs @@ -1,4 +1,4 @@ -use super::genes::panaroo_gene_id_suffix_variants; +use super::genes::cell_gene_id_variants; use super::gff::GffIndex; use super::types::{RescueEvent, SyntenyEvidence}; use crate::parser::PresenceAbsence; @@ -79,7 +79,7 @@ fn cluster_reference_located( .filter(|id| !id.is_empty()) .any(|id| { ref_gff.get_by_id(id).is_some() - || panaroo_gene_id_suffix_variants(id) + || cell_gene_id_variants(id) .into_iter() .skip(1) .any(|variant| ref_gff.get_by_id(variant).is_some()) @@ -309,7 +309,7 @@ fn reference_contexts_for_cluster<'a>( .filter(|id| !id.is_empty()) { let ref_gene = ref_gff.get_by_id(ref_gene_id).or_else(|| { - panaroo_gene_id_suffix_variants(ref_gene_id) + cell_gene_id_variants(ref_gene_id) .into_iter() .skip(1) .find_map(|id| ref_gff.get_by_id(id)) diff --git a/src/rescue/tests.rs b/src/rescue/tests.rs index b18c1ee..1004dc0 100644 --- a/src/rescue/tests.rs +++ b/src/rescue/tests.rs @@ -2626,3 +2626,35 @@ fn pieces_at_both_ends_in_linear_order_are_not_an_origin_wrap() { "{rescue:?}" ); } + +/// A cell's gene ID reaches the GFF whichever tool wrote it: the ID itself +/// first, then Panaroo's suffixes stripped, a PIRATE fission's pieces, and the +/// gene part of a PanTA `sample-contig-gene` ID. +#[test] +fn cell_gene_ids_reach_the_gff_from_every_tool() { + use crate::rescue::cell_gene_id_variants; + assert_eq!( + cell_gene_id_variants("X_dup_2_pseudo"), + vec!["X_dup_2_pseudo", "X_dup_2", "X"] + ); + assert_eq!( + cell_gene_id_variants("(L8_AS00801_00013a:L8_AS00801_00013b)"), + vec![ + "(L8_AS00801_00013a:L8_AS00801_00013b)", + "L8_AS00801_00013a", + "L8_AS00801_00013b" + ] + ); + let panta = cell_gene_id_variants("L1_AS00101-L1_AS00101_contig1-L1_AS00101_00019"); + assert_eq!(panta[0], "L1_AS00101-L1_AS00101_contig1-L1_AS00101_00019"); + assert!(panta.contains(&"L1_AS00101_00019")); + + // In a genome's lookup, a PanTA cell answers for the GFF's own ID. + let pa = crate::parser::PresenceAbsence::new_for_test( + vec!["Rv0019c".into()], + vec!["L1_AS00101".into()], + vec![vec!["L1_AS00101-L1_AS00101_contig1-L1_AS00101_00019".into()]], + ); + let map = build_gene_to_cluster(&pa, GenomeIdx(0)); + assert_eq!(map.get("L1_AS00101_00019"), Some(&ClusterIdx(0))); +} diff --git a/src/search.rs b/src/search.rs index de8cafa..526907e 100644 --- a/src/search.rs +++ b/src/search.rs @@ -82,7 +82,9 @@ const END_BONUS: &str = "40"; /// fails here with its path instead of producing hits on the wrong contig. /// /// `max_secondary` bounds the secondary alignments per query within one -/// genome. Returns hits keyed by (genome, cluster), only for the requested gaps. +/// genome. With `cache_dir`, results are kept and reused between runs (see +/// [`SearchCache`]). Returns hits keyed by (genome, cluster), only for the +/// requested gaps. pub fn batch_search( queries: &HashMap>>, genome_files: &HashMap, @@ -94,6 +96,7 @@ pub fn batch_search( preset: &str, exact_base_alignment: bool, max_secondary: usize, + cache_dir: Option<&std::path::Path>, ) -> Result>> { // One job per genome. A genome listed twice has its gap lists merged, and // the BTreeMap keeps the job order, and so the log, deterministic. @@ -130,12 +133,20 @@ pub fn batch_search( jobs.len() ); - let max_secondary = max_secondary.max(1).to_string(); + let settings = SearchSettings { + min_identity, + preset, + exact_base_alignment, + max_secondary: max_secondary.max(1).to_string(), + }; + let cache = cache_dir + .map(|dir| SearchCache::open(dir, &settings)) + .transpose()?; let pool = rayon::ThreadPoolBuilder::new() .num_threads(threads.max(1)) .build() .context("Failed to build the search thread pool")?; - let results: Vec>> = pool.install(|| { + let results: Vec> = pool.install(|| { use rayon::prelude::*; jobs.par_iter() .map(|(_, path, clusters)| { @@ -144,10 +155,8 @@ pub fn batch_search( path, clusters, cluster_names, - min_identity, - preset, - exact_base_alignment, - &max_secondary, + &settings, + cache.as_ref(), ) }) .collect() @@ -155,19 +164,37 @@ pub fn batch_search( let mut hits: HashMap<(GenomeIdx, ClusterIdx), Vec> = HashMap::default(); let mut kept = 0usize; + let mut cache_use = CacheUse::default(); + let mut genomes_without_minimap2 = 0usize; for ((gi, path, _), result) in jobs.iter().zip(results) { - let found = result.with_context(|| { + let (found, used) = result.with_context(|| { format!( "Search failed for genome {} ({})", genome_names[gi.val()], path.display() ) })?; + if used.searched == 0 && used.reused > 0 { + genomes_without_minimap2 += 1; + } + cache_use.reused += used.reused; + cache_use.searched += used.searched; kept += found.len(); for (ci, hit) in found { hits.entry((*gi, ci)).or_default().push(hit); } } + if let Some(cache) = cache.as_ref() { + info!( + "Search cache {}: {} of {} query sequences reused, {} searched; {} of {} genomes needed no minimap2 run", + cache.dir.display(), + cache_use.reused, + cache_use.reused + cache_use.searched, + cache_use.searched, + genomes_without_minimap2, + jobs.len() + ); + } debug!( "minimap2: {} hits kept across {} (genome, cluster) pairs", kept, @@ -176,46 +203,182 @@ pub fn batch_search( Ok(hits) } +/// The settings that reach minimap2 or filter what it returns. +struct SearchSettings<'a> { + min_identity: f64, + preset: &'a str, + exact_base_alignment: bool, + max_secondary: String, +} + +impl SearchSettings<'_> { + fn minimap2_args(&self) -> Vec<&str> { + let mut args: Vec<&str> = vec!["-x", self.preset, "-t", "1"]; + if self.exact_base_alignment { + args.extend(["-c", "-s", MIN_PEAK_SCORE, "--end-bonus", END_BONUS]); + } + args.extend(["--secondary=yes", "-N", &self.max_secondary, "-p", "0.5"]); + args + } +} + +/// One sequence to align against a genome: a cluster's representative, or one +/// of its alleles. +struct QueryItem<'a> { + name: String, + ci: ClusterIdx, + seq: &'a [u8], +} + +/// A query's hits from its alignment, and from realigning what that alignment +/// left uncovered, each in the order minimap2 reported them. +#[derive(Debug, Clone, Default)] +struct ItemHits { + first: Vec, + remainder: Vec, +} + +impl ItemHits { + /// The same hits under another query's name, for a sequence shared by two + /// queries or read back from the cache. + fn named(mut self, query: &str, cluster: &str) -> Self { + for hit in self.first.iter_mut().chain(self.remainder.iter_mut()) { + hit.query = query.to_string(); + hit.cluster = cluster.to_string(); + } + self + } +} + +/// One genome's hits, by cluster, and how much of its search the cache spared. +type GenomeSearch = (Vec<(ClusterIdx, Hit)>, CacheUse); + +#[derive(Debug, Clone, Copy, Default)] +struct CacheUse { + reused: usize, + searched: usize, +} + /// Align one genome's absent clusters against that genome alone. -#[allow(clippy::too_many_arguments)] fn search_one_genome( queries: &HashMap>>, fasta_path: &std::path::Path, clusters: &[ClusterIdx], cluster_names: &[String], - min_identity: f64, - preset: &str, - exact_base_alignment: bool, - max_secondary: &str, -) -> Result> { + settings: &SearchSettings<'_>, + cache: Option<&SearchCache>, +) -> Result { + let mut items: Vec = Vec::new(); + for &ci in clusters { + let cluster_name = &cluster_names[ci.val()]; + let Some(alleles) = queries.get(cluster_name) else { + continue; + }; + for (ai, seq) in alleles.iter().enumerate() { + let name = if alleles.len() == 1 { + cluster_name.clone() + } else { + format!("{}~~~ALLELE~{}", cluster_name, ai) + }; + items.push(QueryItem { + name, + ci, + seq: seq.as_slice(), + }); + } + } + if items.is_empty() { + return Ok((Vec::new(), CacheUse::default())); + } + + let Some(cache) = cache else { + let all: Vec<&QueryItem> = items.iter().collect(); + let per_item = search_items(&all, fasta_path, cluster_names, settings)?; + return Ok((in_search_order(&items, &per_item), CacheUse::default())); + }; + + let path = cache.file_for(fasta_path)?; + let mut stored = cache.load(&path); + let hashes: Vec = items.iter().map(|item| stable_hash(&[item.seq])).collect(); + // Each sequence not in the cache is searched once, however many queries + // share it: identical sequences align identically. + let mut pending: Vec = Vec::new(); + let mut pending_hashes: ahash::HashSet = ahash::HashSet::default(); + for (k, hash) in hashes.iter().enumerate() { + if !stored.contains_key(hash) && pending_hashes.insert(*hash) { + pending.push(k); + } + } + if !pending.is_empty() { + let to_search: Vec<&QueryItem> = pending.iter().map(|&k| &items[k]).collect(); + let found = search_items(&to_search, fasta_path, cluster_names, settings)?; + for (&k, hits) in pending.iter().zip(found) { + stored.insert(hashes[k], hits); + } + cache.store(&path, &stored)?; + } + let per_item: Vec = items + .iter() + .zip(&hashes) + .map(|(item, hash)| { + stored + .get(hash) + .cloned() + .unwrap_or_default() + .named(&item.name, &cluster_names[item.ci.val()]) + }) + .collect(); + let used = CacheUse { + reused: items.len() - pending.len(), + searched: pending.len(), + }; + Ok((in_search_order(&items, &per_item), used)) +} + +/// Hits in the order a search of every query at once reports them: first +/// alignments in query order, then the realigned remainders by query name. A +/// cached run hands its hits on in this order too, so it decides exactly what +/// an uncached run decides. +fn in_search_order(items: &[QueryItem], per_item: &[ItemHits]) -> Vec<(ClusterIdx, Hit)> { + let mut out = Vec::new(); + for (item, hits) in items.iter().zip(per_item) { + out.extend(hits.first.iter().cloned().map(|hit| (item.ci, hit))); + } + let mut by_name: Vec = (0..items.len()).collect(); + by_name.sort_by(|&a, &b| items[a].name.cmp(&items[b].name)); + for k in by_name { + out.extend( + per_item[k] + .remainder + .iter() + .cloned() + .map(|hit| (items[k].ci, hit)), + ); + } + out +} + +/// Run minimap2 for these queries against one genome. Returns each query's +/// hits, aligned with `items`. +fn search_items( + items: &[&QueryItem], + fasta_path: &std::path::Path, + cluster_names: &[String], + settings: &SearchSettings<'_>, +) -> Result> { + let mut per_item = vec![ItemHits::default(); items.len()]; let mut query_file = NamedTempFile::new().context("Cannot create temp query file")?; - let mut query_to_cluster_idx: HashMap = HashMap::default(); - let mut query_seqs: HashMap = HashMap::default(); + let mut index_of: HashMap<&str, usize> = HashMap::default(); { let mut writer = BufWriter::new(&mut query_file); - for &ci in clusters { - let cluster_name = &cluster_names[ci.val()]; - let Some(alleles) = queries.get(cluster_name) else { - continue; - }; - for (ai, seq) in alleles.iter().enumerate() { - let query_name = if alleles.len() == 1 { - cluster_name.clone() - } else { - format!("{}~~~ALLELE~{}", cluster_name, ai) - }; - writeln!(writer, ">{}", query_name)?; - writer.write_all(seq)?; - writeln!(writer)?; - query_seqs.insert(query_name.clone(), seq.as_slice()); - query_to_cluster_idx.insert(query_name, ci); - } + for (k, item) in items.iter().enumerate() { + writeln!(writer, ">{}", item.name)?; + writer.write_all(item.seq)?; + writeln!(writer)?; + index_of.insert(item.name.as_str(), k); } writer.flush()?; } - if query_to_cluster_idx.is_empty() { - return Ok(Vec::new()); - } let mut target_file = NamedTempFile::new().context("Cannot create temp target file")?; { @@ -228,29 +391,23 @@ fn search_one_genome( .context("Cannot write the temporary search target")?; } - let mut base_args: Vec<&str> = vec!["-x", preset, "-t", "1"]; - if exact_base_alignment { - base_args.extend(["-c", "-s", MIN_PEAK_SCORE, "--end-bonus", END_BONUS]); - } - base_args.extend(["--secondary=yes", "-N", max_secondary, "-p", "0.5"]); + let base_args = settings.minimap2_args(); let target_path = target_file.path().to_string_lossy().to_string(); let paf = run_minimap2(&base_args, &target_path, query_file.path(), fasta_path)?; - let mut hits = Vec::new(); for line in paf.lines() { let Some(fields) = paf_required_fields(line) else { continue; }; - let Some(&ci) = query_to_cluster_idx.get(fields[0]) else { + let Some(&k) = index_of.get(fields[0]) else { continue; }; - let Some(hit) = - hit_from_paf_fields(&fields, fields[0], &cluster_names[ci.val()], fields[5]) - else { + let cluster = &cluster_names[items[k].ci.val()]; + let Some(hit) = hit_from_paf_fields(&fields, fields[0], cluster, fields[5]) else { continue; }; - if hit.identity() >= min_identity { - hits.push((ci, hit)); + if hit.identity() >= settings.min_identity { + per_item[k].first.push(hit); } } @@ -259,36 +416,36 @@ fn search_one_genome( // alignment then stops at the insertion and drops the rest of the chain, // so only one half was reported and the gene could not be seen as // FRAGMENTED. Aligned on its own, the missing stretch finds its half. - if exact_base_alignment { + if settings.exact_base_alignment { let mut remainder_file = NamedTempFile::new().context("Cannot create temp remainder query file")?; - let mut remainders: HashMap = HashMap::default(); + let mut remainders: HashMap = HashMap::default(); { let mut writer = BufWriter::new(&mut remainder_file); // Sorted, so the remainder queries and minimap2's output come in // the same order on every run. - let mut ordered: Vec<(&String, &&[u8])> = query_seqs.iter().collect(); - ordered.sort_unstable_by(|a, b| a.0.cmp(b.0)); - for (query_name, seq) in ordered { - let intervals: Vec<(u64, u64)> = hits + let mut ordered: Vec = (0..items.len()).collect(); + ordered.sort_unstable_by(|&a, &b| items[a].name.cmp(&items[b].name)); + for k in ordered { + let intervals: Vec<(u64, u64)> = per_item[k] + .first .iter() - .filter(|(_, h)| h.query == *query_name) - .map(|(_, h)| (h.query_start, h.query_end)) + .map(|h| (h.query_start, h.query_end)) .collect(); if intervals.is_empty() { continue; } + let seq = items[k].seq; let full_len = seq.len() as u64; for (from, to) in uncovered_stretches(intervals, full_len) { if to - from < REMAINDER_MIN_BP { continue; } - let name = format!("{}~~~REM~{}", query_name, from); + let name = format!("{}~~~REM~{}", items[k].name, from); writeln!(writer, ">{}", name)?; writer.write_all(&seq[from as usize..to as usize])?; writeln!(writer)?; - let ci = query_to_cluster_idx[query_name]; - remainders.insert(name, (query_name.clone(), ci, from, full_len)); + remainders.insert(name, (k, from, full_len)); } } writer.flush()?; @@ -299,24 +456,255 @@ fn search_one_genome( let Some(fields) = paf_required_fields(line) else { continue; }; - let Some((original, ci, offset, full_len)) = remainders.get(fields[0]) else { + let Some(&(k, offset, full_len)) = remainders.get(fields[0]) else { continue; }; + let cluster = &cluster_names[items[k].ci.val()]; let Some(mut hit) = - hit_from_paf_fields(&fields, original, &cluster_names[ci.val()], fields[5]) + hit_from_paf_fields(&fields, &items[k].name, cluster, fields[5]) else { continue; }; hit.query_start += offset; hit.query_end += offset; - hit.query_length = *full_len; - if hit.identity() >= min_identity { - hits.push((*ci, hit)); + hit.query_length = full_len; + if hit.identity() >= settings.min_identity { + per_item[k].remainder.push(hit); } } } } - Ok(hits) + Ok(per_item) +} + +/// minimap2 results kept between runs: one file per genome and search setting, +/// holding the final hits of every query sequence searched in that genome. +/// +/// A search is a function of the genome, the query sequence, the settings that +/// reach minimap2 or filter its output (preset, base-level alignment, the +/// secondary budget, and --min-identity, which also decides what the second +/// pass realigns), and of minimap2 and recore themselves. A file is named after +/// all of those but the query; inside, each query is keyed by its sequence's +/// own hash, and a sequence that found nothing is recorded as searched. A rerun +/// that changes only what happens after the search (the matrix policy, the +/// coverage threshold, synteny) or which gaps are searched (the prevalence +/// floor) aligns only the sequences it has not seen; anything that could change +/// a hit names another file. +pub struct SearchCache { + dir: std::path::PathBuf, + settings: String, +} + +const SEARCH_CACHE_FORMAT: &str = "recore-search-cache-v1"; +const CACHE_FILE_END: &str = "# end"; + +impl SearchCache { + fn open(dir: &std::path::Path, settings: &SearchSettings<'_>) -> Result { + std::fs::create_dir_all(dir) + .with_context(|| format!("Cannot create search cache directory {}", dir.display()))?; + let settings = format!( + "{} recore={} minimap2={} preset={} exact={} secondary={} min_identity_bits={:016x} peak={} end_bonus={} remainder_min={}", + SEARCH_CACHE_FORMAT, + env!("CARGO_PKG_VERSION"), + minimap2_version()?, + settings.preset, + settings.exact_base_alignment, + settings.max_secondary, + settings.min_identity.to_bits(), + MIN_PEAK_SCORE, + END_BONUS, + REMAINDER_MIN_BP + ); + Ok(Self { + dir: dir.to_path_buf(), + settings, + }) + } + + /// The cache file for a genome, named after its bytes and the settings. + fn file_for(&self, fasta_path: &std::path::Path) -> Result { + use std::io::Read; + let mut file = std::fs::File::open(fasta_path) + .with_context(|| format!("Cannot open genome {}", fasta_path.display()))?; + let mut hasher = Fnv128::new(); + let mut buf = vec![0u8; 1 << 20]; + loop { + let n = file + .read(&mut buf) + .with_context(|| format!("Cannot read genome {}", fasta_path.display()))?; + if n == 0 { + break; + } + hasher.write(&buf[..n]); + } + let genome = hasher.finish(); + let name = stable_hash(&[&genome.to_le_bytes(), self.settings.as_bytes()]); + Ok(self.dir.join(format!("{name:032x}.tsv"))) + } + + /// The stored hits by sequence hash. A file that cannot be read back is + /// treated as empty, so its genome is searched again and the file rewritten. + fn load(&self, path: &std::path::Path) -> HashMap { + let Ok(text) = std::fs::read_to_string(path) else { + return HashMap::default(); + }; + match parse_cache_file(&text, &self.settings) { + Some(stored) => stored, + None => { + log::warn!( + "Search cache file {} is unreadable or from other settings; searching that genome again", + path.display() + ); + HashMap::default() + } + } + } + + /// Replace the file in one step, so a run that stops halfway, or another + /// run on the same cache, never sees half a file. + fn store(&self, path: &std::path::Path, stored: &HashMap) -> Result<()> { + let mut file = NamedTempFile::new_in(&self.dir) + .with_context(|| format!("Cannot write to search cache {}", self.dir.display()))?; + { + let mut w = BufWriter::new(&mut file); + writeln!(w, "# {}", self.settings)?; + let mut hashes: Vec<&u128> = stored.keys().collect(); + hashes.sort_unstable(); + for hash in hashes { + writeln!(w, "S\t{hash:032x}")?; + let hits = &stored[hash]; + for (pass, list) in [("1", &hits.first), ("2", &hits.remainder)] { + for h in list { + writeln!( + w, + "{pass}\t{hash:032x}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}", + h.query_length, + h.query_start, + h.query_end, + h.strand, + h.target, + h.target_start, + h.target_end, + h.target_length, + h.matches, + h.alignment_length, + h.mapq + )?; + } + } + } + // A file without this last line was cut short. + writeln!(w, "{CACHE_FILE_END}")?; + w.flush()?; + } + file.persist(path) + .with_context(|| format!("Cannot write search cache file {}", path.display()))?; + Ok(()) + } +} + +/// A cache file's hits by sequence hash, or None when any line is malformed, +/// the file was cut short, or it was written under other settings. Cut short +/// mid-line, a number could otherwise read back as a smaller one. +fn parse_cache_file(text: &str, settings: &str) -> Option> { + let body = text.strip_suffix('\n')?.strip_suffix(CACHE_FILE_END)?; + let mut lines = body.lines(); + if lines.next()?.strip_prefix("# ")? != settings { + return None; + } + let mut stored: HashMap = HashMap::default(); + for line in lines { + let fields: Vec<&str> = line.split('\t').collect(); + let hash = u128::from_str_radix(fields.get(1)?, 16).ok()?; + match fields[0] { + "S" if fields.len() == 2 => { + stored.entry(hash).or_default(); + } + pass @ ("1" | "2") if fields.len() == 13 => { + let mut chars = fields[5].chars(); + let strand = chars.next()?; + if chars.next().is_some() { + return None; + } + let hit = Hit { + query: String::new(), + cluster: String::new(), + query_length: fields[2].parse().ok()?, + query_start: fields[3].parse().ok()?, + query_end: fields[4].parse().ok()?, + strand, + target: fields[6].to_string(), + target_start: fields[7].parse().ok()?, + target_end: fields[8].parse().ok()?, + target_length: fields[9].parse().ok()?, + matches: fields[10].parse().ok()?, + alignment_length: fields[11].parse().ok()?, + mapq: fields[12].parse().ok()?, + }; + let entry = stored.get_mut(&hash)?; + if pass == "1" { + entry.first.push(hit); + } else { + entry.remainder.push(hit); + } + } + _ => return None, + } + } + Some(stored) +} + +/// minimap2's version, asked once per run: its hits can change between versions. +fn minimap2_version() -> Result { + static VERSION: std::sync::OnceLock = std::sync::OnceLock::new(); + if let Some(version) = VERSION.get() { + return Ok(version.clone()); + } + let output = Command::new("minimap2") + .arg("--version") + .stdin(Stdio::null()) + .output() + .context("Failed to run minimap2: is it installed and on PATH?")?; + let version = String::from_utf8_lossy(&output.stdout).trim().to_string(); + if !output.status.success() || version.is_empty() { + anyhow::bail!("minimap2 --version did not report a version"); + } + Ok(VERSION.get_or_init(|| version).clone()) +} + +/// 128-bit FNV-1a: stable across platforms and Rust versions, which the +/// standard library's hasher does not promise, and that is what a file name +/// read back by a later run needs. +struct Fnv128(u128); + +impl Fnv128 { + fn new() -> Self { + Self(0x6c62_272e_07bb_0142_62b8_2175_6295_c58d) + } + + fn write(&mut self, bytes: &[u8]) { + for &byte in bytes { + self.0 ^= u128::from(byte); + self.0 = self + .0 + .wrapping_mul(0x0000_0000_0100_0000_0000_0000_0000_013b); + } + } + + fn finish(&self) -> u128 { + self.0 + } +} + +/// Hash of the parts in order, each followed by its length so that a boundary +/// cannot move between parts. +fn stable_hash(parts: &[&[u8]]) -> u128 { + let mut hasher = Fnv128::new(); + for part in parts { + hasher.write(part); + hasher.write(&(part.len() as u64).to_le_bytes()); + } + hasher.finish() } /// A stretch of a query this short is not worth a second alignment: it is an @@ -531,6 +919,138 @@ mod tests { use super::*; use ahash::HashMapExt; + fn cached_hit(query_start: u64, target: &str) -> Hit { + Hit { + query: String::new(), + cluster: String::new(), + query_length: 900, + query_start, + query_end: 900, + strand: '-', + target: target.into(), + target_start: 10, + target_end: 910, + target_length: 5000, + matches: 899, + alignment_length: 900, + mapq: 60, + } + } + + /// FNV-1a 128 against its published test vectors: a cache file name has to + /// be the same on every platform and in every later build. + #[test] + fn fnv128_matches_the_reference_vectors() { + let fnv = |bytes: &[u8]| { + let mut h = Fnv128::new(); + h.write(bytes); + h.finish() + }; + assert_eq!(fnv(b""), 0x6c62272e07bb014262b821756295c58d); + assert_eq!(fnv(b"a"), 0xd228cb696f1a8caf78912b704e4a8964); + // A boundary between parts cannot move without changing the hash. + assert_ne!(stable_hash(&[b"ab", b"c"]), stable_hash(&[b"a", b"bc"])); + } + + /// A cache file reads back exactly what was stored, sequences with no hit + /// included, and anything cut short or from other settings reads as nothing. + #[test] + fn cache_file_round_trips_and_refuses_damage() { + let dir = tempfile::tempdir().unwrap(); + let cache = SearchCache { + dir: dir.path().to_path_buf(), + settings: "recore-search-cache-v1 test".to_string(), + }; + let mut stored: HashMap = HashMap::new(); + stored.insert( + 7, + ItemHits { + first: vec![cached_hit(0, "ctg1")], + remainder: vec![cached_hit(450, "ctg2")], + }, + ); + stored.insert(9, ItemHits::default()); + let path = dir.path().join("genome.tsv"); + cache.store(&path, &stored).unwrap(); + + let back = cache.load(&path); + assert_eq!(back.len(), 2); + assert!(back[&9].first.is_empty() && back[&9].remainder.is_empty()); + let hits = &back[&7]; + assert_eq!(hits.first.len(), 1); + assert_eq!(hits.remainder.len(), 1); + let (a, b) = (&hits.first[0], &stored[&7].first[0]); + assert_eq!( + ( + a.query_start, + a.strand, + a.target.as_str(), + a.target_length, + a.matches, + a.mapq + ), + ( + b.query_start, + b.strand, + b.target.as_str(), + b.target_length, + b.matches, + b.mapq + ) + ); + assert_eq!(hits.remainder[0].query_start, 450); + + let text = std::fs::read_to_string(&path).unwrap(); + // Cut short anywhere, even on a line boundary, it is refused. + assert!(parse_cache_file(&text[..text.len() / 2], &cache.settings).is_none()); + let without_end = text.trim_end_matches("# end\n"); + assert!(parse_cache_file(without_end, &cache.settings).is_none()); + // Written under other settings, it is not this cache's. + assert!(parse_cache_file(&text, "recore-search-cache-v1 other").is_none()); + // A hit for a sequence never marked as searched is damage too. + let orphan = text.replace("S\t00000000000000000000000000000007\n", ""); + assert_ne!(orphan, text); + assert!(parse_cache_file(&orphan, &cache.settings).is_none()); + } + + /// Cached or not, hits come out in the order one search of every query + /// reports them: first alignments in query order, remainders by name. + #[test] + fn hits_come_out_in_search_order() { + let seq = b"ACGT".to_vec(); + let items = vec![ + QueryItem { + name: "b".into(), + ci: ClusterIdx(0), + seq: &seq, + }, + QueryItem { + name: "a".into(), + ci: ClusterIdx(1), + seq: &seq, + }, + ]; + let named = |start| Hit { + query: String::new(), + ..cached_hit(start, "t") + }; + let per_item = vec![ + ItemHits { + first: vec![named(1)], + remainder: vec![named(2)], + }, + ItemHits { + first: vec![named(3)], + remainder: vec![named(4)], + }, + ]; + let order: Vec<(usize, u64)> = in_search_order(&items, &per_item) + .into_iter() + .map(|(ci, hit)| (ci.val(), hit.query_start)) + .collect(); + assert_eq!(order, vec![(0, 1), (1, 3), (1, 4), (0, 2)]); + } + // ── Hit method tests ── #[test] diff --git a/tests/integration.rs b/tests/integration.rs index 7357f6c..dccfeb6 100644 --- a/tests/integration.rs +++ b/tests/integration.rs @@ -673,3 +673,158 @@ fn test_unsearched_gaps_say_why_instead_of_absent() { assert!(!gc[column("mechanism_detail")].contains("absences_not_searched")); assert!(!gc[column("recommended_wording")].contains("never searched")); } + +/// Every output file except the two that carry a timestamp. +fn same_outputs(a: &Path, b: &Path) { + let mut names: Vec = fs::read_dir(a) + .unwrap() + .map(|e| e.unwrap().file_name().to_string_lossy().to_string()) + .filter(|n| n != "run_manifest.json" && n != "recore_report.html") + .collect(); + names.sort(); + assert!(!names.is_empty()); + for name in names { + let (pa, pb) = (a.join(&name), b.join(&name)); + if pa.is_dir() { + continue; + } + assert_eq!( + fs::read(&pa).unwrap(), + fs::read(&pb).unwrap(), + "{name} differs between {} and {}", + a.display(), + b.display() + ); + } +} + +/// --search-cache keeps minimap2's results: a second run aligns nothing, a run +/// with another matrix policy reuses the same search, and every one of them +/// writes exactly what the same run without a cache writes. +#[test] +fn test_search_cache_reuses_results_without_changing_them() { + if !minimap2_available() { + assert!( + std::env::var_os("RECORE_REQUIRE_MINIMAP2").is_none(), + "RECORE_REQUIRE_MINIMAP2 is set but minimap2 is not in PATH" + ); + eprintln!("SKIPPED: minimap2 not in PATH"); + return; + } + let tmp = tempdir().unwrap(); + let base = tmp.path(); + let mut rng = Lcg(23); + // gB is in the matrix for s1 and s2 only, but s3 carries it; gC is missing + // from s4 for real. Three gaps to search: gB in s3 and s4, gC in s4. + let genomes = ["s1", "s2", "s3", "s4"]; + let clusters = ["gA", "gB", "gC"]; + let in_matrix = |c: &str, g: &str| match c { + "gB" => g == "s1" || g == "s2", + "gC" => g != "s4", + _ => true, + }; + let carried = |c: &str, g: &str| in_matrix(c, g) || (c == "gB" && g == "s3"); + let seqs: Vec = clusters.iter().map(|_| rng.orf(150)).collect(); + let fasta: String = clusters + .iter() + .zip(&seqs) + .map(|(name, seq)| format!(">{name}\n{seq}\n")) + .collect(); + write_file(base, "pan_ref.fa", &fasta); + let mut csv = String::from( + "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,s1,s2,s3,s4\n", + ); + for name in clusters { + let cells: Vec = genomes + .iter() + .map(|g| { + if in_matrix(name, g) { + format!("{name}_{g}") + } else { + String::new() + } + }) + .collect(); + csv.push_str(&format!( + "{name},{name},product,0,0,0,0,0,0,,450,450,450,{}\n", + cells.join(",") + )); + } + write_file(base, "pa.csv", &csv); + let genomes_dir = base.join("genomes"); + for g in genomes { + let mut contig = rng.seq(300); + for (name, seq) in clusters.iter().zip(&seqs) { + if carried(name, g) { + contig.push_str(seq); + contig.push_str(&rng.seq(300)); + } + } + write_file( + &genomes_dir, + &format!("{g}.fasta"), + &format!(">{g}_1\n{contig}\n"), + ); + } + let cache = base.join("cache"); + let run = |out: &str, extra: &[&str]| -> String { + let output = std::process::Command::new(env!("CARGO_BIN_EXE_recore")) + .args([ + "-p", + base.join("pa.csv").to_str().unwrap(), + "-r", + base.join("pan_ref.fa").to_str().unwrap(), + "-g", + genomes_dir.to_str().unwrap(), + "-o", + base.join(out).to_str().unwrap(), + "-t", + "2", + ]) + .args(extra) + .output() + .expect("Failed to run recore binary"); + let log = String::from_utf8_lossy(&output.stderr).to_string(); + assert!(output.status.success(), "{log}"); + log + }; + let cache_arg = cache.to_str().unwrap(); + run("plain", &[]); + let cold = run("cold", &["--search-cache", cache_arg]); + let warm = run("warm", &["--search-cache", cache_arg]); + assert!( + cold.contains("0 of 3 query sequences reused, 3 searched"), + "{cold}" + ); + assert!( + warm.contains("3 of 3 query sequences reused, 0 searched"), + "{warm}" + ); + // The gaps sit in s3 and s4, and neither needed minimap2 the second time. + assert!( + warm.contains("2 of 2 genomes needed no minimap2 run"), + "{warm}" + ); + same_outputs(&base.join("plain"), &base.join("cold")); + same_outputs(&base.join("plain"), &base.join("warm")); + // gB in s3 is the rescue the cache has to hand back intact. + let corrected = + fs::read_to_string(base.join("warm").join("gene_presence_absence_recore.csv")).unwrap(); + assert!(corrected + .lines() + .any(|l| l.starts_with("gB,") && l.contains("RESCUED:"))); + + let strict_plain = run("strict_plain", &["--matrix-policy", "strict"]); + let strict_warm = run( + "strict_warm", + &["--matrix-policy", "strict", "--search-cache", cache_arg], + ); + assert!(!strict_plain.contains("Search cache")); + assert!( + strict_warm.contains("3 of 3 query sequences reused, 0 searched"), + "{strict_warm}" + ); + same_outputs(&base.join("strict_plain"), &base.join("strict_warm")); +} From 25eb06744fec2d8c1104bed69fcd364f8ddc4e86 Mon Sep 17 00:00:00 2001 From: Paula Ruiz-Rodriguez <50167687+Paururo@users.noreply.github.com> Date: Sat, 17 Oct 2026 19:06:41 +0200 Subject: [PATCH 2/2] docs: regenerate the example report from the tutorial run Its list of output files now includes the corrected Rtab. --- docs/assets/example-report.html | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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);}