Butinaモジュールによる類似化合物のクラスタリング

from rdkit import DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit.ML.Cluster import Butina
import numpy as np
import pandas as pd
import time
import matplotlib.pyplot as plt
import collections
#パッケージの読み込み

df = pd.read_pickle("Asinex_Building_Blocks.pkl")
#Asinex_Building_Blocks.pklを読み込み、得られたデータフレームをdfと定義した

mols = df['ROMol']
morgan_fp = [AllChem.GetMorganFingerprintAsBitVect(x, 2, 2048) for x in mols]
#Morganフィンガープリント(半径2, 2048ビット)を作成し、morgan_fpに格納した。

#タニモト距離の三角行列を作成
t1 = time.time()
dis_matrix = []
for i in range(1, len(morgan_fp)):
    similarities = DataStructs.BulkTanimotoSimilarity(morgan_fp[i], morgan_fp[:i], 
                                                      returnDistance = True)
    dis_matrix.extend(similarities)
t2 = time.time()
print('t2 - t1:', t2 - t1)
#タニモト距離の距離行列を作成し、中身をdis_matrixに格納した。
#t2 - t1: 48.877509355545044

cutoff = 0.5    #類似度0.5以上のクラスターをつくる。

t3 = time.time()
clusters = Butina.ClusterData(dis_matrix, len(morgan_fp), cutoff, isDistData = True)
t4 = time.time()
print('t4 - t3:', t4 - t3) 
#Butinaクラスタリングを実施し、結果をclustersに格納した。6490個のクラスターが作成された。
#t4 - t3: 24.776731967926025

#お試しで適当なクラスターの中身を覗いてみる。
num_cluster_random = np.random.randint(0, len(clusters))
cluster_random = [mols[x] for x in clusters[num_cluster_random]]
Draw.MolsToGridImage(cluster_random, molsPerRow = 5)

clusters_count = collections.Counter(len(x) for x in clusters)
#クラスターの要素数ごとに何クラスターあるのか数え、clusters_countに格納した。

#クラスターあたりの要素数をx軸、対応するクラスター数をy軸にプロットする。
left = [x for x in clusters_count.keys()]
height = [x for x in clusters_count.values()]
plt.xlabel('Number of clusters', fontsize=20)
plt.ylabel('Number of molecules', fontsize=20)
plt.bar(left, height)
plt.show()

#クラスター番号をx軸、格納されている分子数をy軸にプロットする。
fig = plt.figure(1, figsize=(10, 4))
plt1 = plt.subplot(111)
plt.axis([0, len(clusters), 0, len(clusters[0])+1])
plt.xlabel('Cluster index', fontsize=20)
plt.ylabel('Number of molecules', fontsize=20)
plt.tick_params(labelsize=16)
plt1.bar(range(1, len(clusters)), [len(c) for c in clusters[:len(clusters)-1]], lw=0)
plt.show()

large_clusters = [x for x in clusters if len(x) >= 100]
#要素数が100以上のクラスターをピックアップした。