在蛋白质计算设计中，ProteinMPNN 凭借其强大的序列生成能力，能够快速给出成百上千个潜在的突变序列。然而，ProteinMPNN 无法直接给出物理意义上的结合自由能变化（$\Delta\Delta G_{\text{bind}}$）。
为了筛选出真正高亲和力的突变体，我们通常需要借助经典力场工具 FoldX。
本文将分享一套成熟的自动化工作流，指导你如何将 ProteinMPNN 的输出自动转化为 FoldX 的突变输入，并批量计算复合物的结合自由能变化 $\Delta\Delta G_{\text{bind}}$。
核心逻辑与工作流设计
计算结合自由能变化（$\Delta\Delta G_{\text{bind}}$）与计算稳定性变化（$\Delta\Delta G_{\text{stability}}$）在 FoldX 中的运行逻辑不同。
- 稳定性 $\Delta\Delta G_{\text{stability}}$：直接通过
BuildModel命令变异后输出的Average_MutantModel.fxout（或Raw_文件）中的总能量差值获取。 - 结合能 $\Delta\Delta G_{\text{bind}}$：必须先使用
BuildModel生成突变体的 3D 结构，然后对野生型（WT）和突变体（MT）的复合物结构分别运行AnalyseComplex命令。两者输出的Interaction Energy之差，才是真正的 $\Delta\Delta G_{\text{bind}}$。
完整管道设计如下：
[WT.pdb] -> FoldX RepairPDB -> [WT_Repair.pdb]
|
[ProteinMPNN FASTA/JSON] --------> [Python 差异比对脚本] -> [individual_list.txt]
|
FoldX BuildModel
|
[Mutant PDBs] & [WT_Repair.pdb]
|
FoldX AnalyseComplex
|
[Python 结果提取] -> [ddG 排序报表]
步骤一：野生型（WT）结构修复
FoldX 对原子碰撞非常敏感。在进行任何突变之前，必须对原始 PDB 进行能量最小化（RepairPDB）。
准备一个文件夹 workspace，将 foldx 可执行文件、WT.pdb 放进去，运行：
# 修复 PDB 结构
./foldx --command=RepairPDB --pdb=WT.pdb
运行后会生成 WT_Repair.pdb，后续所有的突变和能量比对都必须基于这个修复后的结构。
步骤二：解析 ProteinMPNN 输出并生成突变文件
ProteinMPNN 通常会输出一个 .fa (FASTA) 文件，其中第一条是野生型序列，后续是设计出来的突变体序列。
我们需要写一个 Python 脚本，将这些设计序列与野生型序列进行比对，识别出突变位点，并翻译成 FoldX 识别的 individual_list.txt 格式。
警惕陷阱：ProteinMPNN 输出的 FASTA 索引（0, 1, 2...）与 PDB 文件中的真实残基编号（Residue Number，可能从 105 开始，或者中间有缺失）往往不一致。
下面的脚本假设你已经确保了 FASTA 序列与 PDB 链的残基是一一对应的。如果是复杂结构，建议使用Biopython读取 PDB 的真实res_id进行映射。
新建 generate_mutants.py：
import os
import re
def parse_fasta(fasta_path):
"""解析 ProteinMPNN 的 FASTA 文件"""
sequences = []
headers = []
with open(fasta_path, 'r') as f:
current_seq = []
for line in f:
line = line.strip()
if line.startswith('>'):
if current_seq:
sequences.append("".join(current_seq))
current_seq = []
headers.append(line)
else:
current_seq.append(line)
if current_seq:
sequences.append("".join(current_seq))
return headers, sequences
def generate_foldx_mutants(wt_seq, mut_seq, chain_id, pdb_start_num=1):
"""
对比 WT 和 Mutant 序列，生成 FoldX 格式的突变字符串
例如: RA102D (野生Arg, A链, 102号位置, 突变为Asp)
"""
mutations = []
for i, (wt_aa, mut_aa) in enumerate(zip(wt_seq, mut_seq)):
if wt_aa != mut_aa:
pdb_pos = i + pdb_start_num # 映射到 PDB 真实编号
mut_str = f"{wt_aa}{chain_id}{pdb_pos}{mut_aa}"
mutations.append(mut_str)
return mutations
# 配置参数
fasta_file = "mpnn_output.fa" # ProteinMPNN 输出路径
chain_to_mutate = "A" # 发生突变的链 ID
pdb_start_residue = 1 # 该链在 PDB 中的起始残基编号
headers, seqs = parse_fasta(fasta_file)
wt_sequence = seqs[0] # 第一条通常是 WT
# 写入 individual_list.txt 供 FoldX 批量读取
# 格式要求：每一行代表一个突变体，多个突变位点用逗号隔开，分号结尾。
with open("individual_list.txt", "w") as f_out:
for idx, mut_seq in enumerate(seqs[1:], start=1):
muts = generate_foldx_mutants(wt_sequence, mut_seq, chain_to_mutate, pdb_start_residue)
if muts:
mut_line = ",".join(muts) + ";"
f_out.write(f"{mut_line}\n")
print(f"Mutant {idx}: {mut_line}")
else:
# 写入空行或占位符以保持索引一致
f_out.write(";\n")
步骤三：运行 FoldX BuildModel 产生突变体结构
FoldX 的 BuildModel 命令可以根据 individual_list.txt 批量对 WT_Repair.pdb 进行突变。
运行以下命令：
./foldx --command=BuildModel \
--pdb=WT_Repair.pdb \
--mutant-file=individual_list.txt \
--numberOfRuns=1 \
--output-file=mutant_models
注：--numberOfRuns=1 是为了提高运行速度（仅构建 1 个构象）。如果追求极高精度，可以设置为 3 或 5，最后取均值，但计算时间会成倍增加。
运行结束后，你的工作目录下会生成一系列命名形如 WT_Repair_1.pdb, WT_Repair_2.pdb ... 的文件，分别对应 individual_list.txt 中每一行的突变体。
步骤四：批量运行 AnalyseComplex 计算结合能
拿到突变体结构后，我们需要计算复合物中两条链（例如受体 A 和配体 B）之间的结合能。
首先，我们需要测定野生型的结合能：
./foldx --command=AnalyseComplex --pdb=WT_Repair.pdb --analyseComplexChains=A,B
这会生成一个 Summary_WT_Repair_AC.fxout 文件。
接着，我们写一个 Bash 脚本批量对所有生成的突变体运行 AnalyseComplex：
#!/bin/bash
# 定义链的交互界面，例如 A 链和 B 链
CHAINS="A,B"
WT_NAME="WT_Repair"
# 循环处理所有突变体 PDB (假设你有 100 个突变体)
for pdb_file in WT_Repair_[0-9]*.pdb; do
if [ -f "$pdb_file" ]; then
echo "Processing $pdb_file..."
# 运行 AnalyseComplex
./foldx --command=AnalyseComplex --pdb="$pdb_file" --analyseComplexChains=$CHAINS
fi
done
运行该脚本后，每一个突变体 PDB 都会对应输出一个 Summary_WT_Repair_X_AC.fxout 文件（其中 X 是数字）。
步骤五：提取数据并计算 $\Delta\Delta G_{\text{bind}}$
每一个 Summary_*.fxout 文件的最后一行为计算结果。我们需要提取其中的 Interaction Energy 列（通常是第 6 列，视 FoldX 版本可能微调）。
以下 Python 脚本负责提取所有数据，计算 $\Delta\Delta G_{\text{bind}} = \Delta G_{\text{mutant}} - \Delta G_{\text{wt}}$ 并输出 CSV 报表：
import os
import pandas as pd
def get_interaction_energy(fxout_path):
"""从 AnalyseComplex 的 fxout 文件中提取 Interaction Energy"""
if not os.path.exists(fxout_path):
return None
with open(fxout_path, 'r') as f:
lines = f.readlines()
for line in lines:
if line.startswith("WT_Repair") or line.startswith("WT_"):
parts = line.split('\t')
# 交互能通常在第 6 列（索引为 5），请根据你生成的 fxout 实际表头核对
try:
return float(parts[5])
except ValueError:
continue
return None
# 1. 获取野生型结合能
wt_energy = get_interaction_energy("Summary_WT_Repair_AC.fxout")
if wt_energy is None:
raise FileNotFoundError("无法找到野生型的 AnalyseComplex 输出文件，请先运行野生型计算！")
print(f"Wild-Type Binding Energy: {wt_energy} kcal/mol")
# 2. 遍历并提取所有突变体结合能
results = []
# 假设我们有 100 个生成的突变体
for i in range(1, 1000):
fxout_file = f"Summary_WT_Repair_{i}_AC.fxout"
if os.path.exists(fxout_file):
mt_energy = get_interaction_energy(fxout_file)
if mt_energy is not None:
ddg = mt_energy - wt_energy
results.append({
"Mutant_Index": i,
"DG_bind_MT": mt_energy,
"DG_bind_WT": wt_energy,
"ddG_bind": ddg
})
# 3. 输出并排序
df = pd.DataFrame(results)
# 结合能越负，说明结合越紧密。我们需要寻找 ddG_bind < 0 且绝对值尽可能大的突变体
df = df.sort_values(by="ddG_bind")
df.to_csv("protein_mpnn_foldx_ddg_results.csv", index=False)
print("\nTop 5 亲和力提升突变体:")
print(df.head(5))
避坑指南与高级优化
- 残基编号不一致（最经典灾难）
ProteinMPNN 输出的 FASTA 没有 Chain 信息，而 FoldX 严重依赖链 ID（如A）和 PDB 真实编号。如果你的 PDB 存在缺失残基（例如序列第 50 个残基在 PDB 里编号是 120），直接按索引递增会导致 FoldX 在错误的位点引入突变。务必在步骤二中使用Biopython.PDB模块来动态获取残基的真实resseq。 - 结合能变好，但稳定性崩溃了
有时候，AnalyseComplex算出来的 $\Delta\Delta G_{\text{bind}}$ 非常优秀（极负），但这是因为突变体在界面上强行形成了盐桥，而实际上这个突变严重破坏了单体本身的稳定性（导致蛋白根本无法折叠）。
建议：在筛选时，不仅要看 $\Delta\Delta G_{\text{bind}}$，还要看BuildModel输出文件中的总能量变化 $\Delta\Delta G_{\text{stability}}$。只有**两者均不显著变差（甚至都变好）**的突变体才值得进行湿实验验证。 - 多核并行加速
FoldX 本身是单线程运行的。如果你有上千个突变体，直接用 shell 循环会非常慢。建议使用 Bash 中的xargs或GNU Parallel工具进行多进程并行加速：# 使用 GNU Parallel 并行运行 AnalyseComplex，-j 8 表示开 8 个线程 ls WT_Repair_[0-9]*.pdb | parallel -j 8 "./foldx --command=AnalyseComplex --pdb={} --analyseComplexChains=A,B"
