薬理活性情報が紐付いた化合物データを取得する

本記事では、TeachOpenCADDのチュートリアル（T001）を参考に、公共データベースであるChEMBLからComputer-Aided Drug Design（CADD）の勉強/練習用の化合物データを取得します。

実行環境
環境構築
実践
最後に

実行環境

OS: Ubuntu 22.04.3 LTS
CPU: AMD Ryzen 7 5700X 8-Core Processor
GPU: NVIDIA GeForce RTX 3060 Ti (8G)
RAM: 32GB

環境構築

Pythonのライブラリとして提供されているChEMBL webresource client を利用し、API経由でChEMBLからデータを取得します。

必要なライブラリはconda環境を作成してインストールするのが無難です。今回はpython3.9ベースで”chemble”という名前の仮想環境を作成し各種ライブラリをインストールしました。

conda create -n chemble python=3.9
pip install chembl_webresource_client
pip install rdkit
pip install pandas
pip install matplotlib
conda install jupyter

conda環境の作成にはminicondaあるいはanacondaのインストールが必要です。導入がまだの方はどちらかをインストールしましょう。

minicondaのインストールについては以下の記事にまとめてあります。参考になれば幸いです。

実践

チュートリアルと同様、チロシンキナーゼのEGFRをターゲットとする化合物データを取得します。

チュートリアルでは、EGFRのChEMBL ID取得のためにUniport IDによる検索をかけていますが、今回は割愛してChEMBL ID（ChEMBL203）取得後から活性データを取りに行きます。標的を変える場合は、チュートリアルをご参考ください。

以下、jupyter notebook上での利用を想定したPythonコードです。まずは利用するライブラリをインポートします。

# 必要なライブラリのインポート
import math
import pandas as pd
from rdkit.Chem import PandasTools
from chembl_webresource_client.new_client import new_client

次に、APIのエンドポイントを設定します。APIの詳細については、インフォマティクスが主業務でない方は「こんな風に書けばいいんだなぁ〜」くらいの認識で問題ありません。

# APIのエンドポイントを設定
bioactivities_api = new_client.activity
compounds_api = new_client.molecule

APIを経由して活性データを取得します（インターネットへの接続必須）。ChEMBL IDの変数を自力で格納する部分以外は、チュートリアルのコードそのままですがご了承ください。

# 活性データを取得
chembl_id = "CHEMBL203"
bioactivities = bioactivities_api.filter(
    target_chembl_id=chembl_id, type="IC50", relation="=", assay_type="B"
).only(
    "activity_id",
    "assay_chembl_id",
    "assay_description",
    "assay_type",
    "molecule_chembl_id",
    "type",
    "standard_units",
    "relation",
    "standard_value",
    "target_chembl_id",
    "target_organism",
)

print(f"Length and type of bioactivities object: {len(bioactivities)}, {type(bioactivities)}")

2025年8月16日現在で、11,634件のデータが取得できました。

取得したデータはPythonのDict型として保持されていますが、このままだとデータとして扱いにくいためpandasのDataFrame型へと変換します（時間がかかります。私の環境では20分程度でした）。

# pandasのデータフレームへ変換
bioactivities_df = pd.DataFrame.from_dict(bioactivities)
print(f"DataFrame shape: {bioactivities_df.shape}")
bioactivities_df.head()

変換が完了したらデータをキュレーションします。ChEMBLのデータは専門家によってすでにキュレーションているため、ここでの作業は重複処理やデータ型の変換といった比較的簡単なものになります。

# 不必要なカラム削除
bioactivities_df.drop(["units", "value"], axis=1, inplace=True) 

# データ型の変換
bioactivities_df = bioactivities_df.astype({"standard_value": "float64"}) 

# 欠損値の削除
bioactivities_df.dropna(axis=0, how="any", inplace=True)
print(f"DataFrame shape: {bioactivities_df.shape}")

# 単位がnMのデータを抽出
bioactivities_df = bioactivities_df[bioactivities_df["standard_units"] == "nM"]
print(f"DataFrame shape: {bioactivities_df.shape}")

# 重複を削除
bioactivities_df.drop_duplicates("molecule_chembl_id", keep="first", inplace=True)
print(f"DataFrame shape: {bioactivities_df.shape}")

# インデックスのリセット
bioactivities_df.reset_index(drop=True, inplace=True)

# カラム名の変換
bioactivities_df.rename(
    columns={"standard_value": "IC50", "standard_units": "units"}, inplace=True
)
print(f"DataFrame shape: {bioactivities_df.shape}")
bioactivities_df.head()

11,634件が7,334件になりました。以上で活性データの準備は完了です。

CADDのためには化合物の構造情報も必要になります。構造情報もChEMBLから取得可能なので、次は以下のコードで構造データを取得します（こちらも時間がかかります。私の環境では30分程度でした）。

# 構造データの取得
compounds_provider = compounds_api.filter(
    molecule_chembl_id__in=list(bioactivities_df["molecule_chembl_id"])
).only("molecule_chembl_id", "molecule_structures")
compounds = list(compounds_provider)
compounds_df = pd.DataFrame.from_records(
    compounds,
)
print(f"DataFrame shape: {compounds_df.shape}")
compounds_df.head()

データ取得が完了したらこちらもキュレーションを行います。

構造データは様々な形式（SMILES、InChI、InChIKey、etc…）で登録されていますが、今回はチュートリアル同様にSMILESが登録されているデータのみを残します。

# 構造情報が欠損するデータを削除
compounds_df.dropna(axis=0, how="any", inplace=True)
print(f"DataFrame shape: {compounds_df.shape}")

# 超複する化合物を削除
compounds_df.drop_duplicates("molecule_chembl_id", keep="first", inplace=True)
print(f"DataFrame shape: {compounds_df.shape}")

# SMILESが登録されている化合物を抽出
canonical_smiles = []
for i, compounds in compounds_df.iterrows():
    try:
        canonical_smiles.append(compounds["molecule_structures"]["canonical_smiles"])
    except KeyError:
        canonical_smiles.append(None)
compounds_df["smiles"] = canonical_smiles
compounds_df.drop("molecule_structures", axis=1, inplace=True)
print(f"DataFrame shape: {compounds_df.shape}")

# SMILESが欠損するデータを削除
compounds_df.dropna(axis=0, how="any", inplace=True)
print(f"DataFrame shape: {compounds_df.shape}")

構造データの準備ができたので最初に取得した活性データと統合します。

# データの統合
output_df = pd.merge(
    bioactivities_df[["molecule_chembl_id", "IC50", "units"]],
    compounds_df,
    on="molecule_chembl_id",
)

# インデックスのリセット
output_df.reset_index(drop=True, inplace=True)
print(f"Dataset with {output_df.shape[0]} entries.")

以上で、計7,327化合物のデータセットが取得できました。

最後にSDFとして保存するために、各化合物についてRDkitのMolオブジェクトを生成します。また、活性データは対数化して扱うことが多いため、IC50からpIC50へ変換しておきます。

# Molオブジェクトの生成
PandasTools.AddMoleculeColumnToFrame(output_df, smilesCol="smiles")

# pIC50への変換
def convert_ic50_to_pic50(IC50_value):
    pIC50_value = 9 - math.log10(IC50_value)
    return pIC50_value
output_df["pIC50"] = output_df.apply(lambda x: convert_ic50_to_pic50(x.IC50), axis=1)

# pIC50順に並べ替え
output_df.sort_values(by="pIC50", ascending=False, inplace=True)
output_df.reset_index(drop=True, inplace=True)
output_df.head()

以上で最終的なデータセットができました。こちらを任意のディレクトリにSDFとして保存しましょう。

# _Nameプロパティの変更
for id, mol in zip(output_df['molecule_chembl_id'], output_df['ROMol']):
    mol.SetProp('_Name', id)

# データの保存
dataset_dir = os.path.expanduser("~") +'/cadd_training/dataset'
os.makedirs(dataset_dir, exist_ok=True)
PandasTools.WriteSDF(output_df, dataset_dir + '/dataset_EGFR.sdf', properties=list(output_df.columns))