Sanity Checks for Saliency Maps of Molecules
はじめに
これまでQSARの結果を可視化する方法を試してきました。
これらはそれぞれ視覚的になんとなく上手くモデルを説明できてそうに見えます。しかし、使う側からすると"なんとなく"では困ります。しかも前回の記事でも触れたように、そもそも説明可能性の手法自体にも疑問が呈されているのが最近の状況です。このような状況を考えると、できれば定量的に説明可能性を評価したいものです。
最近各種saliency methodが報告されている画像処理の分野でsaliency methodの正当性 (adequacy)を評価する方法が考案・報告されておりますので、今回はこの手法を化合物のQSARで試してみようと思います。
Sanity Checks for Saliency Maps explained in Japanese language - Speaker Deck
www.slideshare.net (機械学習の判断根拠全般の解説。pp.57-61にSanity Checkの解説記載)
"Sanity Checks for Saliency Maps"について
本論文では以下のrandomization testを提唱しています。
- Model parameter randomization test
学習したモデルのsaliency mapと、モデルの重みをランダムに初期化した場合のsaliency mapとを比較。両マップが似ているならば、それは学習したモデルの重みに依存しない (=似ていない方がsaliency methodとして妥当)。
- Data randomization test
学習したモデルのsaliency mapと、教師データの目的変数(ラベル)をランダムに入れ替えて学習した場合のsaliency mapとを比較。両マップが似ているならば、それは教師データに依存しない (=似ていない方がsaliency methodとして妥当)。
"似ている"ことの評価には、saliency mapのSSIMやHOG、Spearman順位相関係数を用いています。
Sanity Check for QSAR Model
今回はNN系ではないのでdata randomization testを適用しました。
学習に用いたトレーニングデータについて可視化を行います。
train_path = os.path.join(RDPaths.RDDocsDir, 'Book/data/solubility.train.sdf')
test_path = os.path.join(RDPaths.RDDocsDir, 'Book/data/solubility.test.sdf')
train_mols = [m for m in Chem.SDMolSupplier(train_path) if m is not None]
test_mols = [m for m in Chem.SDMolSupplier(test_path) if m is not None]
print(len(mols))
# 1282
sol_classes = {'(A) low': 0, '(B) medium': 1, '(C) high': 2}
X_train = np.array([mol2fp(m)[0] for m in train_mols])
y_train = np.array([sol_classes[m.GetProp('SOL_classification')] for m in train_mols], dtype=np.int)
X_test = np.array([mol2fp(m)[0] for m in test_mols])
y_test = np.array([sol_classes[m.GetProp('SOL_classification')] for m in test_mols], dtype=np.int)
ラベルをランダムにシャッフルしたモデルも作成します。
rf_models = []
clf = RandomForestClassifier(random_state=20200119)
clf.fit(X_train, y_train)
print(f'Train: {accuracy_score(y_train, clf.predict(X_train)):.3f}, Test: {accuracy_score(y_test, clf.predict(X_test)):.3f}')
rf_models.append(clf)
for _ in range(10):
y_perm = np.random.permutation(y_train)
_clf = RandomForestClassifier(random_state=20200119)
_clf.fit(X_train, y_perm)
rf_models.append(_clf)
print(f'Shuffuled Data: {accuracy_score(y_perm, _clf.predict(X_train)):.3f}, True Train: {accuracy_score(y_train, _clf.predict(X_train)):.3f}, Test: {accuracy_score(y_test, _clf.predict(X_test)):.3f}')
# Train: 0.992, Test: 0.685
# Shuffuled Data: 0.985, True Train: 0.352, Test: 0.370
# Shuffuled Data: 0.985, True Train: 0.371, Test: 0.401
# Shuffuled Data: 0.985, True Train: 0.363, Test: 0.350
# Shuffuled Data: 0.980, True Train: 0.368, Test: 0.358
# Shuffuled Data: 0.981, True Train: 0.354, Test: 0.389
# Shuffuled Data: 0.983, True Train: 0.365, Test: 0.389
# Shuffuled Data: 0.984, True Train: 0.355, Test: 0.393
# Shuffuled Data: 0.982, True Train: 0.363, Test: 0.444
# Shuffuled Data: 0.986, True Train: 0.378, Test: 0.475
# Shuffuled Data: 0.980, True Train: 0.359, Test: 0.401
論文ではvisual inspection / assessment はあまりよろしくないとの述べられておりましたが、一応ランダムにシャッフルした結果を可視化してみましょう。rdkitのSimilarityMaps.GetAtomicWeightsForModelで描画しました。
正しいデータで学習した場合とシャッフルした場合では、全く違う重みが得られており、saliency methodとして妥当そうに見えます。
次に構造式での描画だけでなく、グラフでも示してみます。正しいデータで学習した時の重みを折れ線グラフ (緑)と右の補助図で示し、シャッフルした時の重みを箱ひげ図で示しました。
正しいデータで学習したモデルの重みとランダムにシャッフルした場合の重みを比較すると、シャッフルした場合の重みが元と全く違う値を示していることが分かります。saliency methodとしては妥当であると言えそうです。
では、シャッフルの前後でsaliencyの結果がどのように変わるか、類似性で評価してみます。原子ごとの重みがarrayで得られていますので、元の重みとシャッフルした場合の重みの相関係数 (Pearson, Spearman)を見てみます。
クラス数が少ないためか、ランダムな試行の中で相関が高く出ている時もありますが、10回の試行全体で見ると相関係数は低く、このチェックでも妥当そうな結果です。
では、データセット全体で重みを計算してみます。RandomForestでのSimilarityMapsに加えて、LightGBMとSHAPも計算してみました。
RandomForest / LightGBM、SimilarityMaps / SHAPいずれも0付近の相関係数を示しており、saliency methodとしては妥当と言えそうです。
RandomForestのトレーニングデータとテストデータでも比較してみました。
トレーニングとテストでAccuracyに差があったので、それが何らか相関係数の差にも反映されないかと思って見てみましたが、ほとんど違いは認められないですね。
おわりに
画像データで提唱されているsanity checkを化合物のQSARにも適用してみました。rdkitのSimilarityMapsもSHAPも、randomization testの結果から妥当なsaliency methodであると言えそうです。一方、定量的にどの程度のばらつきや閾値であれば妥当と言えるのかの感覚は今回の実験だけでは掴めなかったので別のデータセットやアルゴリズムでも試してみます。
[2020/01/29 07:53 typo修正]
[2020/02/01 0:14 原先生のスライドシェアのリンク追加]
Visualization with chainer-saliency
はじめに
以下の2回で補助図を用いた化合物の可視化を試してみました。今回は補助図描画をchainer-saliency の算出結果に適用してみます。
chainer-saliencyについて
以下の論文(1)をchainerで実装したもので、記事(2)にて紹介されています。記事(2)が公開された時点では、chainer-saliencyのレポジトリにありましたが、現在はchainer-chemistry内に入っています。
Saliency計算と可視化
chainer-saliencyにあるexampleのjupyter notebookを参考にし、一部書き方を修正しながら進めます。
notebookでは GraphConvPredictor から定義してますが、chainer-chemistryにはモデルを簡単に定義出来るchainer_chemistry.models.prediction.set_up_predictorが入ってますのでこれを使います。データセットはexampleと同じdelaneyの溶解度データセットです。
# パラメータはexampleのそのまま
n_unit = 32
conv_layers = 4
class_num = 1
predictor = set_up_predictor('ggnn', n_unit=n_unit, conv_layers=conv_layers,
class_num=class_num, postprocess_fn=activation_relu_dropout)
regressor = Regressor(predictor, device=device)
# fitはexampleのまま
fit(regressor, train, batchsize=16, epoch=40, device=device)
学習したモデルを用いてsaliencyを計算します。chainer-chemistry にはVanillaGrad, SmoothGrad, BayesGradの3種類 (文献1,2参照)が実装され、example notebookで紹介されています。
calculator = IntegratedGradientsCalculator(
predictor, steps=5, eval_fun=eval_fun, target_extractor=VariableMonitorLinkHook(predictor.graph_conv.embed, timing='post'),
device=device)
M = 5 # サンプリング数
saliency_samples_vanilla = calculator.compute(
train, M=1, converter=concat_mols)
saliency_samples_smooth = calculator.compute(
train, M=M, converter=concat_mols, noise_sampler=GaussianNoiseSampler())
saliency_samples_bayes = calculator.compute(
train, M=M, converter=concat_mols, train=True)
得られたsaliencyはBaceCalculatorクラスのaggregateで集計します。
method = 'raw'
saliency_agg_vanilla = calculator.aggregate(
saliency_samples_vanilla, ch_axis=3, method=method)
print(f'Before aggregation: {saliency_samples_vanilla.shape} -> After: {saliency_agg_vanilla.shape}')
# Before aggregation: (1, 902, 55, 32) -> After: (902, 55)
saliency_agg_smooth = calculator.aggregate(
saliency_samples_smooth, ch_axis=3, method=method)
print(f'Before aggregation: {saliency_samples_smooth.shape} -> After: {saliency_agg_smooth.shape}')
# Before aggregation: (5, 902, 55, 32) -> After: (902, 55)
saliency_agg_bayes = calculator.aggregate(
saliency_samples_bayes, ch_axis=3, method=method)
print(f'Before aggregation: {saliency_samples_bayes.shape} -> After: {saliency_agg_bayes.shape}')
# Before aggregation: (5, 902, 55, 32) -> After: (902, 55)
これらsaliency_agg_xxxには学習に用いた902化合物ごとに原子の重みが入っており、各化合物のインデックスにアクセスすることで重み/saliencyが取り出せます。このsaliencyを用いてプロットする関数は、これまで用いていた関数を少し変更してます。
重みのscalingにrdkit.Chem.Draw.SimilarityMaps.GetStandardizedWeightsを用いていましたが、chainer-chemistryにも幾つかscalingのメソッドが入っているためこれを指定できる形式にしました。
def plot_mol_saliency(mol, saliency, scaler=None, atoms=['C', 'N', 'O', 'S', 'F', 'Cl', 'P', 'Br']):
n_atom = mol.GetNumAtoms()
symbols = [f'{mol.GetAtomWithIdx(i).GetSymbol()}_{i}' for i in range(mol.GetNumAtoms())]
df = pd.DataFrame(columns=atoms)
saliency = saliency[:n_atom]
num_atoms = mol.GetNumAtoms()
if scaler is not None:
vmax = np.max(np.abs(saliency))
weights = scaler(saliency)
else:
weights, vmax = SimilarityMaps.GetStandardizedWeights(saliency)
arr = np.zeros((num_atoms, len(atoms)))
for i in range(mol.GetNumAtoms()):
_a = mol.GetAtomWithIdx(i).GetSymbol()
arr[i,atoms.index(_a)] = weights[i]
df = pd.DataFrame(arr, index=symbols, columns=atoms)
# 以下同じ
試しに、ID699の化合物を図示してみます。
i = 699 mol = Chem.MolFromSmiles(smiles[i]) plot_mol_saliency(mol, saliency=saliency_agg_bayes[i], scaler=abs_max_scaler, atoms=['C', 'N', 'O', 'S', 'F', 'Cl', 'P', 'Br'])
元々化合物の描画の際のカラーマップ等はchainer-chemistryから拝借していたのでexampleとほぼ同じ描画が出来ました。(今回の記事では補助図は除いてます。)
スケーラーによる描画の違い
前回の記事までは化合物の描画の際、化合物ごとに原子の重み/saliencyはrdkitのGetStandardizedWeightsを用いて、絶対値最大の値でスケーリングしていました。chainer-chemistryでもabs_max_scalerで同じスケーリングを施しています。一方、chainer-chemistryにはabs_max_scaler以外にも、最小値と最大値を用いて0-1にスケーリングするmin_max_scalerが実装されています。(他にもnormalize_scalerがありますが、負の値を含まないsaliencyのためのメソッドなので今回は割愛します。)
min_max_scalerを使って先のID699を可視化してみます。
i = 699 mol = Chem.MolFromSmiles(smiles[i]) plot_mol_saliency(mol, saliency=saliency_agg_bayes[i], scaler=min_max_scaler, atoms=['C', 'N', 'O', 'S', 'F', 'Cl', 'P', 'Br'])
溶解度にプラスに寄与しているとされる原子が濃く描画されているものの、色の濃淡のみで分かりづらくなってしまいました。原子ごとの色はchainer_chemistry.saliency.visualizer.visualizer_utils.red_blue_cmapで決定していましたが、これは元々-1~1にノーマライズされた重み/saliencyに対応するために定義されたものであるため、min_max_scalerとの相性は良くないようです。
では、abs_max_scalerが最適なのか。これを考えるために新しくscalerを考えます。与えた重み/saliencyの絶対値最大ではなく、指定した値でスケーリングするfixed_val_scalerを定義しました。
def _fixed_val_scaler(saliency, max_val, logger=None):
try:
xp = cuda.get_array_module(saliency)
except:
xp = saliency
maxv = max_val
if maxv <= 0:
logger = logger or getLogger(__name__)
logger.info('All saliency value is 0')
return xp.zeros_like(saliency)
else:
return saliency / maxv
def fixed_val_scaler(max_val):
return functools.partial(_fixed_val_scaler, max_val=max_val)
np.unravel_index(np.argmax(np.abs(saliency_agg_bayes)), saliency_agg_bayes.shape), np.max(np.abs(saliency_agg_bayes)) # ((455, 10), 1.1827562)
ID455の化合物がデータセット中で最大の重み (1.1827562)を有することがわかりました。この値をfixed_val_scalerのmax_valとして指定します。ID455の場合は、当然ですがabs_max_scalerとfixed_val_scalerで違いはありません。
次に重み/saliencyの絶対値最大の値が小さい化合物で描画してみます。下記の通り、ID762が2番目に小さい重み (0.114224076)を有するようです。(一番小さい化合物はメタンだったので例には不適でした)
cpid_2nd_min = np.argsort(np.max(np.abs(saliency_agg_bayes), axis=1))[1] cpid_2nd_min, np.max(np.abs(saliency_agg_bayes[cpid_2nd_min])) # (762, 0.114224076)
fixed_val_scalerを用いた方が色が薄くなっています。
このように、abs_max_scalerを用いた場合は、一つの化合物"内"の原子の寄与度を比較する場合には良いですが、化合物"間"で比較する場合には注意が必要そうです。一方、fixed_val_scalerで単純にデータセット内の重み最大値を用いると、各化合物の描画が分かりづらくなってしまいます。下記の通り、個々の化合物の重み/saliencyの和は目的変数とよく相関しており、実際の目的変数/予測値を併記したり、目的変数/予測値で重みを補正する等の対処も可能かもしれません。
scipy.stats.pearsonr(y_train, np.sum(saliency_agg_bayes, axis=1)) # (0.8418738520866434, 2.617349583427535e-243)
おわりに
今回はchainer-chemistryで計算できるsaliencyを元に、スケーリングによる違いを試してみました。説明可能性は機械学習において注目されており、化合物関連でも色々な手法が出てきています。一方、下記のTJOさんの記事でも紹介されていましたが、そもそも機械学習のモデルの説明可能性に対して批判的な意見も出てきているようです。
自分としては、これら出てきた手法をそのまま実務での活用にあたってはまだまだ研究・考察が必要であると言うのが個人的な感想です。単純にケミストのintuitionとの比較だけで良し悪しを論じる前に幾つかのポイントを検証していく必要があると思います。ひとまず思いつくところを書いてみましたが、もっとあるかもしれません。ご意見等いただければ幸いです。
データセット (学習データセットの中に、予測したい化合物の部分構造は十分に含まれているのか。予測したい化合物は適用範囲内か。)
説明可能性の算出アルゴリズム
描画方法 (カラーマップ、重み/saliencyのスケーリング方法、など)
今回のコードはGistにアップしておきました。
新規性が高いけれど合成が難しい部分構造に今後どのように向き合っていくか
遅くなりましたが創薬 (wet) Advent Calendar 2019 第20日目の記事です。とりとめのない雑文になってしまってますがご容赦下さい。
医薬品開発において、化合物の新規性は非常に重要なトピックです。メディシナルケミストは創薬プロジェクトの中で主活性、選択性、薬物動態、等々複数のパラメータを最適化しつつ、さらに特許性を考えながら構造最適化を進めていかなければなりません。
新規性獲得のためのシンプルなアイデアの一つは、珍しい骨格を導入したリガンドを合成することです。では、医薬品においていわゆる"環(ring system)", "骨格 (scaffold)"はどの程度知られているのでしょうか。ここに着目してシステマティックな解析を報告した論文が2009年にPittらによって報告されました (文献1)。この論文ではバーチャルにヒュッケル則を満たす約2.5万のヘテロ環を発生させ、合成可能性の評価、SciFinder/Beilsteinといったデータベースでのチェックを経て、新規性の高い有望なヘテロ環として22種まで絞り込んでます。(なお、文献1にて発生させた約2.5万のヘテロ環はChEMBLからvirtual exploratory heterocyclic library 'VEHICLe'としてダウンロード可能 )
この論文はその後度々リファーされてきましたが、それから10年を経て、改めてこの22種のヘテロ環がどこまで実際に人類によって合成されてきたのかを振り返る論文が発表されました (文献2)。
下図がPittらによって選抜された22のヘテロ環です。この中で9種 (P1, P5–P9, P12, P13, P18)についてこの10年間で実際に合成され、一部は創薬の中で中間体として、あるいは最終検体として利用されました。
文献2では個々の合成例にフォーカスして紹介されています。実際のところ、個別に見ると、誘導化可能な合成経路でない場合も多く、実際に創薬へ利活用するにはまだまだなものも多いのは事実です。しかし、論文で指摘されている通り、Pittらが未踏の骨格を解析・提案したことが医薬/有機合成化学界の興味を惹起したと言う意味で意義深い報告であり、その事後解析であると思います。
では、未だ合成されていないヘテロ環は今後も合成されていくのでしょうか。ケミストはこれらの骨格が新規性が高いことを頭に入れ、積極的に取り組んでいくべきでしょうか。生産性が極めて強く求められている時代です。「そもそも既存の骨格だけでも十分にデザインできる。こんな素性の知れないものをわざわざ合成しようとは思わない」なんて意見もあるでしょうし、それはたしかにその通りだと思います。しかし、新たなヘテロ環を実際に合成・評価することで、新たな知見が得られてきたのも事実です。
もしかしたら、「AIが人の思いつかないような構造をいい感じに出してくれるはずだから、人は出てきた構造を評価・合成すれば良い」と言う意見もあるかもしれません。しかし、以下の観点から、現時点ではその期待通りに事が進むのは難しいのではないかと思うのです。
- 深層生成モデルベース/ルールベースの構造生成の手法は、いずれも基本的には既知の化合物やその部分構造を元にモデル/ルールを作成する方法が多く、今回のような未知のヘテロ環を生成しにくい(と個人的に思ってる)
- 生成したバーチャル構造の評価/スコアリングに「合成難易度」を取り入れていることが多く、今回のような未知のヘテロ環のように合成例が少ない/無い部分構造は低く評価されやすい。
- 新規性が高い部分構造を有する分子の評価を機械学習手法のみで正しく予測することは難しい
- そもそも今の構造生成手法では、いわゆる「特許性」を考慮できていない
2つ目について少し評価してみます。RDKitには合成難易度を評価するSAスコアと言う手法が実装されています。詳細は日本語の紹介記事があるのでそちらを挙げるに留めますが、既知化合物を元にした部分構造の出現頻度と、分子の複雑度を元に評価する手法です。このSAスコアは簡単に計算できるため、深層生成モデルで生成した化合物のスコアリングにも頻出しています。(文献3(a))
比較対象として簡単な分子のスコアも載せます (スコアは1(易)->10(難))。ベンゼンやインドール、アスピリンなどは1-2の範囲です。
一方、P1-P22のスコアは(当然ですが) 2.5-4.5程度と上記よりも高くなっています。
実例も紹介します。Merckが2008年に報告したDAAO阻害剤では単環/2環を有する種々のヘテロ環を有するカルボン酸の変換が検討されており、その中でヘテロ環P13が導入されています (文献4, 構造式はChEMBLより取得)。 上段が文献中のリガンドで最もSAスコアが高い4化合物(SA=3.185-3.083)、下段が最もSAスコアが低い4化合物 (SA=2.384-1.786)です。P13は上段の3番目に入っており、必ずしも他のヘテロ環誘導体よりも抜きん出てはないものの、高いスコアを示す傾向にあることは見て取れます。
SAスコア自体は、大まかな合成難易度の傾向を見る場合には良いものの細かな優先順位付けに使うには注意が必要な手法です。大量の分子を発生させる構造生成手法においてブラックボックス的に用いられると、意図せずに有望な候補分子の優先度を落としてしまうことになりかねません。
今回は部分構造に着目した論文を、簡単な計算と一緒に紹介しました。前述のように、生産性への期待と圧力が高まる昨今では、難しい合成となるデザインをわざわざ考え、チャレンジするのは大変かもしれません。デザインの部分を構造生成技術に頼るとしても、技術の特徴をしっかり把握して活用しなければ期待とは異なる結果に繋がり得ます。生産性ばかりを目指し、新しい化学の裾野を広げていかなければ創薬の土台も広がりません。生産性向上の先に、困難な合成化学に取り組む時間やアイデアを生み出すことにも繋げていければと思うところです。
Heteroaromatic Rings of the Future. Pitt, W. R.; Parry, D. M.; Perry, B. G.; Groom, C. R. J. Med. Chem. 2009, 52, 2952.; ブログでも紹介されている (まだ手つかずの残されたヘテロ環)
‘Heteroaromatic Rings of the Future’: Exploration of Unconquered Chemical Space. Passador, K.; Thorimbert, S.; Botuha, C. Synthesis 2019; 51(02): 384.
(a) RDKitで合成難易度を評価して化合物をスクリーニング; (b)構造生成手法の評価手法全般については Molecular Sets (MOSES): A Benchmarking Platform for Molecular Generation Models.などを参照
The discovery of fused pyrrole carboxylic acids as novel, potent d-amino acid oxidase (DAO) inhibitors. Sparey, T. et al. Bioorg. Med. Chem. Lett. 2008, 18, 3386.
化合物の可視化を機械学習に適用してみた
はじめに
創薬 (dry) Advent Calendar 2019 第15日目の記事です。
前回は化合物を構成する原子ごとの重みを元に、原子ハイライト+ヒートマップ/棒グラフの補助図で可視化する方法を紹介しました。本記事ではその手法を機械学習に適用してみようと思います。
前回の記事では以下の形式で算出した weights を用いれば同様に可視化できるように定義しました。
weights = weight_fn(mol)
本記事ではこの形式に合わせて、以下の3つの手法で重みを計算します。
rdkit.Chem.Draw.SimilarityMaps.GetAtomicWeightsForModel- Feature Importance
- SHAP values
[各手法の概説]
GetAtomicWeightsForModel
RDKitのCookbook (Using scikit-learn with RDKit)でも紹介されている GetSimilarityMapForModelの中で読まれているメソッドです。以下の通り、Fingerprint全体 baseFPを用いて算出した予測値baseProbaと、原子ごとで作成したnewFPを用いて算出したnewProbaの差を以って各原子の重みとしています。
# https://github.com/rdkit/rdkit/blob/master/rdkit/Chem/Draw/SimilarityMaps.py より
def GetAtomicWeightsForModel(probeMol, fpFunction, predictionFunction):
"""
Calculates the atomic weights for the probe molecule based on
a fingerprint function and the prediction function of a ML model.
Parameters:
probeMol -- the probe molecule
fpFunction -- the fingerprint function
predictionFunction -- the prediction function of the ML model
"""
if hasattr(probeMol, '_fpInfo'):
delattr(probeMol, '_fpInfo')
probeFP = fpFunction(probeMol, -1)
baseProba = predictionFunction(probeFP)
# loop over atoms
weights = []
for atomId in range(probeMol.GetNumAtoms()):
newFP = fpFunction(probeMol, atomId)
newProba = predictionFunction(newFP)
weights.append(baseProba - newProba)
if hasattr(probeMol, '_fpInfo'):
delattr(probeMol, '_fpInfo')
return weights
Feature Importance
決定木系アルゴリズムにおいて、各特徴量が判別/回帰にどれぐらい寄与したかを示す指標です。改めて解説するよりはググれば解説記事が見つかるのでそちらを挙げるに留めます。
- <Scikit-learn> ランダムフォレスト回帰のfeature_importances_の定義
- Xgboost : feature_importanceのimportance_type算出方法
- 特徴量重要度にバイアスが生じる状況ご存知ですか？
SHAP（SHapley Additive exPlanations) values
2017年末に出てきた手法ですが、同時に使いやすいライブラリーが出てきたこともあって急速に利活用が進んでいる手法です。ケムインフォマティクスの分野でも早速活用・解析した論文が出ています。よく「ゲーム理論に基づいて...」という説明がなされますが、考え方としてはシンプルに、ある特徴量がない場合の予測値の差分をその特徴量の重要度として計算しているものです。こちらも詳細は先達の記事へ。
- A Unified Approach to Interpreting Model Predictions (元論文)
- 機械学習モデルを解釈する指標SHAPについて (定式の抜粋と利用例も書いてあった)
- slundberg/shap (公式レポ、この通りに動かしたら感覚分かる)
- Interpretation of Compound Activity Predictions from Complex Machine Learning Models Using Local Approximations and Shapley Values (Fingerprint使って得られたSHAPの解析を色々実施)
- ADMET Evaluation in Drug Discovery. 19. Reliable Prediction of Human Cytochrome P450 Inhibition Using Artificial Intelligence Approaches (記述子でSHAPを見ている)
- 機械学習モデルの判断根拠の説明 (SHAPも含めた他の手法が非常に多く取り上げられている)
溶解度データセットを用いた機械学習に適用
まずはRDKitの溶解度データからランダムフォレストでモデルを作ります。
# ref. https://iwatobipen.wordpress.com/2018/11/07/visualize-important-features-of-machine-leaning-rdkit/
sol_classes = {'(A) low': 0, '(B) medium': 1, '(C) high': 2}
X_train = np.array([mol2fp(m)[0] for m in train_mols])
y_train = np.array([sol_classes[m.GetProp('SOL_classification')] for m in train_mols], dtype=np.int)
X_test = np.array([mol2fp(m)[0] for m in test_mols])
y_test = np.array([sol_classes[m.GetProp('SOL_classification')] for m in test_mols], dtype=np.int)
# len(X_train) => 1025, len(X_test) -> 257
clf = RandomForestClassifier(random_state=123)
clf.fit(X_train, y_train)
print(accuracy_score(y_test, clf.predict(X_test)))
# 0.7120622568093385
まずまずのモデルが出来ました。それではGetAtomicWeightsForModelから試してみましょう。最初に幾つか関数を定義します。
# 目的のクラスのprobabilityを返す関数
def get_proba(fp, proba_fn, class_id):
return proba_fn((fp,))[0][class_id]
# デフォルトでGetMorganFingerprintの引数が2048ビットになっているため
def fp_partial(nBits):
return functools.partial(SimilarityMaps.GetMorganFingerprint, nBits=nBits)
# モデルの予測と正解を出力
def show_pred_results(mol, model):
y_pred = model.predict(mol2fp(mol)[0].reshape((1,-1)))
sol_dict = {val: key for key, val in sol_classes.items()}
print(f"True: {mol.GetProp('SOL_classification')} vs Predicted: {sol_dict[y_pred[0]]}")
これらを用いて、以下のように使います。(前回定義したplot_explainable_images を少し変えて、weightsを指定できるようにしています。)
mol = test_mols[141] show_pred_results(mol, clf) # Class 2 (High)に対する値を用いるためget_probaの class_id=2を指定 weights = SimilarityMaps.GetAtomicWeightsForModel(mol, fp_partial(1024), lambda x: get_proba(x, clf.predict_proba, 2)) plot_explainable_images(mol, weight_fn=None, weights=weights, atoms=['C', 'N', 'O', 'S', 'F', 'Cl', 'P', 'Br'])
クラスの予測は外れてますが、全体としてヘテロ原子が赤に、sp2炭素が青になっており、何となくイメージに合っていそうです。
こちらはカルボン酸とメチレン鎖で色分けできていてよりイメージ通りの結果になりました。
続いて、Feature Importanceです。各bitのFeature Importanceはclf.feature_importances_で取得できます。一方、以下で取得できる分子のbitinfoには、各ビットにおける(原子ID, radius)がタプルで入っています。
# https://iwatobipen.wordpress.com/2018/11/07/visualize-important-features-of-machine-leaning-rdkit/
def mol2fp(mol,radius=2, nBits=1024):
bitInfo={}
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=radius, nBits=nBits, bitInfo=bitInfo)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
return arr, bitInfo
fp, bitinfo = mol2fp(test_mols[141])
>>>bitinfo
{0: ((4, 2),),
33: ((0, 0), (5, 0), (10, 0), (6, 2)),
121: ((0, 1), (5, 1), (10, 1)),
179: ((7, 2),),
234: ((1, 2),),
283: ((12, 2),),
314: ((3, 1), (13, 1)),
330: ((11, 2),),
356: ((2, 0), (6, 0), (11, 0), (12, 0)),
378: ((7, 0),),
385: ((9, 1),),
400: ((2, 2),),
416: ((11, 1),),
428: ((7, 1),),
463: ((9, 2),),
493: ((8, 2),),
504: ((12, 1),),
564: ((4, 1), (1, 1)),
650: ((3, 0), (13, 0)),
672: ((6, 1),),
771: ((2, 1),),
849: ((8, 0),),
932: ((8, 1),),
935: ((1, 0), (4, 0), (9, 0))}
各bitのimportance値を、そのbitの中心の原子の値として各原子の重みを計算しました。上記の例のようにbit collisionがあり、その扱いには考慮が必要ですが、今回は単純に各bitに入っている原子数で平均しました。
def weights_from_feat_imp(mol, feature_importance, collided_bits='mean'):
fp, bitinfo = mol2fp(mol)
weights = np.zeros(mol.GetNumAtoms(), )
for bit, infos in bitinfo.items():
for atom_infos in infos:
if collided_bits == 'mean':
# 同じbitに入っている原子数で平均
weights[atom_infos[0]] += feature_importance[bit]/ len(infos)
else:
weights[atom_infos[0]] += feature_importance[bit]
return weights
この方法で計算した重みを可視化したのが下図です。コハク酸のカルボン酸が強くハイライトされているのはイメージに合いますが、全体的に分かりづらいですね。Feature Importanceが非負値であるため、色の濃淡のみでハイライトするこの手法とは相性が悪いようです。Feature Importanceについては、特に重要なbitに絞って(※)可視化した方が視認性が上がるかも知れません。 (※ Visualize important features of machine leaning #RDKit )
最後にSHAPです。公式のレポジトリーを参考にして、以下のようにSHAP値を算出します。後は同様に可視化します。
explainer = shap.TreeExplainer(clf, data=X_train, feature_dependence="independent") shap_values = explainer.shap_values(fp.reshape((1,-1)))[class_id]
SHAPでは正負両方の値が出てくるため、GetAtomicWeightsForModelと同様に感覚的に捉えやすい結果になったと思います。
前回紹介した可視化手法を機械学習と組み合わせ、GetAtomicWeightsForModel / Feature Importance / SHAPという3つの重み計算方法を試してみました。bitinfoを活用すれば原子ごとの重みの計算はできましたが、bit collisionやradiusの考慮等、まだ検討すべき点はありそうです。より大きいサイズのデータセットやモデルチューニングによるパフォーマンス向上の効果も見ておきたいところです。今回は時間が足りず可視化するところまでで終わりましたが、改めて時間を作って追加実験したいと思います。
(2019/12/17 07:30, コードをnbviewerに変更)
(2020/1/13, nbviewer再アップロード)
少し拡張した化合物の可視化を試してみた
はじめに
創薬 (dry) Advent Calendar 2019 第13日目の記事です。
最近では機械学習の結果の解釈性に注目が集まっています。その中でケムインフォマティクスの世界でもモデルの予測解釈性を上げるために、判断の根拠となる部分構造を可視化する手法が検討されています。
[過去の検討例]
- RDKitの類似度マップを用いて原子ごとの寄与を可視化する
- NNの予測根拠可視化をライブラリ化する
- Visualize atom weight of AttentiveFP #DGL #RDKit #Chemoinformatics
- RDKitでフィンガープリントの可視化
- Using the new fingerprint bit rendering code
- Visualize important features of machine leaning #RDKit
- #6 RDKitのDrawMorganBitsはFCFPやunhashed FPでもできるのか確認した
可視化の方向性は大きく2つあります。(A) 1つの化合物に対して重要度の高い原子/判断の根拠となった原子をハイライトする (記事1-3); (B) 重要なfingerprintのbitを抜き出して描画する (記事4-7)。各化合物ごとにどの部分構造が寄与しているかみたい時には(A)の手法が良いのですが、ハイライトした構造式だけではやや定性的過ぎて個人的にはもう少し情報が欲しいなぁと思ってました。
そんな中、構造式だけではなく、原子ごとの重みを別途可視化した図を使っている論文/レポが最近出ていたのでその手法を使って簡単な可視化を試してみました。
[1905.13686] Explainability Techniques for Graph Convolutional Networks
元論文はGCNで予測したモデルの解釈性に関する手法を提案したものですが、詳細は本記事の主題ではないので元論文をご参照頂ければと思います。
個人的に気に入ったのは原子/結合ごとの重みを可視化したヒートマップ+棒グラフで、構造式中で原子をハイライトしただけでは分かりづらい色合いの微妙な違いも分かりやすくなっています。
脂溶性の寄与度可視化への適用
RDkitのドキュメントにも載っているlogPを事例として試しました。元論文とは異なり結合の重みは算出されてないので原子の重みだけを使って少し体裁を変えています。
まず、各原子の重みはドキュメント同様に以下で計算できます。
def calc_crippen_weight(mol):
contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
weights = [x for x,y in contribs]
return weights
ヒートマップに用いる重みのデータは、以下のようにデータフレームとして準備しました。(numpy.arrayのままでも良かった)
symbols = [f'{mol.GetAtomWithIdx(i).GetSymbol()}_{i}' for i in range(mol.GetNumAtoms())]
df = pd.DataFrame(columns=atoms)
contribs = weight_fn(mol) # calc_crippen_weight
num_atoms = mol.GetNumAtoms()
arr = np.zeros((num_atoms, len(atoms)))
for i in range(mol.GetNumAtoms()):
_a = mol.GetAtomWithIdx(i).GetSymbol()
arr[i,atoms.index(_a)] = contribs[i]
df = pd.DataFrame(arr, index=symbols, columns=atoms)
構造式は Chem.Draw.rdMolDraw2D.MolDraw2DCairo で描画します。描画の際の色は原子ごとの重みから計算します。今回は chainer-chemistryで使われているものを拝借しました。
# from chainer-chemistryの visualizer_utils.pyより
def red_blue_cmap(x):
"""Red to Blue color map
Args:
x (float): value between -1 ~ 1, represents normalized saliency score
Returns (tuple): tuple of 3 float values representing R, G, B.
"""
if x > 0:
# Red for positive value
# x=0 -> 1, 1, 1 (white)
# x=1 -> 1, 0, 0 (red)
return 1.0, 1.0 - x, 1.0 - x
else:
# Blue for negative value
x *= -1
return 1.0 - x, 1.0 - x, 1.0
def is_visible(begin, end):
if begin <= 0 or end <= 0:
return 0
elif begin >= 1 or end >= 1:
return 1
else:
return (begin + end) * 0.5
def color_bond(bond, saliency, color_fn):
begin = saliency[bond.GetBeginAtomIdx()]
end = saliency[bond.GetEndAtomIdx()]
return color_fn(is_visible(begin, end))
# 重みを絶対値最大の値でスケーリング
weights, vmax = SimilarityMaps.GetStandardizedWeights(contribs)
vmin = -vmax
# ハイライトする色を算出
atom_colors = {i: red_blue_cmap(e) for i, e in enumerate(weights)}
bondlist = [bond.GetIdx() for bond in mol.GetBonds()]
bond_colors = {i: color_bond(bond, weights, red_blue_cmap) for i, bond in enumerate(mol.GetBonds())}
算出したものを用い、描画します。
drawer = rdMolDraw2D.MolDraw2DCairo(molSize[0],molSize[1])
drawer.drawOptions().useBWAtomPalette()
drawer.drawOptions().padding = .2
drawer.DrawMolecule(
mc,
highlightAtoms=[i for i in range(len(atom_colors))], # 上で計算した原子の色
highlightAtomColors=atom_colors,
highlightBonds=[i for i in range(len(bond_colors))],
highlightBondColors=bond_colors,
highlightAtomRadii={i: .5 for i in range(len(atom_colors))}
)
drawer.FinishDrawing()
これで材料は揃ったので後はレポジトリーにあるコードを適宜変更して描画していくだけです。
RDKitのドキュメントにあるSimilarityMapsを用いた寄与度可視化 (比較用)
今回の手法を適用した可視化です。ヒートマップは左軸に原子IDが振ってあるのですが、構造式中のどれに対応するか分かりづらいので原子IDを埋め込んだ図も合わせて出力するようにしました。
...原子IDを埋め込んでもまだ見づらいですね。このあたりは改善の余地がありそうです。それでも構造式だけよりは寄与度の大小などが見やすくなったのではないでしょうか。
logP以外にもTPSAやGasteiger Chargeでも同様に描画してみました。コードは以下にアップしています。
今回はレポジトリーの手法を原子ごとの重みが計算できる簡単な手法に適用してみました。次は機械学習モデルの結果解釈に本手法を使ってみようと思います。
(2019/12/14 23:00, 誤字修正)
(2019/12/17 07:30, コードをnbviewerに変更)
PD-1/PD-L1タンパクタンパク相互作用阻害剤の特許を調べてみた
この記事は 創薬 Advent Calendar 2017 - Adventar (#souyakuac2017 hashtag on Twitter)の第23日目の記事です．
背景
今年になって抗体以外のPD-1阻害剤のアプローチの特許が徐々に増えてきている様子なので，フリーで可能な範囲で少し調べてみた経緯を紹介します．元々のきっかけは先月末偶然見つけたIncyte社の特許US20170320875A1についてのやりとりです．
PD-1 (Programmed cell death-1)と言えば今や製薬業界では知らぬものはいないほど有名なターゲットで，PD-1またはそのリガンドPD-L1の抗体医薬は画期的新薬として期待と注目を集めています．なのですが，まず非常に高価であることから医療経済面の問題が指摘されています．PD-1抗体の場合，その画期的新薬加算に加えて抗体であるがゆえの高い製造原価が薬価に繋がっています (参考: ポエム: 新薬の薬価算定 – y__sama – Medium)．もう一つ，服薬のし易さから考えると経口剤が望ましいのですが，抗体を含む高分子は経口剤としての開発は難しいと言われています．ではどうするか？最もシンプルなアイデアは「低分子でPD-1阻害剤を見つける」です．しかしこれも言うは易しで，タンパクタンパク相互作用を低分子で狙うのは一般的に難易度が高いことが知られています (もちろん幾つかの成功例はあります)．抗体医薬のニボルマブ (オプジーボ)を持つブリストル・マイヤーズスクイブ (BMS社)も，製造承認を得た2014頃から非抗体の化合物としてマクロサイクル系化合物 (WO2014151634)や低分子化合物 (WO2015034820)を特許開示していました．しかし，某社との共同研究の成果で開発パイプラインに載った話や特許もどんどん出願している等の状況から，どちらかと言うとマクロサイクル系に力を入れているように (筆者には)見え，やはり「低分子は難しいのかな？」と感じていました．そのような背景で見つけたIncyte社特許．見たところBMS社の低分子化合物より分子量の低下に成功している様子です．では，実際のところIncyte社はいつぐらいから研究を始めているのか？BMS社の状況は？他に参入している会社はないのか？以上を特許から探ってみました．
(BMS-M1: WO2014151634, Ex1173, Mw=1880, IC50=2.7 nM [HTRF])
(BMS-1: WO2015034820, Ex202, Mw=419, IC50=18 nM [HTRF])
(BMS-2: WO2015160641, Ex1305, Mw=603, IC50=0.92 nM [HTRF])
(US20170320875A1 [Incyte-5, WO2017192961のファミリー特許], Ex2, Mw=306, IC50<100 nM [HTRF])
特許調査の流れ
以下の手順でpatentscopeを用いてPCT特許を調査しました．
1. BMS社のマクロサイクル系を含む非抗体PD-1化合物を検索
2. Incyte社の化合物を検索
3. 上記2社以外の化合物を検索
細かい手順は割愛しますが，やはり3.の検索の際，フルテキストで"PD-1"を検索すると不要な特許が大量にヒットしてしまいます．そこで今回はヒット数が多くなった場合はFront Page内に"PD-1"を含む特許に限定し，その他，Front Pageにおいて "combination"や "biomarker", "imaging"と言った別の技術を指す単語もANDNOTで除き，BMS社が化合物開示を始めた2014年以降に限定する等で50件以内に絞りました．後はタイトル，Front Pageと最後に本文を見ながらチェックしていきました．
調査結果
| 特許番号 | ID | 出願人 | 代表実施例 | 活性値_nM | Mw | 化合物数 | 発明者数 | 開示日 | 1st_Priority_Data |
|---|---|---|---|---|---|---|---|---|---|
| WO2014151634 | BMS-M1 | BMS | 1173 | 2.714 | 1880 | >900 | 16 | 2014/9/25 | 2013/3/15 |
| WO2015034820 | BMS-1 | BMS | 202 | 18 | 419 | ~300 | 2 | 2015/3/12 | 2013/9/4 |
| WO2015160641 | BMS-2 | BMS | 1305 | 0.92 | 603 | ~400 | 12 | 2015/10/22 | 2014/4/14 |
| WO2016039749 | BMS-M2 | BMS | 10091 | 2.5 | 935 | >1000 | 17 | 2016/3/17 | 2014/9/11 |
| WO2016057624 | BMS-M3 | BMS | 9007 | 5 | 970 | ~200 | 5 | 2016/4/14 | 2015/10/10 |
| WO2016077518 | BMS-M4 | BMS | 11179 | 3 | 1348 | ~400 | 12 | 2016/5/19 | 2014/11/14 |
| WO2016100285 | BMS-M5 | BMS | 3011 | 6.12 | 939 | 51 | 7 | 2016/6/23 | 2014/12/18 |
| WO2016100608 | BMS-M6 | BMS | 9002 | 4.13 | 953 | 31 | 7 | 2016/6/23 | 2014/12/19 |
| WO2016126646 | BMS-M7 | BMS | 1055 | 8 | 906 | 48 | 10 | 2016/8/11 | 2015/2/4 |
| WO2016149351 | BMS-M8 | BMS | 13093 | 8.17 | 956 | 40 | 6 | 2016/9/22 | 2015/3/18 |
| WO2017066227 | BMS-3 | BMS | 1058 | 0.48 | 752 | ~700 | 15 | 2017/4/20 | 2015/10/15 |
| WO2017070089 | Incyte-1 | Incyte | 5 | <100 | 359 | 27 | 7 | 2017/4/27 | 2015/10/19 |
| WO2017087777 | Incyte-2 | Incyte | 8 | <10 | 358 | 40 | 3 | 2017/5/26 | 2015/11/19 |
| WO2017106634 | Incyte-3 | Incyte | 8 | <10 | 372 | 49 | 4 | 2017/6/22 | 2015/11/17 |
| WO2017112730 | Incyte-4 | Incyte | 2 | <10 | 389 | 21 | 3 | 2017/6/29 | 2015/12/22 |
| WO2017118762 | RG_1 | フローニンゲン大 | 60 | <1-1000000 | 572 | 60 | 1 | 2017/7/13 | 2016/1/8 |
| WO2017151830 | BMS-M9 | BMS | 1403 | 7 | 1919 | 30 | 10 | 2017/9/8 | 2016/3/4 |
| WO2017176608 | BMS-M10 | BMS | 1001 | 4 | 981 | 38 | 10 | 2017/10/12 | 2016/4/5 |
| WO2017192961 | Incyte-5 | Incyte | 2 | <100 | 306 | 40 | 3 | 2017/11/9 | 2016/5/6 |
| WO2017202273 | CAMC-1 | 中国医学科学院 | 20 | 2.48 | 621 | 28 | 10 | 2017/11/30 | 2016/5/23 |
| WO2017202274 | CAMC-2 | 中国医学科学院 | 4 | 0.0001 | 598 | 26 | 9 | 2017/11/30 | 2016/5/23 |
| WO2017202275 | CAMC-3 | 中国医学科学院 | 4 | 0.00001 | 647 | 16 | 10 | 2017/11/30 | 2016/5/23 |
| WO2017202276 | CAMC-4 | 中国医学科学院 | 8 | 10.2 | 575 | 17 | 11 | 2017/11/30 | 2016/5/23 |
| WO2017205464 | Incyte-6 | Incyte | 1 | <100 | 405 | 55 | 4 | 2017/11/30 | 2016/5/26 |
...思ったより大変でした．代表実施例番号，活性値，分子量は活性値記載の中で最も高活性のものを選出してます．本当は大体の分子量など物性値の分布も示せれば良かったのですが，時間が足りませんでした (SureChemblからデータは落としたのでいつかやりたい)．化合物数は生物活性テーブル記載の化合物から算出してますが，一部実施例番号が不連続な記載のものがあり概算値として見て下さい．
さて，結果を見るとやはりBMS社が世界に先駆けて研究を進めており ，低分子化合物だけ見ても最も多くの化合物を開示しています (マクロサイクル系BMS-M1~BMS-M10の10報 + 低分子BMS-1~BMS-3の3報)．最初の頃は分子量500を切る化合物も合成してます (BMS-1)が，その後を見ると比較的大きめの分子に注力しているようです (BMS-3)．特許の報数だけ見ると低分子よりマクロサイクル系に力を入れているように見えますが，化合物数を見ると低分子もしっかり合成していますね．印象だけで判断せず実際に調べることが大切でした．BMS社に続いてIncyte社 (Incyte-1~Incyte-6の6報)が続いていますが，それ以外に中国医学科学院 (Chinese Academy of Medical Sciences)も参入していることが分かりました (CAMC-1~CAMC-4)．
少し構造式を見てみましょう．
BMS
BMS社は一貫してビフェニル構造に置換アリールオキシメチル基を導入した誘導体を中心にしており，分子量は大きいものの非常に高活性を達成しています。
(BMS-3: WO2017066227, Ex1058, Mw=752, IC50=0.48 nM [HTRF])
Incyte
一方のIncyte社はビフェニル構造は共通しているもののそこに置換する部位が縮合アリール基やアミドを介したアリール基，ピペリジンなどになっています．正確な活性値は記載がないものの1 nMオーダーの化合物も見出しているようです．
構造意外に目を引くのがIncyte社のPriority Dataです．BMS社の最初の低分子特許 (BMS-1)の開示が2015年3月なのに対して，Incyte-1のPriority Dataが2015年10月19日と約半年後で，もしBMS-1を見てスタートと考えると非常にスピード感ある研究体制が整っていることが推察されます．とは言え，Incyte社はPD-1抗体やその併用などの研究開発も進めており，その次の策として独自に低分子探索を既に始めていたとしても不思議ではありません (Incyte社パイプライン)．構造式を見ると参考にしているようにも見えますが，低分子創薬研究では，以前から独自に進めていて構造活性相関 (SAR)がある程度取得できていた状況に他社から特許が公開され，その導入置換基やSARを参考にすることで活性がブーストされて一気に合成研究が加速されることがしばしばあります．この辺りは残念ながら表には出てこない情報ですが，導入されている置換基の傾向を見ることで類推することはできるでしょう．
また，Incyte-1のPriority Dataから程なくBMS-2が開示になっております (2015/10/22)．BMS-1でもそれなりの化合物数が記載されているので少し間を置いてくるかと思いきや，その半年後に続報が出たのですから，これにはIncyte社研究者も驚いたのではないかと想像します．その後Incyte社は短期間に計3報出願していますが，これはBMS-2を見て早急な対策を講じたのかもしれません．
(Incyte-1: WO2017070089, Ex5, Mw=359, IC50<100 nM [HTRF])
(Incyte-2: WO2017087777, Ex8, Mw=358, IC50<10 nM [HTRF])
(Incyte-3: WO2017106634, Ex8, Mw=372, IC50<10 nM [HTRF])
(Incyte-6: WO2017205464, Ex1, Mw=405, IC50<100 nM [HTRF])
非常に興味深いメディシナルケミストリーを展開しているIncyte社ですが，6報の特許で計11人の発明者がいます．その中で6報全てに名を連ねているのはLiangxing WuとWenqing Yaoの二人です．もう少し彼らについて調べてみます．両者ともPD-1の前はLSD1阻害剤 (WO2015123465 など)やFGFR阻害剤 (WO2014007951 など)を担当していたようです．これらのターゲットについてはIncyte社から論文は出ておりませんが，それぞれ開発パイプラインに入っており，両者が貢献している可能性があります．Wenqing Yaoはそれ以前からPI3K-deltaやc-Met，Pim，H4などの特許にも入っており，PI3K-deltaやc-Metでは臨床開発化合物の論文著者にも入っていました(PI3K-delta阻害剤INCB040093 ; c-Met阻害剤 INCB28060)．それぞれのターゲットにおける関与の仕方は分かりませんが，これだけの臨床開発化合物に関与している研究者ですから，豊富な経験でPD-1の研究においても大きな貢献を果たしていそうです．もしこれでPD-1低分子阻害剤でも臨床まで辿り着いたとしたらもの凄いです．
中国医学科学院
中国医学科学院はBMS同様アリールオキシメチル基が置換したビフェニル構造の化合物ですが，4報ともブロモ基を固定しています．下記には最も高活性でin vivoの試験例も記載されていた化合物を記載しておりますが，HTRFのbinding活性が10-13 Mとナノどころかピコオーダーも超えていると記載されてます．試験系が違うため単純にBMS社やIncyte社との比較はできませんが，これだけ活性が強いと四の五の言わずまずは実際に合成して確かめてみたくなってきます．
(CAMC-2: WO2017202274, Ex4, Mw=598, IC50<10-13 M [HTRF])
(フローニンゲン大学は活性値情報がほとんどなかったため割愛)
実際に調べてみて
BMS社のみが特許を開示していた時は，「薬理的には興味深いターゲットだけど低分子ではまだ難しそう」と思っていましたがIncyte社の化合物を見ている限り，置換基や骨格の許容度もあり十分にチャレンジできそうです．もちろん低分子化すればそれで全てが解決するほど生命科学は単純ではなく，タイムラインでご指摘があったように，逆に選択性が落ちてオフターゲット作用に繋がる可能性もあります．In vitroで十分な活性を示してもIn vivoの薬効になぜか繋がらないこともしばしばあります．しかしこれだけ注目を集めるターゲットですから，近い将来に何らかの研究結果が表に出てくるはずです．今後参入してくる会社の可能性も含めてこれから数年間の動向が注目される状況と言えるでしょう．
創薬化学と特許の関わりについて
この記事は 創薬 Advent Calendar 2017 - Adventar (#souyakuac2017 hashtag on Twitter）の第17日目の記事です．本日から2回に分けて創薬，特にメディシナルケミストと特許について概説したいと思います (知財部の専門ではないので用語の正確性に欠ける点はあるかと思いますがご容赦下さい)．
医薬品開発における特許の重要性や他業種と比べた特殊性については今年良記事が出てましたのでそのリンクを紹介させて頂きます．下記リンクでは特許が製薬企業の売上にどれだけ重要か，特許の種類 (物質特許・用途特許・製剤特許・製法特許について)，特許成立に必要な要件 (産業利用可能性・新規性・進歩性・先願)について紹介されています．ぜひご一読下さい．
医薬品の特許戦略、その重要性とは？『医薬品クライシス』著者・佐藤健太郎氏が語る、医薬研究職の世界 Vol.6
特許を制する者が医薬品を制する 『医薬品クライシス』著者・佐藤健太郎氏が語る、医薬研究職の世界 Vol.7
この記事では，
- 化学者が特許に関わるタイミング
- 特許明細書の構成
- フリーの特許情報源
- 特許を見て分かること (の私見)
について記載します．
化学者が特許の関わるタイミング
研究開始前や創薬初期 (スクリーニングヒット段階)
先行研究調査で文献調査や他社(製薬会社だけでなく大学等も含む)の臨床開発情報と合わせて調査します．自身が狙っている創薬ターゲットで既に化合物が特許公開されていた場合は，その化合物数，構造，生物活性，特許数や出願している製薬企業数，特許公開時期，特許の成立状況を見て，そのターゲットに自社としてどのように取り組むべきか検討します．
合成研究中
化合物のin vitro/in vivo活性やADMEToxを指標に最適化研究を進めますが，特許が取れなければ事業性は一気に低下します．もし首尾よく臨床開発に進め得る化合物プロファイルが得られたとしても公知特許のクレームの範囲に含まれ化合物の新規性が主張できない場合は新規性を獲得できるような構造変換が必要です．メディシナルケミストリー論文では活性やADMEToxに焦点を当てた記述に終始することが多く，新規性の観点からの構造変換は語られないため，製薬企業外から見るとピンと来ないところかもしれません．
他社の化合物は単に眺めて終わりではありません．特許が出ていたらベンチマークとなる化合物を自社でも合成・評価して，プロファイルを知る必要があります．しかしどの化合物を合成すべきか？はそれほど単純な作業ではありません．特許中にはともすれば何百化合物も記載があります．生物活性値 (in vitroすら)記載がないことも多々あり，論文のように最適化経緯が書いてあるわけでもありません．
他社特許中の重要化合物 (key compounds)は，特許明細書の(1)クレーム，(2) 本文 (Description)の活性値や合成量などを見て推定します．(1) クレームでは最初の第1クレームで広くクレームし，その後第2クレーム，第3クレームとどんどん絞り込んでいくように記載されています．そのためより後半のクレーム記載に該当する化合物が重要であるという考え方です．(2) in vivoやin vitro二次アッセイが記載されている時は比較的選択は容易です．一方，in vitroのbiochemicalな活性値のみの場合，どのような目的で重要化合物を推定したいかによってその選び方は変わってきます．単純にin vitro活性を見たいだけなら活性値を基準に選べば良いですが，in vivo活性も検証したいならば薬物動態も考慮して化合物を選ぶ必要があります．実際には薬物動態データの記載がない場合が多いため，脂溶性や極性表面積，分子量などの物性値だけでなく，最後は化学者の経験から選ぶことになります．合成量も重要です．もしある化合物が活性・毒性面で有望なプロファイルを示した場合，その周辺化合物を集中的に合成しているはずで，もしある中間体もしくは最終物が他の化合物よりも多く合成されていた場合はその周辺は有望であろうとの考え方です．有望な化合物周辺は集中的に合成されているだろうとの仮説を元に計算的に重要化合物をピックアップする手法も報告されています (Hattori, K. et al. J. Chem. Inf. Model. 2008, 48 (1), pp 135–142 ; 日本語の記事もあります)．
特許出願とその後
合成研究が進んだら特許出願ですが，出願時期は様々な要素を考えた上で決定されます．出願が遅ければ他社に先んじられる可能性がありますし，一方早すぎれば特許が切れるまでの期間を短くしてしまうことになります．完全に合成研究が終了してから出願することも可能ですし，出願後も1年間は新たな化合物の記載を追加することが可能なため．その期間を見越して早めに出願することも可能です．出願後も合成研究が続く場合は，他社のみならず自社がクレームした範囲も考慮しながら合成を進めていかなければなりません．
特許の構成
多くの特許はPCT出願後 ，1年以内の実施例追加期間を経て，出願から1.5年後に国際知的所有期間 (WIPO)のデータベースPatentscope にて公開されます．日本では毎週木曜日の夕方から金曜日朝が公開タイミングです．PCT出願とはPCT加盟国(日米欧各国や中国も加盟しています)での特許出願手続きを改善するための制度です．特許は公開された時点ではまだ権利化されておらず，同じWO特許の内容で各国の特許法に基づいて審査され，特許が認められると登録特許となります．同じWO特許から派生した各国特許をパテントファミリーと呼びます．また公開した時点では特許番号末尾にAが付きますが，登録特許になるとUSの場合B2が付きます．このように末尾アルファベットで特許の権利化状況が分かります (このような特許の分類を特許種別と言います)．
メルク社のDPP-IV阻害剤sitagliptinを例にすると以下の3つはいずれも同じファミリーの特許です：
- WOの公開特許 WO/2003/004498A1
- USの公開特許 US20030100563A1
- USで権利化された登録特許 US6699871B2
WO特許は「Front Page」「Description」「Claim」「Drawings」「Search Report」から構成されます．
Front Page
特許番号や公開日，発明者などの特許基本情報です．主な箇所だけ取り上げます．
(Front Pageの例, WO2003/004498A1より)
- Pub. No.: 特許番号
- Publication Date: 特許公開日
- IPC: 国際特許分類として，特許がどの分野に関するものかアノテーションしたものです．低分子医薬品特許の場合，多くはA61K (A61K 医薬用，歯科用又は化粧用製剤 )かC07D (C07D 複素環式化合物 )が付与されていることが多いです．テキスト検索でヒット数が多い時の絞り込みなどの参考になります．
- Applicants，Inventors: 出願人および発明者
- Priority Data: 優先権主張日です．各国移行の際にこの日付から権利を主張できます．公開日の1.5年前となっていることが確認できます．
Description
特許本文です．まず特許の背景情報が記載され，次に細かく特許の詳細説明 (クレームする化合物や用語の定義)，製造法 (仮想の合成方法)，実際の合成例 (参考例，実施例)と続き，最後に生物活性評価方法 (とデータ)が記載されています．実際は詳細説明や製造法がかなり長々と記載されているため，前から読むよりは後の生物活性から見ていくと特許の概要が掴みやすいです．
疾患名やターゲットタンパク名などが広く記載されている場合が多いため，単純なテキスト検索やマイニングの対象にする場合は注意が必要です．また、WIPOの[Description]タブから見れるテキストはOCRのため，化合物名が間違いが散見されます．
Claims
特許が主張する内容です．先述の通り最初の第1クレームでは広い範囲でクレームされ，そこから徐々に絞り込んで記載されていきます．また，特許でクレームされる化合物の一般式をマーカッシュ構造 (マルクーシュ構造)と呼びます．
(マーカッシュ構造の例, WO2003/004498A1より)
Drawings
特許内にグラフや図表はこの項目にまとめてある場合があります．
Search Report
WO特許明細書の最後に付いていて，該当特許と関連する先行技術 (論文，特許)とその関連の内容が簡単にまとめてあるため，特許性の参考になります．
情報源
商用特許データベースを使っている方も多いと思いますが，ここではそれらは取り上げずフリーのサイトを幾つか挙げておきます．
PCT Patentscope
WO特許やPCT加盟国の特許が収載されており，特許明細書PDFも入手できます．また最近機械翻訳に力を入れており，英語以外の言語で公開されている特許でもGoogle Chromeの翻訳ツールバーのような機能が使えます (但し，日本語や韓国語など一部の言語以外はOCR)．最近アカウント登録すると構造式検索機能ができるようになりました． 但しやはりOCRであるがゆえの精度の低さや，中間体などの本文中の全てのテキストを対象に拾ってしまうので今時点での性能は期待できません．もし構造式から検索するならばSureChembl の方が使い勝手が良いです．
Espanet
欧州特許庁のサイトです．検索結果の特許のファミリー特許検索が可能で，先行している他社特許の各国移行状況が分かります．
(WO2003/004498の検索結果．左の"patent family"をクリックすると各国のファミリー特許一覧を見ることができる)
Google Patents
US特許に対応してます．WO特許より先にUS特許が公開されることもあるため，Google Scholarなどでアラートをセットしておくと即時性の高い情報が収集でき便利です．WIPOのテキストよりキレイなテキストで検索できるのも大きなメリットです．
(WO2003/004498 [上]と対応するUS20030100563A1[下]の化合物名)
特許を見て分かること (の私見)
上述の通り特許は論文と異なり構造変換の説明も無ければ細かな生物データも無いことが多いです．しかし論文とは異なり企業のビジネスと直結しているため，特許を他の情報と組み合わせて見ていくとその企業の研究戦略や実態がおぼろげながら浮かび上がってきます．
一つの企業から出てくる同一ターゲットの特許を経時的に追っていけばそのターゲットにかけている化学者の人数が推定できます．あるメンバーの名前を経時的に見ていった時，その化学者が別のターゲットの特許で出てきたら，もしかしたら前のターゲットの合成研究に何かの区切りがついたのかもしれません．もしその前のターゲットがいつまでもパイプラインに出てこなかったら,それは経営方針の変更かもしれませんし，前に進められない何かネガティブな要因が見つかったのかもしれません (開発に見合う化合物を創出できなかったのかもしれないし，ターゲット由来の副作用が回避できなかったのかもしれません)．また，特許開示より前にメディシナルケミストリーの論文が出ていたとしたら，その論文の化合物は特許で保護するだけのものではなかったとも考えられます．
他社の特許と見比べていくと，ある会社の化合物を見て別の会社が合成研究を始めたのではないかと思われるような事例も見つかります．そのタイミングやキャッチアップの速さでその企業のマネジメントのスピード感も想像できます．
これらは全て憶測で何が正解かは分かりません．それでも特許を他の情報源と組み合わせて他社の研究動向やその研究者の思考過程に思いを巡らせたり，周囲の方と色々議論してみると色々面白い発見や意見が出てくるはずです．
23日には実際にPD-1の特許を調査して考察した内容について紹介したいと思います．
(更新履歴: 2017/12/23 - 誤字修正)
