def annotate_evidence_network(
bundle,
annotation,
fasta,
output,
*,
query_fasta=None,
fields=("KEGG_ko", "EC", "GOs", "PFAMs"),
max_records=100_000,
):
"""Attach annotations to exact sequences from a recorded study FASTA.
An accession match must also agree in sequence. Alternative names can match
by full protein SHA-256; I/L variants remain distinct. Conflicting annotations
for duplicate query sequences fail instead of selecting an arbitrary row.
Only annotated reporting representatives enter the group-to-term exports.
"""
out = Path(output)
if out.exists() or out.is_symlink():
raise FileExistsError(f"Annotation output exists: {out}; choose a new directory")
fields = tuple(fields)
if not fields or len(set(fields)) != len(fields) or not set(fields) <= set(FUNCTION_FIELDS):
raise ValueError("Choose distinct supported annotation namespaces")
network = load_evidence_network(bundle)
fasta = Path(fasta).resolve(strict=True)
fasta_hash = _sha(fasta)
expected = {network["source_summary"].get("annotation_targets", {}).get("sha256")}
manifest = network.get("input_manifest", [])
if isinstance(manifest, list):
expected.update(
row.get("sha256")
for row in manifest
if isinstance(row, dict)
and str(row.get("path", "")).endswith(
(".fasta", ".fa", ".faa", ".fasta.gz", ".fa.gz", ".faa.gz")
)
)
if fasta_hash not in expected:
raise ValueError("Supply the study's representatives.fasta or a recorded search FASTA")
proteins = read_fasta(fasta, max_records=max_records)
loaded = load_eggnog_annotations(annotation, query_fasta=query_fasta, max_records=max_records)
by_hash = defaultdict(list)
for query, (_, _, digest) in loaded["queries"].items():
by_hash[digest].append(query)
nodes, edges, rows, used_queries = {}, [], [], set()
maps = {field: [] for field in fields}
for node in network["nodes"]:
if node["node_type"] != "reported_protein":
continue
accession = node["label"]
row = dict(
protein_id=accession,
node_id=node["id"],
protein_sha256="",
query="",
match_rule="",
annotation_status="sequence_not_supplied",
selected_representative=node["selected_representative"],
raw_annotation={},
)
if accession in proteins:
digest = proteins[accession][2]
row["protein_sha256"] = digest
query = loaded["aliases"].get(accession)
if query is not None:
if loaded["queries"][query][2] != digest:
raise ValueError(
f"Annotation reuses accession for a different sequence: {accession}"
)
row["match_rule"] = "accession_and_exact_sequence"
else:
candidates = sorted(by_hash.get(digest, []))
signatures = {
json.dumps(
{k: v for k, v in loaded["rows"][q].items() if k != "query"}, sort_keys=True
)
for q in candidates
}
if len(signatures) > 1:
raise ValueError("Duplicate query sequences carry conflicting annotations")
query = candidates[0] if candidates else None
row["match_rule"] = "exact_sequence_alias" if query else ""
row["annotation_status"] = "not_annotated"
if query is not None:
used_queries.add(query)
values = loaded["rows"][query]
if not set(fields) <= set(values):
raise ValueError("eggNOG output lacks a requested annotation namespace")
row.update(
query=query,
annotation_status=values["annotation_status"],
raw_annotation=dict(values),
)
for field in fields:
for term in function_terms(values[field], field):
key = (
"function:" + hashlib.sha256((field + "\0" + term).encode()).hexdigest()
)
nodes[key] = dict(
id=key, node_type="annotated_function", namespace=field, label=term
)
edges.append(
dict(
source=node["id"],
target=key,
relation="orthology_annotation",
namespace=field,
query=query,
protein_sha256=digest,
)
)
if node["selected_representative"]:
maps[field].append({"study_group": accession, "term": term})
rows.append(row)
if not used_queries:
raise ValueError(
"No annotation query matches a graph protein with a supplied exact sequence"
)
if _sha(fasta) != fasta_hash or any(_sha(p) != h for p, h in loaded["source_files"].items()):
raise ValueError("Functional annotation inputs changed during export")
status_counts = dict(Counter(row["annotation_status"] for row in rows))
scope = (
"Orthology-based protein annotations; peptide memberships and reporting assignments "
"remain unchanged. No peptide-to-function inference, activity or flux is asserted."
)
data = dict(
schema_version=1,
graph_id=network["graph_id"],
proteins=rows,
nodes=list(nodes.values()),
edges=edges,
status_counts=status_counts,
source_fasta={"path": str(fasta), "sha256": fasta_hash},
annotation_source_files=loaded["source_files"],
annotation_provenance=loaded["provenance"],
query_linkage=loaded["linkage"],
unused_annotation_queries=sorted(set(loaded["queries"]) - used_queries),
scope=scope,
)
out.mkdir(parents=True, exist_ok=False)
(out / "annotations.json").write_text(json.dumps(data, indent=2) + "\n")
columns = [
"protein_id",
"node_id",
"protein_sha256",
"query",
"match_rule",
"annotation_status",
"selected_representative",
*fields,
]
write_table(
out / "protein_annotations.tsv",
columns,
[
{
**{k: row[k] for k in columns if k not in fields},
**{k: row["raw_annotation"].get(k, "") for k in fields},
}
for row in rows
],
)
write_table(
out / "function_nodes.tsv", ["id", "node_type", "namespace", "label"], data["nodes"]
)
write_table(
out / "function_edges.tsv",
["source", "target", "relation", "namespace", "query", "protein_sha256"],
edges,
)
for field, pairs in maps.items():
write_table(out / f"group_terms_{field}.tsv", ["study_group", "term"], pairs)
files = {p.name: _sha(p) for p in out.iterdir() if p.is_file()}
receipt = dict(
status="PASS",
schema_version=1,
graph_id=network["graph_id"],
protein_nodes=len(rows),
function_nodes=len(nodes),
annotation_edges=len(edges),
status_counts=status_counts,
files=files,
scope=scope,
)
(out / "COMPLETE.json").write_text(json.dumps(receipt, indent=2) + "\n")
return receipt