解读蛋白质物理互作网络

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
import networkx as nx
RECIPE_ID='05-string-physical-network'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'network': '19b6b87b2a66cd2e11f0c8c82ee712659dc7e8a9eaded9e36238604a5223d136'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['network']

def physical_graph(snapshot, threshold=0.7):
    if not 0 <= threshold <= 1:
        raise ValueError("Network threshold must be between 0 and 1")
    graph = nx.Graph()
    for mapping in snapshot["mapping"]:
        graph.add_node(mapping["preferredName"],string_id=mapping["stringId"])
    for row in snapshot["edges"]:
        a, b = row["preferredName_A"], row["preferredName_B"]
        if a == b or row["score"] < threshold:
            continue
        if row["ncbiTaxonId"] != 9606 or a not in graph or b not in graph:
            raise ValueError("Unexpected network species or neighbor")
        attrs = {key:float(row.get(key,0)) for key in ("score","escore","dscore","tscore")}
        if not graph.has_edge(a,b) or graph[a][b]["score"] < attrs["score"]:
            graph.add_edge(a,b,**attrs)
    return graph

print("Method ready: physical_graph")

Method ready: physical_graph

network=snapshot("network")
mapping=pd.DataFrame(network["mapping"])
display(mapping[[c for c in ["queryItem","queryIndex","stringId","preferredName","ncbiTaxonId"] if c in mapping]])
print("Query:",network["query"])
print("Snapshot API threshold:",network["request"]["required_score"],"; added nodes:",network["request"]["add_nodes"])
print("Version:",network["version"])

Query: ['TP53', 'MDM2', 'MDM4', 'ATM', 'ATR', 'CHEK1', 'CHEK2', 'BRCA1', 'BRCA2', 'CDKN1A', 'BAX', 'BCL2']
Snapshot API threshold: 400 ; added nodes: 0
Version: [{'stable_address': 'https://version-12-5.string-db.org', 'string_version': '12.5'}]

THRESHOLD=.7
graph=physical_graph(network,THRESHOLD)
print("Nodes / edges / isolates:",graph.number_of_nodes(),graph.number_of_edges(),list(nx.isolates(graph)))
positions=nx.spring_layout(graph,seed=42,weight="score")
fig,ax=plt.subplots(figsize=(8,6))
nx.draw_networkx_edges(graph,positions,ax=ax,alpha=.5,width=[1+2*graph[a][b]["score"] for a,b in graph.edges],edge_color="#829b8f")
nx.draw_networkx_nodes(graph,positions,ax=ax,node_size=[450+80*graph.degree(n) for n in graph],
    node_color=["#18785f" if n=="TP53" else "#c4ded3" for n in graph],edgecolors="white",linewidths=1.5)
nx.draw_networkx_labels(graph,positions,ax=ax,font_size=10)
ax.set_title("STRING physical network: human query, confidence ≥0.7");ax.axis("off")
plt.tight_layout();plt.show()

Nodes / edges / isolates: 12 32 []

edge_table=pd.DataFrame([dict(protein_a=a,protein_b=b,**attrs) for a,b,attrs in graph.edges(data=True)],
    columns=["protein_a","protein_b","score","escore","dscore","tscore"])
evidence=edge_table.sort_values("score",ascending=False).head(15)
fig,ax=plt.subplots(figsize=(7,5))
if len(evidence):
    im=ax.imshow(evidence[["score","escore","dscore","tscore"]],vmin=0,vmax=1,cmap="YlGnBu",aspect="auto")
    ax.set(xticks=range(4),xticklabels=["Combined","Experiments","Database","Text mining"],
           yticks=range(len(evidence)),yticklabels=evidence.protein_a+" — "+evidence.protein_b)
    fig.colorbar(im,ax=ax,label="Channel score")
else:ax.text(.5,.5,"No edges at this threshold",ha="center",transform=ax.transAxes)
ax.set_title("Evidence channels for the highest-score edges");plt.tight_layout();plt.show()
display(edge_table.head(10))

	protein_a	protein_b	score	escore	dscore	tscore
0	TP53	MDM2	0.999	0.999	0.7	0.927
1	TP53	BAX	0.997	0.995	0.0	0.560
2	TP53	ATM	0.999	0.999	0.0	0.651
3	TP53	ATR	0.999	0.999	0.0	0.113
4	TP53	MDM4	0.999	0.999	0.7	0.760
5	TP53	BRCA2	0.999	0.999	0.0	0.190
6	TP53	CHEK2	0.999	0.999	0.0	0.518
7	TP53	BCL2	0.999	0.999	0.9	0.680
8	TP53	CDKN1A	0.999	0.999	0.0	0.343
9	TP53	CHEK1	0.999	0.999	0.0	0.223

thresholds=[.4,.7,.9,.95,.99]
sensitivity=pd.DataFrame([dict(threshold=t,edges=g.number_of_edges(),isolates=len(list(nx.isolates(g))),
    components=nx.number_connected_components(g)) for t in thresholds for g in [physical_graph(network,t)]])
fig,ax=plt.subplots(figsize=(7,3))
ax.plot(sensitivity.threshold,sensitivity.edges,"o-",label="Edges")
ax.plot(sensitivity.threshold,sensitivity.isolates,"s-",label="Isolates")
ax.set(xlabel="Minimum combined confidence",ylabel="Count",title="Network sensitivity to the confidence threshold")
ax.legend();plt.tight_layout();plt.show();display(sensitivity)

	threshold	edges	isolates	components
0	0.40	38	0	1
1	0.70	32	0	1
2	0.90	27	0	1
3	0.95	23	0	1
4	0.99	20	0	1

centrality=nx.betweenness_centrality(graph,weight=None)
communities=nx.community.louvain_communities(graph,weight="score",seed=42) if graph.number_of_edges() else [{n} for n in graph]
membership={n:i+1 for i,nodes in enumerate(communities) for n in nodes}
node_table=pd.DataFrame([dict(protein=n,degree=graph.degree(n),betweenness=centrality[n],community=membership[n],
    string_id=graph.nodes[n]["string_id"]) for n in graph]).sort_values(["degree","protein"],ascending=[False,True])
display(node_table)

assert graph.number_of_nodes()==12
assert nx.number_of_selfloops(graph)==0
assert all(attrs["score"]>=THRESHOLD for _,_,attrs in graph.edges(data=True))
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
edge_table.to_csv(OUTPUT_DIR/"physical-edges.csv",index=False)
node_table.to_csv(OUTPUT_DIR/"proteins.csv",index=False)
sensitivity.to_csv(OUTPUT_DIR/"threshold-audit.csv",index=False)
nx.write_graphml(graph,OUTPUT_DIR/"physical-network.graphml")
print("Checks passed; GraphML, nodes, edges and threshold audit exported")

Checks passed; GraphML, nodes, edges and threshold audit exported

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://string-db.org/api/json/version',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged