如何在PyMOL中计算参考与目标蛋白比对残基的Delta Phi Psi值

背景 #

我已经用PyMOL将同一家族的多个蛋白结构与参考结构完成比对，现在需要计算参考结构和每个目标结构中每一对比对残基的Phi、Psi角差值（Delta Phi Psi）。我知道可以用phi_psi命令提取Phi/Psi值，也已经有了序列比对结果，但不知道怎么对应每对比对残基来计算差值。

结构比对基础代码（基于Super工具） #

Super工具用于比对两个选择对象，它执行与align不同的、不依赖序列的基于结构的动态规划比对，随后通过一系列优化循环消除高相对变异性的配对以提升拟合效果（和align类似），在低序列相似性蛋白比对中，Super比align更稳健。

涉及的PDB结构：参考结构3HR8，目标结构1XMV、1XP8

import pymol
from pymol import cmd
import time

list_pdbs = ['1xp8' , '3hr8', '1xmv']

pymol.finish_launching()
time.sleep(5)

# 加载所有PDB结构
for _i in list_pdbs :
    try :
        cmd.fetch(_i)
    except :
        print(f'无法加载 {_i}')

# 移除非蛋白原子
cmd.remove('hetatm')

# 基于CA原子用Super做多结构比对
sel = 'name CA'
ref = list_pdbs[1]  # 3HR8作为参考结构
method = 'super'
cmd.do(f"extra_fit {sel} , {ref} , {method}, object = aln_super_all")

# 显示序列比对视图
cmd.do('set seq_view , 1')
cmd.select('aln_super_all')
print(cmd.get_type('aln_super_all'))

计算Delta Phi Psi的解决方案 #

核心是将参考结构和目标结构的比对残基一一对应，再分别提取角度值计算差值，具体步骤如下：

步骤1：编写函数提取单个结构的Phi/Psi值 #

创建函数，输入结构名，返回该结构所有可计算角度残基的Phi/Psi字典（键为「链ID+残基号」，值为(phi, psi)）：

def get_phi_psi(structure):
    phi_psi_dict = {}
    # 遍历结构的每个蛋白残基
    for resi in cmd.get_model(f'{structure} and polymer').residues:
        res_id = f"{resi.chain}{resi.resi}"
        # 提取数值类型的phi/psi角（第二个参数为1）
        phi_val = cmd.phi_psi(f'{structure} and resi {resi.resi} and chain {resi.chain}', 1)[0]
        psi_val = cmd.phi_psi(f'{structure} and resi {resi.resi} and chain {resi.chain}', 1)[1]
        # 跳过无法计算角度的末端残基
        if phi_val is not None and psi_val is not None:
            phi_psi_dict[res_id] = (phi_val, psi_val)
    return phi_psi_dict

步骤2：获取比对残基的映射关系 #

通过PyMOL的get_raw_alignment命令提取比对对象中残基的配对关系：

# 获取比对的残基映射列表，每个元素是(参考残基信息, 目标残基信息)
alignment = cmd.get_raw_alignment('aln_super_all')

步骤3：遍历比对残基计算差值 #

结合残基映射和Phi/Psi字典，逐个计算并输出差值：

# 获取参考结构的Phi/Psi数据
ref_phi_psi = get_phi_psi(ref)

# 遍历所有目标结构
targets = [pdb for pdb in list_pdbs if pdb != ref]
for target in targets:
    print(f"\n=== {ref} 与 {target} 的 Delta Phi Psi 结果 ===")
    target_phi_psi = get_phi_psi(target)
    # 遍历每一对比对残基
    for ref_res, target_res in alignment:
        ref_res_id = f"{ref_res[0]}{ref_res[1]}"
        target_res_id = f"{target_res[0]}{target_res[1]}"
        # 仅处理双方都有有效角度的残基对
        if ref_res_id in ref_phi_psi and target_res_id in target_phi_psi:
            r_phi, r_psi = ref_phi_psi[ref_res_id]
            t_phi, t_psi = target_phi_psi[target_res_id]
            delta_phi = round(t_phi - r_phi, 2)
            delta_psi = round(t_psi - r_psi, 2)
            print(f"残基对: {ref_res_id} ↔ {target_res_id} | ΔPhi: {delta_phi}° | ΔPsi: {delta_psi}°")

完整整合代码 #

将上述计算部分整合到基础代码末尾：

import pymol
from pymol import cmd
import time

list_pdbs = ['1xp8' , '3hr8', '1xmv']

pymol.finish_launching()
time.sleep(5)

# 加载所有PDB结构
for _i in list_pdbs :
    try :
        cmd.fetch(_i)
    except :
        print(f'无法加载 {_i}')

# 移除非蛋白原子
cmd.remove('hetatm')

# 基于CA原子用Super做多结构比对
sel = 'name CA'
ref = list_pdbs[1]  # 3HR8作为参考结构
method = 'super'
cmd.do(f"extra_fit {sel} , {ref} , {method}, object = aln_super_all")

# 显示序列比对视图
cmd.do('set seq_view , 1')
cmd.select('aln_super_all')

# -------------------------- Delta Phi Psi计算部分 --------------------------
def get_phi_psi(structure):
    phi_psi_dict = {}
    for resi in cmd.get_model(f'{structure} and polymer').residues:
        res_id = f"{resi.chain}{resi.resi}"
        phi_val = cmd.phi_psi(f'{structure} and resi {resi.resi} and chain {resi.chain}', 1)[0]
        psi_val = cmd.phi_psi(f'{structure} and resi {resi.resi} and chain {resi.chain}', 1)[1]
        if phi_val is not None and psi_val is not None:
            phi_psi_dict[res_id] = (phi_val, psi_val)
    return phi_psi_dict

# 获取比对残基映射
alignment = cmd.get_raw_alignment('aln_super_all')

# 计算并输出差值
ref_phi_psi = get_phi_psi(ref)
targets = [pdb for pdb in list_pdbs if pdb != ref]

for target in targets:
    print(f"\n=== {ref} 与 {target} 的 Delta Phi Psi 结果 ===")
    target_phi_psi = get_phi_psi(target)
    for ref_res, target_res in alignment:
        ref_res_id = f"{ref_res[0]}{ref_res[1]}"
        target_res_id = f"{target_res[0]}{target_res[1]}"
        if ref_res_id in ref_phi_psi and target_res_id in target_phi_psi:
            r_phi, r_psi = ref_phi_psi[ref_res_id]
            t_phi, t_psi = target_phi_psi[target_res_id]
            delta_phi = round(t_phi - r_phi, 2)
            delta_psi = round(t_psi - r_psi, 2)
            print(f"残基对: {ref_res_id} ↔ {target_res_id} | ΔPhi: {delta_phi}° | ΔPsi: {delta_psi}°")

关键说明 #

• get_raw_alignment返回的每个元素是两个元组，分别对应参考/目标结构的残基信息（格式为(链ID, 残基号, 残基名)）
• phi_psi命令第二个参数设为1时，返回数值类型的角度，方便后续计算
• 代码自动跳过无法计算Phi/Psi角的末端残基

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