Allergen and antigen modules for the ten proteome runs; third aggregate report

AGENTS.md

@@ -90,7 +90,7 @@ All optimizations maintain full backward compatibility with the existing API. No
 
 - Added: modules `gpi`, `disorder`, `aga1_model` and the thresholds table; six calls (`gpi_anchored_protein`, `long_disordered_protein`, `gpi_long_disorder_protein`, `ser_thr_gpi_wall_candidate`, `flocculin_like_candidate`, `model_surface_candidate`); acceptance on 10 proteomes with smoke status files for two calls in two taxa; `cellsurface_sorting_hat_aggregate` and a first report (`docs/reports/data/aga1_aggregate_2026-10-26/`).
 - Later (PR after 161): `analysis/sorting_hat_run/aga1_domain_modules.py` converts the finished Pfam and TMHMM scans and runs repeat02, and a second aggregate report (`docs/reports/data/aga1_aggregate_domains_2026-10-26/`) leaves only the allergen and antigen calls NA.
-- Not done: the hand-check verdicts (owner), per-clade threshold recalibration, orthogroup/diversity statistics, allergen BLAST and antigen lookup on the ten proteomes.
+- Not done: the hand-check verdicts (owner), per-clade threshold recalibration, orthogroup/diversity statistics, allergen BLAST and antigen lookup on the ten proteomes (done afterwards: `analysis/sorting_hat_run/aga1_allergen_antigen_modules.py` and `docs/reports/data/aga1_aggregate_allergen_antigen_2026-10-26/`).
 
 ---
 

analysis/sorting_hat_run/aga1_allergen_antigen_modules.py

@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+"""Add the allergen_homology and antigen_lookup modules to existing sorting-hat work directories and re-run the engine.
+
+Per species folder `<runs-root>/<species>/` (with `p.faa` and `wd/`): BLASTP of the proteome against the IUIS fungal allergen sequences
+(options of scripts/sorting_hat/blast_allergen.sbatch: -evalue 1 -max_target_seqs 200 -seg no), module `allergen`, module `antigen`
+(applicable taxon 246410 only; every other proteome gets `not_applicable`), then the engine into `<species>/out_allergen_antigen`.
+
+The antigen ranking is keyed on C. immitis RS gene IDs (CIMG_nnnnn). The UniProt proteome of C. immitis RS is keyed on UniProt accessions, so a
+protein map (accession -> gene) is made from the `GN=CIMG_...` field of the UniProt FASTA header (`--cimm-uniprot-fasta`); accessions without
+such a field are `not_in_reference`. Run it with `python3 -I`; blastp and makeblastdb must be on PATH. Nothing is overwritten.
+"""
+
+import argparse
+import hashlib
+import re
+import subprocess
+import sys
+import tempfile
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[2]
+sys.path.insert(0, str(ROOT / "src"))
+
+APPLICABLE = "246410"
+RANKING = ROOT / "analysis/cocci_antigens/cocci_antigen_ranking.tsv"
+ISOALLERGEN = ROOT / "analysis/allergen_scoping/iuis_fungal_isoallergens.tsv"
+ALLERGEN = ROOT / "analysis/allergen_scoping/iuis_fungal_allergens.tsv"
+CAP, EVALUE, SEG = 200, "1", "no"
+BLAST_COLUMNS = "6 qseqid sseqid pident length qlen slen bitscore evalue"
+GN = re.compile(r"\bGN=(CIMG_\d+)\b")
+
+
+class StepError(Exception):
+    pass
+
+
+def sha256(path):
+    return hashlib.sha256(Path(path).read_bytes()).hexdigest()
+
+
+def uniprot_protein_map(fasta, out):
+    """accession -> CIMG gene from the UniProt header (`>tr|ACC|NAME ... GN=CIMG_nnnnn ...`); returns (n_mapped, n_proteins)."""
+    n = m = 0
+    with open(fasta) as fh, open(out, "w") as o:
+        o.write("protein_id\tgene_id\n")
+        for line in fh:
+            if not line.startswith(">"):
+                continue
+            n += 1
+            parts = line[1:].split()[0].split("|")
+            acc = parts[1] if len(parts) == 3 else parts[0]
+            g = GN.search(line)
+            if g:
+                o.write(f"{acc}\t{g.group(1)}\n")
+                m += 1
+    return m, n
+
+
+def blast_version():
+    return subprocess.run(
+        ["blastp", "-version"], check=True, capture_output=True, text=True
+    ).stdout.splitlines()[0]
+
+
+def run_species(species, taxon, runs_root, nodes, allergen_fa, pmap, tmp, threads):
+    from cellsurface_sorting_hat.cli import main as engine
+    from cellsurface_sorting_hat.modules.cli import main as module
+
+    base = Path(runs_root) / species
+    if (base / "out_allergen_antigen").exists():
+        raise StepError(f"{base / 'out_allergen_antigen'} exists; nothing is overwritten")
+    fa, wd = base / "p.faa", base / "wd"
+    for name in ("allergen_homology", "antigen_lookup"):
+        if (wd / "modules" / f"{name}.json").exists():
+            raise StepError(f"{wd / 'modules' / name}.json exists; nothing is overwritten")
+    blast_out = tmp / f"{species}.blast.tsv"
+    subprocess.run(
+        ["blastp", "-query", str(fa), "-db", str(tmp / "allergens"), "-evalue", EVALUE, "-max_target_seqs", str(CAP), "-seg", SEG,
+         "-outfmt", BLAST_COLUMNS, "-num_threads", str(threads), "-out", str(blast_out)],
+        check=True,
+    )  # fmt: skip
+    common = ["--fasta", str(fa), "--workdir", str(wd)]
+    steps = [
+        ["allergen", *common, "--blast", str(blast_out), "--allergen-fasta", str(allergen_fa), "--blast-version", blast_version(),
+         "--evalue", EVALUE, "--seg", SEG, "--max-target-seqs", str(CAP)],
+        ["antigen", *common, "--taxon", taxon, "--applicable-taxa", APPLICABLE, "--ranking", str(RANKING), "--protein-map", str(pmap)],
+    ]  # fmt: skip
+    for argv in steps:
+        if module(argv) != 0:
+            raise StepError(f"{species}: module {argv[0]} failed")
+    rc = engine(
+        [
+            "--fasta",
+            str(fa),
+            "--taxon",
+            taxon,
+            "--taxdump",
+            str(nodes),
+            "--workdir",
+            str(wd),
+            "--out",
+            str(base / "out_allergen_antigen"),
+        ]
+    )
+    if rc != 0:
+        raise StepError(f"{species}: the engine failed")
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
+    ap.add_argument("--runs-root", required=True)
+    ap.add_argument("--taxdump", required=True, help="nodes.dmp")
+    ap.add_argument("--taxa", required=True, help="TSV with columns species, taxon_id")
+    ap.add_argument(
+        "--cimm-uniprot-fasta",
+        required=True,
+        help="UniProt FASTA of C. immitis RS (UP000001261) for the accession-to-gene map",
+    )
+    ap.add_argument("--threads", type=int, default=4)
+    a = ap.parse_args(argv)
+    from cellsurface_sorting_hat.modules.allergen import build_allergen_fasta
+
+    taxa = {}
+    for line in Path(a.taxa).read_text().splitlines()[1:]:
+        f = line.split("\t")
+        taxa[f[0]] = f[1]  # species, taxon_id
+    with tempfile.TemporaryDirectory() as t:
+        tmp = Path(t)
+        allergen_fa = tmp / "iuis_fungal_allergens.faa"
+        n, skipped = build_allergen_fasta(ISOALLERGEN, allergen_fa, ALLERGEN)
+        if n >= CAP:
+            raise StepError(
+                f"{n} allergen sequences; the -max_target_seqs cap {CAP} must exceed the database size"
+            )
+        print(
+            f"allergen database: {n} sequences ({len(skipped)} skipped); sha256 {sha256(allergen_fa)}"
+        )
+        subprocess.run(
+            [
+                "makeblastdb",
+                "-in",
+                str(allergen_fa),
+                "-dbtype",
+                "prot",
+                "-out",
+                str(tmp / "allergens"),
+                "-logfile",
+                "/dev/null",
+            ],
+            check=True,
+        )
+        pmap = tmp / "protein_map.tsv"
+        m, total = uniprot_protein_map(a.cimm_uniprot_fasta, pmap)
+        print(
+            f"protein map: {m} of {total} UniProt accessions have GN=CIMG_; sha256 {sha256(pmap)}; ranking sha256 {sha256(RANKING)}"
+        )
+        for species, taxon in taxa.items():
+            run_species(species, taxon, a.runs_root, a.taxdump, allergen_fa, pmap, tmp, a.threads)
+            print(f"{species}: done")
+    return 0
+
+
+if __name__ == "__main__":
+    try:
+        sys.exit(main())
+    except StepError as err:
+        print(f"FATAL: {err}", file=sys.stderr)
+        sys.exit(1)

docs/TECHNICAL-NOTES.md

@@ -154,3 +154,6 @@ context, phylogenetic statistics, a scale-out runner, an interactive page.
 Second report with the Pfam, TMHMM and repeat02 modules added (`analysis/sorting_hat_run/aga1_domain_modules.py`):
 `docs/reports/data/aga1_aggregate_domains_2026-10-26/`. Remaining NA calls: those that need the allergen BLAST and antigen lookup.
 
+Third report with the allergen and antigen modules (`analysis/sorting_hat_run/aga1_allergen_antigen_modules.py`, local BLASTP, no HPCC):
+`docs/reports/data/aga1_aggregate_allergen_antigen_2026-10-26/`. `iuis_allergen_homolog` stays not assessable (no active allergen Pfam family).
+

docs/reports/data/aga1_aggregate_allergen_antigen_2026-10-26/README.md

@@ -0,0 +1,25 @@
+# Aggregate report with allergen and antigen modules: 10 proteomes (2026-10-26)
+
+Third aggregate report. Adds `allergen_homology` and `antigen_lookup` to the work directories of
+`docs/reports/data/aga1_aggregate_domains_2026-10-26/`. Inputs stay outside the repository; no sequences are stored.
+
+Made with `python3 -I analysis/sorting_hat_run/aga1_allergen_antigen_modules.py --runs-root <runs> --taxdump nodes.dmp --taxa analysis/sorting_hat_run/aga1_acceptance_taxa.tsv --cimm-uniprot-fasta <UP000001261.faa>`
+(local, no HPCC job), then `cellsurface_sorting_hat_aggregate --manifest manifest.tsv --taxdump nodes.dmp --names names.dmp --rank class --out <new folder>`.
+
+* **Allergen:** BLASTP (local blastp 2.17.0+; `-evalue 1 -max_target_seqs 200 -seg no`, the options of `scripts/sorting_hat/blast_allergen.sbatch`) of each proteome against the IUIS fungal allergen sequences (111 sequences, 5 skipped; database sha256 `86a559cf...a16a7`).
+* **Antigen:** module `antigen`, applicable taxon 246410 only. Every other proteome is `not_applicable`. The ranking is `analysis/cocci_antigens/cocci_antigen_ranking.tsv` (sha256 `cf04b9d4...318bc`), keyed on *C. immitis* RS gene IDs. The UniProt proteome of *C. immitis* RS has accessions, so the accession-to-gene map was made from the `GN=CIMG_` field of the UniProt headers: 9,442 of 9,712 accessions map (sha256 of the map `452177f4...e977`). The `protein_map.tsv` that made the earlier Cimm record is on HPCC only and was not used; this is a different identifier route, not a re-run of that record.
+* Engine output: `<species>/out_allergen_antigen`; `categories.yaml` sha256 unchanged.
+
+| sha256 | file |
+|---|---|
+| 29baf95dfcf50db3c75d3ae7db91c4e31330a2dc0e112943958444c52f3e67e2 | counts_by_group.tsv |
+| 3b7fa046fe1a79a5c480adb4b16351ad9513dd05f321c91b517bf30f0fe80b06 | counts_by_proteome.tsv |
+| e31a36e20cfe4158056ef231809891e6f273755c8c9d75b43ded732c366c17dd | proteomes.tsv |
+| 16fec78d7b5f5b56e39709c9b0da91e69cdc921b35441256fd20adcc1f717f6d | report.md |
+
+## Reading the numbers
+
+* `cocci_specificity_rank_top15`: *C. immitis* RS 1,334 of 8,679 assessable proteins (153.7 per 1,000, as expected for a top 15% rule). Every other proteome is `not_applicable`, which the report shows as NA (nothing is assessable there), not as 0. `serodiagnostic_marker_candidate`: 135 in *C. immitis*, none elsewhere.
+* `iuis_allergen_similarity` (identity at least 35% over at least 80 residues to an IUIS fungal allergen): 8.8 to 13.9 per 1,000 in every proteome; a similarity call, not an allergen call, and unvalidated.
+* `iuis_allergen_homolog` is still 100% not assessable: `pfam_allergen` is `unavailable` because both allergen families are inactive in the family table (owner decision of 2026-10-17, issue #19).
+* About 11% of *C. immitis* proteins are not assessable for the antigen rank (accessions with no `GN=CIMG_` or no ranking row).