最近要做一個(gè)藥物分子屬性預(yù)測(cè)的課題，在跑別人現(xiàn)成的模型時(shí)，出現(xiàn)了花兩天時(shí)間都解決不了的Bug。這讓我開始反思，無腦套用網(wǎng)上的模型真的好嗎？之前對(duì)“一知半解”嗤之以鼻，覺得自己怎么樣都不會(huì)成為那個(gè)對(duì)知識(shí)對(duì)學(xué)問敷衍的人。可是為了趕進(jìn)度，自己慢慢的也變成了一個(gè)知其然而不知其所以然的人了。

無意中讀到蔡元培先生的北大就職演說里上說的話：平時(shí)則放蕩冶游，考試則熟讀講義，不問學(xué)問之有無，惟爭(zhēng)分?jǐn)?shù)之多寡；試驗(yàn)既終，書籍束之高閣，毫不過問，敷衍三四年，潦草塞責(zé)，文憑到手，即可借此活動(dòng)于社會(huì)，豈非與求學(xué)初衷大相背馳乎？

先生的話簡(jiǎn)直是古往今來都受用，這個(gè)演說拿到現(xiàn)在來說更是醍醐灌頂。于是乎我決定好好的從現(xiàn)成的代碼開始學(xué)習(xí)，追根溯源，把整個(gè)過程走通一遍，不必?zé)o腦解決BUG和跑通代碼收獲得多得多。使用深度學(xué)習(xí)處理分子數(shù)據(jù)最重要的一步就是將分子轉(zhuǎn)換成機(jī)器學(xué)習(xí)算法可以處理的數(shù)據(jù)格式。因此這里先來學(xué)習(xí)一下怎么將SMILE生成分子圖。

注：這篇博客不是教程，僅是我對(duì)代碼里一下不明白的地方做的批注和擴(kuò)展。

---

（1）化學(xué)特征

要對(duì)藥物分子進(jìn)行化學(xué)特征分析和計(jì)算時(shí)，需要先導(dǎo)入一個(gè)特征庫(kù)，創(chuàng)建一個(gè)特征工廠，并通過特征工廠計(jì)算化學(xué)特征。

• 獲取特征庫(kù)：RDConfig.RDDataDir目錄下的’BaseFeatures.fdef’
• 構(gòu)建特征工廠：ChemicalFeatures.BuildFeatureFactory(fdefName)
• fdef_name：特征庫(kù)文件
• 使用特征工廠搜索特征：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)

搜索到的每個(gè)特征都包含了該特征家族（例如供體、受體等）、特征類別、該特征對(duì)應(yīng)的原子、特征對(duì)應(yīng)序號(hào)等信息。

• 特征家族信息：GetFamily()
• 特征類型信息：GetType()
• 特征對(duì)應(yīng)原子：GetAtomIds()
• 特征對(duì)應(yīng)序號(hào)：GetId()
• 如果分子包含坐標(biāo)信息，化學(xué)特征也會(huì)包括原子坐標(biāo)：GetPos()

對(duì)于有些原子的特征，其在原子層面上無法判斷化學(xué)特征，因此需要在分子層面上進(jìn)行提取。例如化學(xué)中的donor（供體）和acceptor（受體）概念。

mol_conformers = mol.GetConformers() assert len(mol_conformers) == 1 # 獲取所有原子坐標(biāo), 該方法返回一個(gè)numpy數(shù)組 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化學(xué)包RDkit的教程_Glimmer_r的博客-CSDN博客_python rdkit

# 獲取所有原子數(shù)目 num_atoms = mol.GetNumAtoms() for u in range(num_atoms):# 訪問單個(gè)原子atom = mol.GetAtomWithIdx(u)# 獲取它的標(biāo)簽、價(jià)電子、原子元素周期編號(hào)symbol = atom.GetSymbol()atom_type = atom.GetAtomicNum()aromatic = atom.GetIsAromatic()# 獲取雜化類型hybridization = atom.GetHybridization()# 獲取與該原子連接的氫原子個(gè)數(shù)num_h = atom.GetTotalNumHs()

（3）化學(xué)鍵

通過GetBondBetweenAtoms()可以返回一個(gè)化學(xué)鍵的對(duì)象rdkit.Chem.rdchem.Bond，這個(gè)方法的輸入?yún)?shù)是兩個(gè)原子在數(shù)組中的編號(hào)，如果化學(xué)鍵不存在，則返回None

# 獲取所有原子數(shù)目 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:# 獲得化學(xué)鍵的類型bond_type = e_uv.GetBondType()

（4）將分子SMILES生成DGLGraph

摘自博文：(6條消息) 將分子SMILES生成DGLGraph_wufeil的博客-CSDN博客_smiles轉(zhuǎn)graph

1. 導(dǎo)入相關(guān)的包

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轉(zhuǎn)化為RDKIT的mol對(duì)象，同時(shí)生成一個(gè)空的dgl圖

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

3. 為DGL圖添加節(jié)點(diǎn)

G.add_nodes(molecule.GetNumAtoms())

4. 提取原子的特征和化學(xué)鍵特征

獲取原子的特征，特征包括：元素種類、隱含價(jià)、價(jià)電子、成鍵、電荷、雜化類型

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)

獲取邊特征，包括：是否為單鍵、雙鍵、三鍵、成環(huán)、芳香環(huán)、共軛

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. 提取每一個(gè)原子的特征（節(jié)點(diǎn)特征和邊特征）生成圖

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添加原子/節(jié)點(diǎn)特征 G.edata['w'] = torch.from_numpy(np.array(edge_features)) #dgl添加鍵/邊特征

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。