TeachOpenCADDで学ぶ創薬ケモインフォマティクス入門（前半）

はじめに

前回までの溶解度予測シリーズでは，分子の物性予測を学びました．
今回は創薬の世界に足を踏み入れてみます．

がんの薬物標的として有名な EGFR（上皮成長因子受容体）の阻害剤を題材に，創薬計算化学（CADD）の基礎を学んでいきます．

教材として，ベルリン Charité 大学の Volkamer 研究室が公開しているTeachOpenCADD（https://projects.volkamerlab.org/teachopencadd/） を参考にしています．

T001〜T035 まで 35 以上のトピックがカバーされた，創薬計算化学（CADD）の体系的なチュートリアル集です．
本シリーズでは，その中から T001–T006 および T022 を取り上げ，データ取得からフィルタリング，類似度検索，QSAR モデル構築までを学びます．
AIにこれを元に学習コースを作成してもらいながら進めました．
創薬に関しては素人なので，誤った知識等あったらご指摘いただけますとありがたいです．

今回は以下を行います．

• ChEMBL データベースから EGFR 阻害剤の活性データを取得
• 活性指標 pIC50 の意味と分布の確認
• Lipinski の Rule of Five による薬物らしさのフィルタリング
• PAINS フィルタによる偽陽性化合物の除外
• フィンガープリントを使った類似度検索

コードに入る前に，今回扱う題材の背景を簡単に整理します．

EGFR（上皮成長因子受容体：Epidermal Growth Factor Receptor）
細胞の表面にあるタンパク質で細胞の増殖や生存をコントロールします．
正常な細胞では成長因子が結合すると活性化し、細胞の分化や増殖、組織の修復を促します．
しかし，がん細胞では EGFR が異常に活性化し，細胞が止まらずに増殖してしまいます．
EGFR 阻害剤はこの機能をブロックする薬で，ゲフィチニブやエルロチニブが有名です．

ChEMBL
欧州バイオインフォマティクス研究所が管理する世界最大の薬物活性データベースです．
数百万の化合物とその生物活性データ（IC50 など）が登録されており，無料で使えます．

IC50（阻害濃度50%）
標的タンパク質の活性を 50% 阻害するのに必要な薬の濃度です．
IC50 が小さいほど，少量で効く強い薬ということになります．

pIC50
IC50 は値が小さいほど良い薬ですが，値の範囲が広い（0.003 nM 〜 数万 nM）ので扱いにくい．
そこで溶解度の logS と同じ発想で対数を取ります．pIC50 が大きいほど強い阻害剤です．

IC50	pIC50	意味
0.003 nM	11.5	非常に強い
1 nM	9.0	強い
1000 nM	6.0	弱い
100000 nM	4.0	ほぼ効かない

ライブラリのインポート

TeachOpenCADD の GitHub から直接 CSV を取得しています．

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit import Chem
from rdkit.Chem import Draw, Descriptors

データの取得と確認

TeachOpenCADD の GitHub から EGFR 阻害剤のデータを取得します．

# TeachOpenCADD の GitHub から EGFR のデータを直接取得
url = "https://raw.githubusercontent.com/volkamerlab/teachopencadd/master/teachopencadd/talktorials/T001_query_chembl/data/EGFR_compounds.csv"

df = pd.read_csv(url)

実際のデータの中身を確認します．

df.head()

実際に使っていくのはsmiles列とpIC50列だけになりそうです．

基本情報を整理します．

print(f"=== データの概要 ===")
print(f"化合物数: {df.shape[0]}")
print(f"列数: {df.shape[1]}")
print(f"\n=== 列名と意味 ===")
print(df.columns.tolist())
print(f"\n=== 欠損値 ===")
print(df.isnull().sum())
print(f"\n=== 統計量 ===")
print(df[["IC50", "pIC50"]].describe())

5568 化合物のデータがあります．
統計量から，IC50では10-3~10+7の幅広い値を取るのに対して，対数をとったpIC50は1.6~11.5と扱いやすい範囲になっていることがわかります．

EDA

pIC50 の分布

まず目的変数である pIC50 の分布を確認します．

plt.figure(figsize=(8, 5))
plt.hist(df["pIC50"], bins=40, edgecolor="black")
plt.xlabel("pIC50")
plt.ylabel("Count")
plt.title(f"Distribution of pIC50 (n={len(df)})")
plt.show()

概ね正規分布に近い分布です．対数スケールなので追加の変換は不要です．

強い/弱い阻害剤の構造

強い阻害剤 Top 5 と弱い阻害剤 Top 5 の構造を見てみます．

# 強い阻害剤 Top5 の構造
top5 = df.nlargest(5, "pIC50")
print("=== 強い阻害剤 Top 5 ===")
for _, row in top5.iterrows():
    print(f"  {row['molecule_chembl_id']}: pIC50 = {row['pIC50']:.2f}, IC50 = {row['IC50']:.3f} nM")

mols_top = [Chem.MolFromSmiles(smi) for smi in top5["smiles"]]
legends_top = [f"{row['molecule_chembl_id']}\npIC50={row['pIC50']:.1f}"
               for _, row in top5.iterrows()]
img = Draw.MolsToGridImage(mols_top, molsPerRow=5, subImgSize=(250, 200),
                            legends=legends_top)
display(img)

強い阻害剤の構造を見ると，2つ以上の縮環含窒素芳香環を持つものが多く，ピリミジン骨格が目立ちます．
これらの構造的特徴は，後のフィンガープリントやMCS（Maximum Common Substructure）解析で定量的に裏付けることができそうです．

# 弱い阻害剤 Top5 も
bottom5 = df.nsmallest(5, "pIC50")
print("\n=== 弱い阻害剤 Top 5 ===")
for _, row in bottom5.iterrows():
    print(f"  {row['molecule_chembl_id']}: pIC50 = {row['pIC50']:.2f}, IC50 = {row['IC50']:.1f} nM")

mols_bottom = [Chem.MolFromSmiles(smi) for smi in bottom5["smiles"]]
legends_bottom = [f"{row['molecule_chembl_id']}\npIC50={row['pIC50']:.1f}"
                  for _, row in bottom5.iterrows()]
img2 = Draw.MolsToGridImage(mols_bottom, molsPerRow=5, subImgSize=(250, 200),
                             legends=legends_bottom)
display(img2)

一方，弱い阻害剤にはフェノールを部分骨格に持つものが多そうです．
AI曰く「フェノール OH は極性が高く，EGFR タンパク質との疎水的な結合ポケットには不利なのかもしれません」とのことです．

Lipinski の Rule of Five

ここからは，取得した化合物を「薬物らしさ」でフィルタリングします．
Lipinski の Rule of Five は，1997年に Pfizer の Christopher Lipinski が提唱した経口投与で効く薬の経験則です．

以下の4条件を全て満たす分子は，経口吸収されやすいとされています．
1. 分子量 ≤ 500
2. LogP ≤ 5
3. 水素結合ドナー数 ≤ 5
4. 水素結合アクセプター数 ≤ 10

化学的に解釈するとこんな感じです．

• 分子量が大きすぎると細胞膜を通れない
• LogP が高すぎると水に溶けず吸収されない
• 水素結合が多すぎると細胞膜（脂質二重層）を通過しにくい

溶解度予測で使った RDKit 記述子がそのまま登場するので，同じ要領で計算します．
まず SMILES から Mol オブジェクトを作成し，変換に失敗した分子がないか確認します．

from rdkit.Chem import Descriptors

df["mol"] = df["smiles"].apply(Chem.MolFromSmiles)
df = df[df["mol"].notnull()].reset_index(drop=True) # True の行だけ残す（None の行を削除）
print(f"有効な化合物数: {len(df)}")

無事に全化合物について変換できています．

4記述子の分布と基準値

次に Lipinski の4記述子を計算し，各記述子の分布と基準値を可視化します．
赤い点線が Lipinski の基準値で，その左側にある分子が条件を満たしています．

# Lipinski の4記述子を計算
df["MolWt"] = df["mol"].apply(Descriptors.MolWt)
df["LogP"] = df["mol"].apply(Descriptors.MolLogP)
df["HBD"] = df["mol"].apply(Descriptors.NumHDonors)
df["HBA"] = df["mol"].apply(Descriptors.NumHAcceptors)

# 各記述子の分布と Lipinski 基準を可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
descriptors = [("MolWt", 500), ("LogP", 5), ("HBD", 5), ("HBA", 10)]

for ax, (col, threshold) in zip(axes.ravel(), descriptors):
    ax.hist(df[col], bins=40, edgecolor="black")
    ax.axvline(x=threshold, color="red", linestyle="--", linewidth=2)
    n_pass = (df[col] <= threshold).sum()
    ax.set_title(f"{col} (pass: {n_pass}/{len(df)}, {n_pass/len(df)*100:.1f}%)")
    ax.set_xlabel(col)
    ax.set_ylabel("Count")
    ax.legend()

plt.tight_layout()
plt.show()

MolWtとLogPは基準以上の分子が多い一方，HBDとHBAはほぼ全分子が通過しています．
EGFR 阻害剤は縮環芳香族が多いとさきほどTop5で確認しましたが，それらはLipinskiルールにより経口吸収しにくいので，そのバランスを実現するのが難しさになりそうです．

実際にLipinskiを違反する薬も存在するそうですが，今回は学習目的なので単純にLipinskiでフィルタリングをかけることにします．

フィルタリング結果

# Lipinski の全条件を満たすかどうか
df["Lipinski_pass"] = (
    (df["MolWt"] <= 500) &
    (df["LogP"] <= 5) &
    (df["HBD"] <= 5) &
    (df["HBA"] <= 10)
)

n_pass = df["Lipinski_pass"].sum()
n_fail = (~df["Lipinski_pass"]).sum()
print(f"通過: {n_pass}件 ({n_pass/len(df)*100:.1f}%)")
print(f"不通過: {n_fail}件 ({n_fail/len(df)*100:.1f}%)")

ヒートマップ（記述子と pIC50 の相関）

次に，「どの記述子が pIC50 と相関が強いか」を確認します．

cols_for_corr = ["MolWt", "LogP", "HBD", "HBA", "pIC50"]
plt.figure(figsize=(8, 6))
sns.heatmap(df[cols_for_corr].corr(), annot=True, cmap="coolwarm", center=0, fmt=".2f")
plt.title("Correlation: Lipinski descriptors vs pIC50")
plt.tight_layout()
plt.show()

全記述子の相関が弱い（最大でもMolWtのr=0.26）ことが分かります．
溶解度予測ではLogPがr=-0.83で圧倒的だったのと対照的です．

薬物活性は「分子全体の物性」だけでは決まらず，「EGFRタンパク質の結合ポケットにどうフィットするか」という立体的な要因が大きいためです．また，Lipinski の4記述子は「薬になりそうか」のフィルタであって，「強い薬かどうか」の予測にはあまり関係ないことが分かります．

記述子同士の相関では，MolWtとHBAが0.68と高いです．
大きい分子ほどアクセプター（O，N）が多い傾向にあるということで，多重共線性の候補です．

散布図（各記述子 vs pIC50）

また，pIC50と各記述子の散布図を作成し，全体的な傾向を把握したいと思います．

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
for ax, col in zip(axes.ravel(), ["MolWt", "LogP", "HBD", "HBA"]):
    ax.scatter(df[col], df["pIC50"], alpha=0.1, s=10)
    ax.set_xlabel(col)
    ax.set_ylabel("pIC50")
    ax.set_title(f"pIC50 vs {col}")
plt.tight_layout()
plt.show()

4つとも明確な傾向が見えず，点が広く散らばっています．
溶解度予測の MolWt vs logS では綺麗な右下がりが見えましたが，今回は r = 0.26 程度の弱い相関なのでこんなものかなと思います．
Lipinskiの4記述子だけではpIC50 を予測するのは難しそうだと分かります．あくまでも経口吸収のしやすさであり，薬としての指標は別であると再認識しました．

KDE プロット（Lipinski Pass vs Fail）

Lipinski を通過した化合物と通過しなかった化合物で，pIC50 の分布に違いがあるかをKDE プロットで確認します．先ほどTop5は多環芳香族が多いことを確認した通り，pIC50が高いがLipinskiは通過しなかった（見落とし）があるかを見てみようと思います．

plt.figure(figsize=(8, 5))
for label, group in df.groupby("Lipinski_pass"):
    name = "Pass" if label else "Fail"
    group["pIC50"].plot.kde(label=name)
plt.xlabel("pIC50")
plt.legend()
plt.show()

Fail（青）の方が pIC50 がやや右（高活性）に分布しています．
「Lipinski に落ちる分子の方が活性が高い」という一見矛盾した結果ですが，
Top5構造で確認した通り，EGFR 阻害剤は縮環芳香族が多く分子量や LogP が高くなりがちで，そういう構造ほど EGFR に強く結合するためだと考えられます．
Lipinski はあくまで「経口吸収できるか」のフィルタであって，「効く薬かどうか」とは別の基準であることがここからも分かります．

PAINS フィルタ

ここからは，PAINS（Pan-Assay INterference compoundS）フィルタを適用し，アッセイで偽陽性を出しやすい部分構造を持つ化合物を取り除きます．
創薬では「活性があるように見えて，実は偽物」を除外するのが重要とのことです．

これらのWebサイトでも勉強させていただきました．

また，PAINS フィルタの仕組みは以下の通りです：

1. 約480個の「問題のある部分構造パターン」が SMARTS 形式で格納されている
2. 各化合物に対して，HasMatch で部分構造一致を判定
3. 1つでもマッチしたら PAINS 該当

from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams

# PAINS フィルタを設定
params = FilterCatalogParams()  # 空のフィルタ設定を作る
params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)  # PAINSパターン集を追加
catalog = FilterCatalog(params)  # フィルタ本体を作る

# 各化合物が PAINS の部分構造パターンにマッチするか判定
df["is_PAINS"] = df["mol"].apply(catalog.HasMatch)
n_pains = df["is_PAINS"].sum()
n_clean = (~df["is_PAINS"]).sum()
print(f"PAINS 該当: {n_pains}件 ({n_pains/len(df)*100:.1f}%)")
print(f"クリーン:   {n_clean}件 ({n_clean/len(df)*100:.1f}%)")

PAINS 該当分子の構造

実際にどんな分子がPAINSに引っかかったのか見てみます．

# PAINS に該当する分子の例を表示
pains_mols = df[df["is_PAINS"]].head(6)
print("=== PAINS 該当分子の例 ===")
for _, row in pains_mols.iterrows():
    print(f"  {row['molecule_chembl_id']}: pIC50 = {row['pIC50']:.2f}")

mols_pains = [Chem.MolFromSmiles(smi) for smi in pains_mols["smiles"]]
legends_pains = [f"{row['molecule_chembl_id']}\npIC50={row['pIC50']:.1f}"
                 for _, row in pains_mols.iterrows()]
img_pains = Draw.MolsToGridImage(mols_pains, molsPerRow=3, subImgSize=(250, 200),
                                  legends=legends_pains)
display(img_pains)

PAINSのどのパターンにマッチしたかも調べてみます．

pains_mols = df[df["is_PAINS"]].head(10)

for _, row in pains_mols.iterrows():
    mol = row["mol"]
    entry = catalog.GetFirstMatch(mol)  # 最初にマッチしたパターンを取得
    if entry:
        print(f"{row['molecule_chembl_id']}: {entry.GetDescription()}")

該当する分子を見ると，以下の構造がPAINSに引っかかっていたことがわかります．

• catechol_A: カテコール部位は金属キレート能があるので，圧制の金属イオンに結合して偽陽性が出る
• anil_di_alk_A: アニリン+ジアルキル置換（今回はピペラジンが直接ベンゼンに結合した部分）は酸化されやすい

PAINSに該当した474件の出現頻度についてもみてみます．

from collections import Counter

pattern_names = []
for _, row in df[df["is_PAINS"]].iterrows():
    entry = catalog.GetFirstMatch(row["mol"])
    if entry:
        pattern_names.append(entry.GetDescription())

counts = Counter(pattern_names)
for name, count in counts.most_common(10):
    print(f"  {name}: {count}件")

これらの構造については，ざっくり以下のような理由でPAINSに登録されています．

• anil_di_alk_A: アニリンのジアルキル置換．酸化されやすく偽陽性を出す
• catechol_A: カテコール（隣接2つのOH基）．金属キレートで偽陽性を出す
• mannich_A: マンニッヒ塩基．求核剤（タンパク質のシステインなど）と共有結合しやすい
• ene_cyano_A: シアノエン（C=C-CN）．マイケル付加でタンパク質と非特異的に反応する
• azo_A: アゾ基（N=N）．代謝で分解されやすく，結果が不安定になる
• ene_one_ene_A: エノン系（C=C-C=O-C=C）．マイケルアクセプターとして非特異的に反応する
• quinone_A: キノン．レドックス活性が高く，酸化還元サイクルで偽陽性を出す
• indol_3yl_alk: インドール3位のアルキル置換．非特異的結合を起こしやすい
• ene_five_het_B: 五員環ヘテロ環に隣接するエン構造．反応性が高い
• hzone_phenol_A: フェノールヒドラゾン．金属キレートと酸化の両方で偽陽性を出す

KDE（PAINS vs Clean）

PAINS に該当する化合物は活性も低いのかを KDE で比較してみます．

# PAINS 該当 vs クリーンで pIC50 の分布を比較
plt.figure(figsize=(8, 5))
for label, name in [(True, "PAINS"), (False, "Clean")]:
    df[df["is_PAINS"] == label]["pIC50"].plot.kde(label=name)
plt.xlabel("pIC50")
plt.title("pIC50 distribution: PAINS vs Clean")
plt.legend()
plt.show()

PAINS（青）はピークが pIC50 ≈ 5.5 で左（低活性）に偏り，Clean（オレンジ）はピークが3つあって右（高活性）寄りです．
Lipinskiでフィルターした時と異なり，これはPAINSは本当は効いていないのに活性があるように見える化合物なので，削れるのは8.5%と少ないですが活性の有無のスクリーニングに効果がありそうです．

Lipinski + PAINS の両方でフィルタリング

これらの両方でフィルタリングをすることで，どこまで化合物を絞れるか試してみます．
「Lipinskiをパス」「PAINSに引っ掛からない」の両方に当てはまるデータフレームを新しく作ります．

# Lipinski + PAINS の両方でフィルタリング
df_filtered = df[(df["Lipinski_pass"]) & (~df["is_PAINS"])].reset_index(drop=True)

print(f"元データ:              {len(df)}件")
print(f"Lipinski 通過:         {df['Lipinski_pass'].sum()}件")
print(f"PAINS クリーン:        {(~df['is_PAINS']).sum()}件")
print(f"両方通過:              {len(df_filtered)}件 ({len(df_filtered)/len(df)*100:.1f}%)")

約半数が落ちました．
Lipinski（特に MolWt と LogP）で45%，PAINSで9%落とせていたので，感覚としても合います．

Lipinski「体に入るか」，PAINS「実験結果が信用できるか」という異なる観点から2つのフィルターをかけることで，元の 5568 化合物から 2797 化合物（50.2%）に絞り，より信頼性の高い，薬物らしいデータセットになったと言えると思います．

ケモインフォマティクスの論文でも，膨大なデータセットの中から，複数の基準（フィルター）を定めて，化合物を絞り込んでいき，そして量子化学計算など別の評価を行なっていくという手順はよくみます．

類似度検索（Tanimoto）

フィルタリング後の 2797 化合物に対して，フィンガープリントによる類似度検索を行います．
前回やった溶解度予測ではフィンガープリントを「ML の特徴量」として使いましたが，ここでは「2分子の類似度を計算するツール」として使います．その判定には，Tanimoto類似度というものを用います．

ゲフィチニブは最も有名な EGFR 阻害剤の一つなので，「この薬に似た化合物を探す」という創薬の典型的なアプローチとして使います．

from rdkit.Chem import AllChem
from rdkit import DataStructs

# 既知の強い EGFR 阻害剤(ゲフィチニブ)を参照分子とする
gefitinib_smiles = "COc1cc2ncnc(Nc3ccc(F)c(Cl)c3)c2cc1OCCCN1CCOCC1"
gefitinib = Chem.MolFromSmiles(gefitinib_smiles)

ゲフィチニブの構造を見てみます．

display(Draw.MolToImage(gefitinib, size=(300, 200)))

Tanimoto 類似度とは

2つのフィンガープリント間の「重なり具合」を 0〜1 で数値化したものです．

Tanimoto = (共通の ON ビット数) / (どちらかが ON のビット数)

• 1.0: 完全に同じ構造
• 0.7 以上: かなり似ている（同じ骨格の類縁体）
• 0.3 以下: ほとんど別物

つまり，フィンガープリントを比較して，共通してONになっている部分が多いほど構造が似ているということです．
ケモインフォマティクスではTanimotoが最も標準的に使われているそうです．

ゲフィチニブのフィンガープリントも作成しておきます．

# 参照分子のフィンガープリント
ref_fp = AllChem.GetMorganFingerprintAsBitVect(gefitinib, radius=2, nBits=2048)

ゲフィチニブとの類似度分布

まず，全化合物とゲフィチニブの類似度を計算し，その分布を作成します．

# 全化合物との Tanimoto 類似度を計算
def calc_tanimoto(mol):
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
    return DataStructs.TanimotoSimilarity(ref_fp, fp)

df_filtered["tanimoto_gefitinib"] = df_filtered["mol"].apply(calc_tanimoto)

# 類似度の分布
plt.figure(figsize=(8, 5))
plt.hist(df_filtered["tanimoto_gefitinib"], bins=50, edgecolor="black")
plt.xlabel("Tanimoto similarity to Gefitinib")
plt.ylabel("Count")
plt.title("Distribution of Tanimoto similarity")
plt.axvline(x=0.7, color="red", linestyle="--", label="similarity = 0.7") #0.7を基準としてみる
plt.legend()
plt.show()

このデータセットはいろんな化合物が含まれているので，類似している分子はかなり少数であることが視覚的にわかります．

類似度 0.7 以上がどれくらいあるか，定量的にも確認してみます．

threshold = 0.7
n_similar = (df_filtered["tanimoto_gefitinib"] >= threshold).sum()
n_total = len(df_filtered)
print(f"Tanimoto ≥ {threshold}: {n_similar}件 ({n_similar/n_total*100:.1f}%)")

類似度 0.7 以上の化合物は 24件（0.9%）と少なく，大半の化合物はゲフィチニブとは異なる骨格を持っていることが分かります．

類似化合物 Top 10

次に，Top10化合物についての情報を見てみます．

# 類似度が高い Top 10
top_similar = df_filtered.nlargest(10, "tanimoto_gefitinib")
for _, row in top_similar.iterrows():
    print(f"  {row['molecule_chembl_id']}: similarity={row['tanimoto_gefitinib']:.3f}, pIC50={row['pIC50']:.2f}")

similarityが1.0となっているCHEMBL939はゲフィチニブ自身で，元からデータセットに含まれていたことが確認できました．Top6に関しては構造も描画します．

# Top 6 の構造を描画
mols_similar = [Chem.MolFromSmiles(smi) for smi in top_similar["smiles"].head(6)]
legends_similar = [f"sim={row['tanimoto_gefitinib']:.2f}\npIC50={row['pIC50']:.1f}"
                   for _, row in top_similar.head(6).iterrows()]
img_similar = Draw.MolsToGridImage(mols_similar, molsPerRow=3, subImgSize=(300, 200),
                                    legends=legends_similar)
display(img_similar)

pIC50だけを見るとゲフィチニブよりも高いものも多いです．実際に薬として使えるかどうかは，他の要素もあり，ゲフィチニブはそれらを総合的にクリアしたものだと考えられます．

類似度と pIC50 の関係

類似度が高い化合物は活性も高いのか，という点に関して散布図と相関係数で確認します．

gef_row = df_filtered[df_filtered["molecule_chembl_id"] == "CHEMBL939"]
gef_pIC50 = gef_row["pIC50"].values[0]

plt.figure(figsize=(8, 6))
plt.scatter(df_filtered["tanimoto_gefitinib"], df_filtered["pIC50"], alpha=0.2, s=10)
plt.scatter(1.0, gef_pIC50, color="red", marker="*", s=300,
            edgecolor="black", zorder=5, label=f"Gefitinib (pIC50={gef_pIC50:.1f})")
plt.axvline(x=0.7, color="red", linestyle="--", alpha=0.5)
plt.xlabel("Tanimoto similarity to Gefitinib")
plt.ylabel("pIC50")
plt.title("Tanimoto similarity vs pIC50")
plt.legend()
plt.show()

# 相関係数
corr = df_filtered["tanimoto_gefitinib"].corr(df_filtered["pIC50"])
print(f"Tanimoto similarity vs pIC50 の相関: r = {corr:+.3f}")

r = 0.469と，「似た構造 → 似た活性」若干の相関はあると言えます．
また，他には以下の点が特徴として言えます．

1. 類似度が低くても高活性の分子がある（左上の点）
Tanimoto = 0.1〜0.2 なのに pIC50 = 9〜10 の化合物が散見されます．
ゲフィチニブとは全く違う骨格で EGFR を阻害している分子だと思われます．

2. 類似度が高くても活性がばらつく（右側の縦の広がり）
Tanimoto = 0.7 以上でも pIC50 = 6〜8 とばらつきがあります．
これは創薬では activity cliff（活性崖）と呼ばれる現象で，構造がほぼ同じでも置換基1つで活性が桁違いに変わるケースです．

activity cliffの例も具体的に見てみます．

# Tanimoto ≥ 0.7 の化合物だけ抽出
similar = df_filtered[df_filtered["tanimoto_gefitinib"] >= 0.7].copy()

# pIC50 が最も高いものと最も低いものを比較
best = similar.nlargest(1, "pIC50").iloc[0]
worst = similar.nsmallest(1, "pIC50").iloc[0]

print(f"最強: {best['molecule_chembl_id']}, pIC50={best['pIC50']:.2f}, sim={best['tanimoto_gefitinib']:.3f}")
print(f"最弱: {worst['molecule_chembl_id']}, pIC50={worst['pIC50']:.2f}, sim={worst['tanimoto_gefitinib']:.3f}")
print(f"pIC50 の差: {best['pIC50'] - worst['pIC50']:.2f}")

# 2つの構造を並べて描画
mols_cliff = [Chem.MolFromSmiles(best["smiles"]), Chem.MolFromSmiles(worst["smiles"])]
legends_cliff = [
    f"Strong\npIC50={best['pIC50']:.1f}\nsim={best['tanimoto_gefitinib']:.2f}",
    f"Weak\npIC50={worst['pIC50']:.1f}\nsim={worst['tanimoto_gefitinib']:.2f}"
]
img = Draw.MolsToGridImage(mols_cliff, molsPerRow=2, subImgSize=(350, 250),
                            legends=legends_cliff)
display(img)

構造は確かに似ていますが，芳香環部分，アミノ基の置換基値，シアノ基の有無，メトキシ基orアミド基，など，結構違うようにも感じます．
tanimoto > 0.7の化合物のうち，もっと構造が近くてpIC50の差が大きな分子を探してみました．

from rdkit.Chem import Draw
from itertools import combinations

similar = df_filtered[df_filtered["tanimoto_gefitinib"] >= 0.7].copy()

best_cliff = {"sim": 0, "diff": 0}
for i, j in combinations(similar.index, 2):
    fp_i = AllChem.GetMorganFingerprintAsBitVect(similar.loc[i, "mol"], radius=2, nBits=2048)
    fp_j = AllChem.GetMorganFingerprintAsBitVect(similar.loc[j, "mol"], radius=2, nBits=2048)
    sim = DataStructs.TanimotoSimilarity(fp_i, fp_j)
    diff = abs(similar.loc[i, "pIC50"] - similar.loc[j, "pIC50"])
    
    if sim >= 0.85 and diff > best_cliff["diff"]:
        best_cliff = {"sim": sim, "diff": diff, "i": i, "j": j}

i, j = best_cliff["i"], best_cliff["j"]
mols_cliff = [similar.loc[i, "mol"], similar.loc[j, "mol"]]
legends_cliff = [
    f"{similar.loc[i, 'molecule_chembl_id']}\npIC50={similar.loc[i, 'pIC50']:.1f}\nsim={best_cliff['sim']:.2f}",
    f"{similar.loc[j, 'molecule_chembl_id']}\npIC50={similar.loc[j, 'pIC50']:.1f}\nsim={best_cliff['sim']:.2f}"
]
img = Draw.MolsToGridImage(mols_cliff, molsPerRow=2, subImgSize=(350, 250),
                            legends=legends_cliff)
display(img)

炭素数がひとつ変わるだけでpIC50の値が8.0から6.3へとかなり変わります．
pIC50 で 1.7 の差は IC50 に換算すると約50倍の差に相当し，わずかな構造変化が活性に大きく影響することが分かります．

溶解度は分子全体の物性で決まるので予測しやすかったのに対し，活性は「タンパク質の結合ポケットにどうフィットするか」という精密な相互作用に依存するため，その予測の方が本質的に難しいということが実感できます．

そのため，創薬での類似度検索は「完璧な予測」ではなく，候補の絞り込みとして使われるそうです．

まとめ

今回は ChEMBL から EGFR 阻害剤 5568 化合物のデータを取得し，Lipinski / PAINS フィルタによるデータクリーニングと，類似度検索を行いました．

フィルタリングの推移

ステップ	化合物数	割合
元データ	5568件	100%
Lipinski 通過	3050件	54.8%
PAINS クリーン	5094件	91.5%
両方通過	2797件	50.2%

Lipinski は「体に入るか」（物性的な薬物らしさ），PAINS は「実験結果が信用できるか」（偽陽性の除外）と，全く異なる観点からデータの質を管理しています．

EDA から分かったこと

• Lipinski の4記述子と pIC50 の相関は全て弱い（最大でも r = 0.26）
• 溶解度予測では LogP が r = -0.83 で圧倒的だったのと対照的
• 薬物活性は分子全体の物性だけでは決まらず，フィンガープリントなど部分構造の情報が重要になると予想される

類似度検索

• ゲフィチニブを参照分子として Tanimoto 類似度で検索
• 「似た構造ほど活性が高い傾向」が見えたが，activity cliff も確認
• 溶解度予測ではフィンガープリントを「ML の特徴量」として使ったが，今回は「類似度計算のツール」として使用

次回は，化学空間の可視化（t-SNE）と最大共通部分構造（MCS）の探索を行い，最後に機械学習によるQSARモデルの構築に進みます．