見出し画像
Photo by ia19200102

機械学習で創薬探索のプロセスをかじってみました（パート１：データの取得と前処理）

leap

バイオインフォマティクスに疎い私ですが、創薬探索のプロセスを経験すべく、機械学習で化合物のデータを用いた分析をしてみました！

参考にしたのは、Data Professor という名前でユーチューバーをされているタイ人でバイオインフォマティクスの元准教授の Chanin Nantasenamat さんのビデオ Bioinformatics Project from Scratch - Drug Discovery Part 1 (Data Collection and Pre-Processing) になります。説明欄にビデオで使用されているコードのリンクも掲載されています。なお、ここで記述したコードは、自らがやりやすいように一部変更を加えていますのでご了承下さい。

開発環境

MacBook Air
Google Colaboratory
Python 3.6.9

目的

薬剤候補化合物による、エストロゲン受容体 α1(estrogen receptor alpha 1)に対する阻害活性を推測するための機械学習モデルを構築する。

概要

  • ChEMBL データベースからエストロゲン受容体 α1(estrogen receptor alpha 1) に対する阻害活性が調べられている化合物のデータを取得する。

  • 取得したデータを前処理し、阻害活性のレベルに従って分類する。

  • データを Google drive に保存する。

データの取得

ChEMBLからデータを取得するためのパッケージのインストール

ChEMBL とは200万以上の化合物の化学構造や阻害活性などのデータが保存されているデータベースです。今回はここから化合物のデータを取得します。

! pip install chembl_webresource_client


必要なライブラリーのインポート

import pandas as pd
from chembl_webresource_client.new_client import new_client



ターゲット分子の探索

今回は estrogen receptor alpha 1 をターゲット分子にしてみました。ビデオでは coronavirus 3C-like proteinase (CHEMBL3927) をターゲット分子にしています。

ターゲット分子 Estrogen receptor alpha 1 の情報を取得

target = new_client.target
target_query = target.search('Estrogen receptor alpha')
targets = pd.DataFrame.from_dict(target_query)
targets
3048 rows × 9 columns

このようにChEMBLのサイトで検索した場合と同様の検索結果が出てきます。

ターゲット分子の選択
次に、この結果から実際にターゲットとして使用したい分子のインデックス番号（ここでは「２」で　Estrogen receptor alpha1、CHEMBL206）を以下のコードに入れると、target_chembl_id カラムに入っている ID番号が出力されます。

selected_target = targets.target_chembl_id[2] 
selected_target

’CHEMBL206'

ターゲット分子に対する阻害活性のデータが存在する化合物の抽出

この Estrogen receptor alpha1（CHEMBL206）の活性への影響が過去に調べられており、かつIC50 値（nM）が阻害活性の単位として使用されている化合物のみを抽出します。

　なお、IC50（half maximal inhibitory concentration）とは、ターゲット分子の50%を阻害するのに必要な濃度のことを言います。IC50 値が低いほど、つまり濃度が低いほど阻害活性が高いということになります。ターゲット分子の半数を邪魔するのに少量で済んでしまう、ということです。

activity = new_client.activity
res =activity.filter(target_chembl_id=selected_target).filter(standard_type="IC50")
df = pd.DataFrame.from_dict(res)
df
4567 rows × 45 columns

このように、計4567の化合物の化学構造、阻害活性等の沢山の情報が出力されました。このデータを csv 形式で保存します。

df.to_csv('estrogenRa_bioactivity_data.csv', index=False)


データの前処理

データが欠損している行の削除

分析に必要なstandard_value 、canonical_smilesカラムのいずれかの値が欠損している行を削除します。この処理により、”4567 rows × 45 columns”から”3543 rows × 45 columns”と、約1000の化合物のデータが削除されました。

# standard_value カラムの値が欠損している行の削除　行数が3551に
df2 = df[df.standard_value.notna()] 
# canonical_smiles カラムの値が欠損している行の削除　行数が3543に
df2 = df2[df2.canonical_smiles.notna()] 
df2

　なお、canonical smiles とはSMILES 記法といって化合物の分子構造を原子記号などを使って一次元配列にしたもので「O[C@]1(C(F)(F)F)CCCC[C@H]1Nc1ccc(F)cc1」などのように表されています。つまり、各化合物の立体構造とその化学的な性質についての情報であるとご理解いただければいいかと思います。これらの情報はターゲット分子を阻害しうるか否かを予測する上で非常に重要です。
SMILES 記法の詳細は SMILES記法 その1 をご覧ください。


阻害活性データの前処理

IC50の値によって阻害活性のレベルを分類します。

  • ≦ 1000 nM: 活性有り（Active）

  • ≧ 10,000 nM: 不活性（Inactive）

  • >1000 nM、<10,000 nM: Active と Inactive の間（Intermediate）

低濃度であるほど、ターゲット分子に対する阻害活性が高いとみます。


bioactivity_class = []
for i in df2.standard_value:
  if float(i) >= 10000:
    bioactivity_class.append("inactive")
  elif float(i) <= 1000:
    bioactivity_class.append("active")
  else:
    bioactivity_class.append("intermediate")


必要な情報のみを抽出したデータフレームを作成

'molecule_chembl_id', 'canonical_smiles', 'standard_value'の３つのカラムのみを取り出し、

selection = ['molecule_chembl_id', 'canonical_smiles', 'standard_value']
df3 = df2[selection]
df3

上述した阻害活性のラベルのデータと結合します。

# bioactivity_class データを上の３つのカラムでできたデータフレームと結合
pd.concat([df3, pd.Series(bioactivity_class, name = 'bioactivity_class')], axis=1)

が、このコードがどうやってもうまく動きませんでした😢（データの中の文字の一部に問題があったのか？）。よって、以下のより手間がかかる方法を用いて結合しました。

# molecule_chembl_id データのリストの作成
mol_cid = []
for i in df2.molecule_chembl_id:
  mol_cid.append(i)

# canonical_smiles データのリストの作成
canonical_smiles = []
for i in df2.canonical_smiles:
  canonical_smiles.append(i)

# standard_value データのリストの作成
standard_value = []
for i in df2.standard_value:
  standard_value.append(i)

# これらのリストと'bioactivity_class'のリストを結合
data_tuples = list(zip(mol_cid, canonical_smiles, bioactivity_class, standard_value))
df3 = pd.DataFrame( data_tuples,  columns=['molecule_chembl_id', 'canonical_smiles', 'bioactivity_class', 'standard_value'])
df3


このデータを csv 形式で保存します。

df3.to_csv('estrogenRa_bioactivity_preprocessed_data.csv', index=False)


Google Drive へのファイルの保存


データファイルを Google Drive にも保存できます！まず、Google colab からGoogle Drive にアクセスできるように以下のコードを実行します。コードの実行後、ブラウザー上で Google account のパスワードを入力するなどの指示がありましたが、その通りに進めるとうまくいきました。

from google.colab import drive
drive.mount('/content/gdrive/', force_remount=True)

次にファイルを保存するための新しいフォルダーを作成します。ここではbiodataというフォルダ名にします。

! mkdir "/content/gdrive/My Drive/Colab Notebooks/biodata"
! cp bioactivity_data.csv "/content/gdrive/My Drive/Colab Notebooks/biodata"

作成したフォルダにファイル名を指定して保存します。

! cp estrogenRa_bioactivity_preprocessed_data.csv "/content/gdrive/My Drive/Colab Notebooks/biodata"

ls コマンドで保存されたか確認します。保存されていたらファイル名が出力されます。

! ls "/content/gdrive/My Drive/Colab Notebooks/biodata"
# ls -l で保存時間などのファイルのより詳細な情報がみれます。

estrogenRa_bioactivity_preprocessed_data.csv


ファイルの最初の数行を確認できます。

! head estrogenRa_bioactivity_preprocessed_data.csv


次回


次回のパート２では、探索的データ分析を行います。化合物の化学的な特徴を調べて、経口薬の候補になり得るかを検討します。また、阻害活性のレベルによる結果の違いが、統計的に有意であるかどうかも検討します。

この記事が気に入ったらサポートをしてみませんか？