在三维空间中分析血红蛋白

旋转真实 1A3N 坐标，观察血红素与邻近残基，将结构链编号与 RefSeq 蛋白序列联系起来。

生物信息学 · 第 06 课🔋 离线核心分析⏱ 60 分钟内嵌真实数据📦 numpy · pandas · scipy · matplotlib · py3Dmol · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python

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='06-hemoglobin-structure'
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

读取内嵌快照

In [2] · Python

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={'structure': '86cdd64dd8465a5f1448470db5067ea1ebccc2b61faad19710f78e9714d56e29', 'sequences': 'a8ae02a43c841e2932efafaf82a0bf78199baa8921a821fcbd774c8d7649e8b7'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['structure', 'sequences']

方法函数

In [3] · Python

def parse_pdb(text):
    """First model; choose highest occupancy altloc, then blank/A/lexical order."""
    rows, first_model, inside = {}, None, True
    for line in text.splitlines():
        record = line[:6].strip()
        if record == "MODEL":
            number = line[10:14].strip()
            if first_model is None:
                first_model = number
            inside = number == first_model
        if record == "ENDMDL" and inside:
            break
        if record not in ("ATOM","HETATM") or not inside:
            continue
        try:
            atom = dict(record=record,serial=int(line[6:11]),atom=line[12:16].strip(),altloc=line[16:17].strip(),
                residue=line[17:20].strip(),chain=line[21:22].strip(),resnum=int(line[22:26]),icode=line[26:27].strip(),
                x=float(line[30:38]),y=float(line[38:46]),z=float(line[46:54]),
                occupancy=float(line[54:60].strip() or 0),bfactor=float(line[60:66].strip() or 0),element=line[76:78].strip())
        except ValueError as exc:
            raise ValueError("Invalid PDB numeric fields") from exc
        key = (record,atom["chain"],atom["resnum"],atom["icode"],atom["residue"],atom["atom"])
        rank = (atom["occupancy"],atom["altloc"]=="",atom["altloc"]=="A")
        old = rows.get(key)
        if old is None or rank > old[0] or (rank == old[0] and atom["altloc"] < old[1]["altloc"]):
            rows[key] = (rank,atom)
    columns = ["record","serial","atom","altloc","residue","chain","resnum","icode","x","y","z","occupancy","bfactor","element"]
    atoms = pd.DataFrame([entry[1] for entry in rows.values()],columns=columns)
    if atoms.empty:
        raise ValueError("No atoms in the first model")
    if not np.isfinite(atoms[["x","y","z","occupancy"]].to_numpy()).all():
        raise ValueError("Nonfinite coordinates or occupancy")
    return atoms

print("Method ready: parse_pdb")

Method ready: parse_pdb

分步分析

1. 查看模型前先读结构信息

1A3N 是人类脱氧血红蛋白的 X 射线实验结构，包含 α、β 蛋白链及血红素配体。分辨率描述实验，不能直接当成每个原子坐标的统一误差范围。旋转模型前先检查实验方法、分辨率及结构编号。本课不执行结构预测。

In [4] · Python

structure=snapshot("structure")
pdb_text=structure["pdb"];entry=structure["entry"]
print("PDB:",structure["pdb_id"])
print("Method:",[item["method"] for item in entry["exptl"]])
print("Resolution (Å):",entry["rcsb_entry_info"].get("resolution_combined"))
print("Title:",entry["struct"]["title"])
print("Coordinate file:",len(pdb_text),"characters")

PDB: 1A3N
Method: ['X-RAY DIFFRACTION']
Resolution (Å): [1.8]
Title: DEOXY HUMAN HEMOGLOBIN
Coordinate file: 461700 characters

2. 解析链、插入码及替代构象

本课小型解析器按 PDB 固定列读取 ATOM/HETATM。选择第一个模型，按 occupancy 为每个原子选取替代构象，同分时使用确定性规则。残基标识同时包含链、残基编号和插入码。审计保留 occupancy 为零的原子，坐标分析则排除它们。本解析器针对 PDB 格式，不适用于任意 mmCIF 文件。

In [5] · Python

atoms=parse_pdb(pdb_text)
usable=atoms.loc[atoms.occupancy>0].copy()
protein=usable.loc[usable.record=="ATOM"]
ca=protein.loc[protein.atom=="CA"].sort_values(["chain","resnum","icode"])
chain_summary=ca.groupby("chain").agg(observed_ca=("atom","size"),first_residue=("resnum","min"),last_residue=("resnum","max"))
display(chain_summary)
display(usable.loc[usable.record=="HETATM"].groupby("residue").size().rename("atoms").to_frame())
print("Parsed atoms:",len(atoms),"; zero occupancy:",int((atoms.occupancy==0).sum()))

	observed_ca	first_residue	last_residue
chain
A	141	1	141
B	145	2	146
C	141	1	141
D	145	2	146

	atoms
residue
HEM	172
HOH	451

Parsed atoms: 4993 ; zero occupancy: 37

3. 旋转模型并观察血红素

预览使用内嵌坐标与应用随附的 3Dmol 库。拖动旋转，双指缩放。cartoon 展示蛋白骨架，sticks 突出血红素；隐藏水分子以减少遮挡。可以将 STYLE 在 cartoon 与 stick（API 使用单数）之间切换并重新运行。颜色和标签用于辅助观察，不构成实验证据。

In [6] · Python

import py3Dmol
STYLE="cartoon"
view=py3Dmol.view(width="100%",height=420)
view.addModel(pdb_text,"pdb")
view.setStyle({"hetflag":False},{STYLE:{"colorscheme":"chain"}})
view.setStyle({"resn":"HEM"},{"stick":{"colorscheme":"orangeCarbon","radius":.18}})
view.setStyle({"resn":"HOH"},{})
view.setBackgroundColor("#f8faf8");view.zoomTo()
display(HTML(view._make_html()))
print("Interactive 1A3N view: drag to rotate; pinch/scroll to zoom")

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Interactive 1A3N view: drag to rotate; pinch/scroll to zoom

4. 将观测骨架投影到坐标空间

静态图使用每个观测残基的一个 Cα 原子，按链着色。它无需 WebGL，能与交互模型互相核对。只连接残基编号连续的观测点；缺失残基不能被连成已经观测到的骨架。

In [7] · Python

fig=plt.figure(figsize=(7,5));ax=fig.add_subplot(111,projection="3d")
for chain,rows in ca.groupby("chain"):
    coords=rows[["x","y","z"]].to_numpy();ax.scatter(*coords.T,s=8,label=chain)
    numbers=rows.resnum.to_numpy()
    for i in range(len(rows)-1):
        if numbers[i+1]-numbers[i]==1:ax.plot(*coords[i:i+2].T,color="#afbdb5",lw=.5,alpha=.5)
ax.set(xlabel="x (Å)",ylabel="y (Å)",zlabel="z (Å)",title="Observed Cα atoms in the deposited coordinates")
ax.legend(title="Chain");plt.tight_layout();plt.show()

5. 连接序列编号与观测链

RefSeq HBB 翻译得到含起始甲硫氨酸的 147 aa 序列，成熟 β 链惯用编号通常省略它。先读取沉积文件的 SEQRES 预期链序列，再与 ATOM 中的残基比较。缺少坐标并不表示蛋白中不存在该残基；将变异映射到结构时必须区分这两种情况。

In [8] · Python

seqres={}
for line in pdb_text.splitlines():
    if line.startswith("SEQRES"):
        seqres.setdefault(line[11],[]).extend(line[19:70].split())
hbb=next(r for r in snapshot("sequences")["records"] if r["gene"]=="HBB")
print("RefSeq HBB amino acids:",len(hbb["translation"]),"; without starting methionine:",len(hbb["translation"][1:]))
for chain,rows in ca.groupby("chain"):
    expected=len(seqres.get(chain,[]));observed=len(rows)
    print("Chain",chain,"SEQRES:",expected,"; usable Cα:",observed,"; difference:",expected-observed)
display(ca.loc[ca.chain=="B",["chain","resnum","icode","residue"]].head(10))

RefSeq HBB amino acids: 147 ; without starting methionine: 146
Chain A SEQRES: 141 ; usable Cα: 141 ; difference: 0
Chain B SEQRES: 146 ; usable Cα: 145 ; difference: 1
Chain C SEQRES: 141 ; usable Cα: 141 ; difference: 0
Chain D SEQRES: 146 ; usable Cα: 145 ; difference: 1

	chain	resnum	residue
1070	B	2	HIS
1080	B	3	LEU
1088	B	4	THR
1095	B	5	PRO
1102	B	6	GLU
1111	B	7	GLU
1120	B	8	LYS
1129	B	9	SER
1135	B	10	ALA
1140	B	11	VAL

6. 测量距离并定义接触图

对 B 链计算以 Å 为单位的 Cα 欧氏距离。接触定义为距离 <8 Å，并排除自身及序列位置相隔不超过 2 的残基。阈值和原子选择共同定义这幅描述性接触图，不能当成通用生化互作标准。距离矩阵必须对称，对角线必须为零。

In [9] · Python

from scipy.spatial.distance import cdist
chain_b=ca.loc[ca.chain=="B"].copy()
xyz=chain_b[["x","y","z"]].to_numpy()
distances=cdist(xyz,xyz)
separation=np.abs(chain_b.resnum.to_numpy()[:,None]-chain_b.resnum.to_numpy()[None,:])
contacts=(distances<8)&(separation>2)
fig,axes=plt.subplots(1,2,figsize=(10,4))
im=axes[0].imshow(distances,cmap="viridis");fig.colorbar(im,ax=axes[0],label="Cα distance (Å)")
axes[1].imshow(contacts,cmap="Greys",vmin=0,vmax=1)
for ax,title in zip(axes,["Chain B pairwise distances","Contacts <8 Å, sequence separation >2"]):
    ax.set(xlabel="Observed residue index",ylabel="Observed residue index",title=title)
plt.tight_layout();plt.show()
assert np.allclose(distances,distances.T) and np.allclose(np.diag(distances),0)
print("Unique contact pairs:",int(np.triu(contacts,1).sum()))

Unique contact pairs: 390

7. 寻找一个血红素附近的残基

选择 B 链血红素，计算每个蛋白残基的非氢原子到该配体的最小距离。保留 4 Å 内的残基，满足条件的其他链残基也保留。邻近关系可用于选择观察目标，但不能计算亲和力或证明催化作用。B-factor 图表示模型位移参数，不能当成 AlphaFold pLDDT。

In [10] · Python

heavy=usable.loc[~usable.element.isin(["H","D"])]
heme=heavy.loc[(heavy.residue=="HEM")&(heavy.chain=="B")]
assert not heme.empty
near=[]
for key,rows in heavy.loc[heavy.record=="ATOM"].groupby(["chain","resnum","icode","residue"]):
    minimum=float(cdist(rows[["x","y","z"]],heme[["x","y","z"]]).min())
    if minimum<=4:near.append(dict(chain=key[0],resnum=key[1],icode=key[2],residue=key[3],minimum_distance_A=minimum))
neighbors=pd.DataFrame(near).sort_values("minimum_distance_A")
display(neighbors)
fig,ax=plt.subplots(figsize=(8,3))
ax.plot(chain_b.resnum,chain_b.bfactor,color="#416aa6")
ax.set(xlabel="Deposited chain B residue number",ylabel="Cα B-factor (Å²)",title="Experimental model displacement parameters")
plt.tight_layout();plt.show()

	chain	resnum	residue	minimum_distance_A
11	B	92	HIS	2.261944
4	B	63	HIS	3.296014
14	B	102	ASN	3.387431
17	B	141	LEU	3.414291
1	B	41	PHE	3.479508
12	B	96	LEU	3.519500
6	B	67	VAL	3.526266
15	B	103	PHE	3.589757
10	B	91	LEU	3.707475
0	B	31	LEU	3.744089
2	B	42	PHE	3.756942
7	B	70	ALA	3.800856
8	B	71	PHE	3.814101
13	B	98	VAL	3.836636
16	B	106	LEU	3.899402
3	B	45	PHE	3.917097
5	B	66	LYS	3.975199
9	B	88	LEU	3.985119

8. 导出坐标、表格与测量结果

保留原始坐标文件及派生摘要，方便在其他预览器中复现。距离注明单位及链选择。导出的是观测原子及残基标识，不会补齐缺失残基或预测其坐标。

In [11] · Python

assert set(ca.chain)=={"A","B","C","D"}
assert len(heme)>0 and len(chain_b)>100
assert np.isfinite(distances).all()
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
(OUTPUT_DIR/"1A3N.pdb").write_text(pdb_text)
chain_summary.to_csv(OUTPUT_DIR/"chains.csv")
neighbors.to_csv(OUTPUT_DIR/"heme-neighbors.csv",index=False)
pd.DataFrame(distances,index=chain_b.resnum,columns=chain_b.resnum).to_csv(OUTPUT_DIR/"chain-B-distances-A.csv")
atoms.to_csv(OUTPUT_DIR/"selected-first-model-atoms.csv",index=False)
print("Checks passed; original PDB and derived measurements exported")

Checks passed; original PDB and derived measurements exported

自主练习与参考答案

为什么 SEQRES 中的残基数可能多于有效 ATOM Cα 记录？
在序列间隔 >2 的条件不变时，比较 6、8、10 Å 的接触数。

参考解释。 SEQRES 描述沉积的预期聚合物序列，部分原子可能未解析或 occupancy 为零。对固定坐标增大距离阈值不会减少接触数，邻近关系仍不能证明亲和力。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

在 py3Dmol 中高亮导出的 B 链血红素邻近残基，使用链、残基编号及插入码构造选择器，仅使用已有内嵌 PDB。

提示词 2

为接触图添加残基标签表，使每个观测索引对应 chain/resnum/icode，缺失残基仍不进入距离矩阵。

人工验收。 检查第一个模型、插入码标识、确定性替代构象选择、距离对称及零对角线，并确认原始 PDB 不变；实验 B-factor 不当作 pLDDT。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [12] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://data.rcsb.org/rest/v1/core/entry/1A3N',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

来源、快照与下一步

RefSeq HBB/HBD — NCBI public sequence data; retrieved 2026-10-09. SHA-256 ea64e38e8e5ffad06f6d8770eaee066ef556e1488c4b9a6b340cbd445abe04ad. GenBank XML CDS uses 1-based inclusive coordinates
PDB 1A3N coordinates — CC0 1.0; retrieved 2026-10-09. SHA-256 a4db006269e9b6f9db0eb5447a5d40e59d67857484a9294b9ca71d642ffe9861. Original response preserved in the embedded snapshot.
RCSB 1A3N entry metadata — CC0 1.0; retrieved 2026-10-09. SHA-256 e8fc0a94cd71b53bfa5736021ad877eaf68588632fb5be2c4f30040e340d11da. Original response preserved in the embedded snapshot.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。