技術資料
研究背景
- 有機合成では効率性や環境面から化学反応の設計に生体触媒(酵素)が利用される機会が増加
- 酵素は化学触媒よりも環境に優しく，特定の反応物に強く作用して化学反応を早く進める
- 酵素は4桁のEC番号で分類されている
- EC番号を機械学習で予測することで，酵素分析(以下)を行う実験コストや時間の削減になる
- 代謝経路の解析で得られた未知の酵素の性質の特定
- 特定の化学反応を触媒する酵素の探索
目次
- 1節：必要な事前知識(特に4節)についての説明
- 2節：環境構築
- 3節：データ加工
- 4節：数値実験の再現
- データセットのEC番号分類を実行
- 必要に応じて，1節を参考
1. 事前知識
この研究で用いる機械学習・特徴量エンジニアリング手法についての説明
多クラス分類について
- 機械学習の分類タスクに該当
- 正解クラス(ラベル)が与えられた学習データを正しいクラスに分類する問題を分類器に学習させ，分類モデルを作成
- テストデータを分類モデルに入力し，正しいクラスに分類できるかを評価する
- クラスが3種以上ある場合が多クラス分類(2種は2クラス分類)
評価指標
- 混同行列などを用い，Precision・Recall・F1-Scoreで評価するケースがよく見られる
- Accuracyでも評価することはあるが，クラスのデータ数に偏りがあると，データの多いクラスの正解率に左右されてしまい，データの少ないクラスの分類精度を適切に評価できない場合がある．
- 特に酵素反応のEC番号データは偏りが非常に大きいため，F1-Scoreが重要となる
Random Forests (RF)について
以下のページ参照
機械学習と特徴選択
交差検証について
RFを選んだ理由
説明可能な機械学習モデルを用いたため
- 有機合成の研究者に将来的に利用してもらうため，「なぜそのような予測をしたのか」の説明ができないとシステムの信頼性に欠ける
- 深層学習などの手法では度々なぜそのような予測をしたのかの説明が難しい(ブラックボックス)場合がある
- RFはいわゆる「人の手で出来ないレベルのIF分岐」を行うモデルで，「決定木のそれぞれのノードで特徴量に応じてデータが分類されている」という説明ができる
- これはあくまで卒論向けの主張であり，現在はRFよりも良い手法があると考えている
- 論文参考文献[5]の手法では，深層学習系の多層パーセプトロンで高精度予測しているだけでなく，Shaplay値系の手法を組み込むことで，化合物のどの部分構造がEC番号分類に大きく寄与しているかを説明(可視化)できるようにしている．
2.1 環境構築
本プログラムはAnaconda上で実装する
Anacondaインストール
以下のページからインストール
https://www.anaconda.com/download#downloads
rdkitチャネル(仮想環境)の作成
仮想環境の補足
- Anacondaは通常"base"と呼ばれるチャネル(仮想環境)を用いるが，新しいチャネル"rdkit"を作る
- バージョンが1.0.2(2021年版)のscikit-learnを用いるため，メイン環境として"base"，本技術資料のの再現環境として"rdkit"を用いることを推奨
手順
WindowsボタンからAnaconda Promptを開く
バージョン取得
python -V
rdkitチャネル作成
# 研究当時のバージョン
conda create -n rdkit python=3.8.8 -y
※baseチャネルからrdkitチャネルに変更(Anaconda Prompt)
conda activate rdkit
rdkitライブラリの導入(pipでは不可)
conda install rdkit -c conda-forge
jupyterカーネル導入
pip install jupyter environment_kernels
jupyter起動
jupyter notebook
Jupyterのプログラムを起動するとき
anaconda promptを起動し，(base)の部分で以下を入力
(base) C:\Users\～～～>activate rdkit
base→rdkitに切り替わったら，Jupyter起動
(rdkit) C:\Users\～～～>jupyter notebook
- インストールしたモジュールはbaseではなく，rdkitチャネル内のディレクトリに保存される
ディレクトリを辿って「master_ECprediction」まで行き，「～.ipynb」を起動する
2.2 プログラム実行時の補足
実行するプログラムファイル
以下2ファイルで構成される
- 機械学習に用いるデータセットの作成を行う「ECprediction_dataset.ipynb」
- データセットから多クラス分類モデルを作成し，EC番号予測を行う「ECprediction_model.ipynb」
プログラムの取り扱い説明
- プログラムはそれぞれセルで分割されており，1セルずつ順次実行していく
- 計算に時間がかかる箇所の状態確認が重要のため，通常のPythonファイルとは異なり，あえてプログラムを分割している
- セルをクリックし，shift + enterでセル内のプログラムを実行できる
- 変数名，もしくはprint(変数名)をセルの下部に記載することで変数の中身を確認可能
- 初回時は最上セルのpipインストールでモジュールを導入
リスタートポイントについて
途中からプログラムを実行できるようにあらかじめ保存したファイルを読み込み再開できる仕組み
- 上部セルのimportを実行したあと，任意のリスタートポイントから実行可能('''のコメントアウトを外す)
- プログラムを完全再現するのに8日程度かかるため，実行時間が長いセルは飛ばし，最寄りのリスタートポイントから次のセルに進むことを推奨する
- 実行時間(Intel Core i7-9700K 3.6GHz の場合)を参考にして次に進む
- 特に各反応式の特性値の計算，記述子選択，パラメータ調整に多大な時間がかかる
- リスタートポイントから読み込まれるファイル名の語尾には「0」を記載し，通常実行で保存されるファイルと区別している
3. モデル構築用データセット作成(ECprediction_dataset.ipynb)
- 各セルでどのような処理を行っているかの大まかな説明はコメントアウトに記載
- 以下では，セル実行時の補足を記載
3.1 プログラム順次実行
モジュールインストール
pip install rdkit
元データセットを1つの化合物SMILESに分解
- 元データセット*1は，左辺と右辺のSMILESを連結している
- 左辺と右辺を隔てる'|X.X.X.X>>'と，第1項～第N項を隔てる'.'に着目し，左辺のSMILESリスト，右辺のSMILESリストに分解
- MetaNetXにおけるEC 2.A.3.1などのデータは削除している
- 2次元配列のleft_sms，right_smsにそれぞれ格納
特性値計算不可の化合物SMILES削除
SMILES形式の化合物から210種の特性値を計算する際，発散値を持つ化合物が存在するため，そのような化合物を除去
リスタートポイント1から5の前までのセル
- 右辺左辺から全ての化合物(51478個)を抽出し，RDKit構造オブジェクトに変換(～初めの3セル)．
- 1化合物に208種類の特性値を計算し，発散値とnan値を持つ削除化合物リストの作成(～リスタートポイント4)
- バージョン依存の関係から2種の記述子を除外して実行
- SMILES反応式の左辺と右辺を走査し，削除化合物リストを含む反応式を除去(～リスタートポイント5手前)
整理された左辺と右辺のSMILES反応式データフレームが出力される(left_smsDF5, right_smsDF5)
210次元特徴ベクトルの作成
リスタートポイント5から最後までのセル
- 左辺と右辺のSMILES反応式をRDKit構造オブジェクトに変換(～初めの3セル)
- 構造オブジェクトでは計算不可の特性値を持つ化合物を検索できないため，検索後もう一度変換
- 左辺，右辺それぞれの化合物に対して，210種類の特性値を計算(～リスタートポイント6)
- 重い処理のため，Jupyter以外のブラウザやプログラムを閉じることを推奨
- 左辺の各特性値の和と右辺の各特性値の和の差を取り，210種の特性値変化量を要素に持つ特徴ベクトルを作成(～最後のセルまで)
4. EC番号分類モデルの作成，テストデータの分類 (ECprediction_model.ipynb)
- 各セルでどのような処理を行っているかの大まかな説明はコメントアウトに記載
- 以下では，セル実行時の補足を記載
4.1. モジュールインストール(多クラス分類，SMOTE用)
pip install scikit-learn==1.0.2
pip install imblearn
4.2. データ削減1
「データ削減」セル1,2
SPS記述子の除去
- 計算不可の構造オブジェクトが多いためdrop()で除外
相関係数1の記述子・データの片方削除(下図参照)
- 全て同じ値になる記述子は過学習の原因のため片方のみ残す
- 多重共線性などの問題で相関係数が高い特徴量を除外するケースは結構存在する
- データも同様に片方のみ残す
- クラスのデータ数を減らさないため，あえて違うクラスとして残して予測している研究もあるが，異なるクラスに誤分類される可能性も高まることから，本研究では対象外
np.float32の最小・最大値を超える要素を持つ特徴ベクトルの削除
全ての特性値で0のデータを削除
4.3. データ削減2
「データ削減」セル最後
データ数が5個以下のものを削除する
- 層化5分割交差検証を実行するためにはデータ6個以上必要
- 6個では学習:テスト = 4:1分割で，学習5個/テスト1個となり，学習データのみに交差検証を行うため
- 10分割などの場合も，どのクラスもテストデータに最低1個は含まれる必要がある
- 【未実施】確率で同じデータを抽出する層化シャッフル抽出の場合は学習データが3個でも可能？
EC番号が3桁目まで存在するクラスのみ抽出
- 現状のデータはEC T.X.X.X (T=1,2, ... ,7)でラベル付けされている．
- 本研究では3桁目の予測のため，3桁まで分かっているデータのみ用いる
- 遺伝子配列の手法では4桁まで予測されており，いずれはフィンガープリント・物理化学特性値でも開発されると予想(or 開発する必要がある)
- 正規表現のOR演算記法(|)を用い，EC X.X.X.X，EC X.X.X.-，EC X.X.X.S(Sは MetaNetX 固有の文字区別)のデータを抽出する
削除するデータ
- EC7のデータ(データ数が少ないため)
- EC X.-.-.-やEC X.X.-.-などのデータ(1桁分類,1～2桁分類になら使用可能)
【懸念事項】ラベル変更によるクラス分布の変動(交差検証)
4桁までのクラスで層化抽出したデータを1桁or1～3桁に変更した場合，各クラスの学習データとテストデータの割合は一定かどうか(層化抽出が崩れていないか)
- 割合を検証したところ，どのクラスも学習75～85%，テスト15～25%で分布していたため，train_test_splitの標準通りと考えられる．
- 少なくとも4桁の時点で最低5個あるので，3桁クラスでは必ず5個以上ある(5分割可能)
4.4. 予備実験2
- 本実験に用いる記述子組合せの決定と各分類モデルの精度検証
- 1桁目(EC1～EC6)，2,3桁目(EC T.X.X (T=1,2, ... ,6))のクラス分類それぞれで，精度を向上させる記述子組合せを選択し，(本実験用に)マージ・重複削除する．
学習データ(Xtrain)・テストデータ(Xtest)の分割
この研究では，データのラベル付与が重要となる
- 初めにEC番号が3桁または4桁まで分かっているデータを学習:テスト=4:1で分割する
- 1桁目分類の際は，学習・テストデータそれぞれでEC番号3,4桁→1桁にラベルを貼り替え
- 「EC X (1桁目)クラス分類」セクションのセル2を参照
- 2,3桁目分類の際は，1桁目がTの学習(テスト)データのみを抽出し，3,4桁→3桁にラベルを貼り替え
- 「EC 1.X.X (EC 1の2,3桁目)クラス分類」セクションのセル2参照
記述子選択を分けて行う理由
本来ならラベルの貼り替えなしで，3,4桁のラベルが付いた状態で，1～3桁目 or 1～4桁目のクラス分類(記述子選択)を行う方が効率的
- 【時間的な制約】Xtrainに対し交差検証で1～3桁目分類する場合，実行時間が75時間以上(その間PCオーバフローの懸念)推定されていた．(37672データ，148クラス)
- 4桁目までの分類となるとクラス数が増加し，さらに時間を要する
- EC X分類(37672データ，6クラス)，EC 1.X.X(6380データ，64クラス)と分けて分割することで，実行負担を減らせる
- 【先行研究と同様】タンパク質配列を用いた手法*2では，1桁目を分類した後，2～4桁目を分類するモデルが作成されており，その方法に倣った．
- 実際，テストデータを1桁目分類モデルに入力し，正解したデータのみ2,3桁目分類モデルに入力して，正解したデータ数を調べることも可能だが，この研究では行っていない
EC X (1桁目)クラス分類
SMOTEなしで記述子選択を行う
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2値分類TP, FPなどの説明
