RDKit基础操作
目錄
- 預讀內容
- 分子的表示
- 分子指紋Fingerpoint
- SMILES
- InChIKey
- Graph
- 藥物與靶標的相互作用預測
- 讀寫分子
- 操作分子
- 修改分子
- 處理2D分子
- 指紋和相似性
RDKit 是一個常用的生物化學信息python工具包。它提供了大量對化學分子2D或3D的計算操作，可生成用于機器學習的分子描述符，以及提供其他更強大的功能；
RDKit的安裝可以使用Conda完成：
conda install -c rdkit rdkit安裝完成后，可以測試：
import rdkitprint(rdkit.__version__) # 2017.09.1預讀內容
分子的表示
對于機器學習算法而言，如果該特征本身是數值變量那么可以使用它本身作為輸入，對于類別變量而言，最直接的方式便是通過one-hot encoding的方式進行表示，那么同樣的，對于一個化合物分子，不管是大分子還是小分子，其均有相應的結構與之依附，那么對這些結構的不同表示方式，也就決定了模型的特征表示方式；
分子指紋Fingerpoint
表示藥物的一種方法是分子指紋。 指紋的最普遍類型是一系列二進制數字（位），代表分子中是否存在特定的子結構。 因此，藥物（小化合物）被描述為0和1的向量（數組）。如下圖所示：
這種表示方式的優點是簡單快速，但是，很明顯，將分子編碼為二進制向量不是一個可逆的過程（這是有損的轉化）。 即，我們可以將一個能夠表示結構信息的分子式編碼成分子指紋，但是卻不可以從分子指紋中推斷出該分子有怎樣的結構；
SMILES
表示分子的另一種方法是將結構編碼為文本。 將圖結構數據轉換為文本內容，并在機器學習模型中使用文本（編碼字符串）作為輸入。 Simplified Molecular-Input Line-Entry System（SMILES）是最受歡迎的表示之一；
簡化分子線性輸入規范（SMILES）是一種用字符串明確描述分子結構的規范。SMILES支持被分子編輯軟件導入并轉換成二維或三維分子結構，其中，轉換成二維圖形可以使用"結構圖生成算法"；SMILES對于每個結構的唯一性依賴于生成它的算法，這被稱為規范SMILES。
舉例關于SMILES編碼的幾個特點：
- 原子的化學符號表示：原子由各自的原子符號表示，省略氫原子，氫原子根據價數進行復原，比如水的表示就是O，乙醇是CCO；
- 雙鍵用=表示，三鍵用#表示：含有雙鍵的二氧化碳為O=C=O，含有三鍵的氰化氫為C#N；
- 碳鏈上的分支用圓括號表示：丙酸表示為CCC(=O)O；
InChIKey
盡管SMILES在化學家和機器學習研究人員中非常受歡迎，但它并不是唯一可用于表示藥物的基于文本的表示形式；
規范SMILES存在無法自由使用的問題，因為它的生成算法是商業性的，于是提出InChI，以開發可自由使用的化合物規范表示法，InChI是以人類可以理解的形式編寫的分子信息。
InChIKey 是對 InChI 運用 SHA-256 算法處理后得到的哈希值，它的出現是為了解決 InChI 長度不定的問題。與 InChI 相比，InChIKey 具有這樣幾個特點：
- 長度固定，永遠是27個字母
- 與 InChI 幾乎一一對應，但有很小的概率（1/10億）出現兩個 InChI 對應同一個InChIKey
- 不可讀，字符串本身沒有意義，必須轉換回 InChI 才能讀
- 在實際使用中，可以用InChIKey 作為關鍵字檢索出對應的 InChI，再做進一步的使用
Graph
直接使用圖結構表示分子，這是目前常用于GNN處理的數據格式，它能夠有效保留分子的結構信息；
藥物與靶標的相互作用預測
藥物發現中的一個重要問題是預測特定藥物是否可以結合特定蛋白質，即藥物-靶標相互作用（DTI）預測任務；
靶標是指體內具有藥效功能并能被藥物作用的生物大分子，如某些蛋白質和核酸等生物大分子
我們可以像下面這樣構造DTI預測任務：
描述：預測化合物與蛋白質結合親和力的二元分類輸入：化合物和蛋白質表示向量輸出：{0,1}或[0-1]中的實數(概率)讀寫分子
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函數file參數支持定義輸出位置；
操作分子
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’>；
根據原子編號獲取鍵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采用的是最小成環原則: 比如兩個環并在一起，從數學上看可以看成是2個小環，一個大環。在rdkit中計算環的時候只考慮小環。如下圖所示，我們構建了一個4元環并3元環的分子 “OC1C2C1CC2”：
我們可以看到該分子是由3元環并4元環組成的，其中原子2和原子3位于3元環和4元環的公共邊上：
其中：
- ‘IsInRing()’: 判斷是否在環上
- ‘IsInRingSize(n)’：判斷是否在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] """通過計算，我們發現示例分子一共有2個環，第一個環由3個原子（原子1,2,3）組成， 第二個環由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形式；
產生2D坐標AllChem.Compute2DCoords(m)；AllChem.Compute2DCoords(m)產生固定的獨一無二的取向，可用于作為繪制分子的模板，如果有多個分子共享一個骨架，我們希望繪圖的時候能夠保證在這些分子中的骨架方向是一致的。 首先把骨架提取成繪圖模板，然后在繪圖上添加基團，舉例子：
我們可以看到這4個結構都含有吡咯并嘧啶（c1nccc2n1ccc2）的子結構， 但是在圖片上他們的子結構取向不一致，不便于比較他們的結構。 我們可以采用模板繪制法，固定公共子結構的取向：
這樣我們就可以很清楚的看出這4個分子取代基團和位置的差異。
指紋和相似性
RDKit 內置了多種分子指紋計算方法，如：
- 拓撲指紋 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个网站，为设计师提供用不完的免费素