RDKitには化合物の類似度に基づいてクラスタリングを行うモジュールが用意されています。
その中の一つにButinaモジュールと呼ばれるクラスタリングアルゴリズムがあり、
化合物間の距離行列を用意することで類似化合物の分類を行うことが可能です。

今回はMorganフィンガープリントを基に化合物間のタニモト距離を計算し、
それをButinaモジュールで解析して化合物のクラスタリングを行う方法を紹介します。

###環境

Windows 11
python 3.8.8
RDKit 2021.09.4

###下準備
今回もAsinex社のビルディングブロックを題材として扱います。
まずは以前の記事中で作成したpickleファイルを読み込みます。

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と定義した

読み込んだ分子群からMorganフィンガープリントを作成します。

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

22,525化合物の距離行列の計算には約1分の時間を要します。
要素数が大きくメモリの消費量も著しいため、以下のButinaモジュール使用後は
dis_matrixの変数を削除してしまってもいいかもしれません。

###Butinaモジュールによるクラスタリング
Butina.ClusterData関数を使用し、クラスタリングを実施します。

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

この関数を使用する際は主に以下の引数を指定します。
（isDistDataについてよく理解できなかったため、参考例に倣ってTrueを指定しています。
　詳しくはrdkit.ML.Cluster.Butinaモジュールの公式ドキュメントをご覧ください。）

クラスタリングに用いる距離行列
クラスタリングを行う分子数
どれぐらいの化合物距離でクラスタリングするか
isDistDataはとりあえずTrueを指定する（デフォルトだとFalse）

今回はcutoff = 0.5を指定したため、タニモト距離が0.5以下のクラスターが形成されました。
この値を例えば0.2にするとタニモト距離が0.2以下の分子群がクラスタリングされます。
値を小さくすればより細分化され類似度の高い分子群を形成でき、逆に大きな値を指定すれば
類似度が下がるもののより要素数の大きなクラスターが形成されることになります。

Butina.ClusterData関数の戻り値は化合物番号が格納されたタプルのタプルとなります。

格納された化合物番号を利用すれば以下のように中身の分子を取り出すことが可能です。
以下の3行を実行し、適当なクラスターの中身を図示してみましょう。
（ただし多くのクラスターは格納されている分子数が1の小さなクラスターです）

#お試しで適当なクラスターの中身を覗いてみる。
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)

(num_cluster_random = 260の場合)

なんとなく似た化合物が集まっていることが目視で分かるかと思います。

###出力結果を用いたデータ閲覧・加工の例

どんな大きさのクラスターが何個ずつあるのか図示して確認してみます。
要素数1のクラスターが大多数であることが確認できます。

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()

要素数が大きなクラスターだけ集めると後々役に立つことがあるかもしれません。
（他のクラスタリング手法の教材にする等）
試しに要素数が100以上のクラスターだけ集めてみましょう。

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

###おわりに
Butinaクラスタリングモジュールを用いたクラスタリングを実施いたしました。
今回Morganフィンガープリントから計算したタニモト距離を基にクラスタリングしましたが、
その他の指標による距離行列を用意すれば同様の手法でクラスタリングが可能です。

類似化合物のクラスタリングができれば、各クラスターから代表分子を選ぶことで
偏りが小さく多様な分子が含まれるライブラリーを構築することが可能です。

###参考

###これまでの記事のシリーズ

RDKit入門①：まずは環境構築
 RDKit入門②：1分子の読み込みと部分構造検索
 RDKit入門③：複数分子の読み込み (前編)と分子群の描画
 RDKit入門④：複数分子の読み込み (後編)とデータフレームの加工
 RDKit入門⑤：データフレーム内の分子群に対する部分構造検索
 RDKit入門⑥：Morganフィンガープリントの作成とそれを用いたタニモト係数の計算による分子類似性評価
 RDKit入門⑦：反応式の取り扱い(前編)
RDKitで化学反応を扱った際の生成物群を整理し図示する方法
 RDKit入門⑧：反応式の取り扱い(後編)
RDKitで化学反応を繰り返しつつ分子の重複数をカウント・記載する
 RDKit入門⑨：データフレーム内の分子群に対する化学反応の適用とExcelファイルへの出力
 RDKit入門⑩：データフレーム内の分子群に対する完全一致検索
 指定した回数だけ化学反応を実施する関数の作成
 BRICSフラグメントを利用した分子構造生成関数の簡略化