硬核原创，转载请注明出处：
https://leytton.blog.csdn.net/article/details/130406553
一、前言
在用Pytorch图神经网络对化学分子进行数据分析的时候，经常使用现有的数据集。看到自动处理完毕的数据结构，里面的特征值让我们一脸懵逼，不知道代表的是什么含义。本文将带大家分析这些数据结构的来龙去脉。
二、数据原始特征
在使用图神经网络（GNN）对化学分子进行水溶性预测的实验中，加载了MoleculeNet的ESOL数据集。我们打开原始的csv文件，结构是这样的（非专业翻译，有误恳请留言纠正）：
| 表头 | 含义 | 示例 |
|---|---|---|
| Compound ID | 化合物ID | 2-pyrrolidone |
| ESOL predicted log solubility in mols per litre | ESOL预测对数溶解度(mol/L) | 0.243 |
| Minimum Degree | 最小度 | 1 |
| Molecular Weight | 分子量 | 85.10600000000001 |
| Number of H-Bond Donors | 氢键供体数 | 1 |
| Number of Rings | 环数 | 1 |
| Number of Rotatable Bonds | 可旋转键数 | 0 |
| Polar Surface Area | 极性表面积 | 29.1 |
| measured log solubility in mols per litre | 测量对数溶解度(mol/L) | 1.07 |
| smiles | 分子SMILES字符串 | O=C1CCCN1 |
从smiles中可以看到，这个分子有O、C、N、H（一般省略）四种元素，除去H有六个原子。
三、分析预处理数据
加载数据
如下面代码，加载ESOL数据集后将其打印出来：
from torch_geometric.datasets import MoleculeNet
data = 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 Draw
molecule = 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):
continue
if 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 Chem
mol = 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
