magattacaのブログ

日付以外誤報

TeachOpenCADD トピック３ 〜部分構造による化合物フィルタリング〜

TeachOpenCADD

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

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

github.com

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

以下、訳。

トークトリアル 3

化合物フィルタリング：好ましくない部分構造

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

Maximilian Driller and Sandra Krüger

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

いくつか私たちのスクリーニングライブラリーに含めたくない部分構造があります。このトークトリアルでは、そのような好ましくない部分構造の様々なタイプを学び、そしてRDKitを使ってそれらの部分構造を見つけ、ハイライトする方法を学びます。

学習の目標

理論

  • 好ましくない部分構造とは何か?
  • 汎分析干渉化合物 Pan Assay Interference Compounds (PAINS)

実践

  • ChEMBLデータベースから化合物のセットを読み込む（トークトリアル 2　で準備したもの）
  • RDKitでの実装を利用し好ましくない部分構造をフィルターにかけ取り除く
  • 独自の好ましくない部分構造のリストを作成し、フィルタリングを実施
  • 部分構造の検索とハイライト

レファレンス

理論

好ましくない部分構造

部分構造には好ましくないものがあります。例えば毒性あるいは反応性があるといった理由や、薬物動態学的特性が好ましくないという理由、あるいは特定のアッセイに干渉する可能性が高いという理由です。 今日の医薬品探索ではよくハイスループットスクリーニング (HTS wikipedia)　を実施します。好ましくない部分構造をフィルタリングすることで、より望ましいスクリーニングライブラリを構築することができます。これにより、スクリーニングの前にライブラリの数を減らすことができ、時間と資源の節約につながります。

Brenkらは顧みられない病気の治療のための化合物の探索に使用するスクリーニングライブラリーをフィルタリングするため、好ましくない部分構造のリストを作成しました(Chem. Med. Chem. (2008), 3,435-444)。好ましくない特徴の例として、ニトロ基（変異原性の問題）、スルホン酸基とリン酸基（薬物動態学的特性が好ましくない可能性が高い）、2-ハロピリジンとチオール基（反応性の問題）、といったものが挙げられます。

好ましくない部分構造のリストは上記の文献に報告されており、このトークトリアルの実践編でも使用します。

汎分析干渉化合物 Pan Assay Interference Compounds (PAINS)

概要

PAINS (PAINS wikipedia) は実際には偽陽性にも関わらず、HTSでしばしばヒットとして見出される化合物です。PAINSは特定の一つのターゲット分子と反応するというよりもむしろ、非特異的に多数のターゲット分子と反応する傾向があります。通常、非特異的な結合もしくはアッセイの構成要素との相互作用により、様々なアッセイで様々なタンパク質に対して見かけ上結合します。

f:id:magattaca:20200418212342j:plain

Figure 1:汎分析干渉化合物（Pan Assay Interference Compounds、PAINS）の観点における特異的結合と非特異的結合。図はwikipediaより引用

Baellらによる文献で用いられたフィルター

Baellらはアッセイのシグナルに干渉する部分構造に焦点を当てました (J. Med. Chem. (2010), 53(7),2719-2740)。そのようなPAINSを見つけるのに役立つ部分構造について記述し、部分構造フィルタリングに利用できるリストを提供しました。

実践

データの読み込みと可視化

まず、必要なライブラリをインポートし、トークトリアル T2 でフィルタリングしたデータセットを読み込み、はじめの分子群を描画します。

import pandas
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import rdFMCS
from rdkit.Chem import AllChem
from rdkit.Chem import Descriptors
from rdkit.Chem import Draw
from rdkit import DataStructs
from rdkit.Chem import PandasTools
import matplotlib.pyplot as plt

filteredData = pandas.read_csv("../data/T2/EGFR_compounds_lipinski.csv", delimiter=";", index_col=0)
filteredData.drop(['HBD','HBA','MW','LogP'], inplace=True, axis=1) # 不必要な情報の除去
print ('Dataframe shape: ', filteredData.shape) # データフレームの次元をプリント
filteredData.head(5)

# Dataframe shape:  (4523, 6)
molecule_chembl_id units IC50 smiles pIC50 rule_of_five_conform
0 CHEMBL63786 nM 0.003 Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 11.522879 yes
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849 yes
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849 yes
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910 yes
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 yes 
PandasTools.AddMoleculeColumnToFrame(filteredData, smilesCol='smiles') # 分子を列に追加
# 最初の２０個の分子を描画
Draw.MolsToGridImage(list(filteredData.ROMol[0:20]), 
                    legends=list(filteredData.molecule_chembl_id[0:20]), 
                    molsPerRow=4)

f:id:magattaca:20200418212826p:plain

RDKitを使用してPAINSをフィルタリング

PAINSフィルターはすでにRDKitに実装されています(RDKit Documentation)。使ってみましょう！

from rdkit.Chem.FilterCatalog import *
params = FilterCatalogParams()
# 全　PAINS (A, B and C)　からカタログを構築
params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
catalog = FilterCatalog(params)

# フィルタリングしたデータを格納するための空のデータフレームを作成
rdkit_highLightFramePAINS = pandas.DataFrame(columns=('CompID', 'CompMol', 'unwantedID'))
rdkit_noPAINS = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles','pIC50'))
rdkit_withPAINS = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles', 'pIC50','unwantedID'))

# フィルタリングしたデータフレームのインデックスと行についてforループを回す
for i,row in filteredData.iterrows():
    curMol = Chem.MolFromSmiles(row.smiles) # 現在の分子
    match = False # Falseにmatchを設定
    rdkit_PAINSList = []
    # 最初のmatchを取得
    entry = catalog.GetFirstMatch(curMol)
    
    if entry!=None:
        # 現在の好ましくない部分構造の名前をリストに追加
        rdkit_PAINSList.append(entry.GetDescription().capitalize())
        # データフレームに関連するマッチング情報を追加
        rdkit_highLightFramePAINS.loc[len(rdkit_highLightFramePAINS)] = [row.molecule_chembl_id, curMol,
        entry.GetDescription().capitalize()]
        match = True
    if not match:
        # PAINSを含まない化合物のデータフレームに追加
        rdkit_noPAINS.loc[len(rdkit_noPAINS)] = [row.molecule_chembl_id, row.smiles, row.pIC50]
    else: 
        # PAINSを含む化合物のデータフレームに追加
        # 好ましくない部分構造を含むデータフレームに関連する情報を加える
        rdkit_withPAINS.loc[len(rdkit_withPAINS)] = [row.molecule_chembl_id, row.smiles, row.pIC50, entry.GetDescription().capitalize()]

print("Number of compounds with PAINS: %i"%(len(rdkit_withPAINS)))
print("Number of compounds without PAINS: %i (=remaining compounds)"%(len(rdkit_noPAINS)))

rdkit_highLightFramePAINS.head(10)

# Number of compounds with PAINS: 394
# Number of compounds without PAINS: 4129 (=remaining compounds)
CompID CompMol unwantedID
0 CHEMBL93032 Mol Catechol_a(92)
1 CHEMBL2029429 Mol Anil_di_alk_a(478)
2 CHEMBL2029428 Mol Anil_di_alk_a(478)
3 CHEMBL4071474 Mol Anil_di_alk_a(478)
4 CHEMBL4126630 Mol Anil_di_alk_a(478)
5 CHEMBL3663926 Mol Anil_di_alk_a(478)
6 CHEMBL3663925 Mol Anil_di_alk_a(478)
7 CHEMBL3663936 Mol Anil_di_alk_a(478)
8 CHEMBL3663923 Mol Anil_di_alk_a(478)
9 CHEMBL3663934 Mol Anil_di_alk_a(478)

訳注(04/2020)

訳者はPythonに詳しくないのでrdkit_highLightFramePAINS.loc[len(rdkit_highLightFramePAINS)] の部分に関する備忘録です。
locはPandas DataFrameで行をindexで指定する方法です。 len(~)でDataFrameの行数を取得することですでにデータが入っている行数がわかります。
DataFrameのindexは0から始まっているので、DataFrameの行数をindexとして指定することで、新しいデータを新しい行（既存の行数+1の行）に追加することになるのだと思います。

また、PAINSの利用例はRDKit Documentationのこちらに同様のコードが見られます。

訳注ここまで

好ましくない/毒性の懸念される部分構造でフィルタリング(Brenkのリスト)

RDKitには、PAINSのような好ましくない部分構造のリストがいくつかすでに実装されています。ですが、実装されていないRDKit外部のリストを使って部分構造の一致検索を行うことも可能です。ここではBrenkらによる文献(Chem. Med. Chem. (2008), 3,435-444)のsupporting informationで提供されているリストを使用します。

注）　データをダウンロードし、dataフォルダにcsvファイルとして保存済みです（形式は name-space-SMARTS です）。まず、データを読み込みます。

unwantedSubs = []
unwantedNames = []

for line in open('../data/T3/unwantedSubstructures.csv', 'r'):
    if not line.startswith("#"): # ヘッダーを無視
        splitted = line.strip().split(" ") # 各行を分割
        m = Chem.MolFromSmarts(splitted[1]) # SMARTSから分子を生成
        name = splitted[0].capitalize() # nameに名前を保存
        unwantedNames.append(name) # 好ましくない部分構造の名前をリストに追加
        unwantedSubs.append(m) # 好ましくない部分構造をリストに追加
print("Number of unwanted substructures in list =", len(unwantedSubs)) # 好ましくない部分構造の数を表示

# Number of unwanted substructures in list = 104

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
「unwantedSubstructures.csv」ファイルは１行目のみ「#name smart」となっています。if not line.starwith(#)の部分で「#」から始まる１行目(header)を飛ばしています。
strip()は文字列の先頭・末尾の余分な文字を削除するメソッドで、空白文字としてスペース、タブ、改行も取り除けるのでcsvの改行文字を除いているのだと思います。
split(" ")はスペースで分割しているので、改行文字を取り除いた各行のname-space-SMARTS形式をspaceで分割しnameとSMARTSにしているのだと思います。
こう考えるとsplitted[0]がname、splitted[1]がSMARTSになっていることがわかりやすい気がします。

訳注ここまで

２、３個部分構造を見てみましょう　(すべてのSMARTSが表示できるわけではありません。一部を選んでいます。)。

Chem.Draw.MolsToGridImage(list(unwantedSubs[2:5]), subImgSize=(200, 300), legends=unwantedNames[2:5])

f:id:magattaca:20200418215706p:plain

これらの好ましくない部分構造とマッチするものがあるか、フィルタリングしたデータフレームの中を検索してみます。

# フィルタリングしたデータを格納するためのデータフレームを作成
highLightFrameUNW = pandas.DataFrame(columns=('CompID', 'CompMol', 'unwantedID', 'unwSubstr'))
noUnwanted = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles','pIC50'))
withUnwanted = pandas.DataFrame(columns=('ChEMBL_ID', 'smiles', 'pIC50','unwantedID'))
molsToDraw = []

# データセットの各化合物に対して
for i,row in filteredData.iterrows(): # フィルタリングしたデータフレームのインデックスと行についてfor ループを回す
    curMol = Chem.MolFromSmiles(row.smiles) # 現在の分子
    match = False # Falseにmatchを設定
    unwantedList = []
    molsToDraw.append(curMol)
    
    # 全ての好ましくない部分構造を検索
    for idx, unwSub in enumerate(unwantedSubs):
        # 部分構造があるかチェック
        if curMol.HasSubstructMatch(unwSub): # 現在の分子が好ましくない部分構造を有する場合
            match = True # matchをTrueに設定
            unwantedList.append(unwantedNames[idx]) # 好ましくない部分構造の名前をリストに追加
            # 関連する情報をデータフレームに追加
            highLightFrameUNW.loc[len(highLightFrameUNW)] = [row.molecule_chembl_id, curMol, unwantedNames[idx], unwSub]
    if not match: # 一致する部分構造が見つからない場合
        noUnwanted.loc[len(noUnwanted)] = [row.molecule_chembl_id, row.smiles, row.pIC50]
        # 欲しい部分構造のデータフレームに関連情報を追加
    else: # 一致が見つかった場合
        withUnwanted.loc[len(withUnwanted)] = [row.molecule_chembl_id, row.smiles, row.pIC50, unwantedList] #　好ましくない部分構造データフレームに関連する情報を追加

print("Number of compounds with unwanted substructures: %i"%(len(withUnwanted)))
print("Number of compounds without unwanted substructures: %i (=remaining compounds)"%(len(noUnwanted)))

# Number of compounds with unwanted substructures: 2356
# Number of compounds without unwanted substructures: 2167 (=remaining compounds)

highLightFrameUNW.head(4)
CompID CompMol unwantedID unwSubstr
0 CHEMBL63786 Mol Polycyclic-aromatic-hydrocarbon Mol
1 CHEMBL1243316 Mol Aliphatic-long-chain Mol
2 CHEMBL1243316 Mol Michael-acceptor Mol
3 CHEMBL1243316 Mol Triple-bond Mol

部分構造を分子の中で直接ハイライトすることもできます。

first_highLightFrameUNW = highLightFrameUNW.head(8) # リストの最初の8エントリーのサブセット

# 分子を描画し、好ましくない部分構造をハイライトする 
Draw.MolsToGridImage(list(first_highLightFrameUNW["CompMol"]), subImgSize=(400,300),
    molsPerRow=2, highlightAtomLists=
    [m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])], 
    legends=list(first_highLightFrameUNW["CompID"]+": "+first_highLightFrameUNW["unwantedID"]))

f:id:magattaca:20200418215814p:plain

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
[m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])]の部分ですが、enumerateはインデックス番号と要素を同時に取得するメソッドで、ここではそれぞれiとmを割り当てています。
要素m分子がもつ好ましくない部分構造を見つけるために、RDKitのGetSubstructMatchを使用しており、検索対象の好ましくない部分構造をunwSubstr列のインデックスiを参照するという形で実行されています。

訳注ここまで

例をSVGファイルとして保存します。

# イメージをファイルに保存
img = Draw.MolsToGridImage(list(first_highLightFrameUNW["CompMol"]), subImgSize=(400,300),
    molsPerRow=3, highlightAtomLists=
    [m.GetSubstructMatch(first_highLightFrameUNW["unwSubstr"][i]) for i,m in enumerate(first_highLightFrameUNW["CompMol"])], 
    legends=list(first_highLightFrameUNW["unwantedID"]), useSVG=True)

# SVGデータの取得
molsvg = img.data

# 不透明な背景を透明に置き換える
molsvg = molsvg.replace("opacity:1.0", "opacity:0.0");
molsvg = molsvg.replace("12px", "24px");

# 変更されたSVGデータをファイルに保存
f = open("../data/T3/substructures.svg", "w")
f.write(molsvg)
f.close()

好ましくない部分構造を持つ化合物と持たない化合物のリストを保存。

# 好ましくない部分構造を有する化合物をcsvファイルに書き出し
withUnwanted.to_csv("../data/T3/EGFR_compounds_lipinski_noPAINS.csv", sep=',') 

# 好ましくない部分構造を持たない化合物をcsvファイルに書き出し
noUnwanted.to_csv("../data/T3/EGFR_compounds_lipinski_noPAINS_noBrenk.csv", sep=',') 

# 好ましくない部分構造を有する化合物のcsvファイルの最初のいくつかを表示
noUnwanted.head()
ChEMBL_ID smiles pIC50
0 CHEMBL53711 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849
1 CHEMBL35820 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849
2 CHEMBL53753 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910
3 CHEMBL66031 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910
4 CHEMBL176582 Cn1cnc2cc3ncnc(Nc4cccc(Br)c4)c3cc21 11.000000

見つかった好ましくない部分構造をさらに解析します。

# 最も頻度の高い化合物の数を数えます。
unwCounts = {}
for ele in highLightFrameUNW.values:
    unw = ele[2] # highLightFrameUNWデータフレームに含まれる好ましくない部分構造のID
    if unwCounts.get(unw, "empty") == "empty": # 好ましくない構造のIDがまだ辞書に含まれていない場合
        unwCounts[unw] = [1, ele[3]] # 1を付与し、辞書の好ましくない構造に追加
    else: # もしkey （好ましくない構造のID）がすでに存在しているなら、出現回数の値を１増やす 
        unwCounts[unw][0] += 1

frequentUNW = []
frequentUNWNames = []

# unwCountsに含まれる構造: IDをkey、出現頻度と分子をvalueとしてもつ辞書
# 例　('acyclic-C=C-O', [7, <rdkit.Chem.rdchem.Mol object at 0x7fa58fc06710>])

# 部分構造の頻度で辞書をソートする
for key, value in sorted(unwCounts.items(), key=lambda kv: kv[1][0], reverse=True):
    frequentUNW.append(value[1]) # 部分構造
    frequentUNWNames.append(key)

訳注(04/2020)

Pythonやプログラミング初心者仲間のための蛇足・・・
for ele in highLightFrameUNW.valuesの部分ですが、valuesはDataFrameの実際のデータの値で、NumPyの配列ndarrayです。 DataFramehighLightFrameUNWを参照していただければわかると思いますが、ele[2]は２列目unwantedID列を、ele[3]は3列目unwSubstr列の要素を指定しています。
unwCounts.get(unw, "empty")の部分ですが、get(key[, default])は辞書型においてkey が辞書にあれば key に対する値を、そうでなければ default を返します。ですので、今回は辞書unwCountsにkeyunw(unwntedID, ele[2])があるかどうかを検索し、なければdefalutとしてempytを返すというようになっています。辞書unwCountsのkeyはunwntedID、valueはリスト[出現回数, 好ましくない部分構造]となっているので、辞書に含まれていない好ましくない部分構造の場合は新しくエントリーを辞書に追加して値を出現回数を1に設定、すでに辞書に含まれている構造の場合は出現回数[unw][0]をインクリメントするという手順になっているのだと思います。

訳注ここまで

頻度の高い部分構造の上位８個を描画します。

# 上位８個の頻度の高い部分構造
Draw.MolsToGridImage(mols=list(frequentUNW[0:8]), 
                     legends=list(frequentUNWNames[0:8]),
                     molsPerRow=4)

f:id:magattaca:20200418215914p:plain

ディスカッション

このトークトリアルでは、好ましくない部分構造の検索を行う２つの方法を学びました。: * RDKitにすでに実装されている FilterCatalogクラス、そして、 * RDKitに含まれていない外部のリストとHasSubstructMatch()関数です。

実際のところ、PAINSを部分構造検索で実装することもできます。また、BrenkらによるリストもすでにRDKitに実装されています。実装済みのその他のリストについては (RDKit Documentation)を参照してください。

これまで、HasSubstructMatch()関数を使用してきましたが、これは一つの化合物あたり一つのマッチ結果しか出しません。GetSubstructMatches()　関数を使用すると、一つの化合物に含まれる全ての部分構造を見つけることができます。 同様にPAINSに関して、一つの分子あたり、最初にマッチしたもの GetFirstMatch() だけをみてきました。もし、全てのPAINSをフィルタリングして除去したいのであればこれで十分です。ですが、ある分子のもつ全ての危険な部分構造を見るために、GetMatches()を使用することもできます。一つの分子あたり全てのマッチする部分構造を考慮しているわけではないので、最後に描画した部分構造が実際に最も頻度の高いものであったと言うことはできません。ですが、頻度が高いということだけは明らかです。

見つかった部分構造は２つの異なる方法で処理できます。 * 部分構造検索をフィルターとして適用し、引っかかった化合物を、資金と時間の節約のため、さらなる試験からは除外する * あるいは警告として使用することもできます。好ましくない部分構造を有する分子にフラグを立てることができます。（例えば、化学者や毒性学者•・・といった）専門家の目から見て、経験に基づき判断することができるかもしれません。もし、各部分構造がそこまでクリティカルなものでなければ、スクリーニングの対象化合物として含めることも可能です。

ここでは、機械学習に使用するために、あまりに多くの化合物を削除することはしたくなかったので、好ましくない部分構造をフィルタリングで除去することはしません。また、部分構造によるフィルタリングは、後ほど実際のスクリーニング実験を行うまえに適用することもできます。警告（アラート）のフラッグを設定することも可能で、そうすれば（PAINSやBrenk等どのようなリストにでも順じて）好ましくない部分構造についての情報を保持し、あとで考慮することが可能です。

クイズ

  • なぜスクリーニングライブラリーから「PAINS」を取り除くことを考えるべきなのでしょうか？これらの化合物の問題点とは何でしょうか？
  • ある好ましくない部分構造を取り除く必要がないような状況を見つけることはできますか？
  • このチュートリアルで使用した部分構造はどうやってエンコードされていましたか？

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。

今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。

PAINSについて:
* 化学の新しいカタチ：RDKitのPAINSフィルターで化合物をスクリーニング

RDKitを用いた部分構造検索
* 化学の新しいカタチ：RDKitを用いた部分構造検索とMCSアルゴリズム

さて、今回の記事では医薬品として好ましくない部分構造を検索し、データセットから取り除く方法が取り上げられています。一方で、好ましくない部分構造が状況によって変わりうるののだ、ということも指摘されています。例えば冒頭で好ましくない部分構造の例としてニトロ基や2-ハロピリジンといったものが指摘されています。ではこれらの構造が「絶対に医薬品としてありえないか？」というとその答えはNOです。

ニトロ基を有する医薬品の例として高血圧・狭心症治療剤のニフェジピン、より最近の承認薬としては大塚製薬の肺結核治療薬デラマニドが挙げられます。

「ニトロ基がなぜ嫌われるのか?」というと生物学者の方にはアッセイの検出を思い浮かべていただくとわかりやすいかもしれません。学生実習の題材でも使われるアルカリホスファターゼの酵素活性の検出では、酵素の機能により放出されるp-nitrophenyl phosphate (pNPP)の蛍光性が検出原理として使われています。光を吸収、呈色するような分子構造、なんとなく飲みたくないですよね。（大学受験を終えたばかりでこんなブログにたどりついてしまった方にはキサントプロテイン反応というと分かりやすいかもしれません。）

また、2-ハロピリジンそのものではありませんが現在話題のファビピラル(アビガン)も類似の部分構造を有しています。ピラジン骨格ですがN横の脱離基（ハロゲン）としては反応性が高い構造かもしれません。「部分構造をどの範囲までとるか？」という基準で議論が分かれそうです。ぜひ専門家の皆様のご意見を伺いたいところです。

TeachOpenCADD トピック2 〜指標による化合物フィルタリング〜

TeachOpenCADD

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

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

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

github.com

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

以下、訳。

トークトリアル 2

化合物フィルタリング: ADME とリードライクネスのクライテリア

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

Michele Ritschel and Mathias Wajnberg

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

ChEMBLから取得した化合物 (talktorial 1)をリードライクネスのクライテリアに基づきフィルタリングします。我々のスクリーニングライブラリーから薬らしくない分子を取り除くのが目的です。

  • 化合物の生物学的利用可能性（バイオベイラビリティ）に関連するパラメータを計算(リピンスキーのルールオブファイブ)
  • ルールオブファイブのクライテリアに基づきChEMBLから集めた化合物をフィルタリング
  • レーダーチャート形式でパラメーターをプロット

学習の目標

理論

  • ADME - 吸収: absorption, 分布: distribution, 代謝: metabolism and 排出: excretion
  • リードライクネスとリピンスキーのルールオブファイブ
  • リードライクネスの文脈におけるレーダーチャートのバリエーションと解釈

実践

  • 例となる化合物について物理化学的パラメータを計算
  • 複数の分子についてそれぞれの物理化学的パラメーターを比較するため棒グラフを作成
  • ルールオブファイブを満たしているかチェックするための関数を定義
  • ChEMBLから取得したデータセット全体にルールオブファイブを適用
  • ルールオブファイブでフィルターをかけたデータセットのレーダーチャートを作成。これにより一つのプロットでルールオブファイブのクライテリアに関するプロパティの可視化が可能

レファレンス

理論

バーチャルスクリーニングでは、ある化合物が特定のターゲット分子に対して結合し、相互作用するかということを予測することができます。ですが、新しい医薬品を見つけ出すためには、その化合物がターゲット分子に到達し、最終的には体内から好ましい経路で取り除かれることが重要です。したがって、化合物が実際に体の中に取り込まれ、ターゲット分子に到達するためのて特定のバリアーを通過することができるかということも考慮しなければいけません。代謝的に安定か？そしてターゲット分子に対して作用しなくなったらどのように排出されるのか？これらの過程は薬物動態学（pharmacokinetics）の分野で研究されています。薬物動力学（pharmacodynamics）（「医薬品は私たちの体に何をするのか？」）とは対照的に、薬物動態学は「私たちの体のなかで医薬品に対して何が起こるのか？」という問いを扱います。

ADME

薬物動態学は主に４つのステップに分けられます。 Absorption：吸収、 Distribution：分布、 Metabolism：代謝、そして Excretion：排出です。 これら４ステップは ADME　と略されます。ADME(T) はToxicology（毒性学）も含むことがあります。 以下ではADMEの過程がより詳細に議論されています。 (ADME wikipedia and Mol Pharm. (2010), 7(5), 1388-1405)

  • 吸収：ある物質が体内に取り込まれる量と取り込みに必要な時間は、各個人とその状態によって異なる様々な要因と、その物質の特性に依存しています。（低い）化合物溶解度と胃内容排出時間（gastric emptying time）、腸通過時間（intestinal transit time）、胃の中での化合物（不）安定性、そして腸壁の（不）透過能、といった要因が全て、医薬品が、例えば経口投与、吸入、そして皮膚への接触のあとでどの程度吸収されるかに影響をあたる可能性があります。

  • 分布：吸収された物質の分布、即ち体内における、血液とそれぞれの組織の間の分布、血液脳関門の透過は、局所的な血流速度や化合物のサイズと極性、そして血清中のタンパク質やトランスポーター酵素への結合の影響を受けます。毒性学における重要な要因としては脂肪組織への疎水性の高い物質の蓄積、あるいは血液脳関門の透過があげられます。

  • 代謝：化合物が体内に入るとすぐに、通常、化合物の代謝が始まります。つまり実際には、この化合物の一部だけがターゲット分子に到達します。主に肝臓と腎臓の酵素が生体異物（xenobiotics：体の外からきた物質）の分解を担っています。毒性のある化合物が取り除かれる場合、吸収された物質の量を減らすことは好ましいことです。その一方で、化学物質の変換は新たな毒性代謝物をうみだす可能性もあります。

  • 排出：化合物と代謝物は排出によって体内から取りのぞかれる必要があります。通常、腎臓（尿）あるいは糞によって排出されます。不完全な排出は結果として、外来物質の蓄積、あるいは通常の代謝を妨げることに繋がります。
Figure 1:ヒト体内におけるADME過程
(figure taken from openclipart.org and adapted)

リードライクネスとリピンスキーのルールオブファイブ

リード化合物 とは有望な特性をもった医薬品開発の候補となる化合物です。望ましい医薬品を見つけ出すために構造修飾を行う出発点となる構造です。生理活性（目標とするターゲット分子に結合する化合物）に加えて、好ましいADME特性が有効な医薬品をデザインする上で重要なクライテリアです。

生物学的利用可能性（バイオアベイラビリティ）は重要なADME特性で、化合物の構造のみに基づいてこの特性を計るために、リピンスキーのルールオブファイブ（Lipinski's rule of five）が発明されました。これは経験則で、化合物の経口投与によるバイオアベイラビリティを見積もるのに役立ちます。

ルールオブファイブによると、ある化合物が次のルールの一つ以上を破っているとき、経口で体内に吸収される可能性が低くなるとされいます。そのルールとは即ち：

  • 分子量が500ダルトン以下
  • 水素結合アクセプターが10以下
  • 水素結合ドナーが５以下
  • LogP (オクタノールー水　分配係数) <= 5

LogP は分配係数あるいはオクタノールー水分配係数とも呼ばれます。通常、疎水性（例：1-オクタノール）と親水性（例：水）の相の間における化合物の分配を測定します。

疎水性分子は水に対する溶解度が低くなる一方で、より親水性の分子（例：水素結合のアクセプターとドナーの数が多いもの）や大きな分子（分子量の大きいもの）は、リン脂質二重膜を透過しづらくなります。

ルールオブファイブに関して、全ての数字が５の倍数になっています。これがこのルールの名称の起源となっています。

(Adv. Drug Deliv. Rev. (1997), 23, 3-25)

レーダーチャート

ルールオブファイブに関する分子の特性を計算した後で、その特性を可視化できると役に立ちます。Ritchieら(Drug. Discov. Today (2011), 16(1-2), 65-72)　はADMEに関する特性のグラフィカルな表現についての概観を示しています。分子の特性を可視化し、創薬化学者による解釈を手助けするために様々な手法があります（例：Craigプロット、フラワープロット、黄金の三角地帯（golden triangle））。

このチュートリアルではPythonの作図ライブラリである matplotlib を使ってレーダーチャートを作成する方法を学びます。その見た目から、レーダーチャート(radar charts wikipedia) はしばしばスパイダープロット、あるいはクモの巣グラフと呼ばれます。３６０度に沿って配置され、それぞれの条件に対して、中心から出発する軸を一つ持っています。各パラメーターの値は軸に沿ってプロットされ直線で結ばれています。影のついたエリアはパラメーターが条件を満たす領域を示します。

Figure 2: 化合物データセットの物理化学的特性を表すレーダーチャート

実践

化合物例を定義し可視化する

ChEMBLから取り出したデータセット全体で実践する前に、４つの化合物を例として取り上げ、その化学的特性を調査したいと思います。必要なライブラリをインポートし、４化合物のSMILESを使って描画します。

from rdkit import Chem
from rdkit.Chem import Descriptors
import pandas as pd
from rdkit.Chem import Draw
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
from math import pi

smiles_1 = 'CCC1C(=O)N(CC(=O)N(C(C(=O)NC(C(=O)N(C(C(=O)NC(C(=O)NC(C(=O)N(C(C(=O)N(C(C(=O)N(C(C(=O)N(C(C(=O)N1)C(C(C)CC=CC)O)C)C(C)C)C)CC(C)C)C)CC(C)C)C)C)C)CC(C)C)C)C(C)C)CC(C)C)C)C' # Cyclosporine
smiles_2 = 'CN1CCN(CC1)C2=C3C=CC=CC3=NC4=C(N2)C=C(C=C4)C' # Clozapine
smiles_3 = 'CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC=CC=C(C)C=CC=C(C)C=CC2=C(CCCC2(C)C)C)C)C' # Beta-carotene
smiles_4 = 'CCCCCC1=CC(=C(C(=C1)O)C2C=C(CCC2C(=C)C)C)O' # Cannabidiol
smiles_list = [smiles_1, smiles_2, smiles_3, smiles_4]
names_list = ['cyclosporine', 'clozapine', 'beta-carotene', 'cannabidiol']
mol_list = [Chem.MolFromSmiles(smiles) for smiles in smiles_list]

Draw.MolsToGridImage(mol_list, legends=names_list, molsPerRow=4)

f:id:magattaca:20200412172944p:plain

ルールオブファイブに関する分子の特性を計算しプロットする

ルールオブファイブに関する化学特性を計算し視覚的に比較します：

MWs = [Descriptors.ExactMolWt(mol) for mol in mol_list]
HBAs = [Descriptors.NumHAcceptors(mol) for mol in mol_list]
HBDs = [Descriptors.NumHDonors(mol) for mol in mol_list]
LogPs = [Descriptors.MolLogP(mol) for mol in mol_list]
parameters = [MWs, HBAs, HBDs, LogPs]
print('Molecular weight of the four compounds:',MWs)
#  Molecular weight of the four compounds: [1201.841367992, 306.184446704, 536.438201792, 314.2245802]
  • 分子ごとの特性を棒グラフでプロット
# 2x2のプロットの枠を作成
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)
axes = [ax1, ax2, ax3, ax4]
x = np.arange(1, len(mol_list)+1)
colors = ['red', 'green', 'blue', 'cyan']

# サブプロットを作成
for index in x-1:
    axes[index].bar(x, parameters[index], color=colors)

# ルールオブファイブの閾値を点線で追加
ax1.axhline(y=500, color="black", linestyle="dashed")
ax1.set_title("molecular weight (Da)")
ax2.axhline(y=10, color="black", linestyle="dashed")
ax2.set_title("# h-bond acceptors")
ax3.axhline(y=5, color="black", linestyle="dashed")
ax3.set_title("# h-bond donors")
ax4.axhline(y=5, color="black", linestyle="dashed")
ax4.set_title("logP")

# 凡例の追加
legend_elements = [mpatches.Patch(color=colors[i], label=names_list[i]) for i in range(len(mol_list))]
legend_elements.append(Line2D([0], [0], color="black", ls="dashed", label="Threshold"))
fig.legend(handles=legend_elements, bbox_to_anchor=(1.25, 0.5))

# サブプロットと凡例を図に追加
plt.tight_layout()

plt.show()

f:id:magattaca:20200412173052p:plain

上の棒グラフでは４例の化合物についてルールオブファイブの特性 （分子量、水素結合ドナー数、水素結合アクセプター数、LogP）を比較しました。４例の医薬品分子が異なる特性を持つことが見てとれます。次のステップでは各化合物のそれぞれについてルールオブファイブを逸脱しているかどうか調べます。

リピンスキーのルールオブファイブに従っているかの調査

ある化合物がルールオブファイブを逸脱しているか否かを調べるための関数を定義し、４つの化合物例に適用します。

def test_rule_of_five(smi):
    m = Chem.MolFromSmiles(smi)
    
    # ルールオブファイブの化学特性を計算
    MW = Descriptors.ExactMolWt(m)
    HBA = Descriptors.NumHAcceptors(m)
    HBD = Descriptors.NumHDonors(m)
    LogP = Descriptors.MolLogP(m)
    
    # ルールオブファイブの条件
    conditions = [MW <= 500, HBA <= 10, HBD <= 5, LogP <= 5]
    # ４つの条件のうち逸脱が一つ以下の時にTrueを返す
    return conditions.count(True) >= 3

for i in range(len(smiles_list)):
    smi=smiles_list[i]
    name=names_list[i]
    print("Rule of five accepted for %s: %s "%(name,test_rule_of_five(smi)))                

# Rule of five accepted for cyclosporine: False 
#  Rule of five accepted for clozapine: True 
#  Rule of five accepted for beta-carotene: False 
#  Rule of five accepted for cannabidiol: True

我々の test_rule_of_five 関数によると４つの化合物例のうち２つはルールオブファイブを満たしていません。ここから、シクロスポリン（cyclosporin）とβ-カロテン（betacarotene）は経口で体内に吸収される可能性が低いと解釈できます。４つとも全て医薬品として市場で手に入るので、それぞれのターゲット分子にどうにかして到達しているはずです。ルールの例外である可能性もありますし、経口投与とは異なる投与ルートで使われている可能性もあります。

EGFRデータセットにルールオブファイブを適用

test_rule_of_five 関数を使って、リピンスキーのルールオブファイブに従っているかを基準にして、目的のメインデータセットにフィルターをかけます。

  • ルールオブファイブに関する化学特性を全て返すように関数を調整
  • メインとなるデータフレームを読み込み (ChEMBL_df)
  • ルールオブファイブの関数をChEMBL_df　に適用
  • ルールを２つ以上逸脱しているかによって ChEMBL_df にフィルターをかける
  • フィルタリングした後のデータフレームを保存
def df_rule_of_five(df):
    
    smi = df['smiles']
    m = Chem.MolFromSmiles(smi)
    
    # ルールオブファイブの化学特性を計算
    MW = Descriptors.ExactMolWt(m)
    HBA = Descriptors.NumHAcceptors(m)
    HBD = Descriptors.NumHDonors(m)
    LogP = Descriptors.MolLogP(m)
    
    # ルールオブファイブの条件
    conditions = [MW <= 500, HBA <= 10, HBD <= 5, LogP <= 5]
    
    # 特性の値とルールオブファイブに従っているか否かについての情報をもつpandasの行を作成
    return pd.Series([MW, HBA, HBD, LogP, 'yes']) if conditions.count(True) >= 3 else pd.Series([MW, HBA, HBD, LogP, 'no'])

ChEMBL_df = pd.read_csv('../data/T1/EGFR_compounds.csv', index_col=0)
print(ChEMBL_df.shape)
ChEMBL_df.head()
# (5428, 5)
molecule_chembl_id units IC50 smiles pIC50
0 CHEMBL63786 nM 0.003 Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 11.522879
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 
# ルールオブファイブの結果を得るために関数を適用 (しばらく時間がかかる可能性があります)
rule5_prop_df = ChEMBL_df.apply(df_rule_of_five, axis=1)

# 条件列に名前をつける
rule5_prop_df.columns= ['MW', 'HBA', 'HBD', 'LogP', 'rule_of_five_conform']

# 計算した値をもつデータセットを連結（concatenate）
ChEMBL_df = ChEMBL_df.join(rule5_prop_df)

# 空の行を削除 --> rule of five
filtered_df = ChEMBL_df[ChEMBL_df['rule_of_five_conform']=='yes']

# データに関する情報
print('# of compounds in unfiltered data set:', len(ChEMBL_df))
print('# of compounds in filtered data set:', len(filtered_df))
print("# of compounds not compliant with Lipinski's rule of five:", (len(ChEMBL_df)-len(filtered_df)))

# フィルターをかけたデータの保存
filtered_df.to_csv('../data/T2/EGFR_compounds_lipinski.csv', sep=';') 
filtered_df.head(10)

#  # of compounds in unfiltered data set: 5428
#  # of compounds in filtered data set: 4523
#  # of compounds not compliant with Lipinski's rule of five: 905
molecule_chembl_id units IC50 smiles pIC50 MW HBA HBD LogP rule_of_five_conform
0 CHEMBL63786 nM 0.003 Brc1cccc(Nc2ncnc3cc4ccccc4cc23)c1 11.522879 349.021459 3 1 5.28910 yes
1 CHEMBL53711 nM 0.006 CN(C)c1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.221849 343.043258 5 1 3.59690 yes
2 CHEMBL35820 nM 0.006 CCOc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OCC 11.221849 387.058239 5 1 4.93330 yes
3 CHEMBL53753 nM 0.008 CNc1cc2c(Nc3cccc(Br)c3)ncnc2cn1 11.096910 329.027607 5 2 3.57260 yes
4 CHEMBL66031 nM 0.008 Brc1cccc(Nc2ncnc3cc4[nH]cnc4cc23)c1 11.096910 339.011957 4 2 4.01220 yes
5 CHEMBL176582 nM 0.010 Cn1cnc2cc3ncnc(Nc4cccc(Br)c4)c3cc21 11.000000 353.027607 5 1 4.02260 yes
6 CHEMBL174426 nM 0.025 Cn1cnc2cc3c(Nc4cccc(Br)c4)ncnc3cc21 10.602060 353.027607 5 1 4.02260 yes
7 CHEMBL29197 nM 0.025 COc1cc2ncnc(Nc3cccc(Br)c3)c2cc1OC 10.602060 359.026939 5 1 4.15310 yes
8 CHEMBL1243316 nM 0.030 C#CCNC/C=C/C(=O)Nc1cc2c(Nc3ccc(F)c(Cl)c3)c(C#N... 10.522879 477.136781 6 3 4.75878 yes
9 CHEMBL3613702 nM 0.037 C=CC(=O)Nc1ccc2ncnc(Nc3cc(F)c(Cl)c(Cl)c3)c2c1 10.431798 376.029395 4 2 4.94380 yes

ルールオブファイブの特性を可視化するためのレーダープロット

まずデータセットの平均と標準偏差を求める関数を定義します。

これらの統計量はあとでデータセットのリピンスキーのルールオブファイブに関するパラメーターをプロットすることに使います。

def get_properties_stats(data_df):
    """
    データセットの物理化学特性の平均と標準偏差を計算する関数
    
    Input: 
    化合物ごとの物理化学特性の値、 HBD, HBA, MW そして　LogPを列として
    （正確にこのままの列名で）含むデータセット
    
    Output:
    物理化学特性（行）それぞれについて平均と標準偏差の値（列）をもつデータフレーム
    """
    properties = ["HBD", "HBA", "MW", "LogP"]
    
    data_stats = []
    
    for i in properties:
        std = data_df[i].std()
        mean = data_df[i].mean()
        da = pd.DataFrame([[round(mean, 2), round(std, 2)]], index=[i], columns=["mean", "std"])
        data_stats.append(da)
    
    data_stats = pd.concat(data_stats)
    
    return data_stats

リピンスキーのルールオブファイブに一致する化合物データセットの統計量を計算します（フィルタリングしたデータセット）。

stats_rof = get_properties_stats(filtered_df)
stats_rof
mean std
HBD 1.88 1.01
HBA 5.97 1.88
MW 413.03 87.93
LogP 4.07 1.19

リピンスキーのルールオブファイブに合致しない化合物データセットの統計量を計算します。

stats_not_rof = get_properties_stats(ChEMBL_df[ChEMBL_df['rule_of_five_conform']=='no'])
stats_not_rof
mean std
HBD 2.30 1.74
HBA 7.99 2.37
MW 588.96 103.05
LogP 5.96 1.44

レーダーチャートで化合物特性を可視化する関数を作成します。stackoverflowのチュートリアルに従って行います。

def plot_radarplot(data_stats, output_path):
    """
    ４つの物理化学特性（HBD、HBA、MW、LogP)の平均と標準偏差に基づくレーダープロットを作成する関数
    
    Input: 
    それぞれの物理化学特性（行）について平均と標準偏差（列）を有するデータフレーム。
    
    Output:
    レーダープロット (ファイルとして保存し、Jupyter notebookにも表示する)。
    """

    # 線を引くためのデータポイントを取得
    std_1 = [data_stats["mean"]["HBD"] + data_stats["std"]["HBD"], 
             (data_stats["mean"]["HBA"]/2) + (data_stats["std"]["HBA"]/2), 
             (data_stats["mean"]["MW"]/100) + (data_stats["std"]["MW"]/100), 
             data_stats["mean"]["LogP"] + data_stats["std"]["LogP"]]
    std_2 = [data_stats["mean"]["HBD"] - data_stats["std"]["HBD"], 
             (data_stats["mean"]["HBA"]/2) - (data_stats["std"]["HBA"]/2), 
             (data_stats["mean"]["MW"]/100) - (data_stats["std"]["MW"]/100), 
             data_stats["mean"]["LogP"] - data_stats["std"]["LogP"]]
    mean_val = [data_stats["mean"]["HBD"], (data_stats["mean"]["HBA"]/2), 
                (data_stats["mean"]["MW"]/100), data_stats["mean"]["LogP"]]

    # （塗りつぶされた）領域のデータポイント（ルールオブファイブ）を取得
    rule_conditions = [5, (10/2), (500/100), 5]
    
    # 特性の名称を定義
    parameters = ['# H-bond donors', '# H-bond acceptors/2', 'Molecular weight (Da)/100', 'LogP']

    # 
    N = len(rule_conditions)

    # フォントサイズを指定
    fontsize = 16

    # 条件の軸の角度
    x_as = [n / float(N) * 2 * pi for n in range(N)]

    # 円型チャートなので各リストの最初の値のコピーを、
    # 各リストの最後にデータとともに追加する必要があります。
    std_1 += std_1[:1]
    std_2 += std_2[:1]
    mean_val += mean_val[:1]
    rule_conditions += rule_conditions[:1]
    x_as += x_as[:1]

    # 図のサイズを指定
    plt.figure(figsize=(8,8))

    # 軸の色を指定
    plt.rc('axes', linewidth=2, edgecolor="#888888")

    # 極座標プロット（polar plot）の作成
    ax = plt.subplot(111, polar=True)

    # 時計回りを指定。つまり：
    ax.set_theta_offset(pi / 2)
    ax.set_theta_direction(-1)

    # y軸のラベル位置を指定
    ax.set_rlabel_position(0)

    # グリッドの色と線のスタイルを指定
    ax.xaxis.grid(True, color="#888888", linestyle='solid', linewidth=2)
    ax.yaxis.grid(True, color="#888888", linestyle='solid', linewidth=2)

    # 放射状の数の指定とラベルの除去
    plt.xticks(x_as[:-1], [])

    # y軸の印（ytick）の指定
    plt.yticks([1, 3, 5, 7], ["1", "3", "5"], size=fontsize)

    # 軸の範囲（limit）の指定
    plt.ylim(0, 7)

    # データをプロット
    # 平均値
    ax.plot(x_as, mean_val, 'b', linewidth=3, linestyle='solid', zorder=3)

    # 標準偏差
    ax.plot(x_as, std_1, 'm', linewidth=2, linestyle='dashed', zorder=3, color='#111111')
    ax.plot(x_as, std_2, 'y', linewidth=2, linestyle='dashed', zorder=3, color='#333333')

    # 領域を塗りつぶす
    ax.fill(x_as, rule_conditions, "#3465a4", alpha=0.2)

    # 正しくフィットしていることを確かめるためにy軸の印のラベルを描く
    for i in range(N):
        angle_rad = i / float(N) * 2 * pi
        if angle_rad == 0:
            ha, distance_ax = "center", 1
        elif 0 < angle_rad < pi:
            ha, distance_ax = "left", 1
        elif angle_rad == pi:
            ha, distance_ax = "center", 1
        else:
            ha, distance_ax = "right", 1
        ax.text(angle_rad, 7 + distance_ax, parameters[i], size=fontsize,
                horizontalalignment=ha, verticalalignment="center")

    # 左上のプロットに合わせて凡例を追加 
        labels = ('Mean', 'Mean + std', 'Mean - std', 'Rule of five area')
        legend = ax.legend(labels, loc=(1.1, .7),
                           labelspacing=0.3, fontsize=fontsize)
    plt.tight_layout()

    # プロットを保存　- テキストボックスを含めるために bbox_inches　を使用：
    # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1
    plt.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True)

    # 極座標プロットを表示
    plt.show()

まず、ルールオブファイブでフィルタリングしたデータセットをプロットします。

plot_radarplot(stats_rof, "../data/T2/radarplot_rof.png")

f:id:magattaca:20200412173555p:plain

上で作られたレーダチャートで、青色の四角は物理化学特性がルールオブファイブの基準の中に収まっている領域を表します。青線はフィルタリングしたデータセットの平均値を結んだもので、点線は標準偏差を結んだものです。平均値はリピンスキーのルールのどれも逸脱していないということがわかります。ですが、標準偏差をみると、いくつかの特性値はまだ閾値よりも大きくなっています。これは許容範囲内です。４つの特性のルールのうち一つは逸脱する可能性があるということを覚えておいてください。

次に、ルールオブファイブを逸脱した化合物を眺めます。

plot_radarplot(stats_not_rof, "../data/T2/radarplot_not_rof.png")

f:id:magattaca:20200412173624p:plain

ルールオブファイブに従っていない化合物の特性のほとんどが、LogPと分子量によってルールを外れていることがわかります。

ディスカッション

リピンスキーのルールオブファイブは経口吸収性（バイオアベイラビリティ）に焦点をあてています。医薬品は経口以外の投与ルートでも投与可能です（例：吸入、経皮、注射）。ルールオブファイブは経口投与における生物学的利用可能性の指標であって、例外はあるのだということに注意してください。バイオアベイラビリティを題材に、いくつかのADME特性の一つを見てきました。

ADME特性の全体像を手にするために利用可能なwebサーバーとプログラムがあります。　例： SwissADME.

クイズ

  • ルールオブファイブで表現された化学特性はADMEに対してどのような影響を与えうるでしょうか？
  • ３つ、あるいは４つのルールを逸脱するような分子を見つけるか、デザインしてください。
  • 上で作成したレーダーチャートに新しい化合物の情報を追加でプロットするにはどうすればよいでしょうか？

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
ところでレーダーチャートの利用と書き方を初めて見ました。
設定項目が多くて面倒ですが一つのプロットで複数のパラメーターが確認できるという点では便利そうですね。

今回のトピックに関連して日本語で読める記事をいくつかご紹介いたします。

Drug likenessの指標について:
* 化学の新しいカタチ：色々な薬らしさの指標
* Druglikenessについてのよもやま話

RDKitによる記述子計算:
* 化学の新しいカタチ：RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ

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

TeachOpenCADD

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

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

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

github.com

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

以下、訳。

トークトリアル１

化合物データの取得 (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サービスのリソースを示し、二つのリソースを繋ぐラインは、それらが共通の特徴を共有していることを示します。矢印の向きはリソースタイプについての一次情報を見つけられる場所を示しています。点線は２つのリソースの関係が異なる形式で振舞うことを示しています。例えば、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

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

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

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

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のため表示できないといったエラーが出ました。コメントアウトし、残りの３つを表示させています。
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のデータを１行ずつ書き込み
    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の条件ではand、or、notの代わりに&、|、~を使うことでエラーを回避できるそうです。参考
訳注ここまで

いくつか行を削除しましたが、あとでインデックスにしたがって反復処理（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のデータを１行ずつ書き込み
    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で出力ごとダウンロードできるの知らなかった！便利！

はじめに ~TeachOpenCADD 環境構築~

TeachOpenCADD

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

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

github.com

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

以下、訳。

TeachOpenCADD

TeachOpenCADD: オープンソースのパッケージとデータを使ったコンピューター支援（CADD）ドラッグデザインの教育用プラットフォーム

Dominique Sydow, Andrea Morger, Maximilian Driller and Andrea Volkamer

In Silico Toxicology, Institute for Physiology, Universitätsmedizin Berlin, Virchowweg 6, 10117 Berlin

DOI DOI Binder

質問がある場合は teachopencadd@charite.de に連絡してください。もしくは、このページにない講義用チュートリアルについての情報、チュートリアルのもととなるアイデアがあればぜひ連絡してください。

このページは　CC BY 4.0 ライセンスののもとで提供されています。ライセンスの詳細は　http://creativecommons.org/licenses/by/4.0/　をご覧ください。

目次

  • この教育用プラットフォームの目的
  • トピック
  • 使用方法
  • 連絡先
  • ライセンス
  • 引用

この教育用プラットフォームの目的

コンピューター支援ドラッグデザイン (CADD)にとって、ケモインフォマティクスと構造バイオインフォマティクスに関するオープンソースのプログラミングパッケージは、パワフルなツールで、モジュール式かつ再現可能、簡単に共有可能なパイプラインを提供します。ところが、ツールの説明、ドキュメンテーションは利用可能な一方で、CADDの活用方法について焦点をあて、背景にあるコンセプトを教える上で、自由に利用可能な具体例はほんの少ししかありません。例を挙げれば TDT initiative と言ったものがあります。 [1]

ここでは、学生による学生のために開発されたCADDの教育用プラットフォームを提供します。使用するオープンソースパッケージはPythonツール、RDKit [2] 、PyPDB [3] 、そして PyMol [4] です。それぞれのトピックについて、インタラクティブなJupyter Notebook [5] を開発しました。Githubで自由に利用可能で、トピックに関する理論的な背景の詳細な説明と、詳しいコメント付きのPythonコードが書いてあります。EGFR キナーゼを例に、公共の化合物データベース ChEMBL [6] からデータを取得する方法、ドラッグライクネスの観点から化合物をフィルタリングする方法、好ましく無い部分構造をもつ化合物を探し出し取り除く方法を議論します。さらに、化合物の類似性を評価する方法を導入し、化合物のクラスター化を行い、多様性のある化合物セットをつくります。また、新規な活性化合物を予測するためのモデルを構築するため、機械学習のアプローチを導入します。最後に、公共のタンパク質データベース PDB [7] から構造を取得し、それらを使ってリガンドベースの3D ファーマコフォアの生成、オフターゲットを予測するために結合サイト比較を実施します。

このプラットフォームを作成したことの目的は、オープンソースのケモインフォマティクスツールに興味を持っている利用者に、いかに簡単に利用できて、どんな利点があるかということを紹介することです。取り上げるトピックは継続的に増やしていきますし、コミュニティからのコントリビューションもお待ちしています。教育目的の利用だけでなく、この Jupiter notebookは各プロジェクトに合わせて修正し拡張するための出発点としても使えます。

文献:

[1] S. Riniker et al., F1000Research, 2017, 6, 1136 [2] G. Landrum, RDKit [3] W. Gilpin, Bioinformatics, 2016, 32, 156-60 [4] The PyMOL Molecular Graphics System, Version 1.8, Schrödinger, LLC [5] T. Kluyver et al., IOS Press, 2016, 87-90. [6] A. Gaulton et al., Nucleic Acid Res., 2017, 40, D1100-7 [7] H. Berman et al., Nucleic Acid Res., 2000, 28, 235-42

トピック

f:id:magattaca:20200405221750j:plain
TeachOpenCADD文献（D. Sydow et al., J Chem, 2019, 11, 29）Figure 1より

TeachOpenCADDのトピックはいわゆる トークトリアル（talkatorial） の形式（理論とコードのミックスで、プレゼンテーション形式での発表にも利用可能な形式）で取り上げられていて、現在、以下の内容で構成されています：

ケモインフォマティクス :

  • Topic 1. 化合物データの取得：ChEMBL
  • Topic 2. 化合物フィルタリング：ADMEとリードライクネスのクライテリア
  • Topic 3. 化合物フィルタリング：好ましく無い部分構造
  • Topic 4. リガンドベーススクリーニング: 化合物類似性
  • Topic 5. 化合物クラスタリング
  • Topic 6. 最大共通部分構造（Maximum common substructures）
  • Topic 7. リガンドベーススクリーニング：機械学習

構造バイオインフォマティクス：

  • Topic 8. タンパク質データの取得：Protein Data Bank (PDB)
  • Topic 9. リガンドベースファーマコフォア
  • Topic 10. 結合サイトの類似性
  • Topic 11. オンライン API/サーバー を使用した構造ベースのCADD
    • 11a. 潜在的なキナーゼ阻害剤を見つけるためのKLIFSとPubChemのクエリ化
    • 11b. 標的に対して候補のドッキングを行う
    • 11c. 結果の可視化と既知のデータとの比較

教育用の素材は次のフォーマットで提供されています：

  • コーディングベースのJupyter Notebook（topics 1-11）はGithubのこのページにあり、いわゆるトークトリアル（talktorial = talk + tutorial）、つまりプレゼンテーションでも使えるチュートリアルの形式です。
  • GUIベースのKNIME ワークフロー (topics 1-8) は KNIME Hub で手に入ります。

使用方法

レポジトリと必要な依存関係にあるものを全てホストするためにBinderを使うことができます。あるいはローカルでトークトリアルを使うこともできます（レポジトリをダウンロートして、依存関係にあるものをインストールしてください）

Binderの使い方 (対象トークトリアル T1-T8)

以下のリンクに従うことでレポジトリと必要な依存関係にあるものを全てホストするためにBinderを使うことができます:

Binder

セットアップには２、３分かかります。

ローカル環境へのインストール (対象トークトリアル　全て)

Linux

  • 1（トークトリアルを含む）TeachOpenCADDレポジトリのコピーをローカルに落としてください

    1. ... zipアーカイブとしてダウンロードしてunzipしてください： f:id:magattaca:20200405221951p:plain

    2. ... もしくは git パッケージを使って自分のコンピュータにクローンしてください：

git clone https://github.com/volkamerlab/TeachOpenCADD.git

  • 2 Anacondaソフトウェアを使ってパッケージのバージョン管理をクリーンにしてください。

    Anaconda2、あるいはAnaconda3をインストールします。理論上はインストールするAnacondaのバージョンの問題です。ですが、Anaconda2でしかテストしていないので、Anaconda3で同じような動作となるかは保証できません。

    https://docs.anaconda.com/anaconda/install/

  • 3 パッケージ管理システムcondaを使ってトーカトリアルの環境を作ります（ teachopencadd という名前をつけています）

    必要なパッケージをすべて含んでいる環境ファイル (yml file) を提供しています。

conda env create -f environment.yml

Note: この環境を手動でつくることもできます。 "別の方法：conda環境を手動で作る" を参照してください。

  • 4 作成したconda環境をアクティベートしてください。
conda activate teachopencadd

これでconda環境の中で作業することができます。

  • 5 作成したconda環境をJupyter notebook にリンクしてください
python -m ipykernel install --user --name teachopencadd

# FYI: このリンクをアンインストールするには以下のコマンドを使用してください:
#jupyter kernelspec uninstall teachopencadd
  • 6 Jupyter notebookを起動してください。
jupyter notebook

  • 7 メニューからJupyterのカーネルを作成したconda環境に変えてください。

    Kernel > Change kernel > Python [conda env:teachopencadd]

    f:id:magattaca:20200405222419p:plain

  • 8 これで最初のトーカトリアル開始です。Enjoy!

別の方法：conda環境を手動で作る

Note: これは上のステップ３で書かれているymlファイルを使用してconda環境をつくる方法の代わりとなる方法です。

# `teachopencadd` という名前の環境をつくりアクティベートします
conda create -n teachopencadd python=3.6
conda activate teachopencadd

# conda経由でパッケージをインストールします
conda install jupyter  # ipykernelもインストールします
conda install -c conda-forge rdkit  # umpy と pandas もインストールします
conda install -c samoturk pymol  # Linux: freeglut と glewもインストールします
conda install -c conda-forge pmw  # Pymolのターミナルウィンドウを立ち上げるのに必要です。
conda install -c conda-forge scikit-learn  # scipyをインストールします
conda install -c conda-forge seaborn  # matplotlibをインストールします
conda install -c conda-forge chembl_webresource_client
conda install -c conda-forge biopandas
conda install -c conda-forge pypdb
#conda install jupyter ipykernel pandas scikit-learn seaborn  # おそらくすでにインストールされていると思います。

# pip経由でパッケージをインストールします (おそらくデフォルトインストールされています）

pip install chembl_webresource_client biopandas pypdb

Windows

大枠では、LinuxとWindowsで同じように問題なくインストールできます。

唯一の違いはPyMolのインストール作業です。

cd C:/Anaconda2
set path=%path%;C:\Anaconda2
pip install pymol‑2.3.0a0‑cp36‑cp36m‑win_amd64.whl
pip install pymol_launcher‑2.1‑cp36‑cp36m‑win_amd64.whl
# 上の文言はダウンロードした pymol_XXXX.whl ファイルに合わせてください

MacOS

Linux環境と同様にインストールできますが、pymol をオープンソース samoturk condaチャンネルからインストールすることはできませんでした。代わりにschrodinger チャンネルを使うことができます。残念ながらPyMOLを実行するにはSchrödingerライセンスが必要です（ライセンスは 教育目的 の使用には無料で提供されています）。

インストールを手動で行うには、次のコマンドを：

conda install -c samoturk pymol # Installs also freeglut and glew

次に置き換えてください

conda install -c schrodinger pymol  # Installs also freeglut and glew

出くわした問題

RDKit
  • from rdkit.Chem.Draw import IPythonConsole がエラーメッセージ ImportError: No module named 'PIL'を投げる: pip install pillow のインストールを試してください(https://github.com/rdkit/rdkit/issues/1179 を確認してみてください)
  • Draw.MolsToGridImage が Serial errorをだす: データのタイプを Draw.MolsToGridImage(list(df.ROMol), useSVG=True) のようにリストしてください。
pip

バージョン　10.0.X 以降のパッケージマネージャーは以前のものと動作が異なっています。上で使われているコマンドのpip install chembl_webresource_clientは importエラー cannot import name mainを返します:

この問題はpython経由でpipを呼び出すことで解決することができます：

python3 -m pip install chembl_webresource_client
python3 -m pip install pypdb

chembl_webresource_client をインストールする際に、 urllib3 の互換性がないという警告が出るかもしれません。アップデートすることで解決することができます：

python3 -m pip install urllib3 --upgrade

解決方法は次からとりました。:

https://stackoverflow.com/questions/28210269/importerror-cannot-import-name-main-when-running-pip-version-command-in-window

連絡先

質問や提案があれば連絡してください！

皆様からのご意見をお待ちしています！

ライセンス

この仕事は Attribution 4.0 International (CC BY 4.0) のもとでライセンスされています。このライセンスのコピーをみるには次のリンクを参照してください http://creativecommons.org/licenses/by/4.0/ 。

引用

TeachOpenCADDプラットフォームの著者らは次の出資者からの公的出資を受け取っています:

  • Bundesministerium für Bildung und Forschung (Grant Number 031A262C)
  • Deutsche Forschungsgemeinschaft (Grant number VO 2353/1-1)
  • HaVo-Stiftung, Ludwigshafen, Germany
  • Stiftung Charité (Einstein BIH Visiting Fellow project)
  • "SUPPORT für die Lehre" program (Förderung innovativer Lehrvorhaben) of the Freie Universität Berlin
  • Open Access Publication Fund of Charité – Universitätsmedizin Berlin

サイエンティフィックな文献で TeachOpenCADD の教材を利用した時は、我々の文献を引用してください：

TeachOpenCADDプラットフォームの有用さを計り、今後の出資をうけるために役立ちます！

@article{TeachOpenCADD2019, author = {Sydow, Dominique and Morger, Andrea and Driller, Maximilian and Volkamer, Andrea}, doi = {10.1186/s13321-019-0351-x}, journal = {J. Cheminform.}, number = {1}, pages = {29}, title = {{TeachOpenCADD: a teaching platform for computer-aided drug design using open source packages and data}}, url = {https://doi.org/10.1186/s13321-019-0351-x}, volume = {11}, year = {2019} }

@article{TeachOpenCADDKNIME2019, author = {Sydow, Dominique and Wichmann, Michele and Rodr{\'{i}}guez-Guerra, Jaime and Goldmann, Daria and Landrum, Gregory and Volkamer, Andrea}, doi = {10.1021/ACS.JCIM.9B00662}, journal = {J. Chem. Inf. Model.}, publisher = {American Chemical Society}, title = {{TeachOpenCADD-KNIME: A Teaching Platform for Computer-Aided Drug Design Using KNIME Workflows}}, url = {https://pubs.acs.org/doi/full/10.1021/acs.jcim.9b00662}, year = {2019} }

訳以上

追記

以上でREADMEの訳となります。環境構築に関して私の場合、MacBook Pro(2017) macOS High Sierraを使用しており、conda環境での構築を試しています。
他のOSでのインストールやyml、Binder（って何かわかりませんでした・・・）を試していないので、不具合の有無はわかりません。
ソフトウェアのバージョンアップはそれぞれまちまちなので、場合によっては互換性がなくなってしまっているものもあるかもしれません。
ぜひぜひ情報をシェアしてください。（…このページについて本家の連絡先には問い合わせないでください）

素人の試訳なので誤訳、間違い等あると思います。ご指摘いただければ幸いです。

・・・いつかGitHubにちゃんとjupyter notebookをあげたい。

オープンソースで医薬品探索のコンセプトを勉強しよう！〜TeachOpenCADD〜

TeachOpenCADD

TeachOpenCADDをご存知でしょうか？ざっくり言うとCADDの自習用教材です。*1

f:id:magattaca:20200405201048p:plain:w300
TeachOpenCADD

CADD（Coputer-Aided Drug Design、コンピュータ支援医薬品設計）というのは、 その名の通り「医薬品の探索、ドラッグデザインにコンピュータを利用する」ということです。（≒ in silico創薬？）

医薬品探索における計算機利用、シミュレーションの重要性とともにその門戸はどんどんと開かれており、
1. 計算機性能の向上
2. 利用可能な（オープンソース）ソフトウェアが増えた
3. 公開データベース
などがその理由としてよく挙げられています。

それぞれ説明書もあって道具は揃った！よし、あとは実際につかうだけ！・・・とはいかないですよね。
とりわけ私のような初学者にとっては「車あるけどカーナビもスマホも無い。どっち行けばいいの？？」状態です。

ツールと実践の間のハードルをクリアするためには、
1. そもそもなぜそのツールが必要とされたのか、背景にあるコンセプトを学ぶこと
2. 実践的な具体例でツールを利用してみること
が必要です。

つまり、あれです。「先達・・・あらまほし・・・」ってやつですね。

で、そんな迷える仁和寺の法師たちのためのプロジェクトがこちらTeachOpenCADD platformです。

TeachOpenCADD?

TeachOpenCADDの目的は「CADDに馴染みのない学生、研究者のために、そのコンセプトと実践を学習するための教材を提供する」ことで、
チュートリアル形式で一歩ずつ学習を進められるように設計された、学生による学生のためのプラットフォームとなっています。

開発の主体となっているのはシャリテ（The Charité、ベルリン大学病院）のVolkamer研究室で、 構造ベースの手法による毒性予測やオフターゲット予測を主な研究テーマとされているようです。

TeachOpenCADDの特徴として、

  • 無料で利用可能なソースを使用
  • Interactiveに学習できるJupyter notebook形式 (GitHubに公開)

となっており、「EGFRキナーゼ阻害剤の探索」を具体例にしてリガンドベース、構造データベースのCADDの手法を順に学ぶことができるようになっています。

f:id:magattaca:20200405202008p:plain
TeachOpenCADDで学べること（J Cheminform(2019)11:29 Fig1.より）

オープン（Open）ソースでCADDを教える(Teach)教材、TeachOpenCADDということのようですね。

各チュートリアルは自習だけでなく講義で説明しながら使えるようにもなっており、 トークトリアル（talk + tutorial = talktorial)とも呼ばれています。

f:id:magattaca:20200405202310p:plain
実際のチュートリアル（J Cheminform (2019)11:29 Fig2）

以上は私の雰囲気翻訳なので正確な情報はオリジナルを参照してください。

  1. TeachOpenCADD Jupyter Notebooks: Talktorials 1-10.
  2. TeachOpenCADD KNIME workflows
  3. GitHub

EGFRキナーゼ阻害剤

さて、TeachOpenCADDではチュートリアル 10回（GitHubの最新版は11回まで）を通してCADDのコンセプトを学びますが、 一貫して取り上げられている題材「EGFRキナーゼ阻害剤（EGFR-TKI）」というのはどういうものなのでしょうか？

まずは敵を知る必要がある！！*2

Wikipediaによると、EGFR（上皮成長因子受容体、Her1）は細胞膜上の受容体で、膜貫通型糖タンパク質です。細胞内にチロシンキナーゼ活性を有するドメインがあり、細胞外からの刺激を細胞内でリン酸化シグナルとして伝達します。このシグナルは細胞の分化・増殖に関わっており、したがって機能変異により発がん、がんの増殖・浸潤・転移に関与するようになります。*3
EGFR遺伝子の変異は肺がんで多く見られるドライバー変異で、日本人の場合非小細胞肺がん患者で多く変異が見られているそうです。*4

変異による過剰な増殖シグナルをとめることで抗がん作用が期待できるということで、細胞外（抗EGFR抗体）および細胞内（EGFRキナーゼ阻害剤）の２つのアプローチで医薬品の開発、上市されています。

EGFRキナーゼ阻害剤といえばやはり一番最初に思い当たるのはゲフィチニブではないしょうか？薬学においてゲフィチニブは様々な点で取り上げられる医薬品です。

  • 初期の分子標的薬
  • 優先審査により世界に先駆けて日本で承認
  • 副作用（イレッサをめぐる問題）
  • 遺伝子変異と薬効の関係

など様々な文脈で名前を聞いたように思います。*5

すでに第3世代までが承認されているEGFRキナーゼ阻害剤ですが、耐性発現、既存の薬剤で効果のでない変異があるなど依然としてアンメットニーズがあるそうです。*6 また、EGFR陽性の肺がんでは免疫チェックポイント阻害薬の効果が出にくいという報告もあるそうで、 EGFR阻害剤最初の承認から20年近くたった今でもまだまだ研究する余地がありそうです。*7

で、私は何がしたいのか？

Q: ぐだぐだと書いているが結局何がしたいのか？
A: 一人で勉強するのが寂しい

ということでのんびり翻訳して晒していこうと思います。 ライセンスCC BY 4.0だから怒られないはず・・・

ではでは！

~~ (追記 2020/09/28) ~~~

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

github.com

*1:文献. TeachOpenCADD Jupyter Notebooks: Talktorials 1-10.

*2:私は医療従事者でも専門家でもないので記述の正確性に問題があることをご承知おきください。

*3:上皮成長因子受容体_Wikipedia

*4:日本内科学会雑誌2014(103)1355

*5:このあたり付け焼き刃の知識で踏み込めないので興味を持たれた方は調べてみてください。

*6:【肺がん】EGFR-TKI 一次治療に相次ぎ新薬…それでも残るアンメットニーズ

*7:名古屋大学プレスリリース

pymolをimportしようとする話の続きの話

PyMOLのjupyter notebookからのimportに失敗したことについて記事を書いたところ、早速twitterで解決方法を教えていただきました。 @kira_kira_worldさん、@Ag_smith先生ありがとうございました！

解決策に入る前に・・・pythonのバージョンについて

先の記事を書いた時点では気づいていなかったのですがどうもpythonのバージョンに原因があったようです。

私の環境はmacOS High Sierraで

  • Homebrew > pyenv > Anaconda（複数の仮想環境）

みたいな感じになっています。（たぶん）

で、全体（pyenv global）はPython 3.70となっています。 この状態でhomebrewからPyMOLをインストールしたためか、importしようとしたPyMOLは以下のディレクトリ 「/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages」 を呼びに行っていました（先の記事でpathを通した所）。
如何にもpython3.7に紐づいていそうな場所です。

また、別の話として、色々とソフトウェアをインストールするためにAnacondaで複数の仮想環境を作成し、 jupyter notebookのカーネルに追加、切り替えつつ使おうとしています。どうやらここに問題がありそうでした。

先の記事で使用した仮想環境（カーネル）のpythonのバージョンを確認したところPython 3.6.10となっていました。 そこで別のPython 3.70の仮想環境に切り替えたところ、エラーは出ずにimport pymolすることができました。（pathを通す必要はありました）

おそらくPythonのバージョンをきちんと合わせておけば問題なく使用できるのではないでしょうか？ また確認不足で誤情報を垂れ流してしまった。。。すみません。

解決策1）シンボリックリンクの作成

ではでは、教えていただいた方法をご紹介したいと思います！
まずは@Ag_smith先生に教えていただいた方法です。

先の記事では怪しげなファイル「_cmd.cpython-37m-darwin.so」をコピペして「_cmd.so」に名前を変えるという荒技を行いました。 これだと思わぬ不具合を起こす可能性があるということで、こういう場合はシンボリックリンクを作成すれば良いそうです。

シンボリックリンクはUNIX系のOS（MacやLinux）で、Windowsのショートカットに相当するものだそうです。*1*2

つまり「_cmd.cpython-37m-darwin.so」のシンボリックリンクとして「_cmd.so」を作成しておけば、 import pymolはシンボリックリンクを中継として本体のファイルに辿りつくことができる、ということのようです。

作成方法は簡単で、ターミナルで以下を実行すればOKです。

# シンボリックリンクを作成したいファイルが入っているディレクトリに移動
cd /usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages/pymol

# 「ln -s 本体のファイル　リンクの名前」でシンボリックリンクを作成
ln -s _cmd.cpython-37m-darwin.so _cmd.so

できました！

f:id:magattaca:20200321225101p:plain

だそうです。「種類：エイリアス」となっていますがエイリアスとシンボリックリンクは微妙に機能が違うらしい・・・。 このあとjupyter notebookで（エラーの出ていたpython 3.6の仮想環境でも）import pymolできることを確認しました。

スマート！

解決策2）インストール方法を変える（Anaconda CloudのSchrödinger社のチャンネルを利用）

さて、もう一つの解決方法は@kira_kira_worldさんに教えていただいたもので、 condaでSchrödinger社のチャンネルからPyMOLをインストールするというものです。

１行でOK!

conda install -c schrodinger pymol

ですが注意点があります。こちらはSchrödinger社が配布しているものだということです。

Wikipediaによると、 もともとPyMOLはWarren L. DeLano氏により開発されたオープンソースのソフトウェアで、 現在はSchrödinger社によって商品化されているとのことです。 ですので、上記の方法で入手可能なPyMOLはライセンスに注意が必要そうです。

私はライセンス等に詳しくなく判断ができないので、以下手に入る情報を共有させていただきます。（←コンプライアンスを意識した発言）

まずはSchrödinger社のPyMOLのページ。
2020年03月20日現在、こちらではPyMOL 2.3のダウンロードが可能です（Python 3.7を含むそう）。 Windows、macOS、Linuxそれぞれのインストーラ、および上記のAnaconda Channelを利用する方法が提供されています。 また、ライセンスの購入方法についてもリンクがあります。

さらにスクロールしていくとOpen-Source Philosophyの文字が・・・
Schrödinger社が商用として維持、開発しているものの、元々のDelano氏の意向を引き続き大枠をオープンソースプロジェクトとしているようです。 したがってほとんどのソースコードはGitHub上で公開されています。ページのリンクはこちら。

Homebrew経由でインストールを行った際にも、ソースコードは上記のGithubから取得しているようです。 @Ag_smith先生の記事に記載されていて、 GitHubのHomebrewのFomulraでも取得先URLを確認できます。

では、Schrödinger社コンパイル済みのものをインストールした場合どうなるのでしょうか？試してみます。

テストのための仮想環境を作成し、condaで取得します。ちなみにAnaconda CloudのPyMolページはこちら。

# pythonのバージョンを指定して仮想環境pymoltestを作成
conda create -n pymoltest python=3.7

# 仮想環境をactivate
conda activate pymoltest

# PyMOLをインストール
conda install -c schrodinger pymol

# PyMOLを起動
pymol

無事PyMOLのウインドウが立ち上がりました！

f:id:magattaca:20200321225337p:plain

ただしHomebrew経由のインストールとは明確な違いが！

同時にActivationのウィンドウが立ち上がります。
「オフィシャル版なのでライセンスを持っているならactivationしてくださいね」ということのようです。 ライセンスがない場合、評価版として30日間使用することができるようです（skip activation）。

評価期間がすぎると「ウィンドウにNo License Fileの表示がされ使えなくなる機能がある」、とのことですが もうすでにNo License Fileの文字が見えますし、30日すぎるとどうなるのかはまだわかりません。

ちなみに学生であればライセンスは無料で取得することができるそうです。が、私はとうに学生ではなくなってしまったのでとりあえずskip activationで先に進めたいと思います。

さて、こちらのPyMOLがjupyter notebookでimportできるか検証したいと思います。 まず準備としてcondaの仮想環境をjupyter notebookに認識させます。

# 以下ターミナルで実行
# jupyter notebookのカーネルに仮想環境を追加
python -m ipykernel install --user --name pymoltest

# 一旦、仮想環境を抜けます
conda deactivate

また、インストールしたPyMOLの場所を確認するため、PyMOLのコマンド枠で以下を打ち込みます。

import sys
print(sys.path)

今回、それっぽいpathは「'/Users/XXXXX/.pyenv/versions/anaconda3-5.3.1/envs/pymoltest/lib/python3.7/site-packages'」でした（XXXXXは実際にはユーザーアカウント名ですが恥ずかしいので隠しています）。 condaの仮想環境pymoltestの中に入っていることがわかりやすいですね。

準備できたのでjupyter notebookを起動、新しいnotebookを作成してChange kernelでkernelをpymoltestに切り替えます。

f:id:magattaca:20200321225534p:plain

前回同様、pathを通してpymolをimportしてみます。

import sys
sys.path.append('/Users/XXXXX/.pyenv/versions/anaconda3-5.3.1/envs/pymoltest/lib/python3.7/site-packages') # XXXXXは実際にはユーザーアカウント名

import pymol

# importしたpymolのpathを確認
print(pymol.__path__)

無事importすることができました!!　PATHも指定したものとなっていました！

ちなみにpython 3.6の仮想環境のkernelで上記を実行したところ、前回記事同様に「モジュールpymol._cmdが見つからない」というエラーがでてimportできませんでした。 やはりバージョンの確認が大事そうですね。

何はともあれconda経由でもPyMolのインストールとCUIの利用が可能であることが確認できました！

目的は達成できたので、とりあえずオープンソース版での利用することとして仮想環境ごと削除しておきます。

# 以下ターミナルで実行
# jupyter notebookのカーネルからアンインストールしておく
jupyter kernelspec uninstall pymoltest

# condaの仮想環境を削除
conda remove -n pymoltest --all

ターミナルを一度落としてから再度立ち上げpymolを実行することでオープンソース版のPyMOLが立ち上がることを確認できました。

_cmdはどこからきたのか?

さて、jupyter notebookからPyMOLをimportする際の不具合は

  • Pythonのversion違いである可能性がある

こと、解決策として

があることまでわかりました。

で、今回importエラーとなった「pymol._cmd」とは何者か？もう少し調べてみたいと思います。

Github Open-Source PyMOL上でオープンソースのソースコードを確認してみましたが、 ディレクトリmodules/pymolの中にcmd.pyはあるものの_cmdに相当するものはありませんでした。

ビルド（？）する際に作成されたのかな？ということでsetup.pyの中身を確認したところ以下のコードに行き当たりました。

from distutils.core import setup, Extension
# ~ 中略 ~
ext_modules += [
    Extension("pymol._cmd",
              get_sources(pymol_src_dirs),
              include_dirs = inc_dirs,
              libraries = libs,
              library_dirs = lib_dirs,
              define_macros = def_macros,
              extra_link_args = ext_link_args,
              extra_compile_args = ext_comp_args,
              extra_objects = ext_objects,
    ),
# ~以下略~

distutils はPythonのモジュール配布のためのユーティリティで *3、Extensionは「C/C++拡張モジュールをセットアップスクリプトで表すために使われ」るそうです*4。

Extensionの説明と引数(argument)の内容からソースコードからビルドに伴って作成されたC/C++関連のファイルということがなんとなくぼんやりとわかりました。 これを踏まえて「_cmd.cpython-37m-darwin.so」というファイル名を再度眺めてみます。

前回邪魔なのではと削除してしまった「cpython-37m-darwin」の部分ですが、こちらはCythonでコンパイルする際に追加される接尾辞(suffix)のようです*5。 CythonはPythonの「処理速度が遅い」という問題を解決するためにC/C++に変換するもの（？？？？）で*6、 また、darwinはmacOS（Mac OS X）のベースとなっているオープンソースのオペレーティングシステムDarwinのことのようです*7。

なんとなくですがMacでpython 3.7でC/C++関連のコンパイルを行うことで作成されたファイルという感じでしょうか？（・・・浅い理解）

python 3.7に紐づいているからpython 3.6の仮想環境ではうまく動かなかったのかもしれませんね！（平然と誤情報を垂れ流すよ！!）

まとめ

以上、前回に引き続きPyMOLのCUIを利用するための方法についてでした。 適当な記事を書くとすぐに解決策が教えていただけるのでオープンソースのコミュニティに感謝感謝です。ありがとうございました！！

今回も適当調べで誤情報多そうなのでご指摘いただければ幸いです。ではでは！

*1:ソフトリンク

*2:シンボリックリンク (symbolic link)

*3:Distutilsの紹介

*4:Extensionクラス

*5:python 3 cython compiled files extension

*6:深入りしないCython入門

*7:macOS_Wikipedia

JupyterNotebookでpymolのimportに失敗した話

PyMOL便利ですよね！でっかいタンパク質の描画も綺麗で動きも速いですし。

便利なPyMOLがCUIでも使えるらしいと知り、初心者の味方 jupyter notebook で使おうとしたところいろいろと躓いたので備忘録を残します。

まず結論！

  • Pathを通して
  • ファイル名を変更「 _cmd.cpython-37m-darwin.so　 → _cmd.so 」

した！という話をだらだら書きます。

03/20追記
ファイル名変更よりも_cmd.cpython-37m-darwin.soのシンボリックリンク_cmd.soを作成する方が安全とのことです。 @Ag_smith先生ありがとうございました！

PyMOLの状態

PyMOLは @Ag_smith 先生の記事macOS/CentOS 7/Ubuntu 18.04へのオープンソース版PyMOLのインストール方法 にしたがってインストールしました。ちなみにパソコンはMacBookProです。

brew tap brewsci/bio
brew install pymol

で、ターミナルでpymolと打ち込むとウィンドウが立ち上がることを確認しています。

f:id:magattaca:20200320162949p:plain

Jupyter notebookでimport error

早速 import ！ そしてエラーが出た！

import pymol

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-f943b406a205> in <module>
----> 1 import pymol

ModuleNotFoundError: No module named 'pymol'

「pymolなんて名前の奴はない」と怒られました。・・・あるよ。。。

PATHを確認

色々みているとPATHが通っていないのがよくないのでは、ということだったので、 まずはPyMOLがどこにインストールされているか確認して見ることにしました。 （コピペでインストールしてるからソフトの居場所も分からない。）

立ち上げたPyMOLのウィンドウで、コマンドを打ち込む枠（「PyMOL>」と書いてあるところ）に以下を貼り付けました。

import sys
print(sys.path)

以下が表示されました。

['', '/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages', '/usr/local/Cellar/python/3.7.5/Frameworks/Python.framework/Versions/3.7/lib/python37.zip', '/usr/local/Cellar/python/3.7.5/Frameworks/Python.framework/Versions/3.7/lib/python3.7', '/usr/local/Cellar/python/3.7.5/Frameworks/Python.framework/Versions/3.7/lib/python3.7/lib-dynload', '/usr/local/lib/python3.7/site-packages']

HomeBrew経由でインストールしたものは「Cellar」というディレクトリの下に入るそうです。 どうも「'/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages'」が怪しそうです（pymolって書いてあるから）。 こちらを使うことにしましょう。

PATHを通すも新たなエラーが！

.bash_profile にPATHを追加すれば良いようですが、自信がないのでjupyter notebookでテストしてみます。

以下を打ち込み・・・

# PyMOLへのPATHを通す
import sys
sys.path.append('/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages')

import pymol

新たなエラーを得た！

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-3-f943b406a205> in <module>
----> 1 import pymol

/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages/pymol/__init__.py in <module>
    581 # _cmd = sys.modules['pymol._cmd']
    582
--> 583 from . import cmd
    584
    585 cmd._COb = None

/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages/pymol/cmd.py in <module>
    102
    103         import re
--> 104         from pymol import _cmd
    105         import string
    106         import threading

ImportError: cannot import name '_cmd'

今回は ModuleNotFoundError ではなかったので、先ほど追加したPATHで当たりだったようです。 PyMOLの居場所はわかったものの、今度は_cmdというのが見つからないらしい。

FinderでPyMOLのディレクトリを確認すると、似た名前の「cmd.py」というファイルはありましたが、「_cmd.py」は見当たりませんでした。 ひょっとして名前を付け間違えたか？と思い、「cmd.py」を複製し、「_cmd.py」にファイル名を変えてみましたが異なるエラーがでただけでした。

そんなアホなことでは直らないですよねー。アンダーバー、なんか意味ありそうですし・・・。

それっぽいファイルの名前を変えたら動いた！

色々と調べてもさっぱりわからず、もう一度PyMOLのディレクトリを見直すと怪しげなファイルがありました。

_cmd.cpython-37m-darwin.so

拡張子が「.py」でなく名前が長いので素通りしていましたが、最初に「_cmd」と書いてあります。 開いてみても文字化けしてさっぱり中身がわからないファイルでしたが、 「.so」というのはMacOSで使われる共有ライブラリの拡張子（"shared object"の略）だそうです。*1

意味がわかりませんでしたが、とりあえず重要っぽいファイルでファイル名が探しているものに似ている。 ・・・名前変えちゃえ！ということで、ファイルを複製して邪魔っぽい部分（.cpython-37m-darwin）を消し、 「_cmd.so」にしてみました。

で、もう一度jupyter notebookに戻りimportすると今度はエラーを出さずに実行できました！ エラーが出ないだけだと不安なのでPATHを確認してみます。

import pymol
pymol.__path__

['/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages/pymol'] と表示されました！ どうやら通したPATHも合っていたようです！

まとめ

そもそもPyMOLのCUIでどのような機能が使えるかわかっていないので、importはできているっぽいものの本当に正しく使えるのか？ ファイル名を変更してしまったけどよかったのか？等々不安が残ります。とりあえずの備忘録として残しておきます。

おそらく間違っているように思います。詳しい方ご教示いただければ幸いです。

結局 cpython-37m-darwin が何かよくわかりませんでした・・・

*1:.soファイル 【 shared object file 】