RDKit基础操作
目錄
- 預讀內(nèi)容
- 分子的表示
- 分子指紋Fingerpoint
- SMILES
- InChIKey
- Graph
- 藥物與靶標的相互作用預測
- 讀寫分子
- 操作分子
- 修改分子
- 處理2D分子
- 指紋和相似性
RDKit 是一個常用的生物化學信息python工具包。它提供了大量對化學分子2D或3D的計算操作，可生成用于機器學習的分子描述符，以及提供其他更強大的功能；
RDKit的安裝可以使用Conda完成：
conda install -c rdkit rdkit安裝完成后，可以測試：
import rdkitprint(rdkit.__version__) # 2017.09.1預讀內(nèi)容
分子的表示
對于機器學習算法而言，如果該特征本身是數(shù)值變量那么可以使用它本身作為輸入，對于類別變量而言，最直接的方式便是通過one-hot encoding的方式進行表示，那么同樣的，對于一個化合物分子，不管是大分子還是小分子，其均有相應的結構與之依附，那么對這些結構的不同表示方式，也就決定了模型的特征表示方式；
分子指紋Fingerpoint
表示藥物的一種方法是分子指紋。 指紋的最普遍類型是一系列二進制數(shù)字（位），代表分子中是否存在特定的子結構。 因此，藥物（小化合物）被描述為0和1的向量（數(shù)組）。如下圖所示：
這種表示方式的優(yōu)點是簡單快速，但是，很明顯，將分子編碼為二進制向量不是一個可逆的過程（這是有損的轉(zhuǎn)化）。 即，我們可以將一個能夠表示結構信息的分子式編碼成分子指紋，但是卻不可以從分子指紋中推斷出該分子有怎樣的結構；
SMILES
表示分子的另一種方法是將結構編碼為文本。 將圖結構數(shù)據(jù)轉(zhuǎn)換為文本內(nèi)容，并在機器學習模型中使用文本（編碼字符串）作為輸入。 Simplified Molecular-Input Line-Entry System（SMILES）是最受歡迎的表示之一；
簡化分子線性輸入規(guī)范（SMILES）是一種用字符串明確描述分子結構的規(guī)范。SMILES支持被分子編輯軟件導入并轉(zhuǎn)換成二維或三維分子結構，其中，轉(zhuǎn)換成二維圖形可以使用"結構圖生成算法"；SMILES對于每個結構的唯一性依賴于生成它的算法，這被稱為規(guī)范SMILES。
舉例關于SMILES編碼的幾個特點：
- 原子的化學符號表示：原子由各自的原子符號表示，省略氫原子，氫原子根據(jù)價數(shù)進行復原，比如水的表示就是O，乙醇是CCO；
- 雙鍵用=表示，三鍵用#表示：含有雙鍵的二氧化碳為O=C=O，含有三鍵的氰化氫為C#N；
- 碳鏈上的分支用圓括號表示：丙酸表示為CCC(=O)O；
InChIKey
盡管SMILES在化學家和機器學習研究人員中非常受歡迎，但它并不是唯一可用于表示藥物的基于文本的表示形式；
規(guī)范SMILES存在無法自由使用的問題，因為它的生成算法是商業(yè)性的，于是提出InChI，以開發(fā)可自由使用的化合物規(guī)范表示法，InChI是以人類可以理解的形式編寫的分子信息。
InChIKey 是對 InChI 運用 SHA-256 算法處理后得到的哈希值，它的出現(xiàn)是為了解決 InChI 長度不定的問題。與 InChI 相比，InChIKey 具有這樣幾個特點：
- 長度固定，永遠是27個字母
- 與 InChI 幾乎一一對應，但有很小的概率（1/10億）出現(xiàn)兩個 InChI 對應同一個InChIKey
- 不可讀，字符串本身沒有意義，必須轉(zhuǎn)換回 InChI 才能讀
- 在實際使用中，可以用InChIKey 作為關鍵字檢索出對應的 InChI，再做進一步的使用
Graph
直接使用圖結構表示分子，這是目前常用于GNN處理的數(shù)據(jù)格式，它能夠有效保留分子的結構信息；
藥物與靶標的相互作用預測
藥物發(fā)現(xiàn)中的一個重要問題是預測特定藥物是否可以結合特定蛋白質(zhì)，即藥物-靶標相互作用（DTI）預測任務；
靶標是指體內(nèi)具有藥效功能并能被藥物作用的生物大分子，如某些蛋白質(zhì)和核酸等生物大分子
我們可以像下面這樣構造DTI預測任務：
描述：預測化合物與蛋白質(zhì)結合親和力的二元分類輸入：化合物和蛋白質(zhì)表示向量輸出：{0,1}或[0-1]中的實數(shù)(概率)讀寫分子
RDKit 支持從Smiles、mol、sdf 文件中讀入分子獲得分子對象。 Smiles、mol 是通常用于保存單個分子；而sdf格式當初是作為分子庫形式設計的。 因此讀入sdf得到的是分子迭代器，讀入Smiles 和mol 文件是分子對象。
讀入smiles：
smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3' from rdkit import Chem from rdkit.Chem import AllChem mol = Chem.MolFromSmiles(smi)print(mol) # <rdkit.Chem.rdchem.Mol object at 0x0000016D1CC557B0> print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>讀入mol文件：
mol = Chem.MolFromMolFile('rd.mol') print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>讀入sdf文件：
mols_suppl = Chem.SDMolSupplier('ZINC53.sdf') print(type(mols_suppl)) # <class 'rdkit.Chem.rdmolfiles.SDMolSupplier'>顯示變量mols_suppl的類型：<class ‘rdkit.Chem.rdmolfiles.SDMolSupplier’>。 mols_suppl 可以看成是mol的列表，支持索引操作和迭代操作：
mol_1= mols_suppl[0] print(type(mol_1)) # <class 'rdkit.Chem.rdchem.Mol'> for mol in mols_suppl:print(type(mol)) # <class 'rdkit.Chem.rdchem.Mol'>讀入多肽字符串：
seq='GGGGG' mol = Chem.MolFromSequence(seq) smi = Chem.MolToSmiles(mol) print("smi",smi) # smi NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)ORDKit 可以把分子對象保存成Smiles、molBlock、mol文件：
smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3' mol = Chem.MolFromSmiles(smi) smi = Chem.MolToSmiles(mol) print(smi) molblock = Chem.MolToMolBlock(mol) print(molblock) print(molblock,file=open('foo.mol','w+'))此處注意：print(*objects, sep=’ ‘, end=’n’, file=sys.stdout, flush=False) ，python的print函數(shù)file參數(shù)支持定義輸出位置；
操作分子
mol對象中有獲取所有原子的方法GetAtoms()：
from rdkit import Chem smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3' mol = Chem.MolFromSmiles(smi) atoms = mol.GetAtoms() print(type(atoms)) # atoms<class 'rdkit.Chem.rdchem._ROAtomSeq'> print(type(atoms[0])) # <class 'rdkit.Chem.rdchem.Atom'>atoms的類型為：<class ‘rdkit.Chem.rdchem._ROAtomSeq’> 可以看成是atom的列表；
atom的類型：<class ‘rdkit.Chem.rdchem.Atom’>；
mol對象中有獲取所有鍵的方法GetBonds()：
smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3' mol = Chem.MolFromSmiles(smi) bonds = mol.GetBonds() print(type(bonds)) print(type(bonds[0]))bonds的類型為<class ‘rdkit.Chem.rdchem._ROBondSeq’>，可以看成是bond的列表；
bond的類型為：<class ‘rdkit.Chem.rdchem.Bond’>；
根據(jù)原子編號獲取鍵GetAtomWithIdx()：
smi='CC(C)OC(=O)C(C)NP(=O)(OCC1C(C(C(O1)N2C=CC(=O)NC2=O)(C)F)O)OC3=CC=CC=C3' mol = Chem.MolFromSmiles(smi) atom0 = mol.GetAtomWithIdx(0)rdkit采用的是最小成環(huán)原則: 比如兩個環(huán)并在一起，從數(shù)學上看可以看成是2個小環(huán)，一個大環(huán)。在rdkit中計算環(huán)的時候只考慮小環(huán)。如下圖所示，我們構建了一個4元環(huán)并3元環(huán)的分子 “OC1C2C1CC2”：
我們可以看到該分子是由3元環(huán)并4元環(huán)組成的，其中原子2和原子3位于3元環(huán)和4元環(huán)的公共邊上：
其中：
- ‘IsInRing()’: 判斷是否在環(huán)上
- ‘IsInRingSize(n)’：判斷是否在n元環(huán)上
獲取分子中所有的環(huán)，以及每個環(huán)對應的原子組成，以上面的分子為例：
m = Chem.MolFromSmiles('OC1C2C1CC2') ssr = Chem.GetSymmSSSR(m) num_ring = len(ssr) print("num of ring",num_ring) for ring in ssr:print("ring consisted of atoms id:",list(ring))""" num of ring 2 ring consisted of atoms id: [1, 2, 3] ring consisted of atoms id: [4, 5, 2, 3] """通過計算，我們發(fā)現(xiàn)示例分子一共有2個環(huán)，第一個環(huán)由3個原子（原子1,2,3）組成， 第二個環(huán)由4個原子組成（原子4,5,2,3）；
修改分子
增加氫原子：Chem.AddHs()；刪除氫原子：Chem.RemoveHs(m2)；
RDKit 中的分子默認采用隱式H原子形式。 RDKit 中提供了Chem.AddHs()方法，顯式添加H原子。
from rdkit import Chem m = Chem.MolFromSmiles('OC1C2C1CC2') m2 = Chem.AddHs(m) print("m Smiles:",Chem.MolToSmiles(m)) print("m2 Smiles:",Chem.MolToSmiles(m2)) print("num ATOMs:",m2.GetNumAtoms()) print("num ATOMs:",m.GetNumAtoms()) """ m Smiles: OC1C2CCC12 m2 Smiles: [H]OC1([H])C2([H])C([H])([H])C([H])([H])C12[H] num ATOMs: 14 num ATOMs: 6 """處理2D分子
Smiles 可以看成分子的1D形式，分子的平面結構可以看成分子的2D形式；
產(chǎn)生2D坐標AllChem.Compute2DCoords(m)；AllChem.Compute2DCoords(m)產(chǎn)生固定的獨一無二的取向，可用于作為繪制分子的模板，如果有多個分子共享一個骨架，我們希望繪圖的時候能夠保證在這些分子中的骨架方向是一致的。 首先把骨架提取成繪圖模板，然后在繪圖上添加基團，舉例子：
我們可以看到這4個結構都含有吡咯并嘧啶（c1nccc2n1ccc2）的子結構， 但是在圖片上他們的子結構取向不一致，不便于比較他們的結構。 我們可以采用模板繪制法，固定公共子結構的取向：
這樣我們就可以很清楚的看出這4個分子取代基團和位置的差異。
指紋和相似性
RDKit 內(nèi)置了多種分子指紋計算方法，如：
- 拓撲指紋 Chem.RDKFingerprint(mol)
- MACCS 指紋
- Atom Pairs
- topological torsions
- 摩根指紋（圓圈指紋）
拓撲指紋 Chem.RDKFingerprint(x)：
ms = [Chem.MolFromSmiles('CCOC'), Chem.MolFromSmiles('CCO'), Chem.MolFromSmiles('COC')] fps = [Chem.RDKFingerprint(x) for x in ms]MACCS 指紋MACCSkeys.GenMACCSKeys(mol)：
from rdkit.Chem import MACCSkeys fps = [MACCSkeys.GenMACCSKeys(x) for x in ms]相似性計算方法有：
- Tanimoto, 默認的方法
- Dice,
- Cosine,
- Sokal,
- Russel,
- Kulczynski,
- McConnaughey, and
- Tversky.
比較下面3個分子的相似性
分子1： CC(=O)CC(C1=CC=C(C=C1)[N+]([O-])=O)C1=C(O)C2=CC=CC=C2OC1=O 分子2： CC(=O)CC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2 分子3： CCC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2基于MACCS 指紋和Dice 相似性方法計算相似性：
from rdkit import DataStructs from rdkit.Chem import MACCSkeys import rdkit from rdkit import Chem from rdkit.Chem import Draw smis=['CC(=O)CC(C1=CC=C(C=C1)[N+]([O-])=O)C1=C(O)C2=CC=CC=C2OC1=O', 'CC(=O)CC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2', 'CCC(C1=CC=CC=C1)C1=C(O)C2=C(OC1=O)C=CC=C2' ] mols =[] for smi in smis:m = Chem.MolFromSmiles(smi)mols.append(m)fps = [MACCSkeys.GenMACCSKeys(x) for x in mols] sm01=DataStructs.FingerprintSimilarity(fps[0],fps[1],metric=DataStructs.DiceSimilarity)sm02=DataStructs.FingerprintSimilarity(fps[0],fps[2],metric=DataStructs.DiceSimilarity)sm12=DataStructs.FingerprintSimilarity(fps[1],fps[2],metric=DataStructs.DiceSimilarity) print("similarity between mol 1 and mol2: %.2f"%sm01) print("similarity between mol 1 and mol3: %.2f"%sm02) print("similarity between mol 2 and mol3: %.2f"%sm12)""" similarity between mol 1 and mol2: 0.78 similarity between mol 1 and mol3: 0.70 similarity between mol 2 and mol3: 0.92 """從相似性，我們可以看出分子2和分子3比較相似。
總結
- 上一篇： 医疗图像配准-点云配准总结
- 下一篇： 这15个网站，为设计师提供用不完的免费素