rdkit | 药物分子进行片段分解

rdkit | 药物分子进行片段分解

化学信息学家在某些药物研发应用场景里，不仅仅关注整体药物分子，有时也需要将所谓的类药性好的片段单独提取出来，提取出类药性质分子片段的共性，用于数据库的建设或者AI训练。
例如阿司匹林，可以被分解为苯环，羧基、乙醛以及一个单独的氧原子，共四个常见类药片段的组合。

接下来的代码采用了rdkit的BRICS算法，BRICS基于常见的反应，选择片段断键的位点，提供了化学合成意义上的可行性。

方案一

目前rdkit 有更简洁的方案一,更新如下，相比于方案二更加简洁：

from rdkit.Chem import BRICS
aspirin= Chem.MolFromSmiles('CC(=O)OC1=CC=CC=C1C(O)=O')
fragments=BRICS.BRICSDecompose(aspirin,allNodes=None, minFragmentSize=1, 
onlyUseReactions=None, silent=True, keepNonLeafNodes=False, singlePass=False, returnMols=False)
print (sorted['fragments'])
output: ['[1*]C(C)=O', '[16*]c1ccccc1[16*]', '[3*]O[3*]', '[6*]C(=O)O']

Arguments说明

• allNodes 需要指明需要包含的节点分子，相对复杂，一般不会使用；
• minFragmentSize，指明最小的片段至少需要包含的重原子个数，该例中定义为2时，'[3*]O[3*]'这个醚片段不会被拆分，而是与 ‘[16*]c1ccccc1[16*]‘的苯环合并成为’[3*]Oc1ccccc1[16*]’；
• onlyUseReactions, BRICS是基于反应的方式决定拆分位点，这里可以定义仅用什么反应拆分，使用情况较少；
• silent，若不关闭，会打印具体使用了什么反应进行拆分的相关信息；
• keepNonLeafNodes，设定为True时，会一同返回未被完全拆分的中间较大片段；
• singlePass，设定为True时返回片段仅最多包括一个断裂位点的结果，例如 ‘[16*]c1ccccc1[16*]‘结果会变为’[16*]c1ccccc1C(=O)O’与’[3*]OC©=O’，避开同一片段被多次反应断键；
• returnMols，设定为True时返回的片段不是SMILES的形式，而是rdkit.Mol的形式。

方案二

相对方案一更加复杂，但可以了解并操作更多细节，

from rdkit import Chem
from rdkit.Chem import BRICS

def fragment_recursive(mol, frags):
    try:
        bonds = list(BRICS.FindBRICSBonds(mol))
        if len(bonds) == 0:
            frags.append(mol_to_smiles(mol))
            return frags
        idxs, labs = list(zip(*bonds))
        bond_idxs = []
        for a1, a2 in idxs:
            bond = mol.GetBondBetweenAtoms(a1, a2)
            bond_idxs.append(bond.GetIdx())
        order = np.argsort(bond_idxs).tolist()
        bond_idxs = [bond_idxs[i] for i in order]
        broken = Chem.FragmentOnBonds(mol,
                                      bondIndices=[bond_idxs[0]],
                                      dummyLabels=[(0, 0)])
        head, tail = Chem.GetMolFrags(broken, asMols=True)
        #print(mol_to_smiles(head), mol_to_smiles(tail))
        frags.append(mol_to_smiles(head))
        return fragment_recursive(tail, frags)
    except Exception as e:
        print (e)
        pass

aspirin= Chem.MolFromSmiles('CC(=O)OC1=CC=CC=C1C(O)=O')
fragments=fragment_recursive(aspirin, [])
print (fragments)

# > output: ['*C(C)=O', '*O*', '*c1ccccc1*', '*C(=O)O']

可以看到，输出的片段保留了阿司匹林被切断的位点，用通配原子符号*表示，可视化出来的效果是。