化学分子降维算法umap、聚类算法hdbscan和kmeans

参考：
https://github.com/iwatobipen/chemo_info/blob/master/chemicalspace2/HDBSCAN_Chemoinfo.ipynb

https://www.cnblogs.com/liulunyang/p/14275612.html

本文章主要方法：
1、通过umap获取分子的降维特征数据
2、再通过hdbscan或kmeans聚类方法聚类（hdbscan不支持指定最后聚类总数量，可能非常多；kmeans可以指定总的聚类数）

化学分子降维算法umap、聚类算法hdbscan和kmeans

代码已经上传github：https://github.com/lonngxiang/molecular-clustering

安装

## 安装
1）pip install umap-learn  -i https://pypi.douban.com/simple2）hdbscan pip安装报错：
ERROR: Failed building wheel for hdbscan
Failed to build hdbscan
ERROR: Could not build wheels for hdbscan, which is required to install pyproject.toml-based projects解决方法conda安装：conda install -c conda-forge hdbscanconda安装后调用报错：
__init__() got an unexpected keyword argument 'cachedir'解决方法joblib降级：
pip install joblib==1.1.0

代码：

##第一步  指纹获取，下来是2048维
from rdkit import Chem
from rdkit.Chem import DataStructs,AllChem, Draw
import numpy as npX = []
for sim in results["hits"]["hits"]:fragment = sim["_source"]["fragment"]mol = Chem.MolFromSmiles(fragment)arr = np.zeros((0,))fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)DataStructs.ConvertToNumpyArray(fp, arr)X.append(arr)
print('{} mols loaded'.format(len(X)))## 第二部 umap进行降维，默认是降维成2维，umap.UMAP里可以修改n_componentsimport umapdef tanimoto_dist(a,b):dotprod = np.dot(a,b)tc = dotprod / (np.sum(a) + np.sum(b) - dotprod)return 1.0-tcumap_X = umap.UMAP( metric=tanimoto_dist).fit_transform(X)## 第三部 对umap降维处理的数据再次用hdbscan进行降维
import hdbscan
from collections import Countercluster_umap = hdbscan.HDBSCAN(min_cluster_size=5, gen_min_span_tree=True)
cluster_umap.fit(umap_X)cluster_umap.labels_      ##（hdbscan降维后的聚类数）
Counter(cluster_umap.labels_).keys(),Counter(cluster_umap.labels_)  ## 查看各聚类中心的数据量 ## 或kmeans进行聚类
from sklearn.cluster import KMeanskmeans=KMeans(n_clusters=10) ##指定10个聚类中心
kmeans.fit(umap_X)
kmeans.labels_

## 最终聚类效果
import seaborn as snspalette = sns.color_palette()
label_prob = zip(cluster_umap.labels_, cluster_umap.probabilities_)
plt.scatter(umap_X.T[0], umap_X.T[1], c = cluster_umap.labels_, cmap='plasma')

化学分子抽取同一聚类的数据查看

这是聚类0中心的数据展示；my_score函数参考这里（sa_score化合物合成难度打分）：https://blog.csdn.net/weixin_42357472/article/details/127546269

cluster_umap_dict ={}
for index, value in enumerate(cluster_umap.labels_.tolist()):print("索引：" ,index, ", 值：" ,value)value =str(value)if value not in cluster_umap_dict:cluster_umap_dict[value] = [index]else:cluster_umap_dict[value].append(index)new_mols2 = []
sa_scores2 = []for ii in cluster_umap_dict["0"]:  ##这是聚类0中心的数据展示sim = results["hits"]["hits"][ii]   ##results["hits"]["hits"]  这是总共分子的原始数据fragment = sim["_source"]["fragment"]if Counter(fragment)["*"] ==1:# 统计词频new_fragment = ("N(C([C@H](N)CCCCN)=O)[C@H](C)C*" + fragment).replace("*","")# print(new_fragment)mol = Chem.MolFromSmiles(new_fragment)new_mols2.append(mol)sa_score = my_score(mol)sa_scores2.append(sa_score)Draw.MolsToGridImage(new_mols2,subImgSize=(300,300), molsPerRow=5,legends=["sa_score:  "+str(ii)for ii in sa_scores2])

这是聚类20中心的数据展示

new_mols2 = []
sa_scores2 = []for ii in cluster_umap_dict["20"]:sim = results["hits"]["hits"][ii]fragment = sim["_source"]["fragment"]if Counter(fragment)["*"] ==1:# 统计词频new_fragment = ("N(C([C@H](N)CCCCN)=O)[C@H](C)C*" + fragment).replace("*","")# print(new_fragment)mol = Chem.MolFromSmiles(new_fragment)new_mols2.append(mol)sa_score = my_score(mol)sa_scores2.append(sa_score)Draw.MolsToGridImage(new_mols2,subImgSize=(300,300), molsPerRow=5,legends=["sa_score:  "+str(ii)for ii in sa_scores2])

hdbscan预测新数据approximate_predict

参考：https://hdbscan.readthedocs.io/en/latest/prediction_tutorial.html

注意： hdbscan.HDBSCAN里prediction_data=True需要设置为True；二是hdbscan.approximate_predict预测里数据格式是array，如果是单条需要放列表里

import hdbscancluster_umap = hdbscan.HDBSCAN(min_cluster_size=5, gen_min_span_tree=True, prediction_data=True)
cluster_umap.fit(umap_X)##测试新数据cluster, prob = hdbscan.approximate_predict(cluster_umap,[[ 7.5436993 ,  5.556767  ],[ 8.059909  ,  1.9189416 ]])
cluster, prob