最近要做一個藥物分子屬性預測的課題，在跑別人現成的模型時，出現了花兩天時間都解決不了的Bug。這讓我開始反思，無腦套用網上的模型真的好嗎？之前對“一知半解”嗤之以鼻，覺得自己怎么樣都不會成為那個對知識對學問敷衍的人?？墒菫榱粟s進度，自己慢慢的也變成了一個知其然而不知其所以然的人了。

無意中讀到蔡元培先生的北大就職演說里上說的話：平時則放蕩冶游，考試則熟讀講義，不問學問之有無，惟爭分數之多寡；試驗既終，書籍束之高閣，毫不過問，敷衍三四年，潦草塞責，文憑到手，即可借此活動于社會，豈非與求學初衷大相背馳乎？

先生的話簡直是古往今來都受用，這個演說拿到現在來說更是醍醐灌頂。于是乎我決定好好的從現成的代碼開始學習，追根溯源，把整個過程走通一遍，不必無腦解決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 osfdef_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] = 1elif 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:continuee_uv = mol.GetBondBetweenAtoms(u, v)if e_uv is None:bond_type = Noneelse:# 獲得化學鍵的類型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添加鍵/邊特征

總結

以上是生活随笔為你收集整理的学习RDKit的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。