Detect pangenome reference graphs from rGFA content (#148)

.gitignore

@@ -65,4 +65,8 @@ data/
 output/
 
 # Classification run outputs (large generated JSON, written at repo root)
-/*_classifications.json
+/*_classifications.json
+
+# Literal $TMPDIR dir, created when a script writes to an unexpanded "$TMPDIR"
+${TMPDIR}/
+$TMPDIR/

Makefile

@@ -1,4 +1,4 @@
-.PHONY: test lint classify classify-hprc classify-and-report download classify-bam classify-vcf classify-fastq classify-fasta classify-headers classify-bed coverage-report validation-report all-reports download-hprc validate-hprc clean help
+.PHONY: test lint classify classify-hprc classify-and-report download classify-bam classify-vcf classify-fastq classify-fasta classify-gfa classify-headers classify-bed coverage-report validation-report all-reports download-hprc validate-hprc clean help
 
 help:
 	@echo "meta-disco — AnVIL file metadata classification"
@@ -13,6 +13,7 @@ help:
 	@echo "  make classify-vcf       Classify VCF files (network required)"
 	@echo "  make classify-fastq     Classify FASTQ files (network required)"
 	@echo "  make classify-fasta     Classify FASTA files (network required)"
+	@echo "  make classify-gfa       Classify GFA/rGFA graph files (network required)"
 	@echo "  make classify-bed       Classify BED files"
 	@echo "  make classify-hprc      Classify HPRC catalog files (network required)"
 	@echo "  make coverage-report    Generate coverage report from latest run"
@@ -41,7 +42,7 @@ classify-and-report: classify classify-hprc all-reports
 download:
 	python scripts/download_anvil_metadata.py
 
-classify-headers: classify-bam classify-vcf classify-fastq classify-fasta
+classify-headers: classify-bam classify-vcf classify-fastq classify-fasta classify-gfa
 
 classify-bam:
 	python scripts/classify_headers.py --type bam -i data/anvil/anvil_files_metadata.json -o output/anvil/bam_classifications.json -w 4
@@ -55,6 +56,9 @@ classify-fastq:
 classify-fasta:
 	python scripts/classify_headers.py --type fasta -i data/anvil/anvil_files_metadata.json -o output/anvil/fasta_classifications.json -w 10
 
+classify-gfa:
+	python scripts/classify_headers.py --type gfa -i data/anvil/anvil_files_metadata.json -o output/anvil/gfa_classifications.json -w 10
+
 classify-bed:
 	python scripts/classify_bed_files.py --metadata data/anvil/anvil_files_metadata.json
 

docs/classification-hierarchy.md

@@ -120,17 +120,22 @@ file_format (extension)
 
 **data_type**: `pangenome`, `pangenome.reference`
 
-`data_type` describes structure: linear sequence (FASTA — `sequence`, refining to `assembly`/`assembly.reference`) versus a sequence graph (nodes/edges/haplotype paths), which is `pangenome`. The reference-graph case — an HPRC minigraph-cactus (`mc`) graph used as the alignment coordinate system — refines to `pangenome.reference`. PGGB and single-sample assembly graphs stay `pangenome`.
+`data_type` describes structure: linear sequence (FASTA — `sequence`, refining to `assembly`/`assembly.reference`) versus a sequence graph (nodes/edges/haplotype paths), which is `pangenome`. The reference-graph case — a graph used as the alignment coordinate system — refines to `pangenome.reference`, detected two ways:
+
+- **From file content (preferred)** — an rGFA whose segments carry stable rank 0 (`SR:i:0`) plus a stable name (`SN:Z:`) defines a reference backbone. `fetch_gfa_segment_tags` reads the S-line tags from the head of the file and `classify_from_gfa_segment_tags` emits `rgfa_stable_rank_reference`. When those tags are in the head the signal is exact; if they are not, no content claim is made and detection falls back to the filename. This catches the HPRC `minigraph` graphs (#148).
+- **From filename** — HPRC minigraph-cactus graphs match the `-mc-` token (tier-2 `pangenome_reference_mc`). Their GFA segments carry no tags, and the PanSN `P`/`W` path lines that identify them sit after every segment line — ~11 GB into the file — so content inspection is not economical.
+
+PGGB and single-sample assembly graphs stay `pangenome`.
 
 **data_modality**: `genomic` ← extension (all graph formats)
 
 **assay_type**: `not_applicable`
 
 **platform**: `not_applicable` (a graph is an assembled product, not raw reads)
 
-**reference_assembly**: `GRCh38/GRCh37/CHM13` ← filename (shared `filename_ref_*` rules; e.g. an `mc-grch38` graph's coordinate system)
+**reference_assembly**: `GRCh38/GRCh37/CHM13` ← filename (shared `filename_ref_*` rules; e.g. an `mc-grch38` graph's coordinate system). Not derived from graph content: no sequence lengths are parsed, so contig-length detection cannot run, and the stable names in the fetched head (`chr1`) are shared between GRCh38 and CHM13.
 
-**Coverage**: Extension-driven, so all listed formats classify. `pangenome.reference` detection keys on the `-mc-` filename token. Auxiliary vg indices (`.min`, `.dist`, `.snarls`, `.gg`) are out of scope (issue #144).
+**Coverage**: Extension-driven, so all listed formats classify. Content inspection covers the text GFA formats (`.gfa`, `.gfa.gz`, `.rgfa`, `.rgfa.gz`); the binary vg/GBWT formats (`.gbz`, `.vg`, `.gbwt`, `.xg`) classify from extension and filename alone. Auxiliary vg indices (`.min`, `.dist`, `.snarls`, `.gg`) are out of scope (issue #144). Single- vs multi-sample graph distinction is deferred (issue #147).
 
 ---
 

scripts/rerun_all_classifications.py

@@ -12,6 +12,42 @@
 from datetime import datetime
 from pathlib import Path
 
+sys.path.insert(0, str(Path(__file__).parent.parent))
+
+from src.meta_disco.file_types import FILE_TYPE_REGISTRY  # noqa: E402
+
+# Phase 1 classifiers that are NOT header-based, so they have their own script
+# rather than a FILE_TYPE_REGISTRY entry.
+NON_HEADER_JOBS = (
+    ("classify_bed_files.py", "bed_classifications.json"),
+    ("classify_images.py", "image_classifications.json"),
+    ("classify_auxiliary_genomic.py", "auxiliary_classifications.json"),
+)
+
+
+def build_parallel_jobs(metadata: Path, output_dir: Path) -> list[tuple]:
+    """Phase 1 jobs: one per header-based file type, plus the non-header scripts.
+
+    The header jobs are derived from FILE_TYPE_REGISTRY rather than hand-listed,
+    so registering a new file type cannot silently skip production. That is what
+    happened to `gfa` in #151: it was added to the registry and to nothing else,
+    so `make classify` never invoked it and graph files fell through to the
+    filename-only Phase 3 catch-all.
+
+    Every output filename here must also appear in output_utils.CLASSIFICATION_FILES
+    or the reports will not read it — pinned by tests/test_orchestration.py.
+    """
+    jobs = [
+        ("classify_headers.py", output_dir / f"{ftype}_classifications.json",
+         ["--type", ftype, "--input", str(metadata)])
+        for ftype in FILE_TYPE_REGISTRY
+    ]
+    jobs += [
+        (script, output_dir / out, ["--metadata", str(metadata)])
+        for script, out in NON_HEADER_JOBS
+    ]
+    return jobs
+
 
 def run_script(script_name: str, output_path: Path, extra_args: list[str] = None):
     """Run a classification script."""
@@ -55,23 +91,7 @@ def main():
     print(f"Re-running classifications with timestamp: {timestamp}")
     print(f"Output directory: {output_dir}")
 
-    # Phase 1: Run independent classifiers in parallel
-    parallel_jobs = [
-        ("classify_headers.py", output_dir / "bam_classifications.json",
-         ["--type", "bam", "--input", str(args.metadata)]),
-        ("classify_headers.py", output_dir / "vcf_classifications.json",
-         ["--type", "vcf", "--input", str(args.metadata)]),
-        ("classify_headers.py", output_dir / "fastq_classifications.json",
-         ["--type", "fastq", "--input", str(args.metadata)]),
-        ("classify_bed_files.py", output_dir / "bed_classifications.json",
-         ["--metadata", str(args.metadata)]),
-        ("classify_images.py", output_dir / "image_classifications.json",
-         ["--metadata", str(args.metadata)]),
-        ("classify_auxiliary_genomic.py", output_dir / "auxiliary_classifications.json",
-         ["--metadata", str(args.metadata)]),
-        ("classify_headers.py", output_dir / "fasta_classifications.json",
-         ["--type", "fasta", "--input", str(args.metadata)]),
-    ]
+    parallel_jobs = build_parallel_jobs(args.metadata, output_dir)
 
     # Track all classification output paths for Phase 3
     all_classification_files = [path for _, path, _ in parallel_jobs]