学习RDKit

最近要做一个药物分子属性预测的课题，在跑别人现成的模型时，出现了花两天时间都解决不了的Bug。这让我开始反思，无脑套用网上的模型真的好吗？之前对“一知半解”嗤之以鼻，觉得自己怎么样都不会成为那个对知识对学问敷衍的人。可是为了赶进度，自己慢慢的也变成了一个知其然而不知其所以然的人了。

无意中读到蔡元培先生的北大就职演说里上说的话：平时则放荡冶游，考试则熟读讲义，不问学问之有无，惟争分数之多寡；试验既终，书籍束之高阁，毫不过问，敷衍三四年，潦草塞责，文凭到手，即可借此活动于社会，岂非与求学初衷大相背驰乎？

先生的话简直是古往今来都受用，这个演说拿到现在来说更是醍醐灌顶。于是乎我决定好好的从现成的代码开始学习，追根溯源，把整个过程走通一遍，不必无脑解决BUG和跑通代码收获得多得多。使用深度学习处理分子数据最重要的一步就是将分子转换成机器学习算法可以处理的数据格式。因此这里先来学习一下怎么将SMILE生成分子图。

注：这篇博客不是教程，仅是我对代码里一下不明白的地方做的批注和扩展。

---

（1）化学特征

要对药物分子进行化学特征分析和计算时，需要先导入一个特征库，创建一个特征工厂，并通过特征工厂计算化学特征。

• 获取特征库：RDConfig.RDDataDir目录下的’BaseFeatures.fdef’
• 构建特征工厂：ChemicalFeatures.BuildFeatureFactory(fdefName)
• fdef_name：特征库文件
• 使用特征工厂搜索特征：GetFeaturesForMol(m)

from rdkit import Chem
from rdkit.Chem import ChemicalFeatures
from rdkit import RDConfig
import os

fdef_name = osp.join(RDConfig.RDDataDir, 'BaseFeatures.fdef')
mol_featurizer = ChemicalFeatures.BuildFeatureFactory(fdef_name)
mol_feats = mol_featurizer.GetFeaturesForMol(mol)

搜索到的每个特征都包含了该特征家族（例如供体、受体等）、特征类别、该特征对应的原子、特征对应序号等信息。

• 特征家族信息：GetFamily()
• 特征类型信息：GetType()
• 特征对应原子：GetAtomIds()
• 特征对应序号：GetId()
• 如果分子包含坐标信息，化学特征也会包括原子坐标：GetPos()

对于有些原子的特征，其在原子层面上无法判断化学特征，因此需要在分子层面上进行提取。例如化学中的donor（供体）和acceptor（受体）概念。

mol_conformers = mol.GetConformers()
assert len(mol_conformers) == 1
# 获取所有原子坐标, 该方法返回一个numpy数组
geom = mol_conformers[0].GetPositions()

for i in range(len(mol_feats)):
    # 特征家族信息,电子供体与受体
    if mol_feats[i].GetFamily() == 'Donor':
        node_list = mol_feats[i].GetAtomIds()
        for u in node_list:
            is_donor[u] = 1
    elif mol_feats[i].GetFamily() == 'Acceptor':
        node_list = mol_feats[i].GetAtomIds()
        for u in node_list:
            is_acceptor[u] = 1

（2）原子属性

参考：(6条消息) 【RDKit】Python化学包RDkit的教程_Glimmer_r的博客-CSDN博客_python rdkit

# 获取所有原子数目
num_atoms = mol.GetNumAtoms()
for u in range(num_atoms):
    # 访问单个原子
    atom = mol.GetAtomWithIdx(u)
    
    # 获取它的标签、价电子、原子元素周期编号
    symbol = atom.GetSymbol()
    atom_type = atom.GetAtomicNum()
    aromatic = atom.GetIsAromatic()
    
    # 获取杂化类型
    hybridization = atom.GetHybridization()
    
    # 获取与该原子连接的氢原子个数
    num_h = atom.GetTotalNumHs()
            

（3）化学键

通过GetBondBetweenAtoms()可以返回一个化学键的对象rdkit.Chem.rdchem.Bond，这个方法的输入参数是两个原子在数组中的编号，如果化学键不存在，则返回None

# 获取所有原子数目
num_atoms = mol.GetNumAtoms()
for u in range(num_atoms):
    for v in range(num_atoms):
        if u == v and not self_loop:
            continue

        e_uv = mol.GetBondBetweenAtoms(u, v)
        if e_uv is None:
            bond_type = None
        else:
            # 获得化学键的类型
            bond_type = e_uv.GetBondType()

（4）将分子SMILES生成DGLGraph

摘自博文：(6条消息) 将分子SMILES生成DGLGraph_wufeil的博客-CSDN博客_smiles转graph

1. 导入相关的包

import numpy as np
import torch
import dgl
from dgl import DGLGraph
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors as rdDesc

2. 将SMILES转化为RDKIT的mol对象，同时生成一个空的dgl图

molecule_smiles='[C@@H](Cl)(F)Br'
G = DGLGraph()
#加载smile生成mol对象
molecule = Chem.MolFromSmiles(molecule_smiles)

3. 为DGL图添加节点

G.add_nodes(molecule.GetNumAtoms())

4. 提取原子的特征和化学键特征

获取原子的特征，特征包括：元素种类、隐含价、价电子、成键、电荷、杂化类型

def get_atom_features(atom):
    possible_atom = ['C', 'N', 'O', 'F', 'P', 'Cl', 'Br', 'I', 'DU'] #DU代表其他原子
    atom_features = one_of_k_encoding_unk(atom.GetSymbol(), possible_atom)
    atom_features += one_of_k_encoding_unk(atom.GetImplicitValence(), [0, 1])
    atom_features += one_of_k_encoding_unk(atom.GetNumRadicalElectrons(), [0, 1])
    atom_features += one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4, 5, 6])
    atom_features += one_of_k_encoding_unk(atom.GetFormalCharge(), [-1, 1])
    atom_features += one_of_k_encoding_unk(atom.GetHybridization(), 
                                           [Chem.rdchem.HybridizationType.SP, 
                                            Chem.rdchem.HybridizationType.SP2,
                                            Chem.rdchem.HybridizationType.SP3, 
                                            Chem.rdchem.HybridizationType.SP3D])       
    return np.array(atom_features) 

获取边特征，包括：是否为单键、双键、三键、成环、芳香环、共轭

def get_bond_features(bond):
    bond_type = bond.GetBondType()
    bond_feats = [
        bond_type == Chem.rdchem.BondType.SINGLE, bond_type == Chem.rdchem.BondType.DOUBLE,
        bond_type == Chem.rdchem.BondType.TRIPLE, bond_type == Chem.rdchem.BondType.AROMATIC,
        bond.GetIsConjugated(),
        bond.IsInRing()
    ]
    return np.array(bond_feats)

5. 提取每一个原子的特征（节点特征和边特征）生成图

node_features = []
edge_features = []

for i in range(molecule.GetNumAtoms()):
    atom_i = molecule.GetAtomWithIdx(i) 
    atom_i_features = get_atom_features(atom_i) 
    node_features.append(atom_i_features)
    
    for j in range(molecule.GetNumAtoms()):
        bond_ij = molecule.GetBondBetweenAtoms(i, j)
        if bond_ij is not None:
            G.add_edges(i,j) 
            bond_features_ij = get_bond_features(bond_ij) 
            edge_features.append(bond_features_ij)

G.ndata['x'] = torch.from_numpy(np.array(node_features)) #dgl添加原子/节点特征
G.edata['w'] = torch.from_numpy(np.array(edge_features)) #dgl添加键/边特征