對分子進行聚類分析，首先必須要考慮的是其描述符的問題，分子描述符通常是非常高維的，必須對其進行降維才好繼續后面的分析，特別分子量特別大的時候。常用的降維手段有PCA,TSNE和UMAP.一說，TSNE用於可視化.
聚類的方法有許多，比如k-means,層次聚類. 但是這兩個一個需要定義k,一個需要定義閾值,這樣需要試錯法合理進行着兩個量的設置，不是很方便. 因而，我選擇使用HDBSCAN,一個基於數據密度的聚類算法，參考文獻如下：
https://link.springer.com/chapter/10.1007%2F978-3-642-37456-2_14
我選擇使用 HDBSCAN.
先導入所需要的各種庫. 后面還會使用到 mlinsght 庫.
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit import RDLogger
from sklearn.model_selection import train_test_split
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.manifold import TSNE
from mlinsights.mlmodel import PredictableTSNE
from hdbscan import HDBSCAN
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()
plot_kwds = {'alpha' : 0.5, 's' : 80, 'linewidths':0}
RDLogger.DisableLog('rdApp.*')
seed = 794
HDBSCAN需要接受一個衡量分子間距離或者相似度的參數，可以自定義，這里我定義了tanimoto-dist函數，使用方法如下面的代碼所示：
使用的時候發現自定義的tanimoto速度太慢了。。。。
oxrs = [("CHEMBL3098111", "Merck" ), ("CHEMBL3867477", "Merck" ), ("CHEMBL2380240", "Rottapharm" ),
("CHEMBL3352684", "Merck" ), ("CHEMBL3769367", "Merck" ), ("CHEMBL3526050", "Actelion" ),
("CHEMBL3112474", "Actelion" ), ("CHEMBL3739366", "Heptares" ), ("CHEMBL3739395", "Actelion" ),
("CHEMBL3351489", "Eisai" )]
fps = []
docs = []
companies = []
mol_list = []
for cid, company in oxrs:
sdf_file = os.path.join("data", cid + ".sdf")
mols = Chem.SDMolSupplier(sdf_file)
for mol in mols:
if mol is not None:
mol_list.append(mol)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
docs.append(cid)
companies.append(company)
fps.append(arr)
fps = np.array(fps)
companies = np.array(companies)
def tanimoto_dist(ar1, ar2): # 接收兩個numpy 數組
a = np.dot(ar1, ar2)
b = ar1 + ar2 - ar1*ar2
return 1 - a/np.sum(b)
# 實例化HDBSCAN類
clusterer = HDBSCAN(algorithm='best',
min_samples=5,
metric='pyfunc',
func=tanimoto_dist)
clusterer.fit(fps)
docs = np.array(docs)
trainIDX, testIDX = train_test_split(range(len(fps)), random_state=seed)
接下來我們來對數據進行可視化，使用到了tsne技術，下面兩圖的分子分別以company作為顏色標簽和HDBSACA聚類結果作為顏色標簽.
tsne = TSNE(random_state=seed)
res = tsne.fit_transform(fps)
plt.clf()
plt.figure(figsize=(12, 6))
sns.scatterplot(res[:,0], res[:,1], hue=companies, **plot_kwds)
plt.clf()
plt.figure(figsize=(12, 6))
palette = sns.color_palette()
cluster_colors = [sns.desaturate(palette[col], sat)
if col >= 0 else (0.5, 0.5, 0.5) for col, sat in zip(clusterer.labels_, clusterer.probabilities_)]
plt.scatter(res[:,0], res[:,1], c=cluster_colors, **plot_kwds)
coloured with company name
coloured with hdbscan label.
圖中灰色的是沒有標簽的數據。結果表明HDBSCAN還是聚的不錯的。
映射新化合物到已定義的化學空間
下一個主題是如何在預定的化學空間中繪制新化合物。 從理論上講，PCA和tSNE不能將新數據添加到預定義的潛在分布空間。 如果用戶想添加新數據，則需要重新計算。 但是我們想將新設計的化合物映射到當前的化學空間。 嗯... 怎么做？
答案是： mlinghts. 這個庫包含 PredictableTSNE 類.從名字上就不難看出它是干啥的了。 這個類接收一個TSNE實例和一個用於預測的估計器.
在下面的例子中，我們使用 random forest 和 GP回歸器.
trainFP = [fps[i] for i in trainIDX]
train_mol = [mol_list[i] for i in trainIDX]
testFP = [fps[i] for i in testIDX]
test_mol = [mol_list[i] for i in testIDX]
allFP = trainFP + testFP
tsne_ref = TSNE(random_state=seed)
res = tsne_ref.fit_transform(allFP)
plt.clf()
plt.figure(figsize=(12, 6))
sns.scatterplot(res[:,0],
res[:,1],
hue=['train' for i in range(len(trainFP))] + ['test' for i in range(len(testFP))])
藍色的是訓練數據，橙色是測試數據.
使用predictiveTSNE 可以預測新化合物在化學空間的所在.下面的例子分別為RF和GP.
rfr = RandomForestRegressor(random_state=seed)
tsne1 = TSNE(random_state=seed)
pred_tsne_rfr = PredictableTSNE(transformer=tsne1,
estimator=rfr,
keep_tsne_outputs=True)
pred_tsne_rfr.fit(trainFP, list(range(len(trainFP))))
pred1 = pred_tsne_rfr.transform(testFP)
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_rfr.tsne_outputs_[:,0],
pred_tsne_rfr.tsne_outputs_[:,1],
c='blue',
alpha=0.5)
plt.scatter(pred1[:,0],
pred1[:,1],
c='red',
alpha=0.5)
gbr = GaussianProcessRegressor(random_state=seed)
tsne2 = TSNE(random_state=seed)
pred_tsne_gbr = PredictableTSNE(transformer=tsne2, estimator=gbr, keep_tsne_outputs=True)
pred_tsne_gbr.fit(trainFP, list(range(len(trainFP)))pred2 = pred_tsne_gbr.transform(testFP)
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_gbr.tsne_outputs_[:,0], pred_tsne_gbr.tsne_outputs_[:,1], c='blue', alpha=0.5)
plt.scatter(pred2[:,0], pred2[:,1], c='red', alpha=0.5))
RF
GP
看起來RF的效果要好些.
預測算法對TSNE的性能有很大的影響。 這對於繪制新數據很有用，但存在設計誤導的風險。
最后，檢查一下聚類標簽和preditiveTSNE 預測數據.
allmol = train_mol + test_mol
fps2 = []
clusterer2 = HDBSCAN(algorithm='best',
min_samples=5,
metric='pyfunc',
func=tanimoto_dist)
# calc fp
for mol in allmol:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps2.append(arr)
fps2 = np.array(fps2)
# cluster
clusterer2.fit(fps2)
# plot
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_rfr.tsne_outputs_[:,0],
pred_tsne_rfr.tsne_outputs_[:,1],
c=clusterer2.labels_[:len(trainFP)],
alpha=0.5,
marker='x')
plt.scatter(pred1[:,0],
pred1[:,1],
c=clusterer2.labels_[len(trainFP):],
alpha=0.5, marker='o')
x是訓練數據, o為測試數據, 顏色 :類別標簽.
不僅訓練數據而且具有相同顏色的測試數據也出現在幾乎相同的區域中。 這表明該方法在這種情況下效果很好。