如何在RDKit的RECAP碎片器中实现碎片与原分子的原子ID映射？

一、BRICSDecompose的原子ID映射实现 #

BRICS.BRICSDecompose直接返回碎片SMILES，无法获取原子映射，但可以结合BRICS.BreakBRICSBonds手动实现带原子ID映射的碎片筛选：

from rdkit import Chem
from rdkit.Chem import BRICS

def brics_decompose_with_mapping(mol, min_fragment_size=3):
    # 断裂BRICS键，得到带虚键的分子
    fragmented_mol = BRICS.BreakBRICSBonds(mol)
    # 获取所有碎片的原子ID列表和分子对象
    frag_atom_ids_list = Chem.GetMolFrags(fragmented_mol, asMols=False)
    frag_mols = Chem.GetMolFrags(fragmented_mol, asMols=True)
    
    # 根据最小碎片大小筛选（统计非氢原子数）
    filtered_mappings = []
    filtered_frags = []
    for atom_ids, frag_mol in zip(frag_atom_ids_list, frag_mols):
        non_h_count = sum(1 for atom in frag_mol.GetAtoms() if atom.GetAtomicNum() != 1)
        if non_h_count >= min_fragment_size:
            filtered_mappings.append(atom_ids)
            filtered_frags.append(frag_mol)
    
    return filtered_frags, filtered_mappings

# 测试示例
mol = Chem.MolFromSmiles("CCOC1=C(C=CC=C1)O[CH]([CH]2OCCNC2)C3=CC=CC=C3")
frags, mappings = brics_decompose_with_mapping(mol, min_fragment_size=4)
# 输出每个碎片对应的原分子原子ID
for idx, (frag, atom_ids) in enumerate(zip(frags, mappings)):
    print(f"碎片 {idx+1} SMILES: {Chem.MolToSmiles(frag)}")
    print(f"对应原分子原子ID: {atom_ids}\n")

二、RECAP算法的原子ID映射实现 #

RDKit的RECAP.RECAPDecompose不返回原子映射，我们可以通过RECAP.GetRECAPBonds找到断裂键，再手动拆分分子获取映射：

from rdkit import Chem
from rdkit.Chem import RECAP

def recap_decompose_with_mapping(mol, min_fragment_size=3):
    # 获取所有符合RECAP规则的断裂键
    recap_bonds = RECAP.GetRECAPBonds(mol)
    # 提取键的索引列表
    bond_indices = [bond.GetIdx() for bond in recap_bonds]
    
    # 断裂指定键，得到带虚键的分子
    fragmented_mol = Chem.BreakOnBonds(mol, bond_indices)
    # 获取所有碎片的原子ID列表和分子对象
    frag_atom_ids_list = Chem.GetMolFrags(fragmented_mol, asMols=False)
    frag_mols = Chem.GetMolFrags(fragmented_mol, asMols=True)
    
    # 根据最小碎片大小筛选
    filtered_mappings = []
    filtered_frags = []
    for atom_ids, frag_mol in zip(frag_atom_ids_list, frag_mols):
        non_h_count = sum(1 for atom in frag_mol.GetAtoms() if atom.GetAtomicNum() != 1)
        if non_h_count >= min_fragment_size:
            filtered_mappings.append(atom_ids)
            filtered_frags.append(frag_mol)
    
    return filtered_frags, filtered_mappings

# 测试示例
mol = Chem.MolFromSmiles("CCOC1=C(C=CC=C1)O[CH]([CH]2OCCNC2)C3=CC=CC=C3")
recap_frags, recap_mappings = recap_decompose_with_mapping(mol, min_fragment_size=4)
for idx, (frag, atom_ids) in enumerate(zip(recap_frags, recap_mappings)):
    print(f"RECAP碎片 {idx+1} SMILES: {Chem.MolToSmiles(frag)}")
    print(f"对应原分子原子ID: {atom_ids}\n")

关键说明 #

• 代码中min_fragment_size默认统计非氢原子数量，若需统计包含氢的总原子数，替换为len(frag_mol.GetAtoms())即可。
• Chem.GetMolFrags返回的原子ID是原分子的原始索引，确保映射准确性。
• 断裂后生成的碎片带虚键（*），若需移除虚键，可手动遍历原子删除同位素标记（RDKit中虚键原子的同位素值为1）。

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