PROTAC DBを掘り起こしてみた
みなさん、PROTACしてますか？流行ってますね。
乗り遅れてる？大丈夫です。そんな簡単に出来ませんw
でも初動は大事だと思うので、公開されてるDBを参考にしながら計画を立てるのもいいんじゃないかなーって思います。
こちらの報告ではPROTAC自体はもとより、POIのリガンドだったりE3リガンドやLinkerの情報が載ってます。
多分、論文ベースとかなので特許や臨床情報などは網羅されてません。とはいえ、眺めてみて損はないと思います。
前述のようにPROTACの各パーツについて情報が記載されてますが、例えばPROTACのデータテーブル(というか多分Rowのこと)で、そのPROTACに使われたBuilding Blockがどれくらいの頻度で登場するか(=人気あるか)が分かりません。このあたりを分析できるようになるといいですよね。
データを読み込んでR分解
まずはダウンロードしてきたcsvファイルを読み込みます。パスは適切に直してください。
私の環境だと low_memory=Falseにしないとエラーが出ました。エラーは出たけど読み込んだデータに問題はなさそうでしたけど、念の為に指定してます。
import pandas as pd # load df df = pd.read_csv('./protac.csv', low_memory=False)
こんな感じのDataFrameになってると思います。
次にこのデータをR分解します。まずはR分解するための関数を定義しておきます。
2つ目の方はR分解したフラグメントをmolファイルのまま保存できます。重たくなるので今回は1つ目のSMILESの方を使ってます。
# smilesからfragmentに分解する関数 def fragment_smiles(smiles): # mol file生成 mol = Chem.MolFromSmiles(smiles) # RECAPによる分解 decomp = Recap.RecapDecompose(mol) # smilesを生成してリストに格納 children = [Chem.MolToSmiles(Chem.MolFromSmiles(x)) for x in list(decomp.children.keys())] return children def fragment_mol(smiles): # mol file生成 mol = Chem.MolFromSmiles(smiles) # RECAPによる分解 decomp = Recap.RecapDecompose(mol) # mol filesを生成してリストに格納 children = [x.mol for x in decomp.children.values()] return children # R分解したフラグメントを抽出してカウントしたdataframeを出力
次にdataframeにR分解したフラグメントのデータを追加していきます。
下記のスクリプトは100行ずつのチャンクにして処理しています。
メリットは少し計算が早くなったことと、計算時間が長いのでどこまで進捗したか把握できるところですかね。
チャンクにしないで .map()のコードだけでもいいと思います。
再度計算すると時間がかかるので、ここでcsvに保存することをお勧めします。
# 100行ずつに分割して処理を行い、df2に連結 df2 = pd.DataFrame() for i in range(0, len(df), 100): df_chunk = df.iloc[i:i+100].copy() # 100行ずつのチャンクを作成 df_chunk['children'] = df_chunk['Smiles'].map(lambda x: fragment_smiles(x)) df2 = pd.concat([df2, df_chunk], ignore_index=True) # df2に連結 print(len(df2)) df2.to_csv('./protac_decomp.csv', index=False)
ちなみにR分解したフラグメントは Datatypeがリストで出力されてるはずですが、保存するとブラケットで囲まれただけの文字列になってます。つまり pd.read_csv()で読み込んでも strとして読み込まれてしまいます。
R分解したdataframeは下記のコードで読み込むと、R分解したリストを復元できます。
# 保存したcsvを読み込む import ast # カスタム関数を定義して、文字列をリストに変換 def str_to_list(s): return ast.literal_eval(s) df2 = pd.read_csv('./protac_decomp.csv', low_memory=False, converters={'children': str_to_list})
分析
R分解が終わったので、次はフラグメントについて分析できるように加工していきます。
まずはどんなE3があるかを .unique()で確認。
df2['E3 ligase'].unique() # array(['VHL', 'CRBN', 'DCAF15', 'MDM2', 'cIAP1', 'XIAP', 'DCAF16', 'AhR', 'IAP', 'RNF4', 'RNF114', 'DCAF11', 'FEM1B', 'FBXO22', 'BRD4', 'Keap1', 'KLHL20', 'DCAF1', 'KEAP1', 'UBR box', 'KLHDC2'], dtype=object)
お馴染みのE3ももちろん含まれていました。E3がそれぞれどれくらい使われているかも確認しておきます。
df2.drop_duplicates(subset='Smiles')['E3 ligase'].value_counts() # E3 ligase # CRBN 3860 # VHL 1914 # cIAP1 94 # MDM2 57 # IAP 48 # XIAP 39 # DCAF16 32 # KEAP1 12 # KLHL20 9 # Keap1 8 # DCAF1 7 # FEM1B 7 # BRD4 5 # RNF114 4 # AhR 4 # DCAF11 3 # KLHDC2 3 # FBXO22 2 # DCAF15 2 # UBR box 2 # RNF4 1 # Name: count, dtype: int64
ぶっちぎりでCRBNとVHLでした。リガンドの合成や入手性も大きいのかなと思います。
次はCRBNでどんなリガンドがよく使われているかが知りたいです。
色々分析するための関数を先に定義しておきます。
# R分解したフラグメントを抽出してカウントしたdataframeを出力 def count(data): df = data # childrenカラムのリストをフラットにする all_children = [child for sublist in df['children'] for child in sublist] # 要素のカウント counter = Counter(all_children) # カウント結果をデータフレームに変換 count_df = pd.DataFrame(counter.items(), columns=['Child', 'Count']) count_df['MW'] = count_df['Child'].map(lambda x : Descriptors.MolWt(Chem.MolFromSmiles(x))) count_df = count_df[count_df['MW']>150] return count_df # count_dfにクラスタリング情報を追加 # https://note.yu9824.com/howto/modern-rdkit-fingerprint/ def clustering(data, split=3): count_df = data # Smilesを分子オブジェクトに変換 mols = [Chem.MolFromSmiles(smi) for smi in count_df['Child']] # 分子フィンガープrintを生成 (Morganフィンガープrintを使用) mfgen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=1024) fps = np.array(mfgen.GetFingerprints(mols, numThreads=1)) # K-meansクラスタリング (クラスタ数は3とする) kmeans = KMeans(n_clusters=split, random_state=0) cluster_labels = kmeans.fit_predict(fps) # クラスタリング結果をDataFrameに追加 count_df['clustering'] = cluster_labels return count_df # count_dfのtop20のフラグメントを描画 def cluster_img(data, num=20): top = data.sort_values(by='Count', ascending=False).head(num) topmol = [Chem.MolFromSmiles(smiles) for smiles in list(top['Child'])] return Draw.MolsToGridImage( topmol, molsPerRow=4, legends=[str(x) for x in top['clustering']], subImgSize=(300,200)) def count_img(data, num=20): top = data.sort_values(by='Count', ascending=False).head(num) topmol = [Chem.MolFromSmiles(smiles) for smiles in list(top['Child'])] return Draw.MolsToGridImage( topmol, molsPerRow=4, legends=[str(x) for x in top['Count']], subImgSize=(300,200))
count()はR分解したフラグメントをカウントしたdataframeを出力します。clustering()はカウントテーブルの中の構造を、フィンガープリント情報からクラスタリングします。第二引数に数値を入れると、その値にクラスタリングしてくれます。cluster_imgやcount_imgはカウントテーブルのトップランクの構造を描画します。デフォルトは20化合物分描画するようにしてあります。前者はクラスタリング番号、後者は出現頻度をlegendに記載しています。とりあえず使ってみましょう。
まずはE3 ligaseをCRBNに指定して、トップランクのフラグメントを確認します。
# E3ligaseを指定 df_crbn = df2[df2['E3 ligase']=='CRBN'] count_df = count(df_crbn) count_img(count_df)
lenalidomide型のリガンドが多いみたいですが、POIリガンドが紛れ込んでます。これをうまく除いてやるためにクラスタリングを使います。
clustering(count_df)
cluster_img(count_df)
# cluster 0がCRBNリガンドっぽい
ぱっと見でcluster 0がCRBNっぽいですね。2はキナーゼ？と有名どころのJQ1が入ってます。
clusterを0に指定して、トップランクの化合物を描画します。
count_df = count_df[count_df['clustering']==0].sort_values(by='Count', ascending=False) count_img(count_df)
分子量が150以上のものを表示するように関数を定義しているので、細かいPEGみたいなのは除いてある状態です。
やはりLenalidomide型のリガンドがよく使われてるみたいです。この位置はCRBNと相互作用を獲得しやすいんじゃなかったかな。
注意したいのは、今回のソースはおそらく論文ベースのデータということです。
TwitterでよくみるPROTACは最近は環状リンカーが多い気がしますが、論文ベースだとみんなPEGから入ってるっていうことなんだろうと想像してます。実際に、Arvinas社のARVシリーズはこのデータに入ってなかった気がします。
今回は分かりやすくE3 ligaseでリガンドを調べてみましたが、標的ごとにどんなリガンドがよく使われているか調べてみてもいいと思います。
気が向いたらリンカーの情報の紐付けもやってみたいと思います。
おしまい。
共結晶構造情報のリガンドが市販かどうかを調べるには
SBDDに限らず、創薬にあたってポジコンやツール化合物があるかどうかは評価系構築の確度に関わる部分です。
ここをしっかりしてもらうことが、我々ケミストの成功の第一歩になります。
若い子「おじさん、すみません。共結晶が報告されてるリガンドを調べてもらうことってできますか？」 おじさん「いいよ。」
数分後...
おじさん「結構な数があるね。」
若い子「これって買うか作るか出来ますか？」
おじさん「うーん調べるならちょっと時間ください」
よく知られてる標的だと結構な数の共結晶が報告されてますよね。
調べるだけで半日(~1日)かかっちゃいそうです。
そんなわけでpythonのスクリプトを書きました (>3日かかった)。
Contents
conda create -n bi -c conda-forge -y python=3.10 pymol-open-source=2.5 BioPython=1.81 pandas=2.0.3 requests=2.31 pypdb=2.3 biopandas=0.4.1 rdkit=2022.09.1 以降のスクリプトは 色んな方法が報告されています。例えば, 下記のスクリプト実行してキーワード検索の結果をリストにすることができます。 jupyterで実行すると 今回は@くろたんくさんの記事を参考にさせてもらいました。 添加剤やイオンなんかも結果に含まれてしまうので、必要に応じてフィルタリングが必要です。 あーやっとスッキリしました。環境構築
miniconda がインストールされているものとします。ちなみにMacでもWinでも動作確認済みです。
今回のスクリプトは下記の仮想環境で動作確認済みです。conda activate biで仮想環境をactivateする必要があります。
すでにinstallされてたら問題ないです。今回は使わないものも入っています。標的のPDBリストを作る
from pypdb import *
q = input('keywords? ')
pdb_ids = Query(q).search()
print('the number of PDB IDs are ' + str(len(pdb_ids)))
inputの入力を求められます。 今回は有名どころでCRBNと入力したところ、2023.9.4現在で81個のPDBがヒットしました。
他にも方法はあるので気が向いた時にでも公開します。PDBに含まれるリガンドをリストアップして、PubChemのSupplierリンクをつける
私の拙いcodingでは例外処理がめんどうで目的を達成できませんでした。感謝です。
早速ズバリ結論を書くと、上述のPDBリストを作った状態で下記を実行します。import requests
import pubchempy as pcp
import pandas as pd
# Ligandの有無を判定してdfを作成する関数
def PDB2L(entry_id):
ligand_list = []
query = '''
{
entry(entry_id:"''' + entry_id + '''") {
rcsb_id
struct {
title
}
nonpolymer_entities {
rcsb_nonpolymer_entity_container_identifiers {
entry_id
}
nonpolymer_comp {
chem_comp {
id
type
}
pdbx_chem_comp_descriptor {
descriptor
type
program
}
}
}
}
}
'''
url = "https://data.rcsb.org/graphql?query=" + query
response = requests.get(url)
json_data = response.json()
# dfの作成
# Ligandがある場合だけrun
if json_data['data']['entry']['nonpolymer_entities'] != None:
# LigandがUNLやUNKではない時はrun
if json_data['data']['entry']['nonpolymer_entities'][0]['nonpolymer_comp']['pdbx_chem_comp_descriptor'] != None:
# 複数のリガンドをdfに代入
for i in json_data.get('data').get('entry').get('nonpolymer_entities'):
entry_id = i.get('rcsb_nonpolymer_entity_container_identifiers')
nonpolymer_comp = i.get('nonpolymer_comp')
ligand_id = nonpolymer_comp.get('chem_comp').get('id')
for data in nonpolymer_comp.get('pdbx_chem_comp_descriptor'):
if (data.get('type') == "InChI") and (data.get('program') == "InChI"):
inchi = [data.get('descriptor')]
type = [data.get('type')]
program = [data.get('program')]
d = {'Ligand_ID':ligand_id, 'InChI':inchi}
ligand_list.append(d)
df = pd.DataFrame(ligand_list)
df['Title'] = json_data['data']['entry']['struct']['title']
return(df)
# ここからfor文
# pdb_ids = ['5UT2', '4FVR'] # for test
cols = ['Ligand_ID', 'SMILES', 'InChI', 'CID', 'MW', 'PDB_ID', 'Title', 'Supplier', 'RCSB_URL', 'Ligand_URL']
Sup_List = pd.DataFrame(columns = cols)
for pdb_id in pdb_ids:
# Ligandがない場合の例外処理
if judge(pdb_id) != None:
if UNL_UNK(pdb_id) != None:
# DataFrameの作成
df = pd.DataFrame(entry_to_ligand(pdb_id))
df['PDB_ID'] = pdb_id
df['CID'] = df['InChI'].map(lambda x : pcp.get_properties('MolecularWeight', x, 'inchi')[0]['CID']).astype(str)
df['MW'] = df['InChI'].map(lambda x : pcp.get_properties('MolecularWeight', x, 'inchi')[0]['MolecularWeight']).astype(float)
df['SMILES'] = df['InChI'].map(lambda x : pcp.get_properties('IsomericSMILES', x, 'inchi')[0]['IsomericSMILES'])
df['Title'] = df['PDB_ID'].map(lambda x : title(x))
# リンクの作成
df['Supplier'] = df['CID'].map(lambda x : 'https://pubchem.ncbi.nlm.nih.gov/compound/' + x + '#section=Chemical-Vendors&fullscreen=true')
df['RCSB_URL'] = df['PDB_ID'].map(lambda x : 'https://www.rcsb.org/structure/' + x)
df['Ligand_URL'] = df['Ligand_ID'].map(lambda x : 'https://www.rcsb.org/ligand/' + x)
# df['Title'] = json_data['data']['entry']['struct']['title']
# df結合
Sup_List = pd.concat([Sup_List, df])
else:
df = pd.DataFrame(columns=cols)
empty = []
empty.append(pdb_id)
df['PDB_ID'] = empty
df['RCSB_URL'].iloc[0] = 'https://www.rcsb.org/structure/' + pdb_id
Sup_List = pd.concat([Sup_List, df])
else:
df = pd.DataFrame(columns=cols)
empty = []
empty.append(pdb_id)
df['PDB_ID'] = empty
df['RCSB_URL'].iloc[0] = 'https://www.rcsb.org/structure/' + pdb_id
Sup_List = pd.concat([Sup_List, df])
Sup_List = Sup_List.reset_index(drop=True)
重複も除いてあげると親切ですね。# MW>150
Sup_List_150 = Sup_List[Sup_List['MW']>150].reset_index(drop=True)
#remove duplicates
Sup_List_Dup = Sup_List[Sup_List['MW']>150].drop_duplicates('Ligand_ID').reset_index(drop=True)
構築したconda環境のメモ
記事の概要
- 基本にしている仮想環境の構築例
- バージョンを指定した仮想環境の共有方法
- ついでにconda command
今までは教材に教わるがままにDockerで環境構築してたんですが、jupyterやOpennMMのようにコンテナ外でブラウザを使うようなパッケージの設定は私のような素人にはまだツラいです。。
というわけで、取り急ぎminicondaでサクッと構築環境を再現するためのメモです。
構築環境
- pymol
conda create -n pymol -c conda-forge python=3.10 pymol-open-source
- openmm (with PyMOL)
conda create -n openmm -c conda-forge python=3.10 openmm openmm-setup pymol-open-source
- openmm (with cudatoolkit, without PyMOL)
conda create -n openmm python=3.10 openmm openmm-setup cudatoolkit=11.2
- chemoinfomatics ref
conda create -n ci -c conda-forge -y python=3.10 numpy pandas scipy matplotlib seaborn scikit-learn boruta_py lightgbm xgboost deap rdkit jupyterlab spyder
- bioinfomatics
conda create -n bi -c conda-forge -y python=3.10 pandas pymol-open-source BioPython
パッケージを色々追加してグチャグチャになっても、とりあえず上記があれば必要最低限を再構築できてます。
pandas=1.4.2のようにバージョン名を指定して追加することもできます。
構築環境の共有方法
上記に加えて色々パッケージを追加した構築環境を他の人に共有したい場合は該当する、仮想環境内で以下のコマンドで仮想環境を出力して共有できます。
conda env export > env_name.yml
出力したymlファイルを元に環境を再構築する場合は下記コマンドです。
conda env create -n new_env -f env_name.yml
この方法を使うとPythonやパッケージのversionまで指定して再構築できるので便利です。
conda commandまとめ
conda環境の一覧
conda info -e
Channelの追加や削除
conda-forge channelの追加
conda config --append channels conda-forge
defaults channelを削除
conda config --remove channels defaults
パッケージ管理
- 現在のconda環境のパッケージ一覧表示
conda list
指定したconda環境のパッケージ一覧表示
conda list -n py38
パッケージの検索
conda search tensorflow
パッケージインストール
- インストール
conda install tensorflow
- バージョン指定してインストール
conda install tensorflow=1.15
- インストール
パッケージの更新
- パッケージを最新に
conda update tensorflow
- すべてのパッケージを最新に
conda update --all
- パッケージを最新に
パッケージの削除
conda remove tensorflow
パッケージリストを出力
conda list --export > package-list.txt
仮想環境作成
pythonバージョンを指定して作成
conda create --name py38 python=3.8
conda create -n py38 python=3.8仮想環境をコピーして作成
例：baseをコピーしてpy38を作る
conda create -n py38 --clone base
- 環境再現のためのリストから環境を再構築
conda create -n py38 --file package-list.txt
conda環境
有効化
conda activate py38
conda環境の出力
conda env export > env_name.yml
ymlファイルから環境を再構築
conda env create -n env_name2 -f env_name.yml
conda環境の終了
conda deactivate
仮想環境の削除
conda remove -n py38 --all
OpenMM用の仮想環境を構築する
OpenMM、というかGROMACSで遊ぼうと思って作った仮想環境の備忘録です。
サンプルコードはこんな感じ。。
conda create -n openmm -c conda-forge python=3.10 openmm openmm-setup pymol-open-source
ちなみに下記のようにcudatoolkit を入れようとすると-c conda-forgeが使えないので、pymol-open-sourceを別に入れた方がいいのかもしれない。
conda create -n openmm -c conda-forge python=3.10 openmm openmm-setup pymol-open-source cudatoolkit=11.2
PyMOLは-c conda-forgeでopen source版をインストールできるので覚えておいて損はないと思います。
膨大なデータから金の粒を見つける
こちらの記事でChEMBLから引っ張ってきたデータを可視化しました。
全体を俯瞰することも大切ですが、キーになりそうな化合物を見つけることも大切です。
というのが今日のお話です。
データのどこに目を向けるのか
論文や特許を見る場合、データのどんなところを気にするでしょうか。
- 活性値
- 選択性
- ADME
- PK
- 薬効
- 開発状況
- 臨床データ
化合物のステージによって記載されている情報が当然違うと思いますが、それぞれのステージのチャンピオン、例えば活性が最も強い化合物や選択性が最も高い化合物が先に進むわけではないので注意が必要です。
今回は全体を俯瞰しているので、化合物あたりのデータ数が多い化合物に注目すればキーになりそうな化合物が見えてくると考えます。
実習
データはこちらでクレンジングしたものを再利用します。
化合物あたりの統計量を調べたいので、Molecule ChEMBL IDでgroupbyします。
文献あたりの記載回数が多ければ、色々データを取るようなキー化合物だろうと想像するのでDocument ChEMBL IDに絞り込みます。
dfi = df.groupby(['Molecule ChEMBL ID'])['Document ChEMBL ID'].describe() dfi[dfi['unique']>2]
- count : データ記載回数
- unique : データが記載された論文数
- top : データ記載頻度の最も高い論文ID
- freq : topの中でデータが記載された回数
という感じだと思います。
1度しか論文に載っていない化合物が多すぎるので、filterで除いてあります。
CHEMBL221959という化合物は25の論文に計30回記載されていますが、論文あたり最も多く記載された回数は2回でした。今回の標的はTYK2というkinaseなので、上市品など有名な化合物で活性と選択性のポテンシャル確認が多かったのかもなと予想します。
どんな化合物でしょうか。
from rdkit.Chem import Draw dft = df[df['Molecule ChEMBL ID']=='CHEMBL221959'] Draw.MolToImage(dft['Molecule_Image'].iloc[0], legend=dft['Molecule Name'].iloc[0])
Tofacitinibという同じく一番最初に上市されたJAKiでした。 ヤッパリネ！
CHEMBL4435170は記載論文数が4報にも関わらず、最多データ記載回数が5回なので色々とデータが報告された化合物かもしれません。
c_list = [
'Molecule Name',
'Document ChEMBL ID',
'Assay Description',
'Standard Type',
'Standard Value',
'Standard Units',
'Murcko_generic_SMILES'
]
dft = df[df['Molecule ChEMBL ID']=='CHEMBL4435170']
dft.loc[:, c_list].sort_values('Document ChEMBL ID')
完全にたまたまなんですが、これまで紹介してきたDeucravacitinibがヒットしました(笑)。
細胞活性やアロステリック活性の報告だったみたいです。
Document IDから文献を調べる方法は次の宿題ということで。。(やり方忘れた。。
もう一度見渡す
Deucravacitinibはやっぱり注目すべき化合物だったんだろうなということで、化合物の起源をもう少し調べたいです。
Murcko_generic_SMILESを使って遡ってみます。
df_deu = df[df['Murcko_generic_SMILES']=='CC(CC1CCCC(CC2CCCC(C3CCCC3)C2)C1)C1CC1'] df_deu.loc[:, ['Molecule Name', 'Compound Key', 'Document ChEMBL ID', 'Document Year']]
あれ？案外報告されてなかったケモタイプなんですね...
じゃあ同じ文献内で分子量の小さいケモタイプを遡って...
みたいにすると化合物起源にたどり着けると思います。
順番が逆になってしまいましたが、Deucravacitinibのストーリーは非常にざっくりとこちらで紹介しています。
ChEMBLのデータは特許が入ってませんし、商用データベースと比べると網羅性は高くありません。
ただ、データとしては重要なものを綺麗にまとめてくれてるので非常に有用だなと思いました。
おしまい。
pandasのカラム番号を検索する小ネタ
こちらです
def c_index(c_name='', data_frame=df): c_zip = zip(list(range(len(data_frame.columns))),data_frame.columns) # c_name = input('column name?') return {key:value for key,value in c_zip if c_name in value}
defaultではdfをDataFrameとして全てのカラムの番号を表示します。
第一引数に指定した文字列を含むカラムとそのカラム番号を表示できます。
第二引数に指定したDataFrameのカラムを取り扱います。
カラムが多すぎるDataFrameで困った
前回までの記事でChEMBLからデータを引っ張ってきました。
DataFrameのカラムが多すぎてチンプンカンプンになるので最後はSummary Tableを作りましたが、実際はカラムの数が多いので把握するのも大変です。
df.columns # Index(['Molecule ChEMBL ID', 'Molecule Name', 'Molecule Max Phase', # 'Molecular Weight', '#RO5 Violations', 'AlogP', 'Compound Key', # 'Smiles', 'Standard Type', 'Standard Relation', 'Standard Value', # 'Standard Units', 'pChEMBL Value', 'Data Validity Comment', 'Comment', # 'Uo Units', 'Ligand Efficiency BEI', 'Ligand Efficiency LE', # 'Ligand Efficiency LLE', 'Ligand Efficiency SEI', 'Potential Duplicate', # 'Assay ChEMBL ID', 'Assay Description', 'Assay Type', 'BAO Format ID', # 'BAO Label', 'Assay Organism', 'Assay Tissue ChEMBL ID', # 'Assay Tissue Name', 'Assay Cell Type', 'Assay Subcellular Fraction', # 'Assay Parameters', 'Assay Variant Accession', 'Assay Variant Mutation', # 'Target ChEMBL ID', 'Target Name', 'Target Organism', 'Target Type', # 'Document ChEMBL ID', 'Source ID', 'Source Description', # 'Document Journal', 'Document Year', 'Cell ChEMBL ID', 'Properties', # 'MW', 'Ligand Efficiendy SEI', 'Ligand Efficiendy BEI', # 'Ligand Efficiendy LE', 'Molecule_Image', 'Murcko_SMILES', # 'Murcko_generic_SMILES', 'Murcko_Mol', 'Murcko_generic_Mol', 'TPSA', # 'Ct HA', 'Ct F', 'HBD', 'HBD_Ro5', 'HBA', 'HBA_Ro5', 'Ct RB', 'Ct Ar', # 'Fsp3', 'QED'], # dtype='object')
カラムの数だけでも65もありました...
大文字と小文字が混ざってたり_が入ってたり#が入ってたりでタイポも増えるしと地味に悩んでいました。
DataFrameからdataを取り出す方法
pandasでは特定のカラムを抽出する方法がいくつか用意されています。
代表的なものは
df['Molecule ChEMBL ID'] #0 CHEMBL2181327 #1 CHEMBL2181315 #2 CHEMBL21156 #3 CHEMBL535 #4 CHEMBL475251 # ... #1806 CHEMBL3968872 #1807 CHEMBL3983594 #1808 CHEMBL3957838 #1809 CHEMBL4563864 #1810 CHEMBL4873733
みたいな方法です。
例えば複数のカラムを抽出したいときは[]の中にさらにカラム名のリストを入れます。
df[['Molecule ChEMBLE ID', 'Molecule Name', 'Smiles']] # Molecule ChEMBL ID Molecule Name Smiles #0 CHEMBL2181327 NaN CN(C)CCN(C)c1ccc2cc1COCC=CCOCc1cccc(c1)-c1ccnc... #1 CHEMBL2181315 NaN C1=CCCOc2cccc(c2)-c2ccnc(n2)Nc2cccc(c2)OCC1 #2 CHEMBL21156 NaN CC(C)(C)c1nc2c3ccc(F)cc3c3c(=O)[nH]ccc3c2[nH]1 #3 CHEMBL535 SUNITINIB CCN(CC)CCNC(=O)c1c(C)[nH]c(/C=C2\C(=O)Nc3ccc(F... #4 CHEMBL475251 R-406 COc1cc(Nc2ncc(F)c(Nc3ccc4c(n3)NC(=O)C(C)(C)O4)... #... ... ... ... #1806 CHEMBL3968872 NaN CN(CC(F)(F)F)C[C@H]1CC[C@H](N2CN(CC#N)C(=O)c3c... #1807 CHEMBL3983594 NaN N#CCN1CN([C@H]2CC[C@H](CNC3(C(F)(F)F)CC3)CC2)c... #1808 CHEMBL3957838 NaN O=C1c2cnc3[nH]ccc3c2N([C@H]2CC[C@H](CN3CCCS3(=... #1809 CHEMBL4563864 NaN COCCn1cc(Nc2cc(NCc3cccc(NC(=O)CC#N)c3)ncn2)cn1 #1810 CHEMBL4873733 NaN Cc1cnc(Nc2ccc3c(c2)CCN3C(=O)[C@H]2CCCN2C(=O)C2...
わたしはいつもカラム名をタイポするので地味にストレスです。
.locメソッドと.ilocメソッド
.loc[]や.iloc[]は(多分)locationとindex locationに由来するメソッドです。
引数はそれぞれ
df.loc[row index, column name] df.iloc[row index number, column index number]
となってます。
DataFrameのindexがrow numberの場合は、どちらも数字で同じ結果が返ってきます。
indexたdatatimeのようにrow number以外を採用しているときはやっかいなので、.iloc[]を使うことが多いのかなと思います。
例えば
df.loc[0, 'Molecule Name'] # nan
のように値を取り出したり代入するときに使います。
第二引数を省略することで行全体も抽出できます。
df.loc[0] # Molecule ChEMBL ID CHEMBL2181327 # Molecule Name NaN # Molecule Max Phase 0 # Molecular Weight 459.59 # #RO5 Violations 0 # ... # HBA_Ro5 7 # Ct RB 4 # Ct Ar 3 # Fsp3 0.333333 # QED 0.577017 # Name: 0, Length: 65, dtype: object
同じ要領で.iloc[]メソッドだと行と列のindex numberを入力するだけで同じことができます。
非常に簡単な反面、正しくindex numberを把握しておく必要があります。
df.iloc[0,1] # nan
df.iloc[0] # Molecule ChEMBL ID CHEMBL2181327 # Molecule Name NaN # Molecule Max Phase 0 # Molecular Weight 459.59 # #RO5 Violations 0 # ... # HBA_Ro5 7 # Ct RB 4 # Ct Ar 3 # Fsp3 0.333333 # QED 0.577017 # Name: 0, Length: 65, dtype: object # このDataFrameはrow index numberを採用しているので.locと同じ結果が返ってくるはず
じゃあカラムのindex numberを簡単に把握できれば少し楽になるかも？？
Column index identifier
というわけでメソッドを作ってみました。
def c_index(c_name='', data_frame=df): c_zip = zip(list(range(len(data_frame.columns))),data_frame.columns) # c_name = input('column name?') return {key:value for key,value in c_zip if c_name in value}
引数なしで使うとカラムにindex numberを割り当てた結果を全て出力します。
c_index() # Output exceeds the size limit. Open the full output data in a text editor # {0: 'Molecule ChEMBL ID', # 1: 'Molecule Name', # 2: 'Molecule Max Phase', # 3: 'Molecular Weight', #... # 63: 'Fsp3', # 64: 'QED'}
第一引数に検索したい文字列を入力すると、文字列を含むカラム名とindex numberの割り当てを出力します。
c_index('As') # {21: 'Assay ChEMBL ID', # 22: 'Assay Description', # 23: 'Assay Type', # 26: 'Assay Organism', # 27: 'Assay Tissue ChEMBL ID', # 28: 'Assay Tissue Name', # 29: 'Assay Cell Type', # 30: 'Assay Subcellular Fraction', # 31: 'Assay Parameters', # 32: 'Assay Variant Accession', # 33: 'Assay Variant Mutation'}
第二引数にDataFrameを入れることも出来ます。面倒なのでデフォルトはdfにしてあります。
これを使うと特定の情報をシンプルに抽出できます。
df.iloc[[0,3,50, 23,11], [0,1,22]] # Molecule ChEMBL ID Molecule Name Assay Description # 0 CHEMBL2181327 NaN Inhibition of recombinant TYK2 # 3 CHEMBL535 SUNITINIB Binding constant for TYK2(JH2domain-pseudokina... # 50 CHEMBL3949572 NaN Enzyme Assay: The inhibitory activity of the c... # 23 CHEMBL1088856 NaN Inhibition of recombinant GST tagged TYK2 prot... # 11 CHEMBL3655081 ABROCITINIB Inhibition of N-terminal GST-tagged human TYK2...
コードの可読性は下がるので必ずしもオススメはしませんが、どんなカラムがあったかも検索できるので便利かもです。
Deucravacitinibの二つのメチル基を観察してみる
最近承認されたTYK2阻害剤のユニークなキャラクターを前々回のこちらの記事で紹介しました。 keetaneblog.hatenablog.com
いつもはPyMOLで遊んだスクショを載せたりしてましたが、最近遊んだpy3Dmolも前回テストしたので、今回はこれを使って構造を紹介します(使ってみたかっただけ)。
まずDeucravacitinibがどんな構造かというと、こんな感じです。
メチル基自体は
- メチルアミド
- メトキシ
- メチルトリアゾール
の３つがあります。
このうち、重要なのは1と2、特に1のメチル基はJAK1~3、TYK2のJH1結合活性を全く示さないことから、血球系の副作用を示さないことがDeucravacitinibと周辺ケモタイプの特徴として報告されています。
これらの特徴がなぜ発現するかを構造生物学的に説明することができます。
残念ながら2023.2現在、DeucravacitinibのPDBは登録されていませんが、類縁体がいくつか報告されています。
こちらはDuecravacitinib類縁体の一つとTYK2 JH2の共結晶構造です。(PDB ID : 6NZH)
左クリックのドラッグで回転、スクロールでズーム、controlキーを押しながらドラッグすると並行移動できます。
ポケットに入り込んだアミドのメチル基周辺を観察してみてください。主鎖と結構近いですよね？ね？
こればっかりは側鎖を表示して重ね合わせないと比べられないですが、JH1ドメインのポケットはどれも側鎖で埋まっているらしいですが、JH2はAlanine pocketと呼ばれる小さな間隙が存在し、そこを埋めることが選択性を発現していると考察しています(と記憶しています。間違ってたらごめんなさい。)
他のJH2ドメインに対する選択性はどうかというと、完全に切ることは難しかったようです。
こちらの文献に類縁体の報告があります。
https://pubs.acs.org/doi/full/10.1021/acs.jmedchem.0c01698
ただし、JH2ドメインにはそもそも機能がないので他のファミリーでは細胞活性が見られず、Deucravacitinibや類縁体はTYK2だけの活性をJH2側からmodulateしているとの報告でした。
selendipityから見出された選択性ですが薬効面と副作用回避の面から非常に魅力的で、大きな市場になることが期待されているようです。
論文だけだと非常にコンパクトですが、これらの仮説を検証していく過程は彼らの非常に高いサイエンスを創造させてくれました。
最後に今回使用したhtmlのスクリプトです。
はてなブログではmarkdown modeに貼り付けるだけで使うことができます。
<script src="https://3Dmol.org/build/3Dmol-min.js"></script> <script src="https://3Dmol.org/build/3Dmol.ui-min.js"></script> <div style="height: 690px; width: 690px; position: relative;" class='viewer_3Dmoljs' data-cid='134821691' data-style='stick'></div>
<script src="https://3Dmol.org/build/3Dmol-min.js"></script> <script src="https://3Dmol.org/build/3Dmol.ui-min.js"></script> <div style="height: 400px; width: 400px; position: relative;" class='viewer_3Dmoljs' data-pdb="6nzh" data-select="chain:A" data-style="cartoon:color=spectrum" data-select1="hetflag:true" data-style1="stick" data-select2="chain:B" data-style2 data-zoomto="chain:A" > </div>
次回(こそ)はデータ解析後編につづく...はず。。。
