この記事は異業種データサイエンス研究会®で行った勉強会の内容＋αです。
proteinBERT
proteinBERT 1は、
- 自己教師あり深層言語モデルの手法をタンパク質に応用した。
- 最先端の先行研究（2022年当時）と比較してモデルが小さく、ほぼ同等の精度。
- 汎用性が高い。
- UniProt データベースの配列と Gene Ontology Database (Go) のデータを用いて事前学習を行っている。
という機械学習モデルです。
論文では、
- 二次構造
- 遠い相同性（ファミリー分類）
- シグナルペプチドの有無
- ニューロペプチド2
- 蛍光
- 安定性
- リン酸化修飾の有無とその位置
などをベンチマークとして性能の評価を行っていました。
ベンチマークの種類によって予測精度にはばらつきがあり、シグナルペプチドや二次構造といった他のプログラムでも高い精度で予測できるものは成績が良く、遠い相同性など先行研究でも難しいとされているタスクに対してはあまり成績が良くないようでした。
Google colab で動くようにしてもらった
proteinBERT のコードは Github で公開されています。
勉強会で Github にある proteinBERTdemo.ipynb のコードを Google Colab に移植して頂きました！
なので、proteinBERT の環境構築ができているところから出発し、私の研究対象である酵素のデータを使って活性（EC番号）を予測できるかやってみたいと思います。
Github には事前学習用のデータセットも用意されていますが、今回は事前学習済みのモデルをダウンロードして使用しました。
古いバージョンでは、'Adam' object has no attribute 'get_weights' というエラーが出ますが、 model_generation.py を修正版に置き換えることで解決しました。
データセットの準備
詳しい ProteinBERT の説明は省略します（というかできません；）が、入力データはこのようになっています。
入力データのアミノ酸配列は、以下のように行列として扱われます。
（徳山さんのスライド）
proteinBERT で実際にベンチマークとして用意されていた入力データは以下のような感じでした。
remote homology（カテゴリー分類）
signalpeptide（二値分類）
fluorescence（連続値）
自前のデータをこの形式にすれば使えそうだったので、私が大学で研究している「糖質加水分解酵素」の配列から活性（EC 番号）を予測できるか試してみることにしました。
糖質関連酵素界隈では、酵素の分類は CAZy database 最もよく参照されています。しかし、CAZy database からアノテーションを含めて配列情報を取得するのがとても大変なので、ここではCUPP のデータベースを利用しました。
探索的データ分析(EDA)
ダウンロードしたテーブル（2023/04/26 取得）は、23816 行でした。proteinBERT の論文では、学習データに、UniRef90 に登録されている、90% 以上の相同性を持たないタンパク質配列を用いていたので、これに倣って cd-HIT を用いて90% 以上の相同性を持つ配列を除きました。その結果、5710 行になりました。
ファミリーとしては、アミラーゼなどが属している GH13 が最も多かったです。
1 つのタンパク質しか登録されていないファミリーは 29 個、2 個だけ登録されているのは 17 ファミリー、3 個だけ登録されているのは 17 ファミリー……でした。（➡ ファミリーによってはデータが1個～しかない）
EC 番号は、ユニークなラベルの数が 298 個でした。EC 3.2.1.4（セルラーゼ）が最も多く、上の通りの結果に。n = 1 の EC 番号は 99 個、そして、3.2.1.* の「加水分解酵素」という情報しかないものが 300 個もありました。ちゃんとやるならこれらは手動でラベルしてやるべきなのでしょうが、とりあえずこのまま行きます。
データセット中の全長の配列の長さの分布は上のようになりました。
酵素は活性ドメイン以外にも様々な付属ドメインが付いていることが多く、それを除いた活性ドメインだけにすると、配列長の分布は上の通りでした。全長 or 活性ドメインのみのどちらを使うのが良いか分からないので、とりあえず両方やってみることにします。
自前のデータで学習（Fine-tuning）させてみる
データを test set に配列を 599 個、validation set に配列を 600 個取り分け、残りを train set としました。
GitHub に用意されていたデモ用コード proteinBERTdemo.ipynb は用意されたベンチマークを for 文で全部実行する形だったのですが、これも一個づつ実行できるようにコードを書き替えて頂きました。
結果
結果はまとめの表と混合行列が出力されますが、以下まとめの表のみ示します。
全長の配列
| Model seq len | # records | Accuracy |
|---|---|---|
| 512 | 285 | 0.807018 |
| 1024 | 279 | 0.422939 |
| 2048 | 34 | 0.235294 |
| 4096 | 1 | 0.000000 |
| All | 599 | 0.594324 |
ドメインのみを抜き出したもの
| Model seq len | # records | Accuracy |
|---|---|---|
| 512 | 538 | 0.773234 |
| 1024 | 59 | 0.305085 |
| 2048 | 2 | 0.000000 |
| All | 599 | 0.724541 |
コードを動かしてみて気づきましたが、開始時のトークンサイズに 512 が指定されている場合、まず配列長が 510 以下のものだけ取り出して、サイズ = 512 で行い、次に 513~1024 、1024~2048 、と順に実行しているようです。
全長の配列を使った場合、約半数が 512 より長かったため、データセットの量は、実質半分程になっていました。
また、配列が長いものは正答率が低かったのですが、これはマルチドメインタンパク質の割合が高いため、難しかったのではないかと思います。
トークンサイズの最小値を 1024 にしたらデータが分断されることなく、全体の精度の値が上がるのではないかと考え、やってみることにしました。
開始時のトークンサイズに 1024 を指定した結果
全長の配列
| Model seq len | # records | Accuracy |
|---|---|---|
| 1024 | 564 | 0.790780 |
| 2048 | 34 | 0.352941 |
| 4096 | 1 | 0.000000 |
| All | 599 | 0.764608 |
ドメインのみを抜き出したもの
| Model seq len | # records | Accuracy |
|---|---|---|
| 1024 | 597 | 0.782245 |
| 2048 | 2 | 0.000000 |
| All | 599 | 0.779633 |
1024 の結果は、512 の結果とほぼ同等のようです。その結果、低い値が除かれて、平均の精度が少し高くなりました。
ドメインのみと全長ではほとんど精度が変わらなかったこと、ドメインのみの方が Fine-tuning にかかる時間が少し短かったことから、このデータでは ドメインのみ × 配列サイズ 1024 が今のところ最適のようです。
一つ一つの結果を見る
自前のデータセットを使った結果、数字の上ではそこそこの精度が出ているように見えました。
しかし、どの配列で誤答し、どれに正解したのかを見るため、Fine-tuning までとテストデータの評価を別々に行うようにコードを書き替えました。
本当は 「配列-正解ラベル-予測ラベル」 という形で出力したかったのですが、私のプログラミング力が足りず、入力時のラベルが内部で異なる数字に置き換えられるところに対応できなかったので、データが一行だけの csv ファイルを入力し、一行だけの混合行列を確認することで無理やり解決しました (- -；
サンプル数が少ないものの正解率
学習データに一回、テストデータに一回しか出てこない EC 番号ラベルを持つものを抽出してみました。
1/7 しか正解できず、学習データが少ないとやはり予測はむずかしいようでした。
多様化（分化）が進んだファミリーの正解率
同一ファミリー（≒全体の配列は比較的似ている）の中でも活性が異なる酵素を多く持つ GH13 ファミリーを見てみます。
GH13 の中では正解率は 0.64 でした。
不正解一覧
正解一覧
予想通りではありますが、GH13サブファミリー24 (GH13_24) など、同じサブファミリーに一つの EC 番号しかないものは正解率が100%で、GH13_20 などサブファミリー内でも異なる EC 番号が付いているものは正解率が低い結果となりました。
GH13 は多様化が進んでいるといっても、マルトース or イソマルトース のようによく似た糖を基質とするため、アミノ酸配列からその違いを見分けるのは難しいのかもしれません。
学習データに存在しないファミリーの配列
やはり、このモデルが全く未知の配列に対してどれだけ機能するかが気になるところです。
酵素のファミリーは現在でもどんどん新設されており、現在は学習データにないファミリーもあります（ファミリーGH167～GH181）。
そこで、学習データには含まれていないこれらの配列の活性を予測させてみる事にしました。
| family | True | Predicted | |
|---|---|---|---|
| GH181 | exo-α-sialidase (EC 3.2.1.18) | 3.2.1.18 | 正解 |
| GH180 | β-glucosidase (EC 3.2.1.21) | 3.2.1.78 (endo-β-1,4-mannase) |
|
| GH179 | β-N-acetylhexosaminidase (EC 3.2.1.52) |
3.2.1.49 (α-N-acetylgalactosaminidase) |
惜しい |
| GH178 | mannan endo-α-1,4-mannosidase (EC 3.2.1.-) |
3.2.1.4 (cellulase) |
|
| GH177 | exo-α-sialidase (EC 3.2.1.18) | 3.2.1.* | |
| GH176 | isoamylase (EC 3.2.1.68) | 3.2.1.* | |
| GH175 | β-glucosidase (EC 3.2.1.21) | 3.2.1.26 (β-fructofuranosidase) |
|
| GH174 | endo-α-1,3-L-fucanase (3.2.1.211) | 3.2.1.39 (endo-1,3-β-glucanase) |
|
| GH173 | β-galactosidase (EC 3.2.1.23) | 3.2.1.35 (hyaluronidase) | |
| GH172 | α-D-Arabinofuranosidase (EC 3.2.1.-) | 3.2.1.39 (endo-1,3-β-glucanase) |
|
| GH171 | peptidoglycan β-N-acetylmuramidase (EC 3.2.1.92) |
3.2.1.52 (β-N-acetylhexosaminidase) |
|
| GH170 | 6-phospho-N-acetylmuramidase (EC 3.2.1.-) |
3.2.1.52 (β-N-acetylhexosaminidase) |
|
| GH168 | endo-α-(1,3)-L-fucanase (EC 3.2.1.211) | 3.2.1.211 | 正解 |
正答率は低いものの、そもそも正解しているファミリーがある、ということに驚きました。
また、基質に N-acetyl 基があることは 3/3 で見抜けている様子です。
proteinBERT は単純な配列相同性だけでなく、ある程度配列の特徴を認識しているのかもしれません。
今後の展望？
- GH以外のファミリーも学習させてみたい。
- 糖質分解酵素とそれ以外の二値分類にしたら、機能未知タンパク質の中から糖質分解酵素を見つけることができないだろうか？
-
Nadav Brandes and others, ProteinBERT: a universal deep-learning model of protein sequence and function, Bioinformatics, Volume 38, Issue 8, March 2022, Pages 2102–2110, https://doi.org/10.1093/bioinformatics/btac020 ↩
-
ニューロペプチド： 神経細胞で作られ、情報伝達を担うペプチド。大きな前駆体タンパク質（最終産物のペプチド配列、シグナルペプチド、スペーサー配列、切断サイトから成る）として合成され、翻訳後に切断される。 ↩