Biopython计算的溶剂可及表面积（SASA）与PyMOL结果不符的原因排查

这种工具间的SASA数值差异其实挺常见的，核心原因基本都集中在参数设置、算法实现细节以及结构处理逻辑这几个方面，咱们一个个拆解：

1. 原子半径参数不匹配 #

Biopython的ShrakeRupley类默认使用的是一套内置的原子半径，而PyMOL默认采用的是Chothia的原子半径参数（或者其自定义的变体）。不同的半径设定会直接改变每个原子的“占位”大小，进而影响可及表面积的计算——比如Gly的α碳原子，两套参数的半径差可能直接导致其暴露面积的计算结果天差地别。

你可以手动指定Biopython使用Chothia半径来对齐PyMOL的设置，示例代码如下：

from Bio.PDB.SASA import ShrakeRupley
from Bio.PDB.MMCIFParser import MMCIFParser

# 定义Chothia标准原子半径（覆盖常见蛋白质原子类型）
chothia_radii = {
    "N": 1.55, "CA": 1.7, "C": 1.7, "O": 1.52,
    "CB": 1.7, "CG": 1.7, "CG1": 1.7, "CG2": 1.7,
    "CD": 1.7, "CD1": 1.7, "CD2": 1.7, "CE": 1.7,
    "CE1": 1.7, "CE2": 1.7, "CE3": 1.7, "CZ": 1.7,
    "CZ2": 1.7, "CZ3": 1.7, "CH2": 1.7, "NH1": 1.55,
    "NH2": 1.55, "OH": 1.52, "SD": 1.85, "SG": 1.85
}

parser = MMCIFParser(QUIET=True)
structure = parser.get_structure("HELLO", "insulin.cif")
sr = ShrakeRupley(radii_dict=chothia_radii)  # 应用自定义半径
sr.compute(structure[0], level="R")

my_list = []
for chain in structure[0]:
    for res in chain:
        my_list.append((res.get_resname(), round(res.sasa, 2)))
print(my_list[0:10])

2. 探针球与采样点的差异 #

• 探针半径：虽然两者默认都用1.4Å（模拟水分子的半径），但偶尔会有工具版本的默认值变更，建议在初始化ShrakeRupley时显式指定：ShrakeRupley(probe_radius=1.4)。
• 采样点密度：Shrake-Rupley算法是通过在原子表面采样点，判断这些点是否被其他原子遮挡来计算SASA的。Biopython默认每个原子采样100个点，而PyMOL的默认采样点数量通常更高（比如1000个）。采样点越少，结果误差可能越大，你可以通过n_points参数调整：ShrakeRupley(n_points=1000)。

3. 结构处理的不一致性 #

你需要确保Biopython和PyMOL加载的是完全相同的结构：

• 氢原子：PyMOL如果加载的是带氢的结构，而Biopython默认不会处理氢原子（因为PDB/MMCIF文件通常不含氢），氢原子的存在会改变表面的可及性。可以用PyMOL导出去除氢原子的结构文件，再用Biopython重新计算。
• 模型/链的选择：你的代码里用了structure[0]（第一个模型），确认PyMOL可视化的是不是同一个模型、同一条链，有没有误选其他链或者多模型结构的情况。
• 缺失原子：检查MMCIF文件中Gly1的原子是否完整，如果Biopython加载时缺失了某些原子，也会导致SASA计算偏小。

4. 残基SASA的计算逻辑核对 #

有时候Biopython的残基级SASA可能和你预期的计算方式有差异，你可以先计算原子级的SASA，手动求和得到残基总SASA，和Biopython给出的res.sasa对比，同时和PyMOL的结果对齐：

sr.compute(structure[0], level="A")  # 切换到原子级计算
for chain in structure[0]:
    for res in chain:
        total_sasa = sum(atom.sasa for atom in res)
        print(f"{res.get_resname()} {res.id[1]}: 手动求和={round(total_sasa,2)}, Biopython残基值={round(res.sasa,2)}")

通过以上几步调整，基本能把Biopython的计算结果和PyMOL的差值缩小到可接受的范围。

内容的提问来源于stack exchange，提问作者Sad_man