RdKit环境中多肽结构优化报错及前药构建技术问询

一、修复多肽"Bad Conformer Id"优化错误 #

Bad Conformer Id报错的核心原因是长多肽（16肽）用默认参数构象嵌入失败，导致MMFF找不到有效构象进行优化。以下是稳定的优化方案：

import os
from rdkit import Chem
from rdkit.Chem import AllChem

aa_seq = "GDYSHASPLRYPEGGG"
pep = Chem.MolFromSequence(aa_seq)
pep = Chem.AddHs(pep)
drug_smiles = "CCCC1SC[C@]2(C)NC(=O)N[C@]12C"
drug = Chem.MolFromSmiles(drug_smiles)
drug = Chem.AddHs(drug)

# 药物构象优化（保留原逻辑）
AllChem.EmbedMolecule(drug, randomSeed=42)
AllChem.MMFFOptimizeMolecule(drug)

# 多肽构象生成与优化修复
# 使用ETKDGv3参数适配多肽构象生成，开启多线程加速
params = AllChem.ETKDGv3()
params.randomSeed = 11
params.numThreads = 4

# 先尝试单次嵌入，失败则生成多个构象筛选
embed_status = AllChem.EmbedMolecule(pep, params)
if embed_status == -1:
    conf_ids = AllChem.EmbedMultipleConfs(pep, numConfs=10, params=params)
    if not conf_ids:
        raise ValueError("多肽构象生成失败，请检查序列合法性")
    # 选取第一个生成的构象作为基础
    pep.SetConformer(pep.GetConformer(conf_ids[0]))

# 先用UFF力场做初步优化（对大分子兼容性更好）
AllChem.UFFOptimizeMolecule(pep)
# 再尝试MMFF优化，兼容则执行，否则保留UFF结果
mmff_props = AllChem.MMFFGetMoleculeProperties(pep)
if mmff_props:
    AllChem.MMFFOptimizeMolecule(pep, mmffProps=mmff_props)
else:
    print("MMFF力场不支持当前多肽结构，已使用UFF完成优化")

二、药物与多肽连接构建前药并导出PDB #

完成优化后，按需求连接分子并导出PDB：

# 1. 定位目标原子（需提前确认原子索引正确性，可通过遍历打印核对）
drug_c0 = drug.GetAtomWithIdx(0)
pep_n0 = pep.GetAtomWithIdx(0)

# 2. 移除指定氢原子
drug.RemoveAtom(14)  # 移除药物H14
pep.RemoveAtom(118)  # 移除多肽H118

# 3. 合并分子并创建共价键
prodrug = Chem.CombineMols(pep, drug)
pep_atom_count = pep.GetNumAtoms()
# 合并后药物原子索引偏移多肽原子总数
prodrug_c0_idx = pep_atom_count + 0
prodrug_n0_idx = 0
prodrug.AddBond(prodrug_n0_idx, prodrug_c0_idx, Chem.BondType.SINGLE)

# 4. 前药构象优化
AllChem.UFFOptimizeMolecule(prodrug)
mmff_props_prodrug = AllChem.MMFFGetMoleculeProperties(prodrug)
if mmff_props_prodrug:
    AllChem.MMFFOptimizeMolecule(prodrug, mmffProps=mmff_props_prodrug)

# 5. 导出PDB文件
Chem.MolToPDBFile(prodrug, "prodrug.pdb")

辅助核对原子索引的代码 #

若不确定原子索引，可执行以下代码遍历打印原子信息：

# 打印药物原子信息
for atom in drug.GetAtoms():
    print(f"索引: {atom.GetIdx()}, 元素: {atom.GetSymbol()}, 关联氢数: {atom.GetTotalNumHs()}")

# 打印多肽原子信息
for atom in pep.GetAtoms():
    print(f"索引: {atom.GetIdx()}, 元素: {atom.GetSymbol()}, 关联氢数: {atom.GetTotalNumHs()}")

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