HETATM对蛋白质Phi-Psi角的影响及脚本计算问题求助

原脚本的问题是默认只处理标准氨基酸残基，遇到HETATM（比如1j73.pdb里第8位的二氨基丁酸）时会直接跳过，同时因为Phi/Psi角的计算依赖前后残基的原子，导致HETATM前后的第7、9位残基角度也无法正常计算。下面是调整后的脚本，能覆盖所有残基的角度计算：

from Bio.PDB import PDBParser, PPBuilder
import sys

def calculate_phi_psi(pdb_file):
    # 初始化PDB解析器，关闭冗余日志
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("protein", pdb_file)
    
    # 准备输出表头和内容行
    output_lines = ["残基ID,链ID,Phi角,Psi角"]
    
    # 遍历结构中的每条链
    for chain in structure.get_chains():
        # 获取链上所有残基（包含HETATM）
        residues = list(chain.get_residues())
        total_res = len(residues)
        
        for idx in range(total_res):
            current_res = residues[idx]
            # 跳过水分子，其他配体可根据需求修改过滤规则
            if current_res.get_resname() in ["HOH", "WAT"]:
                continue
            
            phi_val = None
            psi_val = None
            
            # 计算Phi角：需要前一个残基的C，当前残基的N、CA、C
            if idx > 0:
                prev_res = residues[idx-1]
                try:
                    # 先尝试用内置方法获取
                    phi_val = current_res.get_phi_angle()
                    # 如果内置方法返回None（比如前一个是HETATM），手动计算
                    if phi_val is None:
                        if all(atom in res for res, atom in [(prev_res, "C"), (current_res, "N"), (current_res, "CA"), (current_res, "C")]):
                            phi_val = PPBuilder().calc_phi(prev_res["C"], current_res["N"], current_res["CA"], current_res["C"])
                except Exception:
                    pass
            
            # 计算Psi角：需要当前残基的N、CA、C，后一个残基的N
            if idx < total_res - 1:
                next_res = residues[idx+1]
                try:
                    psi_val = current_res.get_psi_angle()
                    if psi_val is None:
                        if all(atom in res for res, atom in [(current_res, "N"), (current_res, "CA"), (current_res, "C"), (next_res, "N")]):
                            psi_val = PPBuilder().calc_psi(current_res["N"], current_res["CA"], current_res["C"], next_res["N"])
                except Exception:
                    pass
            
            # 格式化残基ID：残序号+插入码（如果有）
            res_id_tuple = current_res.get_id()
            res_num = res_id_tuple[1]
            ins_code = res_id_tuple[2] if res_id_tuple[2] != " " else ""
            full_res_id = f"{res_num}{ins_code}"
            
            # 添加到输出列表
            output_lines.append(f"{full_res_id},{chain.get_id()},{phi_val},{psi_val}")
    
    # 写入结果文件
    with open("phi_psi_all_residues.csv", "w", encoding="utf-8") as f:
        f.write("\n".join(output_lines))
    print("计算完成，结果已保存到 phi_psi_all_residues.csv")

if __name__ == "__main__":
    if len(sys.argv) != 2:
        print("用法: python 脚本名.py <pdb文件路径>")
        sys.exit(1)
    calculate_phi_psi(sys.argv[1])

关键修改说明 #

• 不再过滤HETATM残基，仅跳过无关的水分子，确保修饰氨基酸也被纳入计算
• 针对内置角度计算方法因前后残基是HETATM返回None的情况，增加手动计算逻辑，直接校验必要原子是否存在并计算角度
• 保留所有残基的记录，即使Phi或Psi角无法计算（会显示为None）

针对1j73.pdb的情况：
第8位的二氨基丁酸只要包含N、CA、C原子，就能计算其Phi（依赖第7位残基的C）和Psi（依赖第9位残基的N）；第7位的Psi会用第8位的N计算，第9位的Phi会用第8位的C计算，这三个残基的角度都会出现在输出文件里。

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