def read_evidence_network(groups, *, max_rows=1_000_000, max_edges=1_000_000):
"""Read and reconcile a schema-1 study's complete reported-membership graph.
Each sample/peptide row must occur once. Edge counts refer to distinct
acquisitions reporting that member, not PSM counts, quantities or unique
proteins. Row/edge limits fail explicitly rather than truncate evidence.
"""
max_rows = _integer(max_rows, "max_rows", 1)
max_edges = _integer(max_edges, "max_edges", 1)
root = Path(groups).resolve(strict=True)
names = [
"summary.json",
"study_groups.tsv",
"peptide_to_study_group.tsv",
"peptide_evidence.tsv.gz",
"input_manifest.json",
]
identities = {name: _sha(root / name) for name in names}
summary = json.loads((root / "summary.json").read_text())
source_manifest = json.loads((root / "input_manifest.json").read_text())
if summary.get("schema_version") != 1:
raise ValueError("Network export requires a schema-1 study summary")
samples = summary.get("acquisition_names")
if not isinstance(samples, list) or any(not isinstance(s, str) for s in samples):
raise ValueError("Study acquisition roster must be a list of distinct names")
if len(samples) != len(set(samples)):
raise ValueError("Study acquisition roster must be a list of distinct names")
for sample in samples:
_identity(sample, "acquisition name")
if summary.get("acquisitions") != len(samples):
raise ValueError("Study acquisition count disagrees with its roster")
sample_set = set(samples)
representatives = {}
fields = ["study_group", "selection_rank", "observed_support_peptides", "assigned_peptides"]
for row in _rows(root / "study_groups.tsv", fields, max_rows):
group = _identity(row["study_group"], "study group")
if group in representatives:
raise ValueError(f"Duplicate study group: {group}")
representatives[group] = {key: _integer(row[key], key) for key in fields[1:]}
if summary.get("inferred_study_groups") != len(representatives):
raise ValueError("Study group count disagrees with the dictionary")
if sorted(r["selection_rank"] for r in representatives.values()) != list(
range(len(representatives))
):
raise ValueError("Study selection ranks must cover the dictionary exactly")
assignment = {}
for row in _rows(
root / "peptide_to_study_group.tsv", ["canonical_peptide", "study_group"], max_rows
):
peptide = _identity(row["canonical_peptide"], "canonical peptide")
if (
not peptide.isascii()
or not peptide.isalpha()
or not peptide.isupper()
or "I" in peptide
):
raise ValueError(
f"Peptide is not canonical uppercase I/L-normalized sequence: {peptide}"
)
if peptide in assignment or row["study_group"] not in representatives:
raise ValueError("Duplicate peptide or foreign group in the assignment dictionary")
assignment[peptide] = row["study_group"]
assigned_counts = Counter(assignment.values())
if summary.get("pooled_canonical_peptides") != len(assignment):
raise ValueError("Peptide count disagrees with the assignment dictionary")
if any(assigned_counts[g] != r["assigned_peptides"] for g, r in representatives.items()):
raise ValueError("Assigned peptide counts disagree with the study dictionary")
memberships = defaultdict(set)
peptide_samples = defaultdict(set)
protein_samples = defaultdict(set)
sample_peptides = defaultdict(set)
sample_proteins = defaultdict(set)
declared_degree = {}
fields = [
"sample",
"study_group",
"canonical_peptide",
"reported_proteins",
"n_reported_proteins_sample",
"n_reported_proteins_study",
]
for row in _rows(root / "peptide_evidence.tsv.gz", fields, max_rows):
sample, peptide = row["sample"], row["canonical_peptide"]
if sample not in sample_set or peptide not in assignment:
raise ValueError("Peptide evidence contains a foreign acquisition or peptide")
if row["study_group"] != assignment[peptide]:
raise ValueError("Peptide evidence disagrees with the reporting assignment")
if sample in peptide_samples[peptide]:
raise ValueError(f"Duplicate acquisition/peptide evidence: {sample}, {peptide}")
proteins = row["reported_proteins"].split(";")
for protein in proteins:
_identity(protein, "reported protein accession")
if len(proteins) != len(set(proteins)) or len(proteins) != _integer(
row["n_reported_proteins_sample"], "sample membership", 1
):
raise ValueError("Reported sample membership count is inconsistent")
degree = _integer(row["n_reported_proteins_study"], "study membership", 1)
if declared_degree.setdefault(peptide, degree) != degree:
raise ValueError("Study membership count changes between acquisitions")
peptide_samples[peptide].add(sample)
sample_peptides[sample].add(peptide)
for protein in proteins:
key = (peptide, protein)
if key not in memberships and len(memberships) >= max_edges:
raise ValueError(
f"Network exceeds {max_edges} edges; raise the explicit edge limit"
)
memberships[key].add(sample)
protein_samples[protein].add(sample)
sample_proteins[sample].add(protein)
observed_degree = Counter(peptide for peptide, _ in memberships)
support_degree = Counter(protein for _, protein in memberships)
if set(peptide_samples) != set(assignment):
raise ValueError("Peptide evidence does not cover the complete assignment dictionary")
for peptide, group in assignment.items():
if observed_degree[peptide] != declared_degree[peptide]:
raise ValueError("Study membership total disagrees with the evidence rows")
if (peptide, group) not in memberships:
raise ValueError("Reporting representative lacks the declared study peptide membership")
if any(support_degree[g] != r["observed_support_peptides"] for g, r in representatives.items()):
raise ValueError("Observed support disagrees with the study dictionary")
if summary.get("proteins_with_observed_evidence") != len(protein_samples):
raise ValueError("Observed protein count disagrees with the evidence graph")
if identities != {name: _sha(root / name) for name in names}:
raise ValueError("Study inputs changed during network export")
graph_id = hashlib.sha256(json.dumps(identities, sort_keys=True).encode()).hexdigest()
def node_id(kind, label):
"""Derive a kind-prefixed node ID from the graph identity and node label."""
return kind + ":" + hashlib.sha256((graph_id + "\0" + label).encode()).hexdigest()
nodes = []
for peptide in sorted(assignment):
nodes.append(
dict(
id=node_id("peptide", peptide),
node_type="peptide",
label=peptide,
reported_degree=observed_degree[peptide],
observed_acquisitions=len(peptide_samples[peptide]),
selected_representative=False,
assigned_group=assignment[peptide],
)
)
for protein in sorted(protein_samples):
nodes.append(
dict(
id=node_id("protein", protein),
node_type="reported_protein",
label=protein,
reported_degree=support_degree[protein],
observed_acquisitions=len(protein_samples[protein]),
selected_representative=protein in representatives,
assigned_group="",
)
)
edges = [
dict(
source=node_id("peptide", peptide),
target=node_id("protein", protein),
relation="reported_membership",
observed_acquisitions=len(observed),
assigned_representative=assignment[peptide] == protein,
acquisitions=sorted(observed),
)
for (peptide, protein), observed in sorted(memberships.items())
]
return dict(
schema_version=1,
graph_id=graph_id,
source_files=identities,
source_summary=summary,
input_manifest=source_manifest,
nodes=nodes,
edges=edges,
acquisitions=[
dict(
sample=sample,
canonical_peptides=len(sample_peptides[sample]),
reported_proteins=len(sample_proteins[sample]),
)
for sample in samples
],
scope="Sage-reported peptide/protein memberships and separate study assignments; "
"not reference-wide specificity, protein inference or physical interactions",
)