化學資訊學入門：RDKit 基本操作 part 1 + 分子相似度計算
RDKit是一個可以處理各種化學資訊的開源函式庫。安裝方法在之前的文章有說明過。
另外，在RDKit 的操作中，我們常會使用 SMILES 字串進行分子的操作。關於 SMILES 的部分，可以參考下面的文章。
這次我們會介紹如何使用 RDKit 來讀取分子、畫出結構圖、搜尋子結構…等基本操作。
RDKit 測試
首先，我們先導入 RDKit 最基本函數 Chem，嘗試把一個 SMILES 轉換成 mol 形式。我們嘗試把 “CCCC” 的結構畫出來。
import rdkit
from rdkit import Chem
mol = Chem.MolFromSmiles('CCCC')
mol
canonical SMILES
一個分子通常具有多個同樣有效的 SMILES 字串，像是乙醇可以寫成CCO、OCC或 C(O)C 。這個情況會造成一些應用上的問題，像是利用語言模型來學習 SMILES 時，可能要把乙醇所有的 SMILES 都拿去學習才行。為了解決這個情況，我們可以使用 canonical SMILES（規範 SMILES ），每個分子只有一個 canonical SMILES。
我們使用 RDKit 來觀察一個分子通常具有多個 SMILES 字串的狀況：
smis = ["CN2C(=O)N(C)C(=O)C1=C2N=CN1C",
"CN1C=NC2=C1C(=O)N(C)C(=O)N2C"]
mols = [Chem.MolFromSmiles(smi) for smi in smis]
再來是如何用 RDKit 輸出 canonical SMILES：
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi),canonical=True) for smi in smis]
cans[‘Cn1c(=O)c2c(ncn2C)n(C)c1=O’, ‘Cn1c(=O)c2c(ncn2C)n(C)c1=O’]
只要設定 canonical=True ，RDKit 可以把相同分子不同的 SMILES 轉換成 相同的 canonical SMILES。
Tautomer 互變異構物
互變異構是某些有機化合物的結構在兩種官能基異構物間產生平衡互相轉換的現象。大多數互變異構都涉及氫原子或質子的轉移，以及單鍵向雙鍵的轉變。這個現象會造成同一個分子會有幾個不同化學結構。使用 RDKit 也可以把互變異構物全部繪製出來。首先出入一個擁有互變異構物的分子。
from rdkit import Chem
m = Chem.MolFromSmiles('Oc1c(cccc3)c3nc2ccncc12')
m
像是 canonical SMILES 一樣，如果沒有一個標準化的結構，在使用上會產生問題。所以互變異構物中也是有一個標準結構，我們可以利用 rdMolStandardize 的一些功能輸出標準結構：
from rdkit.Chem.MolStandardize import rdMolStandardize
enumerator = rdMolStandardize.TautomerEnumerator()
enumerator.Canonicalize(m)
rdMolStandardize 也可以把所有互變異構物畫出來：
tauts = enumerator.Enumerate(m)
Draw.MolsToGridImage(tauts)
Grid image
當我們同時想要畫出多個分子的時候，RDKit 提供了 MolsToGridImage 這個函數來繪製多個分子。先取得分子的 mol 形式，然後使用MolsToGridImage 就可以繪製多個分子了，我們也可以指定圖的大小以及每行要畫近幾個分子。
from rdkit.Chem import Draw
smiles = [
'N#CC(OC1OC(COC2OCC(O)C(O)C2O)C(O)C(O)C1O)c1ccccc1',
'c1ccc2c(c1)ccc1c2ccc2c3ccccc3ccc21',
'C=CC1Cc2c(ccc3c2OC2COc4cc(OC)ccc4C2C3=O)O1',
'ClC(Cl)C(c1ccc(Cl)cc1)c1ccc(Cl)cc1'
]
mols = [Chem.MolFromSmiles(smi) for smi in smiles]
Draw.MolsToGridImage(mols, molsPerRow=2, subImgSize=(200, 200))
只提供了最基礎的功能，雖然可以在圖裡加上分子名稱…等追加訊息，但準備上有些麻煩，而且只能輸出圖片的形式。如果需要繪製更精美的分子圖，建議使用 mols2grid，在操作上比較簡單一些，且可以輸出互動式的 html 形式分子圖。
Substructure match
再來要介紹一個化學資訊學中常用的技巧：子結構探索 (Substructure match) ，像是從乙醇的分子內找到羥基(-OH)。HasSubstructMatch 這個函數可以輸出分子中是否有子結構，GetSubstructMatch 和 GetSubstructMatches則會輸出子結構的原子編號。我們用酚的例子來做說明：
- 尋找酚裡的羥基(-OH)
首先，我們先輸入酚 (c1ccccc1O) 的結構，想要尋找的子結構需要用 SMARTS 形式進行輸入，羥基是用O代表，如果不習慣去氫的寫法，也可以使用[OH]來表示：
- 尋找酚裡的芳香族碳
在 SMILES 中，字母大小寫代表了不同的意義，芳香族的碳是用小寫c代表，這裡比較三個相似的函數:HasSubstructMatch: 尋找是否有子結構，輸出是 True / FalseGetSubstructMatch: 尋找子結構的位置，輸出一組原子編號
GetSubstructMatches
Molecular similarity 分子相似度
我們之前介紹了如何使用 Python 計算分子指紋，分子指紋不只是在建模時候的特徵量，另一個重要的用法就是比較分子之間的相似度。
為了比較分子的相似度，我們先載入 ESOL 數據集來取得一些分子。再把分子用 PandasTools 轉換成 mol 形式。
import pandas as pd
from rdkit.Chem import PandasTools
df = pd.read_csv('esol.csv')
PandasTools.AddMoleculeColumnToFrame(esol_data, smilesCol='smiles')選定一個分子當作參考分子，為了方便看出分子相似度，我們把 Benzothiazole 當作參考分子，然後讓 RDKit 讀取這個分子。
from rdkit import Chem
ref_smiles = 'S1C=NC2=C1C=CC=C2'
ref_mol = Chem.MolFromSmiles(ref_smiles)
ref_mol
再來是計算分子指紋，這裡先用 Morgan Fingerprint 來測試，分別計算參考分子的分子指紋 ref_morgan_fps ，還有 ESOL 數據集裡分子的分子指紋 bulk_morgan_fps ，再來是利用 FingerprintSimilarity 來計算分子相似度，最後用 Grid Image 把結果呈現出來。
from rdkit.Chem import AllChem
from rdkit import DataStructs
ref_morgan_fps = AllChem.GetMorganFingerprintAsBitVect(ref_mol,2)
bulk_morgan_fps = [AllChem.GetMorganFingerprintAsBitVect(x,2) for x in esol_data['ROMol']]
similarity_morgan = [DataStructs.FingerprintSimilarity(ref_morgan_fps,x) for x in bulk_morgan_fps]
esol_data['Similarity'] = similarity_morgan
PandasTools.FrameToGridImage(esol_data.head(8), legendsCol="Similarity", molsPerRow=4)
分子相似度是個從 0 到 1 的分數，0 分代表完全不同的分子，1 則是代表完全相同的分子。從上圖可以看到，數據集的分子裡存在 Benzothiazole 和參考分子的結構完全一致，分子相似度等於 1 。至於其他擁有苯環的分子都有大於 0 的相似度。而沒有苯環的分子相似度等於 0。
我們換一個分子指紋試試看，這次嘗試使用 RDKit fingerprint，計算用的函數和 Morgan fingeprint 稍微不同，以下是計算用的程式碼和結果。
ref_rdkit_fps = Chem.RDKFingerprint(ref_mol)
bulk_rdkit_fps = [Chem.RDKFingerprint(x) for x in esol_data['ROMol']]
similarity_morgan = [DataStructs.FingerprintSimilarity(ref_rdkit_fps,x) for x in bulk_rdkit_fps]
esol_data['Similarity'] = similarity_morgan
PandasTools.FrameToGridImage(esol_data.head(8), legendsCol="Similarity", molsPerRow=4)
我們會發現，不同的分子指紋會得到不同的分子相似度，相同的 Benzothiazole 分子相似度還是 1，但其他的分子相似度卻不太相同。所以在比較分子相似度時，分子指紋的選擇也是很重要的因素。
分子相似度是怎麼計算的？
根據前面的說明，在考慮化合物的相似性時，需要使用分子指紋來進行香速度計算。那麼，當給出兩個化合物的分子指紋時，相似度的數值是如下得到的？
相似度的計算方法有很多中，最常用的是所謂的“谷本係數 (Tanimoto coefficient / Tanimoto similarity) ”。要從兩個分子A和B的分子指紋計算 Tanimoto similarity，我們使用以下公式。
- a 是分子 A 的分子指紋中 1 的數量
- b 是分子 B 的分子指紋中 1 的數量
- c 是分子 A 和 B 中共有的 1 的數量
讓我們使用下面的例子來計算 Tanimoto similarity：
其中 a = 5 、 b = 6、c = 3 套用上述的公式後就可以算出 Tanimoto similarity：
結語
這次我們簡單的介紹一些 RDKit 的常用功能，從最基本的 SMILES 讀取、搜尋子結構、繪製大量分子，以及比較分子的相似度。這些都是 2D 結構為主，下一次會介紹 3D 結構的一些功能。
