"""Characterise PAINT's IRD (Inferred from Rapid Divergence) annotations. Questions: 1. What do IRD calls negate (aspect, term, family, clade depth, date)? 2. Does an IRD reach any leaf gene, as a NOT|IBA or otherwise? 3. Are the IRD blocks effective: do clade members still get a positive IBA for the blocked term? 4. Do experimental annotations on clade members contradict (or confirm) the divergence call? 5. Where do our gene reviews meet a gene inside an IRD clade, and what did they do with the blocked term? Inputs (from the fetch scripts in this folder, plus the cached PAINT leaf GAF): ird_nodes.tsv, leaf_iba_loss_sources.tsv, ird_clade_members.tsv.gz, clade_experimental_annotations.tsv.gz, genes/**/*-ai-review.yaml Outputs: ird_experimental_conflicts.tsv experimental positives on clade members to the blocked term or a GO descendant ird_review_overlaps.tsv reviewed genes inside IRD clades, with what the review did with the blocked term RESULTS.md summary tables (regenerated every run) Usage: uv run python projects/IRD_EVIDENCE/analyze_ird.py """ from __future__ import annotations import csv import gzip import re import sys from collections import Counter, defaultdict from pathlib import Path import yaml from oaklib import get_adapter from ai_gene_review.etl.panther_paint import open_text HERE = Path(__file__).parent ROOT = HERE.parent.parent LEAF_GAF = ROOT / ".cache" / "paint" / "gene_association.paint_uniprot.gaf.gz" NODES = HERE / "ird_nodes.tsv" LEAF_TALLY = HERE / "leaf_iba_loss_sources.tsv" MEMBERS = HERE / "ird_clade_members.tsv.gz" EXP = HERE / "clade_experimental_annotations.tsv.gz" OUT_CONFLICTS = HERE / "ird_experimental_conflicts.tsv" OUT_REVIEWS = HERE / "ird_review_overlaps.tsv" OUT_MD = HERE / "RESULTS.md" PREDS = ["rdfs:subClassOf", "BFO:0000050"] ASPECT = {"F": "MF", "P": "BP", "C": "CC"} HTP = {"HTP", "HDA", "HMP", "HGI", "HEP"} def read_tsv(path: Path) -> list[dict[str, str]]: opener = gzip.open if path.suffix == ".gz" else open with opener(path, "rt") as fh: return list(csv.DictReader(fh, delimiter="\t")) def md_table(header: list[str], rows: list[list]) -> str: out = ["| " + " | ".join(header) + " |", "|" + "---|" * len(header)] out += ["| " + " | ".join(str(c) for c in r) + " |" for r in rows] return "\n".join(out) def pct(a: int, b: int) -> str: return f"{100 * a / b:.1f}%" if b else "n/a" def review_index() -> dict[str, list[tuple[Path, dict]]]: """UniProt accession -> [(review path, parsed doc)].""" idx: dict[str, list[tuple[Path, dict]]] = defaultdict(list) id_re = re.compile(r"^id:\s*(\S+)", re.M) for path in sorted(ROOT.glob("genes/*/*/*-ai-review.yaml")): text = path.read_text() m = id_re.search(text) if not m: continue acc = m.group(1).split(":")[-1] idx[acc].append((path, text)) return idx def main() -> int: go = get_adapter("sqlite:obo:go") ird = read_tsv(NODES) members = read_tsv(MEMBERS) exp = read_tsv(EXP) tally = read_tsv(LEAF_TALLY) md: list[str] = ["# IRD evidence: results", "", "Generated by `analyze_ird.py`; do not edit by hand.", ""] # ---- 1. what IRD negates n = len(ird) by_aspect = Counter(ASPECT[r["aspect"]] for r in ird) md += ["## 1. IRD node annotations in PAINT", ""] md.append(f"{n} IRD rows on {len({r['ird_node'] for r in ird})} nodes in {len({r['family'] for r in ird})} PANTHER families. " f"All are negated: {sum(r['qualifier'].startswith('NOT') for r in ird)} " f"({dict(Counter(r['qualifier'] for r in ird))}).") md += ["", md_table(["Aspect", "IRD rows"], [[k, v] for k, v in by_aspect.most_common()]), ""] stale = [r for r in ird if r["ancestral_has_ibd_for_term"] != "true"] both = [r for r in ird if r["ikr_on_same_ancestor_term"] == "true"] md.append(f"- {len(both)} IRD rows share their ancestor and term with an IKR loss elsewhere in the tree.") md.append(f"- {len(stale)} IRD rows point at an ancestral node that no longer has an IBD for the term (stale overrides):") md += ["", md_table(["Family", "IRD node", "Term", "Ancestral node"], [[r["family"], r["ird_node"], f"{r['go_id']} {r['go_label']}", r["ancestral_node"]] for r in stale]), ""] for asp in ("F", "P", "C"): c = Counter((r["go_id"], r["go_label"]) for r in ird if r["aspect"] == asp) md += [f"Top {ASPECT[asp]} terms blocked by IRD:", "", md_table(["Term", "IRD rows"], [[f"{k[0]} {k[1]}", v] for k, v in c.most_common(12)]), ""] years = Counter(r["date"][:4] for r in ird) md += ["IRD rows by annotation year:", "", md_table(["Year", "Rows"], sorted(years.items())), ""] fams = Counter(r["family"] for r in ird) md += ["Families with the most IRD rows:", "", md_table(["Family", "IRD rows"], fams.most_common(10)), ""] # ---- 2. leaf reach md += ["## 2. Do IRD calls reach leaf genes?", "", "Leaf IBA rows in `gene_association.paint_uniprot.gaf.gz`, tallied by which losses sit on the " "ancestor+term the row was projected from:", "", md_table(["Leaf polarity", "Losses on source ancestor+term", "Leaf rows"], [[t["leaf_polarity"], t["losses_on_source_ancestor_term"], t["leaf_rows"]] for t in tally]), ""] ird_only_not = sum(int(t["leaf_rows"]) for t in tally if t["leaf_polarity"] == "NOT" and t["losses_on_source_ancestor_term"] == "IRD") md.append(f"NOT leaf rows whose only source loss is an IRD: **{ird_only_not}**. " "IRD therefore never produces a NOT|IBA; it only stops the ancestral IBD descending.") md.append("") # ---- 3. clades and block effectiveness clade: dict[tuple[str, str], list[dict]] = defaultdict(list) for m in members: clade[(m["ird_node"], m["go_id"])].append(m) sizes = [len(clade.get((r["ird_node"], r["go_id"]), [])) for r in ird] unresolved = sum(1 for s in sizes if s == 0) accs = {m["uniprot"] for m in members if m["uniprot"]} human = {m["uniprot"] for m in members if m["organism"] == "Homo sapiens" and m["uniprot"]} size_bins = Counter("1 (leaf)" if s == 1 else "2-10" if s <= 10 else "11-100" if s <= 100 else ">100" for s in sizes if s) md += ["## 3. IRD clades", "", f"{len(members)} (IRD row, leaf) pairs; {len(accs)} distinct UniProt accessions " f"({len(human)} human). {unresolved} IRD rows could not be resolved to leaves.", "", md_table(["Leaves under the IRD node", "IRD rows"], [[k, size_bins[k]] for k in ("1 (leaf)", "2-10", "11-100", ">100")]), ""] # Leaf positive IBAs on clade members. PAINT can re-annotate a sub-clade below # the IRD node, so a member that still receives the blocked term (or a GO # descendant) by IBA is outside the effective block. desc_cache: dict[str, set[str]] = {} def desc(t: str) -> set[str]: if t not in desc_cache: desc_cache[t] = set(go.descendants(t, predicates=PREDS, reflexive=True)) return desc_cache[t] clade_accs = {m["uniprot"] for m in members if m["uniprot"]} leaf_pos: dict[str, set[str]] = defaultdict(set) with open_text(LEAF_GAF) as fh: for line in fh: if line.startswith("!"): continue c = line.split("\t") if len(c) < 9 or c[1] not in clade_accs or "NOT" in c[3]: continue leaf_pos[c[1]].add(c[4]) reannotated: set[tuple[str, str, str]] = set() for (node, gid), mems in clade.items(): d = desc(gid) for m in mems: if leaf_pos.get(m["uniprot"], set()) & d: reannotated.add((node, gid, m["uniprot"])) eff = {k: [m for m in v if (k[0], k[1], m["uniprot"]) not in reannotated] for k, v in clade.items()} fully = sum(1 for k, v in clade.items() if v and not eff[k]) md.append(f"PAINT re-annotates part of the clade below the IRD node: {len(reannotated)} (IRD row, protein) pairs " f"still receive a positive IBA for the blocked term or a descendant. " f"{sum(1 for k in clade if len(eff[k]) < len(clade[k]))} IRD rows are partly re-annotated, " f"{fully} fully (the IRD then blocks nothing on these leaves).") md.append("") # ---- 4. experimental conflicts exp_by_acc: dict[str, list[dict]] = defaultdict(list) for e in exp: exp_by_acc[e["GENE PRODUCT ID"]].append(e) label_cache: dict[str, str] = {} def label(t: str) -> str: if t not in label_cache: label_cache[t] = go.label(t) or "" return label_cache[t] ird_by_key = {(r["ird_node"], r["go_id"]): r for r in ird} conflicts = [] confirms = [] for (node, gid), mems in clade.items(): d = desc(gid) r = ird_by_key[(node, gid)] for m in mems: for e in exp_by_acc.get(m["uniprot"], []): if e["GO TERM"] not in d: continue row = { "family": r["family"], "ird_node": node, "aspect": ASPECT[r["aspect"]], "blocked_go_id": gid, "blocked_label": r["go_label"], "ird_date": r["date"], "clade_size": len(mems), "uniprot": m["uniprot"], "organism": m["organism"], "symbol": m["symbol"], "qualifier": e["QUALIFIER"], "exp_go_id": e["GO TERM"], "exp_label": label(e["GO TERM"]), "evidence": e["GO EVIDENCE CODE"], "reference": e["REFERENCE"], "assigned_by": e["ASSIGNED BY"], "date": e["DATE"], "paint_reannotated_below": str((node, gid, m["uniprot"]) in reannotated).lower(), } (confirms if e["QUALIFIER"].startswith("NOT") else conflicts).append(row) conflicts.sort(key=lambda x: (x["family"], x["ird_node"], x["blocked_go_id"], x["uniprot"], x["exp_go_id"])) all_conflicts = conflicts if all_conflicts: with OUT_CONFLICTS.open("w", newline="") as fh: w = csv.DictWriter(fh, fieldnames=list(all_conflicts[0]), delimiter="\t", lineterminator="\n") w.writeheader() w.writerows(all_conflicts) n_reann = sum(1 for c in conflicts if c["paint_reannotated_below"] == "true") conflicts = [c for c in conflicts if c["paint_reannotated_below"] == "false"] conf_ird = {(c["ird_node"], c["blocked_go_id"]) for c in conflicts} conf_ird_lt = {(c["ird_node"], c["blocked_go_id"]) for c in conflicts if c["evidence"] not in HTP} conf_pairs = {(c["ird_node"], c["blocked_go_id"], c["uniprot"]) for c in conflicts} md += ["## 4. Experimental annotations inside IRD clades", "", f"Experimental GOA annotations downloaded for clade members: {len(exp)} rows on {len(exp_by_acc)} proteins.", "", f"- Positive experimental annotations to the blocked term or a GO descendant: {len(all_conflicts)} rows; " f"{n_reann} of them are on proteins PAINT re-annotates below the IRD node, and are not conflicts.", f"- Inside the effective block: **{len(conflicts)} rows**, " f"{len(conf_pairs)} (IRD, protein) pairs, touching **{len(conf_ird)} of {n} IRD rows** " f"({pct(len(conf_ird), n)}); {len(conf_ird_lt)} IRD rows have at least one low-throughput conflict.", f"- NOT experimental annotations agreeing with the block: {len(confirms)} rows on " f"{len({(c['ird_node'], c['blocked_go_id']) for c in confirms})} IRD rows.", ""] by_asp = Counter(c["aspect"] for c in conflicts) by_asp_ird = Counter(ASPECT[r["aspect"]] for r in ird if (r["ird_node"], r["go_id"]) in conf_ird) md += [md_table(["Aspect", "Conflicting rows", "IRD rows conflicted", "IRD rows in aspect", "Share"], [[a, by_asp[a], by_asp_ird[a], by_aspect[a], pct(by_asp_ird[a], by_aspect[a])] for a in ("MF", "BP", "CC")]), ""] md += ["Conflicting rows by evidence code:", "", md_table(["Evidence", "Rows"], Counter(c["evidence"] for c in conflicts).most_common()), ""] md += ["Conflicting rows by relation:", "", md_table(["Relation", "Rows"], Counter(c["qualifier"] for c in conflicts).most_common()), ""] # rank IRD rows by distinct conflicting low-throughput proteins per_ird: dict[tuple[str, str], set[str]] = defaultdict(set) for c in conflicts: if c["evidence"] not in HTP: per_ird[(c["ird_node"], c["blocked_go_id"])].add(f"{c['symbol'] or c['uniprot']} ({c['organism']})") ranked = sorted(per_ird.items(), key=lambda kv: (-len(kv[1]), kv[0])) md += ["IRD rows with the most proteins carrying a low-throughput experimental positive (worklist):", "", md_table(["Family", "IRD node", "Aspect", "Blocked term", "Clade size", "Proteins with experimental positive"], [[ird_by_key[k]["family"], k[0], ASPECT[ird_by_key[k]["aspect"]], f"{k[1]} {ird_by_key[k]['go_label']}", len(clade[k]), f"{len(v)}: " + ", ".join(sorted(v)[:6]) + (" ..." if len(v) > 6 else "")] for k, v in ranked[:40]]), ""] # ---- 5. our reviews ridx = review_index() overlaps = [] for (node, gid), mems in clade.items(): r = ird_by_key[(node, gid)] d = desc(gid) for m in mems: for path, text in ridx.get(m["uniprot"], []): hits = sorted({t for t in re.findall(r"GO:\d{7}", text) if t in d}) actions = [] if hits: doc = yaml.safe_load(text) or {} for a in doc.get("existing_annotations") or []: tid = (a.get("term") or {}).get("id") if tid in d: neg = "NOT " if a.get("negated") else "" actions.append(f"{neg}{tid}/{a.get('evidence_type')}:{(a.get('review') or {}).get('action')}") for cf in doc.get("core_functions") or []: cf_terms = [(cf.get("molecular_function") or {}).get("id")] cf_terms += [x.get("id") for x in cf.get("directly_involved_in") or []] cf_terms += [x.get("id") for x in cf.get("locations") or []] for tid in cf_terms: if tid in d: actions.append(f"{tid}:CORE_FUNCTION") overlaps.append({ "review": str(path.relative_to(ROOT)), "uniprot": m["uniprot"], "family": r["family"], "ird_node": node, "aspect": ASPECT[r["aspect"]], "blocked_go_id": gid, "blocked_label": r["go_label"], "clade_size": len(mems), "review_mentions_term_or_descendant": "|".join(hits), "review_use": "; ".join(actions), "experimental_positive_in_goa": str((node, gid, m["uniprot"]) in conf_pairs).lower(), "paint_reannotated_below": str((node, gid, m["uniprot"]) in reannotated).lower(), }) overlaps.sort(key=lambda x: (x["review"], x["blocked_go_id"])) if overlaps: with OUT_REVIEWS.open("w", newline="") as fh: w = csv.DictWriter(fh, fieldnames=list(overlaps[0]), delimiter="\t", lineterminator="\n") w.writeheader() w.writerows(overlaps) used = [o for o in overlaps if o["review_use"]] used_eff = [o for o in used if o["paint_reannotated_below"] == "false"] core = [o for o in used if "CORE_FUNCTION" in o["review_use"]] md += ["## 5. Our gene reviews inside IRD clades", "", f"- {len({o['review'] for o in overlaps})} reviews cover a protein inside an IRD clade " f"({len(overlaps)} (review, blocked term) pairs).", f"- {len(used)} pairs: the review annotates the blocked term or a descendant " f"({len(core)} as a core function); {len(used_eff)} of these are inside the effective block " f"(not re-annotated below the IRD by PAINT).", ""] if used: md += [md_table(["Review", "Blocked term", "Aspect", "Clade", "Review use", "Exp. positive in GOA", "PAINT re-annotates below"], [[o["review"].split("/")[-1].replace("-ai-review.yaml", "") + f" ({o['review'].split('/')[1]})", f"{o['blocked_go_id']} {o['blocked_label']}", o["aspect"], o["clade_size"], o["review_use"], o["experimental_positive_in_goa"], o["paint_reannotated_below"]] for o in used]), ""] OUT_MD.write_text("\n".join(md) + "\n") print(f"wrote {OUT_MD}", file=sys.stderr) return 0 if __name__ == "__main__": sys.exit(main())