RDKit基础操作

RDKit 是一个常用的生物化学信息python工具包。它提供了大量对化学分子2D或3D的计算操作，可生成用于机器学习的分子描述符，以及提供其他更强大的功能；

RDKit的安装可以使用Conda完成：

conda install -c rdkit rdkit

安装完成后，可以测试：

import rdkit

print(rdkit.__version__) # 2017.09.1

预读内容

分子的表示

对于机器学习算法而言，如果该特征本身是数值变量那么可以使用它本身作为输入，对于类别变量而言，最直接的方式便是通过one-hot encoding的方式进行表示，那么同样的，对于一个化合物分子，不管是大分子还是小分子，其均有相应的结构与之依附，那么对这些结构的不同表示方式，也就决定了模型的特征表示方式；

分子指纹Fingerpoint

表示药物的一种方法是分子指纹。 指纹的最普遍类型是一系列二进制数字（位），代表分子中是否存在特定的子结构。 因此，药物（小化合物）被描述为0和1的向量（数组）。如下图所示：

这种表示方式的优点是简单快速，但是，很明显，将分子编码为二进制向量不是一个可逆的过程（这是有损的转化）。 即，我们可以将一个能够表示结构信息的分子式编码成分子指纹，但是却不可以从分子指纹中推断出该分子有怎样的结构；

SMILES

表示分子的另一种方法是将结构编码为文本。 将图结构数据转换为文本内容，并在机器学习模型中使用文本（编码字符串）作为输入。 Simplified Molecular-Input Line-Entry System（SMILES）是最受欢迎的表示之一；

简化分子线性输入规范（SMILES）是一种用字符串明确描述分子结构的规范。SMILES支持被分子编辑软件导入并转换成二维或三维分子结构，其中，转换成二维图形可以使用"结构图生成算法"；SMILES对于每个结构的唯一性依赖于生成它的算法，这被称为规范SMILES。

举例关于SMILES编码的几个特点：

• 原子的化学符号表示：原子由各自的原子符号表示，省略氢原子，氢原子根据价数进行复原，比如水的表示就是O，乙醇是CCO；
• 双键用=表示，三键用#表示：含有双键的二氧化碳为O=C=O，含有三键的氰化氢为C#N；
• 碳链上的分支用圆括号表示：丙酸表示为CCC(=O)O；

InChIKey

尽管SMILES在化学家和机器学习研究人员中非常受欢迎，但它并不是唯一可用于表示药物的基于文本的表示形式；

规范SMILES存在无法自由使用的问题，因为它的生成算法是商业性的，于是提出InChI，以开发可自由使用的化合物规范表示法，InChI是以人类可以理解的形式编写的分子信息。

InChIKey 是对 InChI 运用 SHA-256 算法处理后得到的哈希值，它的出现是为了解决 InChI 长度不定的问题。与 InChI 相比，InChIKey 具有这样几个特点：

• 长度固定，永远是27个字母
• 与 InChI 几乎一一对应，但有很小的概率（1/10亿）出现两个 InChI 对应同一个InChIKey
• 不可读，字符串本身没有意义，必须转换回 InChI 才能读
• 在实际使用中，可以用InChIKey 作为关键字检索出对应的 InChI，再做进一步的使用

Graph

直接使用图结构表示分子，这是目前常用于GNN处理的数据格式，它能够有效保留分子的结构信息；

药物与靶标的相互作用预测

药物发现中的一个重要问题是预测特定药物是否可以结合特定蛋白质，即药物-靶标相互作用（DTI）预测任务；

---

靶标是指体内具有药效功能并能被药物作用的生物大分子，如某些蛋白质和核酸等生物大分子

---

我们可以像下面这样构造DTI预测任务：

描述：预测化合物与蛋白质结合亲和力的二元分类

输入：化合物和蛋白质表示向量

输出：{0,1}或[0-1]中的实数(概率)

读写分子

RDKit 支持从Smiles、mol、sdf 文件中读入分子获得分子对象。 Smiles、mol 是通常用于保存单个分子；而sdf格式当初是作为分子库形式设计的。 因此读入sdf得到的是分子迭代器，读入Smiles 和mol 文件是分子对象。

读入smiles：

smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3'
from rdkit import Chem
from rdkit.Chem import AllChem
mol = Chem.MolFromSmiles(smi)

print(mol) # <rdkit.Chem.rdchem.Mol object at 0x0000016D1CC557B0>
print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>

读入mol文件：

mol = Chem.MolFromMolFile('rd.mol')
print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>

读入sdf文件：

mols_suppl = Chem.SDMolSupplier('ZINC53.sdf')
print(type(mols_suppl)) # <class 'rdkit.Chem.rdmolfiles.SDMolSupplier'>

显示变量mols_suppl的类型：<class ‘rdkit.Chem.rdmolfiles.SDMolSupplier’>。 mols_suppl 可以看成是mol的列表，支持索引操作和迭代操作：

mol_1= mols_suppl[0]
print(type(mol_1)) # <class 'rdkit.Chem.rdchem.Mol'>
for mol in mols_suppl:
   print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>

读入多肽字符串：

seq='GGGGG'
mol = Chem.MolFromSequence(seq)
smi = Chem.MolToSmiles(mol)
print("smi",smi) # smi NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)O

RDKit 可以把分子对象保存成Smiles、molBlock、mol文件：

smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3'
mol = Chem.MolFromSmiles(smi)
smi = Chem.MolToSmiles(mol)
print(smi)
molblock = Chem.MolToMolBlock(mol)
print(molblock)
print(molblock,file=open('foo.mol','w+'))

此处注意：print(*objects, sep=’ ‘, end=’n’, file=sys.stdout, flush=False) ，python的print函数file参数支持定义输出位置；

操作分子

mol对象中有获取所有原子的方法GetAtoms()：

from rdkit import Chem
smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3'
mol = Chem.MolFromSmiles(smi)
atoms = mol.GetAtoms()
print(type(atoms)) # atoms<class 'rdkit.Chem.rdchem._ROAtomSeq'>
print(type(atoms[0])) # <class 'rdkit.Chem.rdchem.Atom'>

atoms的类型为：<class ‘rdkit.Chem.rdchem._ROAtomSeq’> 可以看成是atom的列表；
atom的类型：<class ‘rdkit.Chem.rdchem.Atom’>；

mol对象中有获取所有键的方法GetBonds()：

smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3'
mol = Chem.MolFromSmiles(smi)
bonds = mol.GetBonds()
print(type(bonds))
print(type(bonds[0]))

bonds的类型为<class ‘rdkit.Chem.rdchem._ROBondSeq’>，可以看成是bond的列表；
bond的类型为：<class ‘rdkit.Chem.rdchem.Bond’>；

根据原子编号获取键GetAtomWithIdx()：

smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3'
mol = Chem.MolFromSmiles(smi)
atom0 = mol.GetAtomWithIdx(0)

rdkit采用的是最小成环原则: 比如两个环并在一起，从数学上看可以看成是2个小环，一个大环。在rdkit中计算环的时候只考虑小环。如下图所示，我们构建了一个4元环并3元环的分子 “OC1C2C1CC2”：

我们可以看到该分子是由3元环并4元环组成的，其中原子2和原子3位于3元环和4元环的公共边上：

m = Chem.MolFromSmiles('OC1C2C1CC2')
atom2  = m.GetAtomWithIdx(2)
print("atom2 in ring:",atom2.IsInRing())
print("atom2 in 3-ring:",atom2.IsInRingSize(3))
print("atom2 in 4-ring:",atom2.IsInRingSize(4))
print("atom2 in 5-ring:",atom2.IsInRingSize(5))
"""
atom2 in ring: True
atom2 in 3-ring: True
atom2 in 4-ring: True
atom2 in 5-ring: False
"""

其中：

• ‘IsInRing()’: 判断是否在环上
• ‘IsInRingSize(n)’：判断是否在n元环上

获取分子中所有的环，以及每个环对应的原子组成，以上面的分子为例：

m = Chem.MolFromSmiles('OC1C2C1CC2')
ssr = Chem.GetSymmSSSR(m)
num_ring = len(ssr)
print("num of ring",num_ring)
for ring in ssr:
    print("ring consisted of atoms id:",list(ring))

"""
num of ring 2
ring consisted of atoms id: [1, 2, 3]
ring consisted of atoms id: [4, 5, 2, 3]
"""

通过计算，我们发现示例分子一共有2个环，第一个环由3个原子（原子1,2,3）组成， 第二个环由4个原子组成（原子4,5,2,3）；

修改分子

增加氢原子：Chem.AddHs()；删除氢原子：Chem.RemoveHs(m2)；

RDKit 中的分子默认采用隐式H原子形式。 RDKit 中提供了Chem.AddHs()方法，显式添加H原子。

from rdkit import Chem
m = Chem.MolFromSmiles('OC1C2C1CC2')
m2 = Chem.AddHs(m)
print("m Smiles:",Chem.MolToSmiles(m))
print("m2 Smiles:",Chem.MolToSmiles(m2))
print("num ATOMs:",m2.GetNumAtoms())
print("num ATOMs:",m.GetNumAtoms())
"""
m Smiles: OC1C2CCC12
m2 Smiles: [H]OC1([H])C2([H])C([H])([H])C([H])([H])C12[H]
num ATOMs: 14
num ATOMs: 6
"""

处理2D分子

Smiles 可以看成分子的1D形式，分子的平面结构可以看成分子的2D形式；

产生2D坐标AllChem.Compute2DCoords(m)；AllChem.Compute2DCoords(m)产生固定的独一无二的取向，可用于作为绘制分子的模板，如果有多个分子共享一个骨架，我们希望绘图的时候能够保证在这些分子中的骨架方向是一致的。 首先把骨架提取成绘图模板，然后在绘图上添加基团，举例子：

我们可以看到这4个结构都含有吡咯并嘧啶（c1nccc2n1ccc2）的子结构， 但是在图片上他们的子结构取向不一致，不便于比较他们的结构。 我们可以采用模板绘制法，固定公共子结构的取向：

from rdkit.Chem import Draw
from rdkit import Chem
import matplotlib.pyplot as plt

smis=[
    'COC1=C(C=CC(=C1)NS(=O)(=O)C)C2=CN=CN3C2=CC=C3',
    'C1=CC2=C(C(=C1)C3=CN=CN4C3=CC=C4)ON=C2C5=CC=C(C=C5)F',
    'COC(=O)C1=CC2=CC=CN2C=N1',
    'C1=C2C=C(N=CN2C(=C1)Cl)C(=O)O',
]
template = Chem.MolFromSmiles('c1nccc2n1ccc2')
AllChem.Compute2DCoords(template)
mols=[]
for smi in smis:
    mol = Chem.MolFromSmiles(smi)
    AllChem.GenerateDepictionMatching2DStructure(mol,template)
    mols.append(mol)
img=Draw.MolsToGridImage(mols,molsPerRow=4,subImgSize=(200,200),legends=['' for x in mols])
plt.imshow(img)
plt.show()

这样我们就可以很清楚的看出这4个分子取代基团和位置的差异。

指纹和相似性

RDKit 内置了多种分子指纹计算方法，如：

• 拓扑指纹 Chem.RDKFingerprint(mol)
• MACCS 指纹
• Atom Pairs
• topological torsions
• 摩根指纹（圆圈指纹）

拓扑指纹 Chem.RDKFingerprint(x)：

ms = [Chem.MolFromSmiles('CCOC'), Chem.MolFromSmiles('CCO'), Chem.MolFromSmiles('COC')]
fps = [Chem.RDKFingerprint(x) for x in ms]

MACCS 指纹MACCSkeys.GenMACCSKeys(mol)：

from rdkit.Chem import MACCSkeys
fps = [MACCSkeys.GenMACCSKeys(x) for x in ms]

相似性计算方法有：

• Tanimoto, 默认的方法
• Dice,
• Cosine,
• Sokal,
• Russel,
• Kulczynski,
• McConnaughey, and
• Tversky.

比较下面3个分子的相似性

分子1： CC(=O)CC(C1=CC=C(C=C1)[N+]([O-])=O)C1=C(O)C2=CC=CC=C2OC1=O
分子2： CC(=O)CC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2
分子3： CCC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2

基于MACCS 指纹和Dice 相似性方法计算相似性：

from rdkit import DataStructs
from rdkit.Chem import MACCSkeys
import rdkit
from rdkit import Chem
from rdkit.Chem import Draw
smis=[
    'CC(=O)CC(C1=CC=C(C=C1)[N+]([O-])=O)C1=C(O)C2=CC=CC=C2OC1=O',
'CC(=O)CC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2',
'CCC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2'
]
mols =[]
for smi in smis:
    m = Chem.MolFromSmiles(smi)
    mols.append(m)


fps = [MACCSkeys.GenMACCSKeys(x) for x in mols]
sm01=DataStructs.FingerprintSimilarity(fps[0],fps[1],metric=DataStructs.DiceSimilarity)

sm02=DataStructs.FingerprintSimilarity(fps[0],fps[2],metric=DataStructs.DiceSimilarity)

sm12=DataStructs.FingerprintSimilarity(fps[1],fps[2],metric=DataStructs.DiceSimilarity)
print("similarity between mol 1 and mol2: %.2f"%sm01)
print("similarity between mol 1 and mol3: %.2f"%sm02)
print("similarity between mol 2 and mol3: %.2f"%sm12)

"""
similarity between mol 1 and mol2: 0.78
similarity between mol 1 and mol3: 0.70
similarity between mol 2 and mol3: 0.92
"""

从相似性，我们可以看出分子2和分子3比较相似。