探索真实的配对 RNA-seq 实验

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
RECIPE_ID='03-airway-expression'
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'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['airway']

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

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)

library_sizes=counts.sum(axis=0)
fig,ax=plt.subplots(figsize=(8,3))
ax.bar(range(len(counts.columns)),library_sizes/1e6,color=["#416aa6" if c=="control" else "#18785f" for c in metadata.condition])
ax.set(xticks=range(8),xticklabels=counts.columns, ylabel="Assigned gene counts (millions)",title="Sequencing depth differs between samples")
ax.tick_params(axis="x",rotation=45);plt.tight_layout();plt.show()
assert library_sizes.astype(int).tolist()==airway["library_totals"]
print("Original library totals verified")

Original library totals verified

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

correlation=log_counts.corr()
fig,ax=plt.subplots(figsize=(6,5))
im=ax.imshow(correlation,vmin=.9,vmax=1,cmap="YlGnBu")
labels=[f"{r.donor} / {r.condition}" for r in metadata.itertuples()]
ax.set(xticks=range(8),yticks=range(8),xticklabels=labels,yticklabels=labels,title="Sample correlation after depth adjustment")
plt.setp(ax.get_xticklabels(),rotation=60,ha="right")
fig.colorbar(im,ax=ax,label="Pearson r");plt.tight_layout();plt.show()

variable_genes=log_counts.var(axis=1).nlargest(min(2000,len(log_counts))).index
matrix=log_counts.loc[variable_genes].T.to_numpy()
matrix=matrix-matrix.mean(axis=0)
u,s,vt=np.linalg.svd(matrix,full_matrices=False)
coordinates=u[:,:2]*s[:2];explained=s**2/np.sum(s**2)
fig,ax=plt.subplots(figsize=(7,5))
for donor,rows in metadata.groupby("donor"):
    indices=[metadata.index.get_loc(sample) for sample in rows.index]
    ax.plot(coordinates[indices,0],coordinates[indices,1],color="#bac6bf",lw=1)
for i,row in enumerate(metadata.itertuples()):
    ax.scatter(*coordinates[i],s=85,color="#416aa6" if row.condition=="control" else "#18785f")
    ax.annotate(row.donor+(" C" if row.condition=="control" else " T"),coordinates[i],xytext=(5,4),textcoords="offset points",fontsize=9)
ax.set(xlabel=f"PC1 ({explained[0]:.1%})",ylabel=f"PC2 ({explained[1]:.1%})",title="Donor pairs: blue control → green treated")
plt.tight_layout();plt.show()
print("PCA uses",len(variable_genes),"genes; explained variance:",explained[:2])

PCA uses 2000 genes; explained variance: [0.3400221 0.276398 ]

paired_changes,response=paired_response(log_counts,metadata,min_effect=1,min_consistent=3)
display(response.direction.value_counts().rename("genes").to_frame())
top=response.mean_log2_change.abs().nlargest(16).index
heat=log_counts.loc[top];heat=heat.sub(heat.mean(axis=1),axis=0)
fig,ax=plt.subplots(figsize=(8,6))
im=ax.imshow(heat,aspect="auto",cmap="RdBu_r",vmin=-3,vmax=3)
ax.set(xticks=range(8),xticklabels=labels,yticks=range(len(top)),yticklabels=top,
       title="Largest observed paired changes: row-centered log counts")
plt.setp(ax.get_xticklabels(),rotation=60,ha="right");fig.colorbar(im,ax=ax,label="Centered log2 normalized counts")
plt.tight_layout();plt.show()

	genes
direction
other	15083
down	529
up	527

GENE="ENSG00000103196"
assert GENE in normalized.index
fig,ax=plt.subplots(figsize=(5,4))
for donor,rows in metadata.groupby("donor"):
    ordered=rows.set_index("condition").loc[["control","treated"],:]
    values=[normalized.loc[GENE,rows.index[rows.condition==condition][0]] for condition in ["control","treated"]]
    ax.plot([0,1],values,marker="o",label=donor)
ax.set(xticks=[0,1],xticklabels=["Control","Treated"],ylabel="Normalized counts",title="CRISPLD2: paired donor measurements")
ax.legend(title="Donor",fontsize=8);plt.tight_layout();plt.show()

assert np.isclose(np.exp(np.log(size_factors).mean()),1)
assert np.isclose(explained.sum(),1)
assert np.isfinite(log_counts.to_numpy()).all()
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
metadata.to_csv(OUTPUT_DIR/"samples.csv")
size_factors.to_csv(OUTPUT_DIR/"size-factors.csv")
response.to_csv(OUTPUT_DIR/"response-candidates.csv")
paired_changes.to_csv(OUTPUT_DIR/"paired-log2-changes.csv")
pd.DataFrame(coordinates,index=metadata.index,columns=["PC1","PC2"]).to_csv(OUTPUT_DIR/"pca.csv")
print("Checks passed. No differential-expression p-values were estimated.")

Checks passed. No differential-expression p-values were estimated.

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://raw.githubusercontent.com/bioconductor-source/airway/596678815f4ade04a9711997150949b9782aeddc/DESCRIPTION',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