magattacaのブログ

日付以外誤報

TeachOpenCADD トピック1 〜化合物データの取得 (ChEMBL)~

こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリトークトリアル1をもととしておりCC BY 4.0のライセンスに従っています。
GitHubはてなMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを!

取り上げられている内容の概観のため目次を貼っておきます。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル1

化合物データの取得 (ChEMBL)

Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin

Paula Junge and Svetlana Leng

このトークトリアルの目的

ChEMBLからデータを抽出する方法の学習:

  • ある特定の標的に対して評価済みのリガンドを見つける
  • 取得可能な生理活性データでフィルタリング
  • pIC50値の計算
  • データフレームを結合し、取り出した分子を描画

学習の目標

理論

実践

目標:与えられた標的に対する生理活性データのある化合物のリストを取得すること

  • ChEMBLデータベースへの接続
  • 標的(ターゲット分子)に関するデータの取得 (EGFR kinase)
  • 生理活性データ
    • 生理活性をダウンロードしフィルタリングする
    • データの整型と変換
  • 化合物データ
    • 化合物リストの取得
    • 出力データの準備
  • 出力
    • 最もpIC50値の大きい分子を描画
    • アウトプットファイルの書き出し

レファレンス


理論

ChEMBLデータベース

  • 大規模な生理活性の公開データベース
  • 現在のデータ内容 (10.2018時点):
    • 1.8 百万以上のユニークな構造の化合物
    • 1 百万アッセイに基づく 15 百万以上の活性値
    • ∼12 000 のターゲット分子に紐づけられたアッセイ
  • データソース に含まれているもの:科学文献、PubChemバイオアッセイ、顧みられない病気のための創薬イニシアティブ(Drugs for Neglected Diseases Initiative (DNDi))、BindingDB データベース、などなど
  • ChEMBLのデータは Webインターフェースあるいは EBI-RDFプラットフォームChEMBL Webサービス からアクセスすることができます。

ChEMBL Webサービス

f:id:magattaca:20200411143208j:plain

Figure 1: "ChEMBL webサービス スキーマダイアグラム。楕円型はChEMBL webサービスのリソースを示し、二つのリソースを繋ぐラインは、それらが共通の特徴を共有していることを示します。矢印の向きはリソースタイプについての一次情報を見つけられる場所を示しています。点線は2つのリソースの関係が異なる形式で振舞うことを示しています。例えば、Image リソースは Molecule のグラフィカルな描写を提供します。Figureと説明は次の文献からとりました: Nucleic Acids Res. (2015), 43, 612-620.

ChEMBL webリソースクライアント

  • ChEMBLのデータにアクセスするためのPythonクライアントライブラリ
  • HTTPSプロトコルとのやり取りを行う
  • 結果の遅延評価(lazy evaluation) -> ネットワークリクエストを減らす

化合物活性評価

IC50

Figure 2: IC50値を算出する方法を視覚的に表した図: 阻害を縦軸、log(濃度)を横軸にデータを配置し、阻害の最大と最小を決定する。曲線が50%阻害レベルと交差する点の濃度がIC50

pIC50

  • IC50値の比較を容易にするため、pIC50値をlogスケールで次のように定義しています
    $pIC{50} = -log{10}(IC{50}) $ ($ IC{50}$の単位はM)
  • pIC50が大きくなるほど指数関数的に薬効が強くなっていることを示します
  • pIC50はモル濃度の指標で与えられます (mol/L or M)
    • IC50はM単位に換算したあとで、pIC50に変換する必要があります。
    • nMの場合は: $pIC{50} = -log{10}(IC{50}*10^{-9})= 9-log{10}(IC_{50}) $

IC50と、pIC50に加えて、その他の生理活性の指標として、平衡定数 KI と50%効果濃度(half maximal effective concentration (EC50 ))といったものが使われます。

実践

それでは、我々の興味のあるターゲット分子であるEGFRキナーゼに対して、評価が実施された全ての分子をダウンロードしたいと思います。

ChEMBLデータベースへの接続

まず、ChEMBL webリソースクライアントと他のpythonライブラリーをインポートします。

from chembl_webresource_client.new_client import new_client
import pandas as pd
import math
from rdkit.Chem import PandasTools

訳注(04.2020時点) from chembl_webresource_client.new_client import new_client を実行すると以下のようなでエラーが出ました。

Exception: svg+xml is not an available format (['xml', 'json', 'yaml', 'svg'])

こちらの問題はclientのバージョンを上げることで解決できます。イシュー #71を参照してください。アップデートは

pip install -U chembl-webresource-client

で実行可能でした。 「-U」を忘れずに!
訳注ここまで

APIアクセスのためのリソースオブジェクトを作成します

targets = new_client.target
compounds = new_client.molecule
bioactivities = new_client.activity

ターゲット分子のデータ

  • UniProt Webサイト (https://www.uniprot.org/) から目的とするターゲット分子(EGFRキナーゼ)の UniProt-ID (http://www.uniprot.org/uniprot/P00533) を取得します。
  • UniProt-ID を使用してターゲット分子の情報を取得します。
  • 他のターゲット分子に興味があるときは、ターゲットに合わせて異なるUniProt-IDを選択してください
uniprot_id = 'P00533'
# ChEMBLからターゲット分子の情報を取得しますが、
# 取得する項目を特定のものに絞って取得します。
target_P00533 = targets.get(target_components__accession=uniprot_id) \
                       .only('target_chembl_id', 'organism', 'pref_name', 'target_type')
print(type(target_P00533))
pd.DataFrame.from_records(target_P00533)

# <class 'chembl_webresource_client.query_set.QuerySet'>
cross_references organism pref_name species_group_flag target_chembl_id target_components target_type tax_id
0 [{'xref_id': 'O43451', 'xref_name': None, 'xre... Homo sapiens Maltase-glucoamylase False CHEMBL2074 [{'accession': 'O43451', 'component_descriptio... SINGLE PROTEIN 9606.0
1 NaN Homo sapiens Epidermal growth factor receptor erbB1 NaN CHEMBL203 NaN SINGLE PROTEIN NaN
2 NaN Homo sapiens Epidermal growth factor receptor and ErbB2 (HE... NaN CHEMBL2111431 NaN PROTEIN FAMILY NaN
3 NaN Homo sapiens Epidermal growth factor receptor NaN CHEMBL2363049 NaN PROTEIN FAMILY NaN
4 NaN Homo sapiens MER intracellular domain/EGFR extracellular do... NaN CHEMBL3137284 NaN CHIMERIC PROTEIN NaN

訳注(04.2020時点)
出力したtypeはQuerySetとなっています。Django QuerySetインターフェースに基づいたデザインとなっているそうです。詳しくはchemble_webresource_clientのGitHubを参照してください 出力結果はオリジナルのJupyter Notebookと異なっています。

  • only(~) で項目を絞っているはずですが、データフレームに表示された項目は絞られていません。
  • また、chembl_webresouce_clientのonlyオペレーターの使用方法では引数をリストにすることが推奨されています。
  • index 0 のエントリーが「CHEMBL2074」 となっています。こちらはEGFRではなさそうです。エラーでしょうか?

参考までに、データの取得と絞り込みを2段階に分けた方法を下に示します。データフレームのカラムを選択しているだけです。

target_P00533 = targets.get(target_components__accession=uniprot_id) 
selected_target_values = ['organism', 'pref_name', 'target_chembl_id', 'target_type']
df_trarget_tmp = pd.DataFrame.from_records(target_P00533)[selected_target_values]
df_trarget_tmp
organism pref_name target_chembl_id target_type
0 Homo sapiens Maltase-glucoamylase CHEMBL2074 SINGLE PROTEIN
1 Homo sapiens Epidermal growth factor receptor erbB1 CHEMBL203 SINGLE PROTEIN
2 Homo sapiens Epidermal growth factor receptor and ErbB2 (HE... CHEMBL2111431 PROTEIN FAMILY
3 Homo sapiens Epidermal growth factor receptor CHEMBL2363049 PROTEIN FAMILY
4 Homo sapiens MER intracellular domain/EGFR extracellular do... CHEMBL3137284 CHIMERIC PROTEIN

訳注ここまで

ヒットしたエントリーを確認したのち、最初のエントリーを目的のターゲット分子として選択

CHEMBL203: 単一のタンパク質(single protein)で human Epidermal growth factor receptor (EGFR, 別名 erbB1) を表しています。

target = target_P00533[0]
target
#    {'organism': 'Homo sapiens',
#    'pref_name': 'Epidermal growth factor receptor erbB1',
#    'target_chembl_id': 'CHEMBL203',
#    'target_type': 'SINGLE PROTEIN'}

訳注(04.2020時点)

DataFrameのindex 0 と同一の結果(「CHEMBL2074の出力」)となるかと思いましたが、オリジナルのGitHub のノートブック同様に「CHEMBL203」のが出力されました。QuerySetをDataFrameに変換するタイミングでエラーがあるのでしょうか?

訳注ここまで

選択したChEMBL-IDを保存します。

chembl_id = target['target_chembl_id']
chembl_id
#  'CHEMBL203'

生理活性データ

それでは、目的とするターゲット分子の生理活性データをクエリとして投げたいと思います。

ターゲット分子の生理活性をダウンロードしてフィルタリング

この工程では、生理活性データをダウンロードし、フィルターをかけます。次の項目のみ考慮します。

  • ヒトタンパク質(human proteins)
  • 生理活性タイプ(bioactivity type IC50)
  • 精確な測定結果(関係式が '=' となっているもの)
  • 結合データ (アッセイタイプ 'B')
bioact = bioactivities.filter(target_chembl_id = chembl_id) \
                      .filter(type = 'IC50') \
                      .filter(relation = '=') \
                      .filter(assay_type = 'B') \
                      .only('activity_id','assay_chembl_id', 'assay_description', 'assay_type', \
                            'molecule_chembl_id', 'type', 'units', 'relation', 'value', \
                            'target_chembl_id', 'target_organism')
# len(bioact), len(bioact[0]), type(bioact), type(bioact[0])
len(bioact[0]), type(bioact), type(bioact[0])

#  (43, chembl_webresource_client.query_set.QuerySet, dict)

訳注(04.2020時点)

len(bioact) はNonTypeのため表示できないといったエラーが出ました。コメントアウトし、残りの3つを表示させています。
len(bioact) に相当するデータ(ヒットした化合物の数)は以降のDataFrame化、サイズ確認の段階で確認できます。

訳注ここまで

ChEMBLデータベースにクエリを投げるのが難しければ、先のセルのクエリの結果(11 April 2019)を含むファイルを提供しています。Pythonのパッケージ pickleを使いました。pickelはPythonのオブジェクトをシリアライズし、ファイルに保存、あとでプログラムに読み込めるようにします。(オブジェクトのシリアライズについて、さらに学びたいときには次を参照してください DataCamp )

pickle化した化合物を読み込むには次のセルのコメントアウトを外して実行してください。

#import pickle
#bioact = pickle.load(open("../data/T1/EGFR_compounds_from_chembl_query_20190411.p", "rb"))

生理活性データを整型し変換する

データは辞書のリストととして保存されています。

bioact[0]
{'activity_comment': None,
 'activity_id': 31863,
 'activity_properties': [],
 'assay_chembl_id': 'CHEMBL663853',
 'assay_description': 'Inhibitory concentration against human DNA topoisomerase II, alpha mediated relaxation of pBR322; no measurable activity',
 'assay_type': 'B',
 'bao_endpoint': 'BAO_0000190',
 'bao_format': 'BAO_0000357',
 'bao_label': 'single protein format',
 'canonical_smiles': 'c1ccc(-c2nc3c(-c4nc5ccccc5o4)cccc3o2)cc1',
 'data_validity_comment': None,
 'data_validity_description': None,
 'document_chembl_id': 'CHEMBL1137930',
 'document_journal': 'Bioorg. Med. Chem. Lett.',
 'document_year': 2004,
 'ligand_efficiency': None,
 'molecule_chembl_id': 'CHEMBL113081',
 'molecule_pref_name': None,
 'parent_molecule_chembl_id': 'CHEMBL113081',
 'pchembl_value': None,
 'potential_duplicate': False,
 'qudt_units': 'http://www.openphacts.org/units/Nanomolar',
 'record_id': 206172,
 'relation': '>',
 'src_id': 1,
 'standard_flag': True,
 'standard_relation': '>',
 'standard_text_value': None,
 'standard_type': 'IC50',
 'standard_units': 'nM',
 'standard_upper_value': None,
 'standard_value': '100000.0',
 'target_chembl_id': 'CHEMBL1806',
 'target_organism': 'Homo sapiens',
 'target_pref_name': 'DNA topoisomerase II alpha',
 'target_tax_id': '9606',
 'text_value': None,
 'toid': None,
 'type': 'IC50',
 'units': 'uM',
 'uo_units': 'UO_0000065',
 'upper_value': None,
 'value': '100.0'}

andasのDataFrameに変換します。(数分かかるかもしれません)

訳注(04.2020時点)
上の出力でもクエリでのデータ取得の段階での項目の絞り込みがうまく機能していません。
オリジナルのノートブックでは次にpandas.DataFrame.from_recordsによりDataFrameに変換していますが、項目数が多すぎるためかうまく変換できませんでした。
そこで冗長ですが、以下のセルでは一度csvファイルに出力し、csvをDataFrameに読み込むという作業を行なっています。
csv出力の段階で項目の絞り込みを実施ししています。(クエリのonly の代替)。
また、DataFrame変換後、項目の値の条件によるフィルタリングを念のため再度実施しています。
訳注ここまで

import csv

# 出力したいフィールド (csvのheaderにする)
fields = ['activity_id','assay_chembl_id', 'assay_description', 'assay_type', 'molecule_chembl_id', 'type', 'units', 'relation', 'value', 'target_chembl_id', 'target_organism']

# bioact.csvというファイルを作成し出力する
with open('../data/T1/bioact.csv', 'w') as f:
    writer = csv.writer(f)
    
    # ヘッダーの書き込み
    writer.writerow(fields)
    
    # activityのデータを1行ずつ書き込み
    for act_data in bioact:
        row = []
        
        for field in fields:
            row.append(act_data[field])
        
        writer.writerow(row)
bioact_df_tmp = pd.read_csv('../data/T1/bioact.csv')

bioact_df = bioact_df_tmp.query('type == "IC50" and assay_type == "B" and relation == "="') # 条件に合致する行を残す
bioact_df.head(10)
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 uM = 0.041 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 uM = 0.170 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 uM = 9.300 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 uM = 500.000 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 uM = 3000.000 CHEMBL203 Homo sapiens
5 32780 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL77085 IC50 uM = 96.000 CHEMBL203 Homo sapiens
6 33406 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL443268 IC50 uM = 5.310 CHEMBL203 Homo sapiens
7 34039 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76979 IC50 uM = 264.000 CHEMBL203 Homo sapiens
8 34041 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76589 IC50 uM = 0.125 CHEMBL203 Homo sapiens
9 34049 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL76904 IC50 uM = 35.000 CHEMBL203 Homo sapiens
bioact_df.shape
# (7177, 11)

欠損値のあるエントリーを削除します。

bioact_df = bioact_df.dropna(axis=0, how = 'any')
bioact_df.shape
# (7176, 11)

重複を削除します。同じ化合物(molecule_chembl_id) について複数回の評価データが登録されていることが時々あります。今回は最初の結果だけを採用します。

bioact_df = bioact_df.drop_duplicates('molecule_chembl_id', keep = 'first')
bioact_df.shape
#  (5505, 11)

生理活性値がモル濃度の単位で評価されているものだけを残したいと思います。次のprint文は、どのような単位が含まれているかを確認し、行をいくつか削除したあとでどのようなエントリーが残すか調整するのに役立ちます。

print(bioact_df.units.unique())
bioact_df = bioact_df.drop(bioact_df.index[~bioact_df.units.str.contains('M')])
print(bioact_df.units.unique())
bioact_df.shape

# ['uM' 'nM' 'M' "10'1 ug/ml" 'ug ml-1' "10'-1microM" "10'1 uM"
#   "10'-1 ug/ml" "10'-2 ug/ml" "10'2 uM" '/uM' "10'-6g/ml" 'mM' 'umol/L'
#   'nmol/L']
#  ['uM' 'nM' 'M' "10'-1microM" "10'1 uM" "10'2 uM" '/uM' 'mM']
# (5435, 11)

訳注(04.2020時点)
上のセルの条件 ~bioact_df.units.str.contains('M')~notを意味します。
pandasの条件ではandornotの代わりに&|~を使うことでエラーを回避できるそうです。参考
訳注ここまで

いくつか行を削除しましたが、あとでインデックスにしたがって反復処理(iteration)を行いたいので、インデックスを連続した値にリセットします。

bioact_df = bioact_df.reset_index(drop=True)  # drop=Trueで元の引数を削除
bioact_df.head()
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 uM = 0.041 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 uM = 0.170 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 uM = 9.300 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 uM = 500.000 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 uM = 3000.000 CHEMBL203 Homo sapiens

IC50値をさらに比較可能とするために、全ての単位をnMに変換します。まず、ヘルパー関数を書きます。この関数を次工程でデータフレーム全体に適用します。

def convert_to_NM(unit, bioactivity):
#     c=0
# for i, unit in enumerate(bioact_df.units):
    if unit != "nM":        
        if unit == "pM":
            value = float(bioactivity)/1000
        elif unit == "10'-11M":
            value = float(bioactivity)/100
        elif unit == "10'-10M":
            value = float(bioactivity)/10
        elif unit == "10'-8M":
            value = float(bioactivity)*10
        elif unit == "10'-1microM" or unit == "10'-7M":
            value = float(bioactivity)*100
        elif unit == "uM" or unit == "/uM" or unit == "10'-6M":
            value = float(bioactivity)*1000
        elif unit == "10'1 uM":
            value = float(bioactivity)*10000
        elif unit == "10'2 uM":
            value = float(bioactivity)*100000
        elif unit == "mM":
            value = float(bioactivity)*1000000
        elif unit == "M":
            value = float(bioactivity)*1000000000
        else:
            print ('unit not recognized...', unit)
        return value
    else: return bioactivity
bioactivity_nM = []
for i, row in bioact_df.iterrows():
    bioact_nM = convert_to_NM(row['units'], row['value'])
    bioactivity_nM.append(bioact_nM)
bioact_df['value'] = bioactivity_nM
bioact_df['units'] = 'nM'
bioact_df.head()
activity_id assay_chembl_id assay_description assay_type molecule_chembl_id type units relation value target_chembl_id target_organism
0 32260 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL68920 IC50 nM = 41.0 CHEMBL203 Homo sapiens
1 32267 CHEMBL674637 Inhibitory activity towards tyrosine phosphory... B CHEMBL69960 IC50 nM = 170.0 CHEMBL203 Homo sapiens
2 32680 CHEMBL677833 In vitro inhibition of Epidermal growth factor... B CHEMBL137635 IC50 nM = 9300.0 CHEMBL203 Homo sapiens
3 32770 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL306988 IC50 nM = 500000.0 CHEMBL203 Homo sapiens
4 32772 CHEMBL674643 Inhibitory concentration of EGF dependent auto... B CHEMBL66879 IC50 nM = 3000000.0 CHEMBL203 Homo sapiens

化合物データ

EGFRに対して評価を実施された全ての化合物(とそれぞれの評価方法)を含むデータフレームが手に入りました。それでは、それぞれのChEMBL IDの裏に隠れた化合物を取得したいと思います。

化合物のリストを取得する

ChEMBLから抽出し生理活性データを明確にした化合物を見てみましょう。まず、望みの生理活性データを有する化合物のChEMBL IDと構造を取り出します。

cmpd_id_list = list(bioact_df['molecule_chembl_id'])
compound_list = compounds.filter(molecule_chembl_id__in = cmpd_id_list)\
                              .only('molecule_chembl_id','molecule_structures')

次に、リストをPandasのDataFrameに変換し、重複を削除します。(くどいですが、padas from_records 関数の実行は時間がかかるかもしれません)。

訳注(04.2020時点)

ここでもクエリのデータ取得の段階での項目の絞り込み(only)がうまく機能しておらず、DataFrameへの変換が非常に遅いので、csvの出力ののちDataFrameに変換するという方法をとりました。出力したcsvは17MBあり1時間弱かかりました。

訳注ここまで

# 出力したいフィールド (csvのheaderにする)
cmpd_fields = ['molecule_chembl_id','molecule_structures']

# cmpd.csvというファイルを作成し出力する
with open('../data/T1/cmpd.csv', 'w') as f:
    writer = csv.writer(f)
    
    # ヘッダーの書き込み
    writer.writerow(cmpd_fields)
    
    # activityのデータを1行ずつ書き込み
    for cmpd in compound_list:
        row = []
        
        for field in cmpd_fields:
            row.append(cmpd[field])
        
        writer.writerow(row)
# compound_df = pd.DataFrame.from_records(compound_list)
compound_df = pd.read_csv('../data/T1/cmpd.csv')
compound_df = compound_df.drop_duplicates('molecule_chembl_id', keep = 'first')
print(compound_df.shape)
print(bioact_df.shape)
compound_df.head()
# (5435, 2)
# (5435, 11)
molecule_chembl_id molecule_structures
0 CHEMBL6246 {'canonical_smiles': 'O=c1oc2c(O)c(O)cc3c(=O)o...
1 CHEMBL10 {'canonical_smiles': 'C[S+]([O-])c1ccc(-c2nc(-...
2 CHEMBL6976 {'canonical_smiles': 'COc1cc2c(cc1OC)Nc1ncn(C)...
3 CHEMBL7002 {'canonical_smiles': 'CC1(COc2ccc(CC3SC(=O)NC3...
4 CHEMBL414013 {'canonical_smiles': 'COc1cc2c(cc1OC)Nc1ncnc(O...

この段階では、構造の情報について複数の異なる構造表現を保持しています。この中からカノニカルSMILESだけを残したいと思います。

訳注(04.2020時点)
オリジナルの方法とは異なりcsvからDataFrameにしたことで、カラム「molecule_structures」の要素は見た目はdict型ですがstr型となっています。
このままではcanonical_smilesをvalueとして取り出すことができません。
以下ではastモジュールのliteral_eval関数を適用してstr型からdict型に変換する処理を最初に実行しています。
また、forループで処理する際の、欠損値行のスキップがもともとのif条件ではうまくいかなかったため、float型の場合スキップするという条件に変更しています。
訳注ここまで

# mapで型変換をデータフレーム全体に適用
# 欠損値でエラーが出るためmapの引数 na_actionを利用

import ast
compound_df['molecule_structures'] = compound_df['molecule_structures'].map(ast.literal_eval, na_action='ignore') 
for i, cmpd in compound_df.iterrows():
    # オリジナルのif条件
    # if compound_df.loc[i]['molecule_structures'] != None:
    
    if type(compound_df.loc[i]['molecule_structures']) != float: # Nanがfloat型であることを利用したif条件に書き換えた
        compound_df.loc[i]['molecule_structures'] = cmpd['molecule_structures']['canonical_smiles']

print (compound_df.shape)
# (5435, 2)

出力データの準備

知りたい値をChEMBL IDに基づいて一つのデータフレームに統合します。 * ChEMBL-IDs * SMILES * 単位 * IC50

output_df = pd.merge(bioact_df[['molecule_chembl_id','units','value']], compound_df, on='molecule_chembl_id')
print(output_df.shape)
output_df.head()
# (5435, 4)
molecule_chembl_id units value molecule_structures
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)...
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N...
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1

カラム(列)名を、それぞれIC50とSMILESに変更します。

output_df = output_df.rename(columns= {'molecule_structures':'smiles', 'value':'IC50'})
output_df.shape
# (5435, 4)

SMILES表記が手に入らない化合物については、以降のトークトリアルで使用することができません。ですので、SMILES列に値のない化合物を削除します。

output_df = output_df[~output_df['smiles'].isnull()]
print(output_df.shape)
output_df.head()
# (5428, 4)
molecule_chembl_id units IC50 smiles
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)...
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N...
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1

次のセルを実行すると、小さなIC50値が読みづらいことがわかると思います。そこで、IC50値をpIC50値に変換したいと思います。

output_df = output_df.reset_index(drop=True)
ic50 = output_df.IC50.astype(float) 
print(len(ic50))
print(ic50.head(10))
5428
0         41.0
1        170.0
2       9300.0
3     500000.0
4    3000000.0
5      96000.0
6       5310.0
7     264000.0
8        125.0
9      35000.0
Name: IC50, dtype: float64
# IC50をpIC50に変換し、pIC50の列を付け足す
pIC50 = pd.Series() 
i = 0
while i < len(output_df.IC50):
    value = 9 - math.log10(ic50[i]) # pIC50=-log10(IC50 mol/l) --> for nM: -log10(IC50*10**-9)= 9-log10(IC50)
    if value < 0:
        print("Negative pIC50 value at index"+str(i))
    pIC50.at[i] = value
    i += 1
    
output_df['pIC50'] = pIC50
output_df.head()
molecule_chembl_id units IC50 smiles pIC50
0 CHEMBL68920 nM 41.0 Cc1cc(C)c(/C=C2\C(=O)Nc3ncnc(Nc4ccc(F)c(Cl)c4)... 7.387216
1 CHEMBL69960 nM 170.0 Cc1cc(C(=O)N2CCOCC2)[nH]c1/C=C1\C(=O)Nc2ncnc(N... 6.769551
2 CHEMBL137635 nM 9300.0 CN(c1ccccc1)c1ncnc2ccc(N/N=N/Cc3ccccn3)cc12 5.031517
3 CHEMBL306988 nM 500000.0 CC(=C(C#N)C#N)c1ccc(NC(=O)CCC(=O)O)cc1 3.301030
4 CHEMBL66879 nM 3000000.0 O=C(O)/C=C/c1ccc(O)cc1 2.522879

EGFRを標的として集めた生理活性データ

集めたデータセットを眺めてみましょう。

化合物を描画する

次の工程ではデータフレームに化合物を追加し、もっとも大きいpIC50値をもつ化合物の構造を眺めます。

PandasTools.AddMoleculeColumnToFrame(output_df, smilesCol='smiles')

pIC50値で化合物をソートします。

output_df.sort_values(by="pIC50", ascending=False, inplace=True)
output_df.reset_index(drop=True, inplace=True)

もっとも活性な分子を、もっとも大きなpIC50値の分子として描画します。

output_df.drop("smiles", axis=1).head()
molecule_chembl_id units IC50 pIC50 ROMol
0 CHEMBL63786 nM 0.003 11.522879 Mol
1 CHEMBL53711 nM 0.006 11.221849 Mol
2 CHEMBL35820 nM 0.006 11.221849 Mol
3 CHEMBL53753 nM 0.008 11.096910 Mol
4 CHEMBL66031 nM 0.008 11.096910 Mol

訳注(04.2020時点)

ROMolに構造が描画されず、(SMILESでもない)文字列になっている場合はRDKitとPandasのバージョンの互換性がよくない可能性があります。Pandasを version 0.25.0 以前のものにすることで解決する可能性があります。 RDKit イシュー を参照してください。

訳注ここまで

出力ファイルの書き出し

以降のトークトリアルでデータを使うため、csvファイルとしてデータを保存します。データを保存するときには化合物の列(分子の構造図のみをもつ列)を削除することをお勧めします。

output_df.drop("ROMol", axis=1).to_csv("../data/T1/EGFR_compounds.csv")

ディスカッション

このチュートリアルではChEMBLのデータベースから、私たちの目的とするターゲット分子に対して取得可能な全ての生理活性データを集めました。データセットをIC50あるいはpIC50の測定値をもつものだけにフィルタリングしました。

ChEMBLのデータはさまざまなデータソースに由来するものであることに注意してください。化合物のデータは世界中の様々なラボの様々な人々によって生み出されたものです。したがって、このデータセットを使って予測を行うときには注意する必要があります。結果を解釈し、自分たちの予測にどれだけの信頼を置けるか決める際には、データのソースとデータの生み出されたアッセイの一貫性について、常に慎重に考えることが重要です。

次のチュートリアルでは、手に入れたデータをリピンスキーのルールオブファイブでフィルタリングし、そして好ましくない部分構造でフィルタリングします。もう一つな他の重要なステップとして、データを綺麗に整型し、重複を取り除くことが挙げられます。(まだ、)私たちのトークトリアルでは取り上げられていないので、このタスクに使えるツールとして、標準化ライブラリ (github Francis Atkinson) あるいは MolVS というものがあることに言及しておきたいと思います。

クイズ

  • このトークトリアルではChEMBLから化合物と生理活性データをダウンロードしました。ChEMBLデータベースは他にどのように利用することができるしょうか?
  • IC50とEC50の違いは何でしょうか?
  • ChEMBLから取り出したデータはどのような目的で使いことができるでしょうか?

訳以上

追記

私はCADD、プログラミングの専門性がないので翻訳、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。 また環境の違いのためか、GitHubのコードそのままでは実行できない箇所が散見されました。行き当たりばったりでごり押ししているので誤り等多々あると思います。 オススメの解決策等あればご教示いただければ幸いです。

Jupyter notebookからmarkdwonで出力ごとダウンロードできるの知らなかった!便利!