"""Inventory high-throughput (HTP-family) GO annotations: all of GOA, and this repo's reviews. Covers the four high-throughput evidence codes introduced with the GO consortium's HTP framework (PMID:30715275) that this project studies: HDA ECO:0007005 high throughput direct assay HTP ECO:0006056 high throughput experiment (unspecified) HMP ECO:0007001 high throughput mutant phenotype HGI ECO:0007003 high throughput genetic interaction HEP (high throughput expression pattern) is left to projects/IEP.md. Steps (downloads cached under data/, so reruns are offline; --refresh re-fetches): 1. all HTP-family rows from QuickGO -> data/htp_goa_all.tsv.gz 2. PubMed titles for every PMID reference -> data/htp_reference_titles.yaml 3. GO term labels -> data/htp_term_labels.yaml 4. reviewed rows from genes/*/*/*-ai-review.yaml -> data/htp_reviewed_rows.yaml 5. summary tables on stdout (saved as data/htp_inventory_report.txt by the caller) Predicates: * GOA rows: QuickGO annotation rows, exact ECO codes above (no descendants: these four codes have none in use). * Reviewed rows: existing_annotations whose evidence_type is HDA/HTP/HMP/HGI. Aspect is taken from the GOA download (term id -> aspect); terms absent from it fall back to the gene's own *-goa.tsv, else "?". * "Use" of a row: aspect x polarity -- location (C), function (F), process (P), NOT. * Cross-code papers: references cited by both an HTP-family row and an RCA row in GOA (RCA rows from ../rca/data/rca_goa_all.tsv, built by ../rca/rca_source_catalog.py). Usage: python3 projects/OMICS_EVIDENCE/htp/htp_inventory.py [--refresh] """ from __future__ import annotations import argparse import csv import glob import gzip import json import os import sys import time import urllib.parse import urllib.request from collections import Counter, defaultdict import yaml HERE = os.path.dirname(os.path.abspath(__file__)) ROOT = os.path.abspath(os.path.join(HERE, "..", "..", "..")) DATA = os.path.join(HERE, "data") GOA = os.path.join(DATA, "htp_goa_all.tsv.gz") TITLES = os.path.join(DATA, "htp_reference_titles.yaml") LABELS = os.path.join(DATA, "htp_term_labels.yaml") REVIEWED = os.path.join(DATA, "htp_reviewed_rows.yaml") RCA_GOA = os.path.join(HERE, "..", "rca", "data", "rca_goa_all.tsv") CODES = {"HDA": "ECO:0007005", "HTP": "ECO:0006056", "HMP": "ECO:0007001", "HGI": "ECO:0007003"} DL = "https://www.ebi.ac.uk/QuickGO/services/annotation/downloadSearch" TERMS_API = "https://www.ebi.ac.uk/QuickGO/services/ontology/go/terms/" ESUMMARY = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi" LOADER = getattr(yaml, "CSafeLoader", yaml.SafeLoader) def fetch(url: str, accept: str | None = None, data: bytes | None = None) -> bytes: for attempt in range(4): try: req = urllib.request.Request(url, data=data, headers={"Accept": accept} if accept else {}) with urllib.request.urlopen(req, timeout=600) as resp: return resp.read() except Exception: # noqa: BLE001 - network retry if attempt == 3: raise time.sleep(2 ** (attempt + 1)) raise RuntimeError("unreachable") def load_goa(refresh: bool) -> list[dict]: if refresh or not os.path.exists(GOA): q = urllib.parse.urlencode( { "evidenceCode": ",".join(CODES.values()), "evidenceCodeUsage": "exact", "downloadLimit": 200000, "selectedFields": "geneProductId,symbol,qualifier,goId,goAspect,evidenceCode,goEvidence," "reference,taxonId,assignedBy,date", } ) raw = fetch(f"{DL}?{q}", accept="text/tsv") os.makedirs(DATA, exist_ok=True) with gzip.open(GOA, "wb") as fh: fh.write(raw) with gzip.open(GOA, "rt", newline="") as fh: return list(csv.DictReader(fh, delimiter="\t")) def load_titles(refs: list[str], refresh: bool) -> dict[str, dict]: cache = {} if refresh or not os.path.exists(TITLES) else (yaml.safe_load(open(TITLES)) or {}) todo = [r.split(":", 1)[1] for r in refs if r.startswith("PMID:") and r not in cache] for i in range(0, len(todo), 150): batch = todo[i : i + 150] body = urllib.parse.urlencode({"db": "pubmed", "id": ",".join(batch), "retmode": "json"}) res = json.loads(fetch(ESUMMARY, data=body.encode()))["result"] for pmid in batch: rec = res.get(pmid) or {} cache[f"PMID:{pmid}"] = {"title": rec.get("title", ""), "year": (rec.get("pubdate") or "")[:4]} time.sleep(0.4) with open(TITLES, "w") as fh: fh.write("# Generated by htp_inventory.py (NCBI esummary) -- do not edit\n") yaml.safe_dump(dict(sorted(cache.items())), fh, sort_keys=False, width=120, allow_unicode=True) return cache def load_labels(ids: list[str], refresh: bool) -> dict[str, dict]: cache = {} if refresh or not os.path.exists(LABELS) else (yaml.safe_load(open(LABELS)) or {}) todo = [i for i in ids if i not in cache and i.startswith("GO:")] for i in range(0, len(todo), 100): res = json.loads(fetch(TERMS_API + ",".join(todo[i : i + 100]), accept="application/json")) for t in res.get("results", []): cache[t["id"]] = {"name": t.get("name", ""), "aspect": {"molecular_function": "F", "biological_process": "P", "cellular_component": "C"}.get(t.get("aspect"), "?")} with open(LABELS, "w") as fh: fh.write("# Generated by htp_inventory.py (QuickGO ontology API) -- do not edit\n") yaml.safe_dump(dict(sorted(cache.items())), fh, sort_keys=False, width=120) return cache def reviewed_rows() -> list[dict]: out = [] for path in sorted(glob.glob(os.path.join(ROOT, "genes", "*", "*", "*-ai-review.yaml"))): text = open(path).read() if not any(f"evidence_type: {c}" in text for c in CODES): continue try: doc = yaml.load(text, Loader=LOADER) or {} except yaml.YAMLError: continue org, gene = path.split(os.sep)[-3], path.split(os.sep)[-2] for a in doc.get("existing_annotations") or []: if a.get("evidence_type") in CODES: rv = a.get("review") or {} out.append({ "organism": org, "gene": gene, "code": a["evidence_type"], "term_id": (a.get("term") or {}).get("id"), "term_label": (a.get("term") or {}).get("label"), "negated": bool(a.get("negated")), "reference": a.get("original_reference_id"), "action": rv.get("action", "MISSING"), }) return out def pct(n: int, d: int) -> str: return f"{100 * n / d:5.1f}%" if d else " n/a" def main() -> int: ap = argparse.ArgumentParser() ap.add_argument("--refresh", action="store_true") args = ap.parse_args() rows = load_goa(args.refresh) total = len(rows) refs = sorted({r["REFERENCE"] for r in rows}) titles = load_titles(refs, args.refresh) rev = reviewed_rows() labels = load_labels(sorted({r["GO TERM"] for r in rows} | {r["term_id"] for r in rev if r["term_id"]}), args.refresh) with open(REVIEWED, "w") as fh: fh.write("# Generated by htp_inventory.py -- do not edit\n") yaml.safe_dump(rev, fh, sort_keys=False, width=120) def lab(t: str) -> str: return (labels.get(t) or {}).get("name", "") def title(ref: str) -> str: return (titles.get(ref) or {}).get("title", "") print(f"# HTP-family rows in GOA: {total} ({len(refs)} references)") # Rows overstate claims: one gene x term is repeated once per annotation extension # (e.g. a DUB x each substrate, a TF x each target promoter), so report distinct # gene-product x term pairs alongside rows everywhere. pairs_all = {(r["GENE PRODUCT ID"], r["GO TERM"]) for r in rows} print(f"# distinct gene-product x term pairs: {len(pairs_all)}") print("\n## code x aspect x polarity (rows, distinct gene x term pairs)") c = Counter() cp: dict[tuple, set] = defaultdict(set) for r in rows: k = (r["GO EVIDENCE CODE"], r["GO ASPECT"], "NOT" if r["QUALIFIER"].startswith("NOT") else "pos") c[k] += 1 cp[k].add((r["GENE PRODUCT ID"], r["GO TERM"])) for (code, asp, pol), n in c.most_common(): m = len(cp[(code, asp, pol)]) print(f" {n:6d} {pct(n, total)} pairs={m:6d} {pct(m, len(pairs_all))} {code} {asp} {pol}") print("\n## assigned_by (rows, references, years)") g_rows, g_refs, g_years = Counter(), defaultdict(set), defaultdict(set) for r in rows: g = r["ASSIGNED BY"] g_rows[g] += 1 g_refs[g].add(r["REFERENCE"]) g_years[g].add(r["DATE"][:4]) for g, n in g_rows.most_common(): print(f" {n:6d} {pct(n, total)} {g:18s} refs={len(g_refs[g]):3d} {min(g_years[g])}-{max(g_years[g])}") by_ref: dict[str, list[dict]] = defaultdict(list) for r in rows: by_ref[r["REFERENCE"]].append(r) conc = sorted((len(v) for v in by_ref.values()), reverse=True) print(f"\n## concentration: top 10 refs = {pct(sum(conc[:10]), total)}, top 50 = {pct(sum(conc[:50]), total)}") print("\n## top 40 references") for ref, rs in sorted(by_ref.items(), key=lambda kv: -len(kv[1]))[:40]: codes = Counter(r["GO EVIDENCE CODE"] + "/" + r["GO ASPECT"] for r in rs) terms = Counter(r["GO TERM"] for r in rs).most_common(2) print(f" {len(rs):5d} genes={len({r['GENE PRODUCT ID'] for r in rs}):5d} {ref:16s} " f"{','.join(sorted({r['ASSIGNED BY'] for r in rs}))} " f"{dict(codes)}\n {title(ref)[:110]}\n " + "; ".join(f"{t} {lab(t)} ({n})" for t, n in terms)) for asp, name in [("F", "molecular function"), ("P", "biological process")]: print(f"\n## {name} rows: top terms (code, rows, references)") tc = Counter((r["GO EVIDENCE CODE"], r["GO TERM"]) for r in rows if r["GO ASPECT"] == asp) for (code, t), n in tc.most_common(15): sub = [r for r in rows if r["GO TERM"] == t and r["GO EVIDENCE CODE"] == code] print(f" {n:5d} genes={len({r['GENE PRODUCT ID'] for r in sub}):5d} {code} {t} {lab(t)} " f"refs={len({r['REFERENCE'] for r in sub})}") print(f" -- references behind {name} rows (>=40 rows)") rc = Counter(r["REFERENCE"] for r in rows if r["GO ASPECT"] == asp) for ref, n in rc.most_common(): if n < 40: break print(f" {n:5d} {ref:16s} {title(ref)[:100]}") print("\n## NOT-qualified rows") for r in rows: if r["QUALIFIER"].startswith("NOT"): print(f" {r['SYMBOL']} {r['QUALIFIER']} {r['GO TERM']} {lab(r['GO TERM'])} {r['REFERENCE']} {r['ASSIGNED BY']}") if os.path.exists(RCA_GOA): with open(RCA_GOA, newline="") as fh: rca = list(csv.DictReader(fh, delimiter="\t")) rca_refs = Counter(r["REFERENCE"] for r in rca) both = sorted((set(by_ref) & set(rca_refs)), key=lambda x: -(len(by_ref[x]) + rca_refs[x])) print(f"\n## papers cited by both HTP-family and RCA rows: {len(both)}") for ref in both: h = Counter(f"{r['GO EVIDENCE CODE']}/{r['GO ASPECT']}" for r in by_ref[ref]) rr = Counter(f"RCA/{r['GO ASPECT']}{'-NOT' if r['QUALIFIER'].startswith('NOT') else ''}" for r in rca if r["REFERENCE"] == ref) print(f" {ref:16s} {dict(h)} {dict(rr)} {title(ref)[:80]}") # ---- reviewed corpus ---- n = len(rev) print(f"\n# Reviewed HTP-family rows: {n} in {len({(r['organism'], r['gene']) for r in rev})} gene folders") for r in rev: r["aspect"] = (labels.get(r["term_id"]) or {}).get("aspect", "?") print("\n## code x aspect x action") ca: dict[tuple, Counter] = defaultdict(Counter) for r in rev: ca[(r["code"], r["aspect"] + ("-NOT" if r["negated"] else ""))][r["action"]] += 1 for k, cc in sorted(ca.items(), key=lambda kv: -sum(kv[1].values())): m = sum(cc.values()) print(f" {m:5d} {k[0]} {k[1]:6s} ACCEPT={pct(cc['ACCEPT'], m)} " + ", ".join(f"{a}={v}" for a, v in cc.most_common())) # Location rows split by what kind of location they assert. Vesicle-type terms come # from secreted-vesicle / body-fluid proteomes, whose cargo is a sample of the cell. vesicle = {"GO:0070062", "GO:1903561", "GO:0072562", "GO:0031982", "GO:0065010", "GO:0043230"} generic = {"GO:0016020"} print("\n## reviewed location (C) rows, HDA+HTP: vesicle-type / generic membrane / other") groups = {"vesicle-type (exosome, EV, microparticle, vesicle)": vesicle, "generic membrane": generic} for name, ids in list(groups.items()) + [("all other locations", None)]: sel = [r for r in rev if r["code"] in ("HDA", "HTP") and r["aspect"] == "C" and not r["negated"] and ((r["term_id"] in ids) if ids else (r["term_id"] not in vesicle | generic))] cc = Counter(r["action"] for r in sel) m = len(sel) rej = cc["MARK_AS_OVER_ANNOTATED"] + cc["REMOVE"] + cc["MODIFY"] print(f" {m:5d} {name:52s} ACCEPT={pct(cc['ACCEPT'], m)} KEEP_NON_CORE={pct(cc['KEEP_AS_NON_CORE'], m)} " f"over/remove/modify={pct(rej, m)}") print("\n## organism folders") for k, v in Counter(r["organism"] for r in rev).most_common(12): print(f" {v:5d} {k}") print("\n## reviewed references (>=12 rows): action split") rr2: dict[str, Counter] = defaultdict(Counter) for r in rev: rr2[r["reference"]][r["action"]] += 1 for ref, cc in sorted(rr2.items(), key=lambda kv: -sum(kv[1].values())): m = sum(cc.values()) if m < 12: break print(f" {m:4d} {ref:16s} ACCEPT={pct(cc['ACCEPT'], m)} " + ", ".join(f"{a}={v}" for a, v in cc.most_common() if a != "ACCEPT") + f"\n {title(ref)[:100]}") # Same papers, two codes: how did reviews treat HDA (location) vs RCA (function) rows? rca_rev_path = os.path.join(HERE, "..", "rca", "data", "rca_reviewed_rows.yaml") if os.path.exists(RCA_GOA) and os.path.exists(rca_rev_path): rca_rev = yaml.safe_load(open(rca_rev_path)) or [] shared = set(by_ref) & {r["REFERENCE"] for r in rca} h = Counter(r["action"] for r in rev if r["reference"] in shared) k = Counter(r["action"] for r in rca_rev if r.get("reference") in shared) print(f"\n## papers cited by both codes ({len(shared)}): reviewed actions") for name, cc in [("HTP-family rows", h), ("RCA rows", k)]: m = sum(cc.values()) print(f" {name:16s} n={m:4d} ACCEPT={pct(cc['ACCEPT'], m)} " + ", ".join(f"{a}={v}" for a, v in cc.most_common())) print("\n## reviewed terms with the most non-ACCEPT rows (top 20)") tr: dict[str, Counter] = defaultdict(Counter) for r in rev: tr[f"{r['code']} {r['term_id']} {r['term_label']}"][r["action"]] += 1 for t, cc in sorted(tr.items(), key=lambda kv: -(sum(kv[1].values()) - kv[1]["ACCEPT"]))[:20]: m = sum(cc.values()) print(f" {m - cc['ACCEPT']:4d}/{m:<4d} {t} " + ", ".join(f"{a}={v}" for a, v in cc.most_common())) return 0 if __name__ == "__main__": sys.exit(main())