如何将化学分子SMILES字符串转化为Pytorch图数据结构——ESOL分子水溶性数据集解析

硬核原创，转载请注明出处：
https://leytton.blog.csdn.net/article/details/130406553

一、前言

在用Pytorch图神经网络对化学分子进行数据分析的时候，经常使用现有的数据集。看到自动处理完毕的数据结构，里面的特征值让我们一脸懵逼，不知道代表的是什么含义。本文将带大家分析这些数据结构的来龙去脉。

二、数据原始特征

在使用图神经网络（GNN）对化学分子进行水溶性预测的实验中，加载了MoleculeNet的ESOL数据集。我们打开原始的csv文件，结构是这样的（非专业翻译，有误恳请留言纠正）：

从smiles中可以看到，这个分子有O、C、N、H（一般省略）四种元素，除去H有六个原子。

三、分析预处理数据

加载数据

如下面代码，加载ESOL数据集后将其打印出来：

from torch_geometric.datasets import MoleculeNetdata = MoleculeNet(root="data", name="ESOL")print("Dataset Size:", len(data))
print("Dataset classes:", data.num_classes)
print("Dataset features:", data.num_features)

Dataset Size: 1128
Dataset classes: 734
Dataset features: 9

从结果可以看到，有1128个分子样本，734种类型，每个分子有9个特征。

分析数据

我们选择第11个分子（smiles比较短）进行分析：

print("Sample:", data[10])
print("Sample y:", data[10].y)

Sample: Data(x=[6, 9], edge_index=[2, 12], edge_attr=[12, 3], smiles='O=C1CCCN1', y=[1, 1])
Sample y: tensor([[1.0700]])

可以看到x、edge_index、edge_attr 是二维数组，y可以看成一个值（水溶性）。

关于水溶性参考《理化性质|(logSw和logP)小分子化合物水溶性和脂溶性指标》

画出分子图

根据SMILES字符串，将其分子图画出来：

from rdkit import Chem
from rdkit.Chem import Drawmolecule = Chem.MolFromSmiles(data[10]["smiles"])
# Draw.MolToFile(molecule, "mol.png")
Draw.MolToImage(molecule)

edge_index数据分析

把edge_index数组打印出来：

print(data[10].edge_index.T)

tensor([[0, 1],[1, 0],[1, 2],[1, 5],[2, 1],[2, 3],[3, 2],[3, 4],[4, 3],[4, 5],[5, 1],[5, 4]])

这是O、C、N六个原子的连接关系。

x数据分析

把x数组打印出来：

print(data[10].x.shape)
print(data[10].x)

torch.Size([6, 9])
tensor([[8, 0, 1, 5, 0, 0, 3, 0, 0],[6, 0, 3, 5, 0, 0, 3, 0, 1],[6, 0, 4, 5, 2, 0, 4, 0, 1],[6, 0, 4, 5, 2, 0, 4, 0, 1],[6, 0, 4, 5, 2, 0, 4, 0, 1],[7, 0, 3, 5, 1, 0, 3, 0, 1]])

这就不太看得懂了，看起来像是描述6个原子9个特征的二维数组。

四、真相：从SMILES字符串得到

作者查阅资料无果，那么久只能去分析MoleculeNet中的代码了，到底对原始数据进行了怎样的处理，x中的数据是怎样来的。

点进去看到一个process函数应该是处理数据的，我对其进行了注释：

# Format: name: [display_name, url_name, csv_name, smiles_idx, y_idx]names = {'esol': ['ESOL', 'delaney-processed.csv', 'delaney-processed', -1, -2],'freesolv': ['FreeSolv', 'SAMPL.csv', 'SAMPL', 1, 2],'lipo': ['Lipophilicity', 'Lipophilicity.csv', 'Lipophilicity', 2, 1],'pcba': ['PCBA', 'pcba.csv.gz', 'pcba', -1,slice(0, 128)],'muv': ['MUV', 'muv.csv.gz', 'muv', -1,slice(0, 17)],'hiv': ['HIV', 'HIV.csv', 'HIV', 0, -1],'bace': ['BACE', 'bace.csv', 'bace', 0, 2],'bbbp': ['BBPB', 'BBBP.csv', 'BBBP', -1, -2],'tox21': ['Tox21', 'tox21.csv.gz', 'tox21', -1,slice(0, 12)],'toxcast':['ToxCast', 'toxcast_data.csv.gz', 'toxcast_data', 0,slice(1, 618)],'sider': ['SIDER', 'sider.csv.gz', 'sider', 0,slice(1, 28)],'clintox': ['ClinTox', 'clintox.csv.gz', 'clintox', 0,slice(1, 3)],}def process(self):with open(self.raw_paths[0], 'r') as f: #读取原始数据文件dataset = f.read().split('\n')[1:-1] #按行分割，并去掉第一行dataset = [x for x in dataset if len(x) > 0]  # 去掉空行data_list = []for line in dataset:  #遍历每行line = re.sub(r'\".*\"', '', line)  # 去掉".*"字符串line = line.split(',') #逗号分隔smiles = line[self.names[self.name][3]] #获取到smiles字符串ys = line[self.names[self.name][4]] #获取到y值ys = ys if isinstance(ys, list) else [ys] #将y值统一成数组形式ys = [float(y) if len(y) > 0 else float('NaN') for y in ys] #将y转成float类型y = torch.tensor(ys, dtype=torch.float).view(1, -1) #将y转成torch.float类型# 重点：获取x、edge_index、edge_attr数据，需要查看from_smiles函数data = from_smiles(smiles)  data.y = y  #y处理完毕if self.pre_filter is not None and not self.pre_filter(data):continueif self.pre_transform is not None:data = self.pre_transform(data)data_list.append(data)torch.save(self.collate(data_list), self.processed_paths[0])

从上面分析可以知道，原来x、edge_index、edge_attr数据都是通过将smile字符串传递到from_smiles函数获取到的！

from_smiles函数如下：

def from_smiles(smiles: str, with_hydrogen: bool = False,kekulize: bool = False) -> 'torch_geometric.data.Data':# 太多了省略。。。return Data(x=x, edge_index=edge_index, edge_attr=edge_attr, smiles=smiles)

这下可以参考这个函数的代码进一步分析了。

我们直接指定smiles进行分析：

smiles='O=C1CCCN1'
from rdkit import Chemmol = Chem.MolFromSmiles(smiles)
for atom in mol.GetAtoms():print(f'原子序号:{atom.GetAtomicNum()}, 手性信息:{atom.GetChiralTag()}, 度:{atom.GetTotalDegree()}, 电荷:{atom.GetFormalCharge()}, 连接氢原子数:{atom.GetTotalNumHs()}, 自由基:{atom.GetNumRadicalElectrons()}, 杂化类型:{atom.GetHybridization()}, 芳香性:{atom.GetIsAromatic()}, 是否在环上:{atom.IsInRing()}')

原子序号:8, 手性信息:CHI_UNSPECIFIED, 度:1, 电荷:0, 连接氢原子数:0, 自由基:0, 杂化类型:SP2, 芳香性:False, 是否在环上:False
原子序号:6, 手性信息:CHI_UNSPECIFIED, 度:3, 电荷:0, 连接氢原子数:0, 自由基:0, 杂化类型:SP2, 芳香性:False, 是否在环上:True
原子序号:6, 手性信息:CHI_UNSPECIFIED, 度:4, 电荷:0, 连接氢原子数:2, 自由基:0, 杂化类型:SP3, 芳香性:False, 是否在环上:True
原子序号:6, 手性信息:CHI_UNSPECIFIED, 度:4, 电荷:0, 连接氢原子数:2, 自由基:0, 杂化类型:SP3, 芳香性:False, 是否在环上:True
原子序号:6, 手性信息:CHI_UNSPECIFIED, 度:4, 电荷:0, 连接氢原子数:2, 自由基:0, 杂化类型:SP3, 芳香性:False, 是否在环上:True
原子序号:7, 手性信息:CHI_UNSPECIFIED, 度:3, 电荷:0, 连接氢原子数:1, 自由基:0, 杂化类型:SP2, 芳香性:False, 是否在环上:True

如上所示，这9个特征就是x变量中每个原子的含义，对其进行一些编码变换就构造成了x变量。具体的原子更多的属性，可以参考 RDKit 文档

接下来我们分析edge_attr和edge_index变量含义：

for bond in mol.GetBonds(): #便利所有的键i = bond.GetBeginAtomIdx()j = bond.GetEndAtomIdx()print(f'连接：{i,j},{j,i}')print(f'键的类型:{bond.GetBondType()}, Stereo:{bond.GetStereo()}, 是否共轭:{bond.GetIsConjugated()}')

连接：(0, 1),(1, 0)
键的类型:DOUBLE, Stereo:STEREONONE, 是否共轭:True
连接：(1, 2),(2, 1)
键的类型:SINGLE, Stereo:STEREONONE, 是否共轭:False
连接：(2, 3),(3, 2)
键的类型:SINGLE, Stereo:STEREONONE, 是否共轭:False
连接：(3, 4),(4, 3)
键的类型:SINGLE, Stereo:STEREONONE, 是否共轭:False
连接：(4, 5),(5, 4)
键的类型:SINGLE, Stereo:STEREONONE, 是否共轭:False
连接：(5, 1),(1, 5)
键的类型:SINGLE, Stereo:STEREONONE, 是否共轭:True

这就是分子SMILES字符串转化成图数据结构的过程，可以看到只用到了原始数据里的SMILES字符串和水溶性结果。

在Pytorch官网找了半天没找到数据集的说明资料，等我分析完后，才发现，这里已经有大佬发表了相关文章。不过，如果不知道图结构数据是从SMILES字符串分析得到，很难通过关键字找到这些资料。
How to turn a SMILES string into a molecular graph for Pytorch Geometric