检验响应候选基因的通路富集

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
from scipy.stats import hypergeom
import networkx as nx
RECIPE_ID='04-reactome-enrichment'
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={'airway': '5440d5e8dc3a3da15735ad2592cf908bcfc7f4bdadc89b1404aa48d728bdeaab', 'reactome': '7ac8938695d294fb8c8b56072fca97d7420fafb2f436c6437a8ef72b7c3bfbe9'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['airway', 'reactome']

def validate_counts(counts, metadata):
    if not counts.index.is_unique or not counts.columns.is_unique:
        raise ValueError("Gene and sample identifiers must be unique")
    if not metadata.index.is_unique or list(counts.columns) != list(metadata.index):
        raise ValueError("Sample columns must match metadata order exactly")
    values = counts.to_numpy()
    if not np.isfinite(values).all() or (values < 0).any() or not np.equal(values,np.floor(values)).all():
        raise ValueError("Counts must be finite nonnegative integers")
    if (values.sum(axis=0) == 0).any():
        raise ValueError("A sample has no reads")
    if set(metadata.condition) != {"control","treated"}:
        raise ValueError("Expected control and treated conditions")
    for donor, rows in metadata.groupby("donor", sort=True):
        if len(rows) != 2 or set(rows.condition) != {"control","treated"}:
            raise ValueError("Each donor needs one control and one treated sample: " + str(donor))

print("Method ready: validate_counts")

Method ready: validate_counts

def normalize_counts(counts, min_count=10, min_samples=4):
    values = counts.to_numpy(dtype=float)
    keep = (values >= min_count).sum(axis=1) >= min_samples
    filtered = counts.loc[keep].astype(float)
    positive = (filtered.to_numpy() > 0).all(axis=1)
    if positive.sum() == 0:
        raise ValueError("Median-ratio method needs genes positive in every sample")
    reference_values = filtered.to_numpy()[positive]
    geometric_means = np.exp(np.log(reference_values).mean(axis=1))
    factors = np.median(reference_values / geometric_means[:,None], axis=0)
    if not np.isfinite(factors).all() or (factors <= 0).any():
        raise ValueError("Invalid size factors")
    factors /= np.exp(np.log(factors).mean())
    normalized = filtered / factors
    return filtered, pd.Series(factors,index=counts.columns,name="size_factor"), normalized

print("Method ready: normalize_counts")

Method ready: normalize_counts

def paired_response(log_counts, metadata, min_effect=1.0, min_consistent=3):
    differences = {}
    for donor, rows in metadata.groupby("donor",sort=True):
        control = rows.index[rows.condition == "control"][0]
        treated = rows.index[rows.condition == "treated"][0]
        differences[donor] = log_counts[treated] - log_counts[control]
    differences = pd.DataFrame(differences)
    response = pd.DataFrame(dict(mean_log2_change=differences.mean(axis=1),
                                 positive_pairs=(differences>0).sum(axis=1),
                                 negative_pairs=(differences<0).sum(axis=1)))
    response["direction"] = "other"
    response.loc[(response.mean_log2_change>=min_effect)&(response.positive_pairs>=min_consistent),"direction"]="up"
    response.loc[(response.mean_log2_change<=-min_effect)&(response.negative_pairs>=min_consistent),"direction"]="down"
    return differences, response

print("Method ready: paired_response")

Method ready: paired_response

def bh_adjust(pvalues):
    p = np.asarray(pvalues,dtype=float)
    if not len(p):
        return p.copy()
    if not np.isfinite(p).all() or (p<0).any() or (p>1).any():
        raise ValueError("p-values must be finite and in [0,1]")
    order = np.argsort(p,kind="stable")
    adjusted = np.minimum.accumulate((p[order]*len(p)/np.arange(1,len(p)+1))[::-1])[::-1]
    result = np.empty_like(p); result[order] = np.minimum(adjusted,1)
    return result

print("Method ready: bh_adjust")

Method ready: bh_adjust

def overrepresentation(up, down, expressed, pathways, min_size=10, max_size=500):
    annotated = set().union(*(set(row["genes"]) for row in pathways)) if pathways else set()
    background = set(expressed) & annotated
    columns = ["direction","pathway","name","background_size","query_size","term_size","overlap","pvalue","genes","qvalue"]
    rows = []
    for direction, supplied in (("up",up),("down",down)):
        query = set(supplied) & background
        for term in pathways:
            members = set(term["genes"]) & background
            if not min_size <= len(members) <= max_size:
                continue
            hits = sorted(query & members)
            p = float(hypergeom.sf(len(hits)-1,len(background),len(members),len(query))) if query else 1.0
            rows.append(dict(direction=direction,pathway=term["id"],name=term["name"],
                background_size=len(background),query_size=len(query),term_size=len(members),
                overlap=len(hits),pvalue=p,genes=";".join(hits)))
    result = pd.DataFrame(rows)
    if len(result):
        result["qvalue"] = bh_adjust(result.pvalue.to_numpy())
        result = result.sort_values(["qvalue","pvalue","pathway","direction"],kind="stable")
    else:
        result = pd.DataFrame(columns=columns)
    return background,result

print("Method ready: overrepresentation")

Method ready: overrepresentation

airway = snapshot("airway")
counts = pd.DataFrame(airway["counts"], index=airway["genes"], columns=airway["samples"])
metadata = pd.DataFrame(airway["metadata"]).set_index("sample").loc[counts.columns]
validate_counts(counts, metadata)
print("Complete count matrix:", counts.shape)
def display_scroll_table(frame):
    # Keep identifiers and numeric values intact on narrow device screens.
    table = frame.to_html(border=0, classes="vibeit-scroll-table")
    style = """<style>
    table.vibeit-scroll-table {table-layout:auto !important;width:max-content !important;max-width:none !important;}
    table.vibeit-scroll-table th, table.vibeit-scroll-table td {
      white-space:nowrap !important;word-break:normal !important;padding:6px 8px;}
    </style>"""
    display(HTML(style+'<div style="max-width:100%;overflow-x:auto;-webkit-overflow-scrolling:touch">'+table+'</div>'))
display_scroll_table(metadata)
display_scroll_table(counts.iloc[:6])

Complete count matrix: (63677, 8)

MIN_COUNT, MIN_SAMPLES = 10, 4
filtered, size_factors, normalized = normalize_counts(counts, MIN_COUNT, MIN_SAMPLES)
log_counts = np.log2(normalized + 1)
print("Retained genes:", len(filtered), "of", len(counts))
print("Genes positive in every sample:", int((filtered > 0).all(axis=1).sum()))
display(size_factors.to_frame())

Retained genes: 16139 of 63677
Genes positive in every sample: 16072

paired_changes,response=paired_response(log_counts,metadata,min_effect=1,min_consistent=3)
up=set(response.index[response.direction=="up"])
down=set(response.index[response.direction=="down"])
fig,ax=plt.subplots(figsize=(7,4))
for direction,color in [("other","#bdc5c0"),("up","#bc6541"),("down","#416aa6")]:
    selected=response.direction==direction
    ax.scatter(normalized.mean(axis=1)[selected],response.mean_log2_change[selected],s=8,alpha=.45,color=color,label=direction)
ax.set_xscale("log");ax.axhline(1,color="#888",ls="--",lw=.7);ax.axhline(-1,color="#888",ls="--",lw=.7)
ax.set(xlabel="Mean normalized count",ylabel="Mean paired log2 change",title="Effect screening before pathway testing")
ax.legend();plt.tight_layout();plt.show()
print("Up / down / retained:",len(up),len(down),len(filtered))

Up / down / retained: 527 529 16139

reactome=snapshot("reactome")
pathways=reactome["pathways"]
annotated=set().union(*(set(term["genes"]) for term in pathways))
background=set(filtered.index)&annotated
print("Reactome version:",reactome["version"])
audit=pd.DataFrame([dict(stage="Measured",genes=len(counts)),dict(stage="Expression filter",genes=len(filtered)),
 dict(stage="Annotated background",genes=len(background)),dict(stage="Up candidates mapped",genes=len(up&background)),
 dict(stage="Down candidates mapped",genes=len(down&background))])
display(audit)
print("Unmapped up / down candidates:",len(up-background),len(down-background))

Reactome version: 97

Unmapped up / down candidates: 239 257

background,enrichment=overrepresentation(up,down,filtered.index,pathways,min_size=10,max_size=500)
print("Complete test family:",len(enrichment))
print("Terms with q < 0.05:",int((enrichment.qvalue<.05).sum()) if len(enrichment) else 0)
display(enrichment.drop(columns="genes").head(12))

Complete test family: 2086
Terms with q < 0.05: 5

plot_terms=enrichment.loc[enrichment.overlap>0].head(12).copy() if len(enrichment) else enrichment.copy()
fig,ax=plt.subplots(figsize=(9,5))
if len(plot_terms):
    colors=-np.log10(plot_terms.qvalue.clip(lower=1e-300))
    dots=ax.scatter(plot_terms.overlap/plot_terms.query_size,np.arange(len(plot_terms)),
        s=30+18*plot_terms.overlap,c=colors,cmap="YlGnBu")
    names=[r.direction+": "+r.name[:62] for r in plot_terms.itertuples()]
    ax.set(yticks=range(len(names)),yticklabels=names);ax.invert_yaxis()
    fig.colorbar(dots,ax=ax,label="−log10(BH q)")
else:
    ax.text(.5,.5,"No eligible terms with overlapping candidates",ha="center",va="center",transform=ax.transAxes)
ax.set(xlabel="Overlap / mapped candidate count",title="Reactome overrepresentation: effect-screened candidates")
plt.tight_layout();plt.show()

term_graph=nx.Graph()
sets={}
for r in plot_terms.itertuples():
    key=r.direction+":"+r.pathway;sets[key]=set(r.genes.split(";"))
    term_graph.add_node(key,label=r.name[:28],direction=r.direction)
keys=list(sets)
for i,x in enumerate(keys):
    for y in keys[i+1:]:
        similarity=len(sets[x]&sets[y])/len(sets[x]|sets[y]) if sets[x]|sets[y] else 0
        if similarity>=.2:term_graph.add_edge(x,y,jaccard=similarity)
fig,ax=plt.subplots(figsize=(8,5))
if term_graph:
    pos=nx.spring_layout(term_graph,seed=42)
    nx.draw_networkx(term_graph,pos,ax=ax,labels=nx.get_node_attributes(term_graph,"label"),
        node_color=["#bde3d1" if term_graph.nodes[k]["direction"]=="up" else "#c5d6ee" for k in term_graph],
        node_size=650,font_size=8,edge_color="#a4b4ac")
else:ax.text(.5,.5,"No terms to compare",ha="center",transform=ax.transAxes)
ax.set_title("Shared candidate genes between displayed pathways");ax.axis("off")
plt.tight_layout();plt.show()

if len(enrichment):
    assert enrichment.qvalue.between(0,1).all()
    assert (enrichment.qvalue+1e-15>=enrichment.pvalue).all()
    assert (enrichment.overlap<=enrichment.term_size).all()
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
enrichment.to_csv(OUTPUT_DIR/"reactome-all-tests.csv",index=False)
response.to_csv(OUTPUT_DIR/"response-screen.csv")
pd.Series(sorted(background),name="ensembl_gene").to_csv(OUTPUT_DIR/"background.csv",index=False)
audit.to_csv(OUTPUT_DIR/"universe-audit.csv",index=False)
print("Checks passed; full background and test family exported")

Checks passed; full background and test family exported

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://reactome.org/ContentService/data/database/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