分子数据的获取、解析与结构绘制（RDKit）

在化学信息学中，分子数据的获取与解析是最基础的部分。对于已知的化合物，我们可以从各种数据库中获取；对于未知的化合物，我们也可以构造该分子。在RDKit中，无论是从数据库中获取的分子文件还是自己构建的分子，在解析前通常将其转换成mol对象，并用mol对象进行分子绘制，获取原子、键、原子序号、原子名称各种信息以及获取分子描述符等一系列操作。本文大致分为三个部分：不同来源的分子数据转成mol对象方法、分子图像的绘制与分子结构信息的解析。

一：分子数据转换mol对象

（1）SMILE编码（MolFromSmiles）

'SMILE'是Structure,Meaning,Language,Effect的缩写，它是以字符串的形式表示分子结构。并且每个smiles编码与化学结构是一一对应的。在smiles编码字符串中，主要使用C,H,O,N,S等字母表示相应的原子(不包含H)；相邻两个原子表示单键连接，原子间'='表示双键，'#'表示三键；若该原子含有支链，则用括号表示支链；若有环结构，首位相连的原子后用1，2等表示(通常芳香结构的原子用小写字母)。对于双键，“/”符号表示顺式，“\”表示反式；对于手性原子，“@”表示原子顺时针排列，“@@”为逆时针。

举例：

1:'CCC'丙烷，2:'C=C'乙烯，3:'C#C'乙炔，4:'C(=O)C'丙酮，5:'C1=CC=CC=C1'苯环，6:'C1(C=CC=C2)=C2C=CC=C1'萘。

from rdkit import Chem
s='NC(C)C(=O)O' #丙氨酸smiles编码
smile=Chem.MolFromSmiles(s) #从smiles编码中获取结构信息
Chem.Draw.ShowMol(smile) #画分子图

使用Chem.MolFromSmiles()函数，可以从smiles编码中获取分子结构信息，进而画分子图。

（2）PDB文件、FASTA文件(MolFromPDBFile、MolFromFSATA)

使用函数MolFromPDBFile获取蛋白质文件分子结构，PDB文件数据可以从NCBI、PDB等数据库获取。函数MolFromFASTA获取FASTA文件分子结构，文件可以为核酸序列或多肽序列，函数中参数flavor表示文件是哪类物质序列以及其它信息，详见（3）。

from rdkit import Chem
pdb=Chem.MolFromPDBFile('2rbg.pdb') #读取该文件夹下2rgb.pdb文件
Chem.Draw.ShowMol(pdb)

pdb文件可以与此py文件在同一文件下，此时使用相对路径；若不在同一文件，字符串为pdb文件的绝对路径。

（3）氨基酸、DNA、RNA序列(MolFromSeqence)

该函数可以利用氨基酸序列(单字符)、DNA序列、RNA序列来绘制分子结构。例如'ATTCGC'为一个DNA序列，'AUCCCGU'为一个RNA序列，'EGGYCCCDS'为氨基酸序列。

参数flavor表示该序列属于哪类物质的序列。0表示L-氨基酸组成的肽链；1表示D-氨基酸组成的肽链；2,3,4,5表示RNA，6,7,8,9表示DNA，其中区别在于核酸两端是否有cap：2,6 no cap；3,7 为5’ cap；4,8为3' cap；5,9 both cap。

from rdkit import Chem
peptide=Chem.MolFromSeqence('EGGYCCCDS',flavor=0)
Chem.Draw.ShowMol(peptide)

（4）MOL、MOL2文件(MolFromMolFile、MolFromMol2File)

这两个函数分别可以读取MOL、MOL2文件，进而画出分子结构。

除了以上方法，Chem中还有许多MolFrom系列函数，用于获取分子mol对象。

二：分子图绘制

（1）基础绘制

from rdkit import Chem
import matplotlib.pyplot as plt
s='C=CC(=O)'
mol=Chem.MolFromSmiles(s)
img=Chem.Draw.MolToImage(mol,size=(500,500))
plt.imshow(img)
plt.show()

 绘制分子图时，虽然可以使用ShowMol()来展示，但是最好将其转换成图片的形式，这样有利于图片的保存与后续操作。其中的参数size控制图片大小(准确来说是像素)，size越大，图片越清晰。

（2）多分子绘制

import matplotlib.pyplot as plt
from rdkit import Chem
s=['C=CC(=O)N1CCC(CC1)C2CCNC3=C(C(=NN23)C4=CC=C(C=C4)OC5=CC=CC=C5)C(=O)N',
    'C=CC=CC=CN(C1=CC=CC=C1)C',
    'NC(C)C(=O)O']
smiles=[]
for i in range(3):
    smiles.append(Chem.MolFromSmiles(s[i]))
img=Chem.Draw.MolsToGridImage(mols=smiles,molsPerRow=3,subImgSize=(200,200),legends=['molecule1','molecule2','molecule3'])
plt.imshow(img)
plt.show()

绘制多个分子可以使用MolFromSmiles等系列函数获取结构信息并放置在一个列表中，之后把这个列表当作整体，使用MolsToGridImage等函数绘制分子图片。

在MolsToGridImage()中，参数mols为保存mol的列表，molsPerRow表示每行的分子数，subImgSize为每个分子的图片大小，legends为分子标签，与mols列表中位置对应。

当然，使用MolsToImage()也是可以的，只是此时所有分子都在一行，与MolsToGridImage()相比没有参数molsPerRow。用这个函数绘制多行分子也是可行的，利用plt.subplot()绘制子图的方法即可。

三：分子结构信息解析

（1）mol.GetAtoms()

通过GetAtoms()可以获取该分子所有原子对象。

atoms=mol.GetAtoms()
for atom in atoms:
    print(atom)

1.atom.GetIdx()：获取原子索引号

2.atom.GetSymbol()：获取原子名。

3.atom.GetAtomicNum()：获取该原子的原子序数。

其它：GetDegree()：原子连接数；GetTotalDegree()：原子总连接数；GetFormalCharge()：原子形式电荷；GetHybirdization()：原子杂化方式；GetExplicitValence()原子隐式化合价；GetTotalValence()：原子总化合价。

（2）mol.GetBonds()

通过GetBonds()可以获取该分子所有化学键对象。

bonds=mol.GetBonds()
for bond in bonds:
    print(bond)

1.bond.GetIdx()：获取化学键索引。

2.bond.GetBondType()：获取化学键类型(键数)，如‘SINGLE'、'DOUBLE'等。

3.bond.GetBondTypeAsDouble()：与2相同，此时返回的是double浮点数，如1.0对应SINGLE，表示单键。

4.bond.GetBeginAtom()、bond.GetEndAtom()：获取化学键两端原子对象。若进一步获取原子信息，参照（1）中方法即可。

其它：GetIsAromatic()：是否为芳香键；GetIsConjugated()是否为共轭键；IsInRing()是否在环中；IsInRingSize(n)：是否在n元环中；

（3）mol.GetNumAtoms()

该方法可获取分子中原子数。

n=mol.GetNumAtoms()
print(n)

（4）mol.GetNumBonds()

该方法用于获取分子中化学键个数。

n=mol.GetNumBonds()
print(n)

（5）Chem.FindMolChiralCenters(mol)

该方法用于获取分子中手性原子与其顺反构型。

s='N[C@@](Br)(Cl)OCC[C@@](O)(N)C'
mol=Chem.MolFromSmiles(s)
p=Chem.FindMolChiralCenters(mol)
print(p)
'''-----------------------------------'''
>>>[(1, 'S'), (7, 'R')]

（6）mol.GetConformer().GetAtomPosition()

该方法用于获取分子坐标对象(如果mol包含分子坐标信息)。GetAtomPositon()中参数为原子索引数，最终返回该原子坐标对象，使用list()进行转换可以得到一个包含该原子三维坐标的列表，故也可以用索引的方式分别获取x、y、z坐标。

mol=Chem.MolFromPDBFile('2rbg.pdb')
atoms=mol.GetAtoms()
conformer=mol.GetConformer()
for atom in atoms:
    p=conformer.GetAtomPosition(atom.GetIdx())
    print(list(p))

（7）其它

mol.GetAtomWithIdx(n)：获取索引为n的原子对象；mol.GetBondWithIdx(n)：获取索引为n的键对象；mol.GetAromaticAtoms()：获取芳香环中的原子；mol.GetNumHeavyAtoms()：获取除了氢原子的其它所有原子个数；mol.GetBondBetweenAtoms()：获取两个原子之间的键。

总结：

RDKit作为化学信息学的分析工具，其功能的非常多的，具体用法可见于Python API Reference — The RDKit 2022.03.1 documentation。但实际上，从未来的发展趋势考虑，传统的机器学习方法、分子动力学模拟以及量子化学等方法很可能被深度学习所替代(如费米神经网络、AlphaFold等)，并且更多的、更为准确的分子表示方法也会逐渐取代以往的各种分子描述方法(分子描述符)。这也意味着，RDKit中的很多模块可以不必深入学习，只需要掌握其基本的分子数据读取、结构解析等部分即可。因为这些处理后的数据，可以根据我们的想法进行进一步的构建处理，并用于神经网络模型的训练。