ケモインフォマティクスにおける逆解析では、何か有用そうな化合物構造が得られたとしても、>その化合物が実際に合成できるかどうか<はまた別の話。
仮に化合物提案ができたとしても合成が出来なかったり、そもそもそのような構造が存在し得ない可能性もある。
SA score を用いて化合物の合成のしやすさの見積もることで、有望かつ合成可能性の高い化合物をスクリーニングすることが出来る。
SA スコアとは
- 化合物の合成のしやすさを評価する指標(Synthetic Accessibility score の略)
- 分子構造の複雑さを基準として、スコアを計算している。
- 1~10の値をとり、低いほど合成が容易、高いほど合成困難
SAスコアの元文献は以下を参照
より詳細な説明は、以下のブログに載っています。
RDKitで合成難易度を評価して化合物をスクリーニング(新しい化学のカタチ)
詳細な計算方法
SA score は以下の式で計算される。
(実際には-4〜2.5の値になるものを、-1を掛けて1~10にスケーリングしているとのこと)
$$ SAscore = fragmentScore - complexityPenaty $$
fragmentScoreは、「一般的な化合物によく見られる部分構造は合成しやすい」という過程のもと作成された値。
- Pubchemから取得した100万個の化合物のECFP4フィンガープリントの出現頻度をもとに、各フィンガープリントの重み付けを行い、
- SA scoreを調べたい分子構造に含まれる全てのフラグメントの重みの和を計算し、それをフラグメントの数で割っている。
(間違ってたらすみません...)
The score is calculated as a
sum of contributions of all fragments in the molecule
divided by the number of fragments in this molecule
例)Pubchemの100万化合物の中で最も一般的だった構造(The most common fragments)
「Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions」, J. Cheminf., 2009, 1, 8.
complexityPenaltyは、単純に環状構造, 立体, 大員環, 分子サイズなどの要素を考慮したもの。
計算式は以下とのこと。(詳細は文献参照)
$$ ringComplexityScore = log(nRingBridgeAtoms +1) + log(nSpiroAtoms +1) $$
$$ stereoComplexityScore = log(nStereoCenters+1) $$
$$ macrocyclePenalty = log(nMacrocycles + 1) $$
$$ sizePenalty = natoms**1.005 -natoms $$
使い方
前準備
RDKitのコアモジュールには含まれて無いので、GiHubのContrib/SA_Scoreから必要なスクリプト(sascorer.py)とデータベース(fpsocres.pkl.gz)を入手しておく。
データ
サンプルデータとして化合物可視化ライブラリ「chemplot」に準備されている化合物の水溶解度データセットに含まれている1290化合物を用いる。(「chemplot」についてはこちら)
モジュールインポート〜データの準備まで
import pandas as pd
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.Chem import Draw
from chemplot import load_data
from SA_Score import sascorer
df = load_data("LOGS")
df.head()
smiles | target | |
---|---|---|
0 | CCCCC | -3.18 |
1 | C1CCCC1 | -2.64 |
2 | CCCCCC | -3.84 |
3 | CCCC(C)C | -3.74 |
4 | CCC(C)(C)C | -3.55 |
SA score の算出
# データの準備
smiles_list = df['smiles'].values.tolist()
#RDkitのMOLオブジェクトに変換
mols_list = [Chem.MolFromSmiles(smiles) for smiles in smiles_list]
# 合成難易度を評価する
# 一分子だけ計算する
print(sascorer.calculateScore(mols_list[0]))
sa_score_list = [sascorer.calculateScore(mol) for mol in mols_list]
df['sa_score'] = sa_score_list
"""output
1.699621281696647
"""
# 10化合物だけ見てみる
print(sa_score_list[:10])
"""output
[1.699621281696647, 1.0, 1.2090305482861616, 1.8933213175169303, 1.9944749269843864, 1.0, 1.5726834894626336, 1.2607134650665888, 1.5486420480625398, 1.304762875769466]
"""
## SA scoreの分布を見てみる。
fig = plt.figure()
plt.hist(sa_score_list, bins=20)
plt.title('SA score',fontsize=15)
plt.xlabel('SA socre',fontsize=15)
plt.ylabel('count',fontsize=15)
plt.show()
SA_scoreの高い化合物、低い化合物をチェックしてみる。
# SA_scoreの高い化合物、低い化合物をチェックしてみる。
df.sort_values('sa_score',inplace=True)
# スコアの抽出
low_sascore = df['sa_score'][:5].astype(str).values.tolist()
high_sascore = df['sa_score'][-5:].astype(str).values.tolist()
# 構造情報の抽出
low_sascore_smiles = df['smiles'][:5]
high_sascore_smiles = df['smiles'][-5:]
low_sascore_mols = [Chem.MolFromSmiles(smiles) for smiles in low_sascore_smiles]
high_sascore_mols = [Chem.MolFromSmiles(smiles) for smiles in high_sascore_smiles]
# SA_scoreの高い10化合物、低い10化合物をチェックする。
img = Draw.MolsToGridImage(low_sascore_mols + high_sascore_mols,
molsPerRow=5,
subImgSize=(200,200),
legends=low_sascore+high_sascore
)
img
一応、SA score の値と構造の複雑さ(合成の困難さ)の感覚は合いそう。
注意点
分子に水素を付加する(Chem.AddHs()など)と、原子数のカウントに水素が含まれてしまい、SAスコアのペナルティ項のサイズペナルティが高くなり、SAscoreが高くなってしまうとのこと。
(参考|【RDKit】SAscoreの計算時に注意するポイント)
試しに水素を付加して計算してみる。
mols_list_h = [Chem.AddHs(mol) for mol in mols_list]
sa_score_list_h = [sascorer.calculateScore(mol) for mol in mols_list_h]
# 水素あり、水素なしのヒストグラム
fig = plt.figure()
plt.hist(sa_score_list, alpha=0.7, bins=20,label='no H atom')
plt.hist(sa_score_list_h, alpha=0.7,bins=20,label='add H atom')
plt.title('SA score',fontsize=15)
plt.xlabel('SA score',fontsize=15)
plt.ylabel('count',fontsize=15)
plt.legend(loc='best')
plt.show()
確かに、水素を付加すると総じてSA scoreが高くなっている
以上
参考
参考HP
関連書籍