HOOOS

白嫖 Meta 算力:无显卡如何在 Colab 快速部署 ESMFold 并搞定单点突变分析

0 11 生信搬砖工 ESMFold单点突变
Apple

做结构生物学和计算生物学的同学,或多或少都经历过被显卡支配的恐惧。想跑个 AlphaFold2,光是配环境和下载那几个 TB 的数据库就能让人崩溃,更别提本地那块瑟瑟发抖的 RTX 3060 显卡了。

其实,如果你只是想针对某个靶点蛋白做**单点突变(Single-point Mutation)**分析,预测一下突变前后的结构变化和置信度(pLDDT)差异,完全不需要动用本地显卡。

利用 Google Colab 搭配 Meta 开源的 ESMFold,我们有两种极速部署方案:

  1. 轻量 API 方案(真正零显卡,Colab CPU 就能跑,1秒出结果)
  2. 本地 Transformers 方案(利用 Colab 免费 T4 显卡,适合私密或批量序列)

下面直接上干货,手把手带你跑通“突变序列准备 -> 结构预测 -> 结构叠合比对 -> 差异分析”的全流程。


方案一:极速白嫖流(基于 Meta ESMFold 官方 API)

如果你的蛋白序列不涉及保密协议,且长度在 400 aa 以内,强烈建议使用 Meta 提供的免费 API。这种方式完全不需要显卡,在 Colab 的普通 CPU 节点上即可运行,单条序列预测只需几秒钟。

1. 环境准备

在 Colab 中新建一个 Notebook,先安装必要的辅助库(用于解析 PDB 和计算 RMSD):

!pip install -q biopython py3Dmol biotite

2. Python 预测脚本

我们编写一个函数,直接向 Meta 的 ESMFold API 发送 POST 请求,获取预测好的 PDB 文本。

import requests
import os

def query_esmfold(sequence: str, save_path: str):
    """
    通过 Meta ESMFold API 预测结构并保存为 PDB 文件
    """
    url = "https://api.esmatlas.com/fold/v1/pdb"
    headers = {"Content-Type": "text/plain"}
    
    # 过滤掉序列中的空格和换行
    clean_seq = "".join(sequence.split())
    
    print(f"正在预测序列 (长度: {len(clean_seq)} aa)...")
    response = requests.post(url, data=clean_seq, headers=headers)
    
    if response.status_code == 200:
        with open(save_path, "w") as f:
            f.write(response.text)
        print(f"保存成功: {save_path}")
    else:
        print(f"预测失败,状态码: {response.status_code}")
        print(response.text)

# 以 EGFR 激酶结构域为例
# 野生型 (Wild Type, WT)
wt_seq = "LSTAVGNIKWLDGLEEIRKIVEE"
# 突变型 (Mutant, MUT) - 假设我们将第 10 位的 I 突变为 D (I10D)
mut_seq = "LSTAVGNIKWLDGLEEDRKIVEE"

# 运行预测
os.makedirs("results", exist_ok=True)
query_esmfold(wt_seq, "results/wt.pdb")
query_esmfold(mut_seq, "results/mut.pdb")

方案二:本地部署流(基于 Colab 免费 T4 GPU + Hugging Face)

如果你的序列涉及保密,不希望上传到第三方 API,或者 API 服务偶尔不稳定,可以使用 Colab 的免费 T4 GPU 节点进行本地推理。Meta 的 ESMFold 已经集成到了 Hugging Face 的 transformers 库中,调用非常方便。

1. 切换 Colab 运行时

在 Colab 菜单栏选择:代码执行程序 -> 更改运行时类型 -> 硬件加速器选择 T4 GPU

2. 安装依赖库

由于要在 GPU 上跑模型,需要安装 transformersaccelerate

!pip install -q transformers accelerate biotite py3Dmol

3. 本地推理代码

import torch
from transformers import AutoTokenizer, EsmForProteinFolding

# 1. 加载模型(第一次加载需要下载大约 8GB 的权重,Colab 内网下载极快)
device = "cuda" if torch.cuda.is_available() else "cpu"
model_path = "facebook/esmfold_v1"

tokenizer = AutoTokenizer.from_pretrained(model_path)
model = EsmForProteinFolding.from_pretrained(model_path)
model = model.to(device)

# 优化内存占用
model.esm = model.esm.half()
torch.backends.cuda.matmul.allow_tf32 = True

def predict_local(sequence: str, save_path: str):
    clean_seq = "".join(sequence.split())
    inputs = tokenizer([clean_seq], return_tensors="pt", add_special_tokens=False)
    inputs = {k: v.to(device) for k, v in inputs.items()}
    
    with torch.no_grad():
        outputs = model(**inputs)
    
    # 提取 PDB 字符串
    pdb_string = model.output_to_pdb(outputs)[0]
    
    with open(save_path, "w") as f:
        f.write(pdb_string)
    print(f"本地预测完成,保存至: {save_path}")

# 预测 WT 和 MUT
predict_local(wt_seq, "results/wt_local.pdb")
predict_local(mut_seq, "results/mut_local.pdb")

核心骨架:如何分析突变前后的结构差异?

拿到了 wt.pdbmut.pdb 之后,我们如何定量评估单点突变对结构的影响?主要看两个指标:

  1. pLDDT 的局部变化:判断突变是否导致了局部结构的解折叠或柔性增加。
  2. RMSD(均方根偏差):定量评估突变引起的整体或局部空间位移。

1. 提取并对比 pLDDT 分数

pLDDT 是 ESMFold 对每个氨基酸预测置信度的评分(0-100)。在 PDB 文件中,这个分数被存储在 B-factor 列中。

from Bio.PDB import PDBParser
import numpy as np

def extract_plddt(pdb_file):
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("protein", pdb_file)
    plddts = []
    for model in structure:
        for chain in model:
            for residue in chain:
                # 获取 B-factor (即 pLDDT)
                for atom in residue:
                    plddts.append(atom.get_bfactor())
                    break # 每个残基取一个原子的 B-factor 即可代表该残基
    return np.array(plddts)

wt_plddt = extract_plddt("results/wt.pdb")
mut_plddt = extract_plddt("results/mut.pdb")

print(f"WT 平均 pLDDT: {wt_plddt.mean():.2f}")
print(f"MUT 平均 pLDDT: {mut_plddt.mean():.2f}")

# 打印突变位置(假设是第 17 位残基,索引 16)
mut_pos = 16 
print(f"突变位点 WT pLDDT: {wt_plddt[mut_pos]:.2f}")
print(f"突变位点 MUT pLDDT: {mut_plddt[mut_pos]:.2f}")

2. 计算结构叠合与 RMSD

我们使用 biotite 库对两个结构进行超叠加(Superposition),并计算 Cα 原子的 RMSD,以此来看突变造成的扰动。

import biotite.structure as stripe
import biotite.structure.io.pdb as pdb

# 读取 PDB
wt_struct = pdb.PDBFile.read("results/wt.pdb").get_structure(model=1)
mut_struct = pdb.PDBFile.read("results/mut.pdb").get_structure(model=1)

# 筛选 C-alpha 原子进行比对
wt_ca = wt_struct[wt_struct.atom_name == "CA"]
mut_ca = mut_struct[mut_struct.atom_name == "CA"]

# 结构叠合
fitted_mut, transformation = stripe.superimpose(wt_ca, mut_ca)

# 计算整体 RMSD
rmsd = stripe.rmsd(wt_ca, fitted_mut)
print(f"WT 与 MUT 的 Cα RMSD: {rmsd:.4f} Å")

经验法则

  • 如果 RMSD < 1.0 Å:说明突变没有引起显着的骨架变化,属于保守性突变。
  • 如果 RMSD > 2.0 Å:说明局部或整体骨架发生了明显的构象转变,需要重点关注。

3. 在 Colab 中 3D 可视化对比

使用 py3Dmol 可以在 Colab 的 Output 框里直接交互式查看两个结构的对齐情况:

import py3Dmol

view = py3Dmol.view(width=800, height=600)

# 读取 WT 并设为绿色
with open("results/wt.pdb", "r") as f:
    view.addModel(f.read(), "pdb")
view.setStyle({"model": 0}, {"cartoon": {"color": "green"}})

# 读取 MUT 并设为红色
with open("results/mut.pdb", "r") as f:
    view.addModel(f.read(), "pdb")
view.setStyle({"model": 1}, {"cartoon": {"color": "red"}})

view.zoomTo()
view.show()

避坑指南与使用建议

  1. 序列长度限制:ESMFold 适合预测 400 aa 以下的单链蛋白。当序列长度超过 800 aa 时,Colab 免费 T4 显卡的 16GB 显存极易发生 Out of Memory (OOM) 报错。
  2. 多聚体预测:ESMFold 本身对多聚体支持较弱。如果是多聚体或复合物,建议使用 : 符号连接序列(例如 SEQ1:SEQ2),但效果仍逊色于 AlphaFold-Multimer。
  3. 突变评估局限性:ESMFold 是基于单序列演化信息的,它对“活性口袋”内部微小侧链调整的灵敏度,可能不如基于物理力场的动力学模拟(MD)。建议将 ESMFold 的预测结果作为初步筛选(High-throughput screening)手段,再针对有显著结构变化的突变体进行分子动力学模拟验证。

点评评价

captcha
健康