第三阶段 3.3：蛋白质与结构方向

对应大纲第三阶段 3.3 蛋白质与结构方向目标：理解蛋白质结构预测的实用方法，掌握 UniProt/PDB/ChEMBL 三大数据库的检索和 API 用法预计完成时间：2-3 周

3.3-protein-structure-tutorial.md 结构

蛋白质结构基础：四个层次、20 种氨基酸分类及对药物发现的意义

AlphaFold2 三种使用方式：

ColabFold 在线操作步骤 + 输出文件解读（PAE 图含义）
本地部署命令
ESMFold API（Python 一键预测）
PDB 文件格式逐字段解析 + BioPython 解析代码
pLDDT 颜色编码含义，以及"低 pLDDT = 无序区而非预测失败"的关键认知

UniProt API：搜索、下载条目 JSON、提取功能区域和疾病突变、批量下载蛋白家族序列

PDB RCSB API：按基因名搜索结构（按分辨率排序）、下载 PDB 文件、结构 RMSD 比对

ChEMBL API：靶点搜索、活性数据获取、pIC50 计算、构建 AI 训练数据集

三库联动：从基因名自动找最适合药物对接的 PDB 结构（有配体 + 最高分辨率）

py3Dmol 可视化（Jupyter 内嵌 3D 结构查看，按 pLDDT 着色）

AlphaFold2 vs AlphaFold3 对比，重点解释 AF3 对药物发现的意义

先建立直觉：为什么蛋白质结构很重要？

基因组告诉你"哪个基因有突变"，但真正起作用的是蛋白质。蛋白质的三维结构决定了它的功能——包括它和哪些分子结合、是否可以被药物抑制。

基因序列 → 氨基酸序列（一维）→ 三维折叠结构 → 功能
"一串文字"    "一串字母"             "立体雕塑"      "作用"

传统方法（X 射线晶体学、冷冻电镜）解析一个蛋白质结构需要数年和数百万美元。2021 年 AlphaFold2 的出现让这个过程缩短到几分钟，被称为生命科学的"GPT 时刻"。

3.3.1 蛋白质结构基础

四个层次的结构

理解 AlphaFold 的输出前，需要知道蛋白质结构的四个层次：

层次	含义	类比
一级结构	氨基酸序列（线性）	字符串
二级结构	局部折叠（α 螺旋、β 折叠）	局部的代码模式
三级结构	整条链的三维折叠	函数的完整实现
四级结构	多条链的组装（蛋白质复合物）	微服务之间的协作

二级结构的两种主要元素：

α 螺旋（Alpha Helix）：像弹簧，氢键在链内部形成
β 折叠（Beta Sheet）：像折叠纸，氢键在链之间形成

蛋白质 = 这些元素 + 连接它们的环状区域（Loop）的组合

氨基酸的 20 种分类

蛋白质由 20 种氨基酸组成，它们的化学性质决定了蛋白质的折叠和功能：

类别	代表氨基酸	性质	在结构中的作用
疏水性	Leu, Val, Ile, Phe	不喜水	埋在蛋白质核心，驱动折叠
极性	Ser, Thr, Asn, Gln	可形成氢键	蛋白质表面和结合位点
带正电	Lys, Arg, His	正电	与 DNA/RNA 结合
带负电	Asp, Glu	负电	活性位点，金属离子配位
特殊	Cys（二硫键）, Pro（刚性），Gly（最灵活）		结构稳定性

为什么这对生物信息学重要：

氨基酸替换的致病性评估（Glu→Val 这类改变疏水性的替换更危险）
药物结合位点通常在疏水口袋里
ESM2 这类蛋白质语言模型，embedding 捕获的正是这些化学性质信息

3.3.2 AlphaFold2：原理与实际使用

AlphaFold2 做了什么？

AlphaFold2 的核心思路：从进化信息预测结构。

如果一个位置的氨基酸在不同物种中总是同步变化（协进化），说明它们在三维空间上可能相互接触。AlphaFold2 用多序列比对（MSA）挖掘这种协进化信号，结合 Transformer 架构，直接预测每个原子的三维坐标。

输入：一段蛋白质的氨基酸序列（如 MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNG...）
  ↓ 在数据库中搜索同源序列（MSA）
  ↓ 多头注意力机制分析氨基酸对之间的协进化关系
  ↓ 结构模块迭代精化三维坐标
输出：每个原子的 x,y,z 坐标 + 置信度分数（pLDDT）

pLDDT 置信度分数解读

AlphaFold2 输出的 PDB 文件中，B-factor 列被替换为 pLDDT（predicted Local Distance Difference Test）分数：

pLDDT 范围	颜色（官方配色）	解读
> 90	深蓝	高置信，结构非常可靠
70-90	浅蓝	可信，主链正确
50-70	黄色	低置信，需谨慎
< 50	橙色	非常低，可能是天然无序区域

重要认知：低 pLDDT ≠ 预测失败

pLDDT < 50 的区域很多时候是真实存在的内在无序区域（Intrinsically Disordered Regions, IDR）——它们在生物学中有功能，只是没有固定的三维结构。AlphaFold2 在这里的低置信度实际上是准确的——它正确地识别了"这里没有稳定结构"。

import numpy as np

def analyze_plddt(pdb_file):
    """从 PDB 文件提取 pLDDT 分数（B-factor 列）"""
    plddt_scores = []
    residue_ids  = []
    current_res  = None
    
    with open(pdb_file) as f:
        for line in f:
            if not line.startswith("ATOM"):
                continue
            res_id   = int(line[22:26].strip())
            atom     = line[12:16].strip()
            b_factor = float(line[60:66].strip())  # pLDDT 在 B-factor 列
            
            if res_id != current_res and atom == "CA":  # 每个残基只取 Cα 原子
                plddt_scores.append(b_factor)
                residue_ids.append(res_id)
                current_res = res_id
    
    scores = np.array(plddt_scores)
    print(f"总残基数：{len(scores)}")
    print(f"平均 pLDDT：{scores.mean():.1f}")
    print(f"高置信（>90）比例：{(scores > 90).mean():.1%}")
    print(f"可信（70-90）比例：{((scores >= 70) & (scores <= 90)).mean():.1%}")
    print(f"低置信（<50）比例：{(scores < 50).mean():.1%}")
    return residue_ids, scores.tolist()

方式一：ColabFold（在线，推荐入门）

ColabFold 是 AlphaFold2 的优化版本，速度比原版快 40-60 倍，可直接在 Google Colab 运行：

使用步骤：

打开 ColabFold Notebook（搜索 "ColabFold AlphaFold2"，在 Google Colab 打开）
在 query_sequence 填入氨基酸序列：
```
MAEGEITTFTALTEKFNLPPGNYKKPKLLYCSNG...
```
设置参数：
- num_relax：能量最小化次数，1 或 5 均可
- template_mode：none（不用模板，纯预测）或 pdb100（用 PDB 模板辅助）
- num_models：通常 5（跑 5 个模型，取最好的）
点击 Runtime → Run all
约 10-30 分钟后下载结果（.pdb 文件 + 置信度图）

ColabFold 输出文件：

protein_unrelaxed_rank_001_alphafold2_ptm_model_1_seed_000.pdb  ← 最佳模型
protein_coverage.png     ← MSA 覆盖度（越丰富预测越准）
protein_plddt.png        ← 每个残基的 pLDDT 分数图
protein_PAE.png          ← Predicted Aligned Error（衡量两个残基间距的不确定性）

PAE（Predicted Aligned Error）图的解读：

PAE 是一个 N×N 的矩阵，PAE[i][j] 表示"以残基 i 为参考系时，残基 j 的位置误差预期"：

低 PAE（深色）→ i 和 j 的相对位置可信 → 它们可能在同一刚性域
高 PAE（浅色）→ 相对位置不确定 → 它们可能在不同结构域之间，有灵活连接

解读方法：PAE 图的深色方块 = 蛋白质的结构域（Domain）

方式二：本地部署 AlphaFold2

需要：NVIDIA GPU（≥ 40GB 显存，如 A100）、约 2TB 存储空间、Ubuntu

# 克隆官方仓库
git clone https://github.com/google-deepmind/alphafold.git
cd alphafold

# 下载数据库（约 2TB，需要专业网络环境）
bash scripts/download_all_data.sh /data/alphafold_data/

# 运行预测
python run_alphafold.py \
    --fasta_paths=protein.fasta \
    --max_template_date=2023-01-01 \
    --model_preset=monomer \           # monomer（单链）或 multimer（多链复合物）
    --data_dir=/data/alphafold_data/ \
    --output_dir=./results/

方式三：ESMFold（快速在线 API）

Meta 的 ESMFold 不需要 MSA，预测速度极快（秒级），适合批量预测：

import requests

def esmfold_predict(sequence):
    """调用 ESMFold API 预测蛋白质结构"""
    url = "https://api.esmatlas.com/foldSequence/v1/pdb/"
    headers = {"Content-Type": "application/x-www-form-urlencoded"}
    response = requests.post(url, data=sequence, headers=headers, timeout=300)
    
    if response.status_code == 200:
        return response.text  # PDB 格式的字符串
    else:
        raise RuntimeError(f"API 错误: {response.status_code}")

# 预测并保存
sequence = "MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQDNLSGAEKAVQVKVKALPDAQFEVVHSLAKWKRQTLGQHDFSAGEGLYTHMKALRPDEDRLSPLHSVYVDQWDWERVMGDGERQFSTLKSTVEAIWAGIKATEAAVSEEFGLAPFLPDQIHFVHSQELLSRYPDLDAKGRERAIAKDLGAVFLVGIGGKLSDGHRHDVRAPDYDDWSTPSELGHAGLNGDILVWNPVLEDAFELSSMGIRVDADTLKHQLALTGDEDRLELEWHQALLRGEMPQTIGGGIGQSRLTMLLLQLPHIGQVQAGVWPAAVRESVPSLL"

pdb_content = esmfold_predict(sequence)
with open("esmfold_result.pdb", "w") as f:
    f.write(pdb_content)
print("结构已保存为 esmfold_result.pdb")

PDB 文件格式解析

预测或下载的结构文件都是 PDB 格式：

HEADER    HYDROLASE                               01-JUN-21   7L6R
REMARK   1 ALPHAFOLD MONOMER V2.0 MODEL
SEQRES   1 A  263  MET ALA GLU GLY GLU ILE THR THR PHE THR...
ATOM      1  N   MET A   1      27.340  24.430   2.614  1.00 92.77
ATOM      2  CA  MET A   1      26.266  25.413   2.842  1.00 92.77
ATOM      3  C   MET A   1      26.913  26.639   3.531  1.00 92.77
ATOM      4  O   MET A   1      27.886  26.463   4.243  1.00 92.77
...

每个 ATOM 行的字段：

列	内容	示例
1-6	记录类型	ATOM
7-11	原子序号	1
13-16	原子名称	N, CA, C, O（骨架）
18-20	残基名称	MET（蛋氨酸）
22	链 ID	A
23-26	残基序号	1
31-38	X 坐标（Å）	27.340
39-46	Y 坐标（Å）	24.430
47-54	Z 坐标（Å）	2.614
55-60	占有率	1.00
61-66	B-factor	92.77（AlphaFold 中是 pLDDT）

# 用 BioPython 解析 PDB 文件
from Bio import PDB
import numpy as np

parser = PDB.PDBParser(QUIET=True)
structure = parser.get_structure("protein", "alphafold_result.pdb")

# 遍历残基
model = structure[0]  # AlphaFold 只有一个 model
chain = model["A"]    # 链 A

residues_info = []
for residue in chain.get_residues():
    if PDB.is_aa(residue):  # 只取氨基酸（排除水分子等）
        res_id   = residue.get_id()[1]
        res_name = residue.get_resname()
        plddt    = residue["CA"].get_bfactor()  # Cα 原子的 B-factor = pLDDT
        ca_coord = residue["CA"].get_vector()    # Cα 的三维坐标
        
        residues_info.append({
            "residue_id":   res_id,
            "residue_name": res_name,
            "plddt":        plddt,
            "x": ca_coord[0], "y": ca_coord[1], "z": ca_coord[2]
        })

import pandas as pd
df = pd.DataFrame(residues_info)
print(df.head())
print(f"\n平均 pLDDT: {df['plddt'].mean():.1f}")

# 计算两个残基之间的距离
def ca_distance(chain, res1_id, res2_id):
    diff = chain[res1_id]["CA"].get_vector() - chain[res2_id]["CA"].get_vector()
    return diff.norm()

print(f"\n残基 10 和 50 的 Cα 距离: {ca_distance(chain, 10, 50):.2f} Å")

3.3.3 UniProt：蛋白质序列与功能数据库

数据库分层

UniProt 是目前最权威的蛋白质序列数据库，分为三个层次：

UniProtKB（UniProt Knowledgebase）
├── Swiss-Prot：人工审核条目，约 570,000 个，质量极高
│   命名规则：P53_HUMAN（蛋白质名_物种）
│   典型 AC：P04637（TP53 的 UniProt ID）
└── TrEMBL：自动注释条目，约 2.5 亿个，未经人工验证
    
UniRef：序列聚类（UniRef100/90/50，按序列相似度聚类）
UniParc：序列存档（不同数据库的交叉引用）

核心原则：找已知功能的人类蛋白质用 Swiss-Prot；大规模同源搜索用 UniRef90。

一条 Swiss-Prot 记录包含什么

以 TP53（p53 肿瘤抑制蛋白）为例（AC：P04637）：

基本信息：
  Entry name:   P53_HUMAN
  Accession:    P04637（主要 AC，永久不变）
  Protein:      Cellular tumor antigen p53
  Gene:         TP53
  Organism:     Homo sapiens
  Length:       393 AA
  
功能注释：
  Function:     Acts as a tumor suppressor...
  Subcellular location: Nucleus
  Tissue specificity:   Ubiquitous
  
序列特征（Feature annotations）：
  Region:       1-42     Interaction with MDM2
  Region:       94-292   DNA-binding domain（核心功能区）
  Region:       364-393  Tetramerization domain
  Natural variant: 175 R→H (in LFS)   ← 已知致病突变
  
数据库交叉引用：
  OMIM:         191170   ← 遗传病数据库
  PDB:          2OCJ, 3KMD...  ← 已有的实验结构
  Reactome:     R-HSA-5633007  ← 通路数据库
  DrugBank:     DB00975         ← 靶向该蛋白的药物

UniProt API

import requests

BASE = "https://rest.uniprot.org"

def get_uniprot_entry(accession):
    """获取 UniProt 条目的 JSON 数据"""
    url = f"{BASE}/uniprotkb/{accession}.json"
    r = requests.get(url)
    r.raise_for_status()
    return r.json()

def search_uniprot(query, organism="9606", reviewed=True, limit=10):
    """搜索 UniProt（9606 = 人类，reviewed=True 表示 Swiss-Prot）"""
    params = {
        "query":  f"({query}) AND (organism_id:{organism})" + (" AND (reviewed:true)" if reviewed else ""),
        "format": "json",
        "size":   limit,
        "fields": "accession,id,gene_names,protein_name,length,sequence"
    }
    r = requests.get(f"{BASE}/uniprotkb/search", params=params)
    r.raise_for_status()
    return r.json()["results"]

# 获取 TP53 的完整信息
tp53 = get_uniprot_entry("P04637")
print(f"蛋白质名：{tp53['proteinDescription']['recommendedName']['fullName']['value']}")
print(f"基因名：{tp53['genes'][0]['geneName']['value']}")
print(f"序列长度：{tp53['sequence']['length']} AA")
print(f"功能描述：{tp53['comments'][0]['texts'][0]['value'][:200]}...")

# 提取序列
sequence = tp53["sequence"]["value"]
print(f"\nTP53 序列（前50个氨基酸）：{sequence[:50]}")

# 提取已知的功能区域
features = tp53.get("features", [])
for feat in features:
    if feat.get("type") in ["Region", "Domain", "Active site"]:
        loc = feat["location"]
        start = loc.get("start", {}).get("value", "?")
        end   = loc.get("end",   {}).get("value", "?")
        desc  = feat.get("description", "")
        print(f"  [{start}-{end}] {feat['type']}: {desc}")

# 批量下载同一基因家族的序列（用于后续 MSA 或 ESM2 分析）
def download_family_sequences(gene_name, limit=50):
    results = search_uniprot(f"gene:{gene_name}", reviewed=True, limit=limit)
    
    fasta_lines = []
    for entry in results:
        ac   = entry["primaryAccession"]
        name = entry.get("genes", [{}])[0].get("geneName", {}).get("value", "unknown")
        org  = entry.get("organism", {}).get("scientificName", "unknown")
        seq  = entry.get("sequence", {}).get("value", "")
        fasta_lines.append(f">{ac}|{name}|{org}\n{seq}")
    
    output = "\n".join(fasta_lines)
    with open(f"{gene_name}_family.fasta", "w") as f:
        f.write(output)
    return output

download_family_sequences("BRCA1", limit=20)

3.3.4 PDB：实验蛋白质三维结构数据库

数据库概况

PDB（Protein Data Bank）收录了约 220,000 个实验解析的蛋白质和核酸结构（截至 2024 年）。

实验技术分布：

方法	占比	特点
X 射线晶体学	~85%	分辨率高（1-2 Å），需要结晶
冷冻电镜（cryo-EM）	~10%，快速增长	可解析大复合物，不需要结晶
核磁共振（NMR）	~5%	可看溶液中的动态运动

分辨率含义：

1 Å（埃）= 0.1 nm = 0.0000000001 m

< 2 Å：原子分辨率，能看清每个原子位置
2-3 Å：侧链构象可靠，通常足够药物设计
3-4 Å：主链可靠，侧链不确定
> 4 Å：只能看到整体形状

PDB ID 和文件格式

每个结构有一个 4 字符的 PDB ID（如 2OCJ）。下载后有两种主要格式：

PDB 格式（传统，逐渐被取代）：

ATOM      1  N   MET A   1      27.340  24.430   2.614  1.00 30.00           N

mmCIF 格式（现代标准，支持更大的结构）：

_atom_site.group_PDB
_atom_site.id
_atom_site.type_symbol
ATOM  1  N  N  . MET A 1 1  ? 27.340 24.430 2.614 1.00 30.00 ? 1 MET A N 1

RCSB PDB API

import requests

PDB_API = "https://data.rcsb.org/rest/v1"

def get_pdb_info(pdb_id):
    """获取 PDB 条目的元数据"""
    url = f"{PDB_API}/core/entry/{pdb_id}"
    r = requests.get(url)
    r.raise_for_status()
    return r.json()

def search_pdb(gene_name, organism="Homo sapiens", max_results=10):
    """按基因名搜索 PDB 结构"""
    query = {
        "query": {
            "type": "group",
            "logical_operator": "and",
            "nodes": [
                {
                    "type": "terminal",
                    "service": "text",
                    "parameters": {
                        "attribute": "rcsb_entity_source_organism.taxonomy_lineage.name",
                        "operator": "contains_phrase",
                        "value": organism
                    }
                },
                {
                    "type": "terminal",
                    "service": "text",
                    "parameters": {
                        "attribute": "rcsb_polymer_entity.rcsb_gene_name.value",
                        "operator": "exact_match",
                        "value": gene_name
                    }
                }
            ]
        },
        "return_type": "entry",
        "request_options": {
            "results_content_type": ["experimental"],
            "sort": [{"sort_by": "rcsb_entry_info.resolution_combined", "direction": "asc"}],
            "paginate": {"start": 0, "rows": max_results}
        }
    }
    r = requests.post("https://search.rcsb.org/rcsbsearch/v2/query", json=query)
    r.raise_for_status()
    return [hit["identifier"] for hit in r.json().get("result_set", [])]

def download_pdb(pdb_id, format="pdb"):
    """下载 PDB 结构文件"""
    if format == "pdb":
        url = f"https://files.rcsb.org/download/{pdb_id}.pdb"
    else:  # mmcif
        url = f"https://files.rcsb.org/download/{pdb_id}.cif"
    r = requests.get(url)
    r.raise_for_status()
    filename = f"{pdb_id}.{format}"
    with open(filename, "w") as f:
        f.write(r.text)
    return filename

# 搜索 TP53 的所有实验结构，按分辨率排序
tp53_pdbs = search_pdb("TP53")
print(f"TP53 的 PDB 结构（分辨率从高到低）：{tp53_pdbs}")

# 获取最高分辨率结构的详细信息
best_pdb = tp53_pdbs[0]
info = get_pdb_info(best_pdb)
print(f"\nPDB ID: {best_pdb}")
print(f"分辨率: {info['rcsb_entry_info'].get('resolution_combined', 'N/A')} Å")
print(f"实验方法: {info['rcsb_entry_info']['experimental_method']}")
print(f"发布年份: {info['rcsb_accession_info']['deposit_date'][:4]}")

# 下载结构
download_pdb(best_pdb)

结构比对：实验结构 vs AlphaFold 预测

from Bio import PDB
from Bio.PDB.Superimposer import Superimposer

def compare_structures(pdb_file1, pdb_file2, chain_id="A"):
    """计算两个结构的 RMSD（均方根偏差）"""
    parser = PDB.PDBParser(QUIET=True)
    s1 = parser.get_structure("s1", pdb_file1)
    s2 = parser.get_structure("s2", pdb_file2)
    
    # 提取 Cα 原子（对齐用）
    atoms1 = [res["CA"] for res in s1[0][chain_id] if PDB.is_aa(res) and "CA" in res]
    atoms2 = [res["CA"] for res in s2[0][chain_id] if PDB.is_aa(res) and "CA" in res]
    
    # 取共同长度
    min_len = min(len(atoms1), len(atoms2))
    atoms1, atoms2 = atoms1[:min_len], atoms2[:min_len]
    
    sup = Superimposer()
    sup.set_atoms(atoms1, atoms2)  # 计算最优叠加矩阵
    
    rmsd = sup.rms
    print(f"比对残基数：{min_len}")
    print(f"RMSD：{rmsd:.3f} Å")
    print(f"解读：{'结构高度相似' if rmsd < 1 else '主链一致' if rmsd < 2 else '总体折叠相似' if rmsd < 5 else '结构差异较大'}")
    return rmsd

# 比较实验结构和 AlphaFold 预测
rmsd = compare_structures("2OCJ.pdb", "alphafold_tp53.pdb")

3.3.5 ChEMBL：药物与靶点活性数据库

数据库定位

ChEMBL 是 EMBL-EBI 维护的大型生物活性数据库，收录了：

约 230 万个小分子化合物
约 1.5 万个蛋白质靶点
约 1900 万条活性数据（化合物对靶点的结合/抑制数据）

这是药物发现和虚拟筛选的核心数据源。

核心概念

概念	含义	示例
ChEMBL ID	化合物/靶点唯一标识	CHEMBL941（Imatinib，格列卫）
Target	蛋白质靶点	BCR-ABL（CHEMBL2366516）
Assay	实验活性测定	IC50 测定实验
IC50	半数抑制浓度（nM）	格列卫对 BCR-ABL：< 10 nM
Ki	抑制常数	越小越强效
pIC50	-log10(IC50)，越大越强	pIC50 > 6 (即 IC50 < 1µM) 算活性化合物

ChEMBL API

import requests
from chembl_webresource_client.new_client import new_client

# 安装：pip install chembl_webresource_client

# ── 搜索靶点 ──
target_api = new_client.target

# 按基因名搜索靶点
targets = target_api.filter(
    target_synonym__icontains="EGFR"
).filter(
    organism="Homo sapiens"
)[:5]

for t in targets:
    print(f"ChEMBL ID: {t['target_chembl_id']}")
    print(f"  名称: {t['pref_name']}")
    print(f"  类型: {t['target_type']}")
    print()

# ── 获取某靶点的所有活性化合物 ──
activity_api = new_client.activity

egfr_id = "CHEMBL203"  # EGFR 的 ChEMBL ID
activities = activity_api.filter(
    target_chembl_id = egfr_id,
    standard_type    = "IC50",       # 只要 IC50 数据
    standard_relation= "=",          # 精确值，排除 > 或 <
    assay_type       = "B"           # B=结合实验（Binding），不用 F=功能实验
).only(["molecule_chembl_id", "canonical_smiles", "standard_value", "standard_units"])[:100]

import pandas as pd
df = pd.DataFrame(list(activities))
df["standard_value"] = pd.to_numeric(df["standard_value"], errors="coerce")
df = df.dropna(subset=["standard_value"])
df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)  # nM 转 M，再取 pIC50

print(f"\nEGFR IC50 数据条数：{len(df)}")
print(f"活性化合物（pIC50 > 6）：{(df['pIC50'] > 6).sum()}")
print(df.sort_values("pIC50", ascending=False).head())

# ── 获取化合物结构和属性 ──
molecule_api = new_client.molecule

# Imatinib（格列卫，第一个靶向抗癌药）
imatinib = molecule_api.get("CHEMBL941")
print(f"化合物名: {imatinib['pref_name']}")
print(f"分子量: {imatinib['molecule_properties']['full_mwt']}")
print(f"SMILES: {imatinib['molecule_structures']['canonical_smiles']}")

# Lipinski 五规则相关属性
props = imatinib["molecule_properties"]
print(f"\nLipinski 五规则（口服药物）：")
print(f"  分子量: {props['full_mwt']} (≤500)")
print(f"  logP: {props['alogp']} (≤5)")
print(f"  氢键供体: {props['hbd']} (≤5)")
print(f"  氢键受体: {props['hba']} (≤10)")

# ── 构建靶点的活性数据集（用于 AI 药物发现模型训练）──
def build_activity_dataset(target_chembl_id, activity_type="IC50", min_count=100):
    """
    构建指定靶点的活性数据集
    返回 DataFrame，含 SMILES 和活性标签（pIC50）
    """
    acts = activity_api.filter(
        target_chembl_id = target_chembl_id,
        standard_type    = activity_type,
        standard_relation= "=",
        assay_type       = "B"
    ).only([
        "molecule_chembl_id",
        "canonical_smiles",
        "standard_value",
        "standard_units"
    ])
    
    df = pd.DataFrame(list(acts))
    df["standard_value"] = pd.to_numeric(df["standard_value"], errors="coerce")
    df = df.dropna(subset=["canonical_smiles", "standard_value"])
    df = df[df["standard_units"] == "nM"]
    df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)
    df = df[df["pIC50"].between(0, 12)]  # 过滤异常值
    
    if len(df) < min_count:
        print(f"警告：数据量不足（{len(df)} 条），模型质量可能受限")
    
    return df[["molecule_chembl_id", "canonical_smiles", "pIC50"]].drop_duplicates()

# 构建 EGFR 抑制剂数据集（用于第四阶段 AI 项目）
egfr_dataset = build_activity_dataset("CHEMBL203")
egfr_dataset.to_csv("egfr_activity_dataset.csv", index=False)
print(f"EGFR 数据集：{len(egfr_dataset)} 个化合物")

3.3.6 三个数据库的联动使用

实际工作中，这三个数据库经常需要串联使用：

场景：找某个靶点的最佳结构用于药物对接

import requests
from chembl_webresource_client.new_client import new_client
from Bio import PDB

def find_best_structure_for_docking(gene_name):
    """
    给定基因名，找最适合药物对接的 PDB 结构
    优先：有配体的高分辨率结构
    """
    # 第一步：从 UniProt 获取该蛋白的 PDB 交叉引用
    results = search_uniprot(f"gene:{gene_name}", reviewed=True, limit=1)
    if not results:
        return None
    
    ac = results[0]["primaryAccession"]
    entry = get_uniprot_entry(ac)
    
    pdb_ids = [
        xref["id"]
        for xref in entry.get("uniProtKBCrossReferences", [])
        if xref["database"] == "PDB"
    ]
    print(f"{gene_name}（{ac}）关联的 PDB 结构：{pdb_ids[:10]}...")
    
    # 第二步：从 RCSB 获取每个结构的分辨率，找最高分辨率且有配体的
    best = {"pdb_id": None, "resolution": 99}
    for pdb_id in pdb_ids[:20]:  # 最多检查 20 个
        try:
            info = get_pdb_info(pdb_id)
            res  = info["rcsb_entry_info"].get("resolution_combined", 99)
            has_ligand = info["rcsb_entry_info"].get("nonpolymer_entity_count", 0) > 0
            
            if has_ligand and res < best["resolution"]:
                best = {"pdb_id": pdb_id, "resolution": res}
        except:
            continue
    
    if best["pdb_id"]:
        print(f"\n推荐用于对接的结构：{best['pdb_id']}（分辨率 {best['resolution']:.2f} Å，含配体）")
        download_pdb(best["pdb_id"])
    
    return best

find_best_structure_for_docking("EGFR")

场景：验证 AlphaFold 预测区域的可信度

def validate_prediction_region(pdb_id, start_res, end_res):
    """
    给定 AlphaFold 预测的 PDB 文件和残基范围
    报告该区域的 pLDDT 统计，判断预测是否可靠
    """
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure("af", pdb_id)
    chain = structure[0]["A"]
    
    region_plddts = []
    for res in chain.get_residues():
        if not PDB.is_aa(res):
            continue
        res_id = res.get_id()[1]
        if start_res <= res_id <= end_res:
            plddt = res["CA"].get_bfactor()
            region_plddts.append(plddt)
    
    if not region_plddts:
        print(f"残基 {start_res}-{end_res} 未找到")
        return
    
    avg = np.mean(region_plddts)
    print(f"残基 {start_res}-{end_res} 的 pLDDT 统计：")
    print(f"  平均值：{avg:.1f}")
    print(f"  最低：{min(region_plddts):.1f}")
    
    if avg > 70:
        print(f"  结论：预测可信，可用于功能分析或药物对接")
    else:
        print(f"  结论：预测置信度低，该区域可能为无序区，慎用")

3.3.7 可视化：用 py3Dmol 在 Jupyter 中查看结构

# pip install py3Dmol
import py3Dmol

def view_structure(pdb_id_or_file, style="cartoon", color_by="plddt"):
    """在 Jupyter Notebook 中可视化蛋白质结构"""
    view = py3Dmol.view(width=800, height=600)
    
    if len(pdb_id_or_file) == 4:  # 是 PDB ID
        view.addModel(requests.get(f"https://files.rcsb.org/download/{pdb_id_or_file}.pdb").text, "pdb")
    else:  # 是本地文件路径
        with open(pdb_id_or_file) as f:
            view.addModel(f.read(), "pdb")
    
    if color_by == "plddt":
        # 按 pLDDT 分数着色（AlphaFold 结构）
        view.setStyle({"cartoon": {
            "colorscheme": {
                "prop": "b",
                "gradient": "roygb",  # 红橙黄绿蓝
                "min": 50,
                "max": 90
            }
        }})
    elif style == "cartoon":
        view.setStyle({"cartoon": {"color": "spectrum"}})
    elif style == "surface":
        view.addSurface(py3Dmol.VDW, {"opacity": 0.7, "color": "white"})
        view.setStyle({"cartoon": {"color": "spectrum"}})
    
    # 显示配体（如果有）
    view.setStyle({"hetflag": True}, {"stick": {"colorscheme": "greenCarbon", "radius": 0.3}})
    
    view.zoomTo()
    view.show()

# 查看 EGFR 和 Erlotinib 的结合结构（1IVO）
view_structure("1IVO", color_by="chain")

# 查看 AlphaFold 预测（按 pLDDT 着色）
view_structure("alphafold_result.pdb", color_by="plddt")

3.3.8 AlphaFold2 vs AlphaFold3 的区别

2024 年 Google DeepMind 发布了 AlphaFold3：

维度	AlphaFold2	AlphaFold3
能预测的分子类型	蛋白质（单体/复合物）	蛋白质 + DNA + RNA + 小分子 + 修饰
架构	Evoformer + Structure Module	Diffusion（扩散模型）
复合物预测	需要 AlphaFold-Multimer	直接支持
小分子对接	不支持	支持（直接预测蛋白-配体复合物）
核酸结构	不支持	支持
开源程度	完全开源	代码开源，但权重使用有限制
在线服务	ColabFold	AlphaFold Server（alphafoldserver.com）

对药物发现的意义： AlphaFold3 可以直接预测蛋白质和小分子的结合构象，从原来的"预测靶点结构 → 单独做分子对接"变成了一步完成。这对 AI 药物发现的影响相当于 DALL-E 出现对图像生成的影响。

本阶段完成标志

能解释 pLDDT 分数的含义，知道什么情况下低 pLDDT 是正常的
用 ColabFold 预测过至少一个蛋白质的结构，并能解读 PAE 图
能从 UniProt 查到一个蛋白质的功能区域、已知疾病突变和关联的 PDB 结构
能用 API 下载 ChEMBL 中某靶点的活性数据，构建 pIC50 数据集
能用 BioPython 解析 PDB 文件，提取每个残基的坐标和 pLDDT 分数
知道 AlphaFold2 和 AlphaFold3 的核心区别，以及 AlphaFold3 对药物发现的意义

完成后，进入 第四阶段：AI/ML 在生物信息中的应用——用 ESM2 蛋白质语言模型做迁移学习，以及图神经网络在药物发现中的应用。

第三阶段 3.3：蛋白质与结构方向

AlphaFold2 vs AlphaFold3 对比，重点解释 AF3 对药物发现的意义​

先建立直觉：为什么蛋白质结构很重要？​

3.3.1 蛋白质结构基础​

四个层次的结构​

氨基酸的 20 种分类​

3.3.2 AlphaFold2：原理与实际使用​

AlphaFold2 做了什么？​

pLDDT 置信度分数解读​

方式一：ColabFold（在线，推荐入门）​

方式二：本地部署 AlphaFold2​

方式三：ESMFold（快速在线 API）​

PDB 文件格式解析​

3.3.3 UniProt：蛋白质序列与功能数据库​

数据库分层​

一条 Swiss-Prot 记录包含什么​

UniProt API​

3.3.4 PDB：实验蛋白质三维结构数据库​

数据库概况​

PDB ID 和文件格式​

RCSB PDB API​

结构比对：实验结构 vs AlphaFold 预测​

3.3.5 ChEMBL：药物与靶点活性数据库​

数据库定位​

核心概念​

ChEMBL API​

3.3.6 三个数据库的联动使用​

3.3.7 可视化：用 py3Dmol 在 Jupyter 中查看结构​

3.3.8 AlphaFold2 vs AlphaFold3 的区别​

本阶段完成标志​