对于分子聚类后的结果,或者是某个化合物数据集,有时希望可视化后有个直观的对比来确认不同来源的分子集合所占据的化学空间。 但对于化合物这种高维度数据,直接可视化是不行的,需要先降到2维或3维才能在平面或者立体空间内展示。降维方式多种多样,这里演示最常用的t-SNE和PCA的降到二维平面的做法。
注意一点是,t-SNE与PCA降维效果见仁见智,虽然通认t-SNE在降维中应用更广,但个人认为对于化合物而言没有明确的优劣之分,但一般来说t-SNE相比PCA降维的速度慢数倍,处理分子规模较大时,可以优先考虑PCA方式。
一、tSNE
t-Distributed Stochastic Neighbor Embedding (t-SNE),创建了一个缩小的特征空间,相似的样本由附近的点建模,不相似的样本由高概率的远点建模。在高水平上,t-SNE为高维样本构建了一个概率分布,相似的样本被选中的可能性很高,而不同的点被选中的可能性极小。然后,t-SNE为低维嵌入中的点定义了相似的分布。最后,t-SNE最小化了两个分布之间关于嵌入点位置的Kullback-Leibler(KL)散度。
我们先将化合物转成分子指纹后,再实现可视化
from rdkit import DataStructs
from rdkit.Chem import AllChem
import numpy as np
class FP:
"""
建立一个FP的类方便后续的分子指纹处理
"""
def __init__(self, fp, names):
self.fp = fp
self.names = names
def __str__(self):
return "%d bit FP" % len(self.fp)
def __len__(self):
return len(self.fp)
def get_cfps(mol, radius=2, nBits=1024, useFeatures=False, counts=False, dtype=np.float32):
arr = np.zeros((1,), dtype)
if counts is True:
info = {}
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=nBits, useFeatures=useFeatures, bitInfo=info)
DataStructs.ConvertToNumpyArray(fp, arr)
arr = np.array([len(info[x]) if x in info else 0 for x in range(nBits)], dtype)
else:
DataStructs.ConvertToNumpyArray(AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=nBits, useFeatures=useFeatures), arr)
return FP(arr, range(nBits))
def calFP(smi):
try:
mol = Chem.MolFromSmiles(smi)
fp = get_cfps(mol)
return fp
except Exception as e:
return None
以一个pandas 表格为例, 表格最初仅包含SMILES与所属的类别两列,
df['FP'] = df['canonSMILES'].apply(calFP)
df.head(3)
SMILES cluster FP
Clc1ccc(-c2c[nH]c(-c3cccnc3)n2)cc1 1 1024 bit FP
C=CCOC(=O)C1=C(C)N=C(C)C(C(=O)OCC=C)C1c1cn(-c2... 1 1024 bit FP
O=C(CCc1c[nH]c2ccccc12)NC(Cc1ccccc1)C(=O)O 2 1024 bit FP
然后对FP这一列进行tSNE降维
from sklearn.manifold import TSNE
tsne_model = TSNE(n_components=2, random_state=42, perplexity=30, n_iter=5000)
# 将表格中的分子指纹转为numpy array的形式
X = np.array([x.fp for x in df['FP']])
tsne_result = tsne_model.fit_transform(X)
# 将降维结果加到表格中去
df['tSNE_1'] = tsne_result.T[0]
df['tSNE_2'] = tsne_result.T[1]
df.head(3)
SMILES cluster FP tSNE_1 tSNE_2
Clc1ccc(-c2c[nH]c(-c3cccnc3)n2)cc1 1 1024 bit FP 122.361221 -21.181770
C=CCOC(=O)C1=C(C)N=C(C)C(C(=O)OCC=C)C1c1cn(-c2... 1 1024 bit FP -64.470451 -21.792231
O=C(CCc1c[nH]c2ccccc12)NC(Cc1ccccc1)C(=O)O 2 1024 bit FP 4.770935 64.870430
这时我们可以可视化刚刚降维的结果了。
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
sns.scatterplot(x='tSNE_1', y='tSNE_2', hue='cluster', s = 15, data=df, alpha=0.2, edgecolors='none')
plt.show()
二、PCA降维
PCA降维核心在于提取数据的主要特征分量,使得降维后的各主成分的方差值尽可能保持降维前的总方差。
类似的,我们仅需在计算出FP的表格中,应用PCA的方式,代替tSNE,即可实现降维与可视化
from sklearn.decomposition import PCA
pca_model = PCA(n_components=2, random_state=42)
X = np.array([x.fp for x in df['FP']])
pca_result = pca_model.fit_transform(X)
df['PCA_1'] = pca_result.T[0]
df['PCA_2'] = pca_result.T[1]
sns.set()
sns.scatterplot(x='PCA_1', y='PCA_2', hue='cluster', s = 15, data=df, alpha=0.2, edgecolors='none')
plt.show()