目录
- 一、简介
- 二、读取分子
- 1.读取SMILES/SMARTS
- 2.读取.sdf
- 3.读取.mol
- 4.读取.mol2
- 5.读取其他格式
- 三、输出分子
- 1.输出SMILES/SMARTS
- 2.输出.sdf
- 3.输出.mol
- 4.读取其他格式
- 四、分子可视化
- 1.单个展示
- 2.批量展示
- 3.3D展示
一、简介
让计算机识别分子结构是计算化学码农的必备技能，也是对分子进行后续操作的基础。本文整理和总结了rdkit进行读取、输出和可视化的一些方法，包含对SMILES、SDF、MOL、MOL2、CSV等文件的处理，以及分子的结构展示。
二、读取分子
2.1.读SMILES/SMARTS
2.1.1.直接读字符串
-
从SMILES/SMARTS直接读取
很简单了，不必多说
>>> from rdkit import Chem
>>> m = Chem.MolFromSmiles('C[C@H](O)c1ccccc1')
>>> m = Chem.MolFromSmarts('Cc1ccccc1')
2.1.2.文件批量读取
文件类似这样（格式化的就行）：
| SMILES | Name |
|---|---|
| C1=CC=CC=CC=C1 | 0 |
| c1ccccc1 | 1 |
| c1cocc1 | 2 |
-
从.smi批量读取：SmilesMolSupplier(data, delimiter, smilesColumn, nameColumn, titleLine, sanitize)
data：数据文件
delimiter：分隔符，默认为' '
smilesColumn：SMILES所在列，默认为0
nameColumn：SMILES名称所在列，默认为1
titleLine：是否含有标题行，默认True
sanitize：是否检查正确性，默认True
>>> suppl = Chem.SmilesMolSupplier('data/batch_smiles.smi', delimiter='\t')
>>> mols = [Chem.MolToSmiles(mol) for mol in suppl]
>>> print(mols)
['C1=CC=CC=CC=C1', 'c1ccccc1', 'c1ccoc1']
2.1.3.文本批量读取
-
从文本批量读取：SmilesMolSupplierFromText()
参数基本同上
>>> with open('data/batch_smiles.smi', 'r') as f:
>>> mols_text = f.read()
>>> suppl = Chem.SmilesMolSupplierFromText(mols_text, delimiter='\t')
>>> mols = [Chem.MolToSmiles(mol) for mol in suppl]
>>> print(mols)
['C1=CC=CC=CC=C1', 'c1ccccc1', 'c1ccoc1']
2.1.4.DataFrame批量读取
-
读取DataFrame中的SMILES：AddMoleculeColumnToFrame(frame, smilesCol, molCol, includeFingerprints)
frame：DataFrame对象
smilesCol：SMILES所在列
molCol：新列名，将存放产生的rdkit mol对象
includeFingerprints：是否生成指纹 - 顺便计算下分子量：Descriptors.MolWt()
>>> import pandas as pd
>>> from rdkit.Chem import Descriptors
>>> from rdkit.Chem import PandasTools
>>> df = pd.read_csv('data/smiles_df.csv')
>>> PandasTools.AddMoleculeColumnToFrame(df,'SMILES','mol',includeFingerprints=True)
>>> df['MW'] = df['mol'].apply(Descriptors.MolWt)
>>> df.head(2)
Name SMILES mol MW
0 Lanreotide c1(c2c(cccc2)[nH... <img data-content="rd... 1096.347
1 Lansoprazole Cc1c(OCC(F)(F)F)... <img data-content="rd... 369.368
2.2.读.sdf
2.2.1.文件批量读取
-
从.sdf里批量读取：SDMolSupplier(fileName, sanitize, removeHs, strictParsing)
fileName：文件名
sanitize：检查化合价，计算芳香性、共轭、杂化、kekule，默认True
removeHs：是否隐藏氢原子，默认True
strictParsing：是否使用严格模式进行解析，默认True
>>> suppl = Chem.SDMolSupplier('data/batch.sdf')
>>> mols = [Chem.MolToSmiles(mol) for mol in suppl if mol]
>>> print(mols)
['C1=C\\C=C/C=C\\C=C/1', 'c1ccccc1', 'c1ccoc1']
2.2.2.压缩包批量读取
- 从file object/.gz里读取
>>> import gzip
>>> gz_file = gzip.open('data/batch.sdf.gz', 'r')
>>> suppl = Chem.ForwardSDMolSupplier(gz_file)
>>> mols = [Chem.MolToSmiles(mol) for mol in suppl if mol]
>>> print(mols)
>>> f.close()
['C1=C\\C=C/C=C\\C=C/1', 'c1ccccc1', 'c1ccoc1']
2.3.读.mol
-
从.mol里读取：MolFromMolFile(fileName, sanitize, removeHs, strictParsing)
参数同上
>>> m = Chem.MolFromMolFile('data/output.mol')
>>> print(Chem.MolToSmiles(mol))
c1cocc1
2.4.读.mol2
-
不推荐，容易出bug：MolFromMol2File(...)
参数同上
>>> m = Chem.MolFromMol2File('data/output.mol2')
>>> print(Chem.MolToSmiles(mol))
c1cocc1
2.5.读其他格式：pdb, fasta, peptide, ...
- 其他格式大同小异，不再赘述了，方法如下，感兴趣可自己尝试
# PDB
>>> Chem.MolFromPDBFile()
>>> Chem.MolFromPDBBlock()
# FASTA
>>> Chem.MolFromFASTA()
# peptide
>>> Chem.MolFromSequence()
三、输出分子
3.1.输出SMILES/SMARTS
3.1.1.输出默认式
-
输出SMILES：MolToSmiles(mol, isomericSmiles, kekuleSmiles, canonical, ...)
kekuleSmiles：默认False，不使用kekule时：脂肪族碳用"C"表示（大写），芳香族用"c"表示（小写）
isomericSmiles：默认True，区分同分异构体（"@"表示手性，"\"和"/"表示顺反异构）
canonical：默认True，输出标准SMILES
>>> m1 = Chem.MolFromSmiles('C1=CC=CC=CC=C1')
>>> m2 = Chem.MolFromSmiles('C1=CC=CC=C1')
>>> m3 = Chem.MolFromSmiles('C1=COC=C1')
>>> mols = [m1, m2, m3]
>>> print([Chem.MolToSmiles(mol) for mol in mols])
['C1=CC=CC=CC=C1', 'c1ccccc1', 'c1ccoc1']
3.1.2.输出kekule式
-
输出kekule形式
kekule形式：在符合4N+2规则的芳香体系中，通过使用双键代替小写的碳原子来表示芳香性
4N+2规则：也叫Hueckel规则，在闭环共轭体系中，当π电子数为4n+2时，才具有芳香性
>>> for mol in mols:
>>> Chem.Kekulize(mol)
>>> print([Chem.MolToSmiles(mol, kekuleSmiles=True) for mol in mols])
['C1=CC=CC=CC=C1', 'C1=CC=CC=C1', 'C1=COC=C1']
- 注：m1有共轭结构，但不属于芳香系统。m3中氧提供了2个π电子，碳各提供1个，总数为6，属于芳香系统
3.1.3.设置立体参数
-
不区分同分异构体
通过isomericSmiles控制
>>> m4 = Chem.MolFromSmiles('C[C@H](O)c1ccccc1')
>>> print(Chem.MolToSmiles(m4))
C[C@H](O)c1ccccc1
>>> print(Chem.MolToSmiles(m4, isomericSmiles=False))
CC(O)c1ccccc1
3.1.4.批量输出SMILES
-
批量输出SMILES：SmilesWriter(fileName, delimiter, includeHeader, nameHeader, isomericSmiles, kekuleSmiles)
fileName：输出文件名
delimiter：分隔符，默认为空格' '
includeHeader：是否写入表头，默认True
nameHeader：分子名一列的列名，默认'Name'
isomericSmiles：立体信息，默认True
kekuleSmiles：kekule形式，默认False
>>> writer = Chem.SmilesWriter('data/batch.smi', delimiter='\t')
>>> for i, mol in enumerate(mols):
>>> writer.write(mol)
>>> writer.close()
- 输出结果：就是2.1.2.表格中的样子
3.1.5.批量输出SMILES和属性
批量输出SMILES及属性，通过以下函数进行操作：
mol.GetPropNames()，查看分子属性列表
mol.GetProp()，获取相应属性
mol.SetProp(key, val)，新增属性名key、对应属性值val
writer.SetProps()，设置哪些属性要输出以输出分子量和LogP为例
使用Descriptors计算属性，并添加
>>> writer = Chem.SmilesWriter('data/batch_smiles.smi', delimiter='\t', nameHeader='mol_id')
>>> writer.SetProps(['LOGP', 'MW'])
>>> for i, mol in enumerate(mols):
>>> mw = Descriptors.ExactMolWt(mol)
>>> logp = Descriptors.MolLogP(mol)
>>> mol.SetProp('MW', '%.2f' %(mw))
>>> mol.SetProp('LOGP', '%.2f' %(logp))
>>> mol.SetProp('_Name', 'No_%s' %(i))
>>> writer.write(mol)
>>> writer.close()
>>> print('number of mols:', writer.NumMols())
number of mols: 3
>>> print('mol properties:', [i for i in mol.GetPropNames()])
mol properties: ['MW', 'LOGP']
- 输出结果：在2.1.2.表格中，多了“MW”和“LOGP“两列，不在这里展示了，想要代码和源文件的可以看这里。
3.1.6.输出SMARTS
-
输出SMARTS：MolToSmarts()
这个也不多说了
>>> Chem.MolToSmarts(m3, isomericSmiles=True)
'[#6]1:[#6]:[#8]:[#6]:[#6]:1'
3.2.输出.sdf
3.2.1.批量输出到.sdf
-
批量输出到文件：SDWriter()
使用方法类似于SMILES的批量输出
可以像3.1.5.一样自定义属性信息，并记录在.sdf文件中
>>> writer = Chem.SDWriter('data/batch.sdf')
>>> writer.SetProps(['LOGP', 'MW'])
>>> for i, mol in enumerate(mols):
>>> mw = Descriptors.ExactMolWt(mol)
>>> logp = Descriptors.MolLogP(mol)
>>> mol.SetProp('MW', '%.2f' %(mw))
>>> mol.SetProp('LOGP', '%.2f' %(logp))
>>> mol.SetProp('_Name', 'No_%s' %(i))
>>> writer.write(mol)
>>> writer.close()
- 输出结果：比较长，不展示了，需要的可以自取。
3.2.2.批量输出到.gz
- 批量输出到.gz
>>> outf = gzip.open('data/batch.sdf.gz','wt+')
>>> writer = Chem.SDWriter(outf)
>>> for mol in mols:
>>> writer.write(mol)
>>> writer.close()
>>> outf.close()
3.3.输出.mol
3.3.1输出连接表
- 直接输出：MolToMolBlock()
>>> m1 = Chem.MolFromSmiles('C1CCC1')
>>> print(Chem.MolToMolBlock(m1))
RDKit 2D
4 4 0 0 0 0 0 0 0 0999 V2000
1.0607 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
...
4 1 1 0
M END
3.3.2.输出到.mol
-
输出到文件：MolToMolFile(mol, filename, includeStereo, ...)
mol：mol对象
filename：文件名
includeStereo：立体信息，默认True
>>> m1.SetProp('_Name', 'cyclobutane')
>>> Chem.MolToMolFile(m1, 'data/output.mol')
3.4.输出其他格式:pdb, fasta, xyz...
- 其他格式大同小异，不再赘述了，方法如下，感兴趣可自己尝试
# PDB
>>> Chem.MolToPDBBlock()
>>> Chem.MolToPDBFile()
>>> Chem.PDBWriter()
# FASTA
>>> Chem.MolToFASTA()
# XYZ
>>> Chem.MolToXYZBlock()
>>> Chem.MolToXYZFile()
四、分子可视化
该部分只介绍方法，不贴图了，代码源文件可自行查看。
4.1.单个展示
-
从mol对象到图片：MolToImage(mol, size, kekulize, wedgeBonds, fitImage, ...)
mol：mol对象
size：图片尺寸，默认(300, 300)
kekulize：是否展示kekule形式，默认True（True：芳香系统用实线表示，False：虚线表示）
wedgeBonds：是否展示楔形键，即立体构型，默认True
>>> from rdkit.Chem import Draw
>>> mol = Chem.MolFromSmiles('C[C@H](O)c1ccccc1')
>>> Draw.MolToImage(mol, size=(150,150), kekulize=True)
-
在新窗口中展示图片：ShowMol()
参数基本同上
>>> Draw.ShowMol(mol, size=(150,150), kekulize=False)
-
保存图片MolToFile(mol, filename, size, kekulize, wedgeBonds, ...)
参数基本同上
>>> Draw.MolToFile(mol, 'data/output.png', size=(150, 150))
4.2.批量展示
4.2.1.从DataFrame中展示
-
从df中展示：FrameToGridImage(frame, column, molsPerRow, subImgSize, legendsCol, ...)
frame：DataFrame对象
column：rdkit mol对象所在列
molsPerRow,：每行显示的分子数
subImgSize：图片大小
legendsCol：标题所在列
>>> df = pd.read_csv('data/smiles_df.csv')
>>> PandasTools.AddMoleculeColumnToFrame(df,'SMILES','mol',includeFingerprints=True)
>>> PandasTools.FrameToGridImage(df, column='mol', molsPerRow=5, subImgSize=(200,200), legendsCol='Name')
4.2.2.从mol列表中展示
-
从列表生成分子结构：MolsToGridImage(mols, maxMols, molsPerRow, subImgSize, legends, ...)
部分参数和上面的一致
mols：mol对象列表
maxMols：最多显示的分子数
molsPerRow,：每行显示的分子数
subImgSize：图片大小
legends：图题
>>> mols = df.mol.tolist()
>>> legends = df.Name.tolist()
>>> Draw.MolsToGridImage(mols, maxMols=2, molsPerRow=2, subImgSize=(300,300), legends=legends)
4.3. 3D展示
- 转换3D时，为了得到靠谱的三维构象，一般先加氢：AddHs(mol)
-
通过距离几何算法计算3D坐标：EmbedMolecule(mol, randomSeed, ...)
mol：mol对象
randomSeed：随机种子 - 转换完后再进行一步力场优化，比如MMFF94：MMFFOptimizeMolecule(mol)
>>> m3d = Chem.MolFromSmiles('CNC(=O)N(N(CCCl)S(C)(=O)=O)S(C)(=O)=O')
>>> m3d = Chem.AddHs(m3d)
>>> AllChem.EmbedMolecule(m3d, randomSeed=3)
>>> AllChem.MMFFOptimizeMolecule(m3d)
>>> Draw.MolToImage(m3d, size=(250,250))
