magattacaのブログ

日付以外誤報

Self Distillationって何?AlphaFold2では何をしているの?

先日AlphaFold2についてのPodCastを聞きました。知らないことばかりで「あーなるほどそういう話だったのかー」と興味深かったです。専門の先生方の議論を拝聴できるのはすばらしいですね。

AlphaFold2はディープランニングの専門の方から見ても、面白い技術がたくさん使われているそうですが、ど素人にはそもそもどこが生物学で、どこが深層学習的な話なのかわからないです。

というわけで、今回は深層学習の用語らしい「Self distillation」について調べてみました。ついでにAlphaFold2での使用例もちょっと見たいと思います。

www.nature.com

1. Self distillation

Slef distillationはネットワークの学習の際に用いられる工夫だそうです。

「①Konwledge distillation(知識蒸留)」という異なるモデル間での学習結果の移植を、自分自身(同じモデル間)で行ってパフォーマンスの向上を行ったのが「②Self distillation(自己蒸留)」です。
さらにSelf-Training(自己訓練)で、学習時にノイズを追加したら精度と頑健性が上がったというのが「③Noisy Student Training」というもののようです。

  1. 文献①: Hinton, G., Vinyals, O. & Dean Jeff. Distilling the Knowledge in a Neural Network (arxiv:1503.02531v1)
  2. 文献②: Zhang, L., Song, J., Gao, A., Chen, J., Bao, C. & Ma, K. Be your Own Teacher: Improve the Performance of Convolutional Networks via Self Distilation (arxiv:1905.08094v1)
  3. 文献③: Zie, Q., Luong, M.-T., Hovy, E. & Le, Q.V. Self-training with Noisy Student improves ImageNet classification (arxiv:1911.04252v4)

AlphaFold2では「SI:1.3 Self-distillation dataset」で文献③が引用されていますが、折角なので順番に見ていきます。

1-1. 知識蒸留はモデル圧縮の技術

知識蒸留はGeoffrey Hinton*1らにより提案されたモデル圧縮(Model compression)の手法です(文献①)*2

ディープラーニングのネットワーク構造は精度を高めるために多層・複雑になっています。モデルが複雑になると特徴の学習がうまくいくようになりますが、学習したモデルを使って推論を行うコストも大きくなってしまいます。アプリやサービスで実際に運用(デプロイ)するには推論の速さ計算負荷も重要になるので、モデルの大きさがネックになります。

そこで、「別に学習の段階と推論の段階が同じモデルである必要ないよね!それぞれの段階に合わせたモデル使えばいいよね!」ってのがこの論文の主張です。Introductionでは「昆虫だってエネルギーを蓄える幼虫の段階と、移動して繁殖する成虫では全く体の構造が違うでしょ。機械学習だってトレーニングとデプロイで、それぞれの段階の目的に合わせた最適な構造にしようぜ!」っていってます。・・・たぶん。

f:id:magattaca:20210904141618p:plain

具体的には「学習用の大きなモデル」と「推論用の小さなモデル」を用意します。予め大きなモデルで高い精度が出るように学習しておき、次に、小さなモデルを大きなモデルと同じ推論ができるように学習させます。前者を教師モデルTeacher Model)、後者を生徒モデルStudent Model)と呼ぶことも多いようです。*3

小さなモデルであれば、性能が弱い端末上でも動かすことができるので、さまざまな現場で利用可能性が広がります(ex. エッジコンピューティング?)。*4

1-2. 大事な知識はソフトな知識

知識蒸留は学習結果のモデル間の転移ということがわかりましたが、どのような知識を伝えるのでしょうか? ここで大事なのは「曖昧さ(誤りの可能性)も含めたソフトな知識」となるようです。

クラス分類のタスクを考えた場合、モデルが予想するクラス(ラベル)にはHard labelSoft labelがあります。例えば、ある画像の動物を判断する場合、「猫である(1) or それ以外(0)」で予測を出すのがHard labelです。対して「猫である可能性(0.8)、犬(0.1)、ネズミ(0.04)...」といった確率分布として予測されるのがSoft labelです。

知識蒸留では不正解に対する知識も大事な知識と考えます*5。「猫だけど、猫の中ではネズミよりも犬っぽい猫」みたいな情報にも学習の結果がつまってるよね!ってことでしょうか?

で、Soft Labelにはこの不正解情報があるので、「教師モデルによるSoft Labelの予測と同じ予測を出せるように生徒モデルを学習させれば、知識が転移できる」という感じのようです。

f:id:magattaca:20210904142002p:plain

ところでニューラルネットワークのクラス分類の出力はSoftmax関数を使って出されることが多いそうです。 知識蒸留で使うのは温度付きSoftmax関数と呼ばれるもので、温度(T)が加わっています。 Tの値を大きく(温度を上げる)することで、重要視している不正解ラベルの確率を強調して学習に利用することができるそうです。「温度を上げて情報を一方から他方に移す」から蒸留なんですねー。

f:id:magattaca:20210904142026p:plain

1-3. 自己蒸留 ~Be Your Own Teacher~

知識蒸留においてやりたいことが大体わかりました。モデルを小さくできて便利な知識蒸留ですが、以下のような課題が指摘されています。

  1. 知識転移の効率が悪い:生徒モデルが教師モデルの全ての知識を受け継げるわけではなく、精度で劣りがちになる
  2. 適切な教師モデルの設定が大変:目的に合った教師モデルを別途デザインし、学習させるのが試行錯誤と時間がかかり大変

どちらも教師ー生徒で別々のモデルを使っていることで生じていそうな課題です。

「・・・それならどっちも同じモデル(自分自身)にしちゃえば良いんじゃない?」
「やってみたら精度も上がっちゃったよ!」

っていうのが自己蒸留self distillation)のようです(文献②)

下図の通り、左から右に向かって層が深くなっていく元のモデル(ResNet)を浅いセクションに分割します。浅いセクションを生徒、もっとも深いものを教師ととらえます。右から左へと戻る3色の矢印で自己(教師)の出力結果を、(生徒の)学習にフィードバックさせている様子がわかります。

f:id:magattaca:20210904142057p:plain

モデル一つで済むのは便利ですね。

1-4. どうして自己蒸留で精度が上がるの?

ところで、なぜ知識蒸留を同じモデルで行うと精度が上がるのはなぜでしょうか?

文献②中では、①平坦最小解(flat minima)、②勾配消失(vanishing gradients)の軽減、③識別力の高い特徴量(discriminating feature)の抽出、といった観点で考察が行われています。

①平坦最小解は耳慣れない単語ですが、直感的には以下の図のようになります。

f:id:magattaca:20210904142122p:plain

浅いニューラルネットワーク尖った最小解に陥りやすい一方で、深くパラメータ数の多いモデルは平坦な最小解に収束する可能性が高くなります。尖った最小解のモデルはデータのバイアスに弱く、トレーニングデータに対する精度が良くても、テストデータに対しては上手く機能しない可能性があります。一方で、平坦な最小解となっているモデルはデータセット間の多少のズレに対しても頑健性があり、よりロバストなモデルになっている可能性が高くなります。*6

著者らは「Self distillationでもこの平坦最小解をみつけられているのではないか?」と仮説をたて実験を行なっています。

以下の図の通り、Self distillationありなしのモデルを比較すると、ありではパラメータにノイズを加えても精度が下がりにくく損失が増えにくいという結果になりました。ノイズに強いのでより平坦(flat)なモデルになっているのでは?ということのようです。*7

f:id:magattaca:20210904142149p:plain

1-5. どうして自己蒸留で精度が上がるの? part 2

以上のように元文献中でも考察は行われていますが、Self distillationで精度があがるというのはやっぱり結構な謎らしく、今年の初めにはThree mysteries in deep learning: Ensemble, knowledge distillation, and self-distillation」(Microsoft Research Blog)という記事が出ています。これは以下の文献の著者自身による解説です。

「複数のモデルを組み合わせたアンサンブルで精度が上がるのは何故か?」という問いから出発し、Self distillationの精度上昇について議論されています。

著者らは画像認識のタスクにおいてデータ構造が持つ特徴に着目し、「multi-view」と名付けています。例えば、ネコを認識する場合、ニューラルネットワーク複数の視覚的特徴(ex. 耳、目、ヒゲ、肉球・・・)を見て、ネコとクラス分類をおこなっている、といった感じです。

ですが、データセットの画像全てで、分類に使われる視覚的特徴が全部写っているわけではなく、一部の特徴しか写っていないもの(ex. 顔のアップ、後ろ向き・・・)もたくさんあります。

著者らの理論によると、このようなmulti-viewデータで学習を行ったニューラルネットワークについて、「学習過程のランダムさ(ex. 初期値, random seed)に従って、視覚的特徴のサブセットを素早く学習する」性質があるそうです。*8

この性質から、異なるランダムな初期化を行ったモデルは、それぞれ異なる視覚的特徴のサブセットを学習していることが想定されます。なので、複数のモデルを組み合わせれば(ensemble)、全体として認識できる特徴が増えるのでテストでも精度が上がる、というわけです。*9

f:id:magattaca:20210904142236p:plain

以上がアンサンブルで精度が上がることの概略です。

では、モデルが自己だけのSelf distillationではどうなるでしょうか?

Self distillationの場合、異なる初期値の同じモデル(Teacher、Student)で、それぞれ異なる視覚的特徴サブセットを学習します。さらに、StudentはTeacherの学習した知識を受け継ぐ(蒸留)ように学習を行っているので、モデル2つ分の個々のサブセットのアンサンブルを学習できたかのようになります*10。したがってSelf distillationを行うことで、同じモデルを使っていてもアンサンブルのように精度が高くなる、ということのようです。

f:id:magattaca:20210904142258p:plain

1-6. かわいい生徒には負荷(ノイズ)をかけよ? ~Noisy Student~

自己蒸留がどこから来たのか大体わかりました。続いて文献③ Noisy Studentについてですが、こちらは以下の日本語記事がわかりやすいです。

ai-scholar.tech

ラベル付きデータで学習した教師モデルをつかって、ラベルなしデータの予測を行い「擬似ラベル」を付与します。 元々のラベル付きデータと、擬似ラベルのデータの両方を使って生徒モデルを訓練するのが自己訓練(Self Training)と呼ばれるモデル学習方法だそうです。

Noisy Studentでは、この生徒モデルの学習の際にさらに「ノイズをたくさんかけて」学習させ、得られたモデルを新たな教師モデルとして利用して、同じことを繰り返していくそうです。

f:id:magattaca:20210904142357p:plain

このように負荷をかけて訓練した生徒モデルは、より精度が高く、ロバスト性のあるモデルとなるようです。

Noisy Studentが通常の知識蒸留と異なる点は2つ。

  1. 知識蒸留では(通常)ノイズを使わない
  2. 知識蒸留の生徒モデルは教師モデルよりも小さい(速さを上げるため)が、Noisy studentは同じより大きいモデルを使う  

生徒モデルのキャパシティを上げつつ、ノイズを追加することで学習を難しくすることで、教師モデルよりも良いモデルを作ろうとする点で、Noisy Studentは「知識拡張knowledge expantion)」の手法と捉えることもできる、と文献では指摘されています。

2. AlphaFold2でのSelf distillation

だいたいやりたいことがわかったのでAlphaFold2でSelf distillationがどのようにつかわれているのか見てみましょう!

まずは、ずばり用語がでてくるSupplementary Information 「1.3 Self-distillation dataset」です。

2-1. Self-distillation dataset

ここでは、Self-distillationのデータセットを準備した手順が書かれています。

AlphaFold2は配列を入力として立体構造を出力するので、ラベルを実際のタンパク質の構造とすることができそうです。 従って、ラベル付きデータとして「PDBデータセットの配列と構造の組み合わせ」を使うことができ、ラベルなしデータとしては「Uniprotから抽出した配列」を使用しています。

ラベルなしデータはより正確には、UniClust30データベース(UniProtKBの配列をペアワイズ配列相同性30%レベルでクラスター化したもの)にフィルタリングをいくつかかけたものを用意しています(計:355,993配列)。例えば、UniClust30の各クラスターについて多重配列アラインメント(MSA)を計算した後、他の配列のMSAにも出てくる配列を除くといった操作を行なっています。*11

次に、ラベル付きデータをモデルで学習し教師モデルとしています。学習のこの段階ではモデルを1つだけ使用しており、「single (no re-ranking) "undistilled" model」と呼んでいます。最終的なモデルは複数のモデルのアンサンブルで、各モデルの予測をスコアでランク付けしてから出力しているようなので、それとの違いからネーミングしているようです。

この教師モデルで、先に用意したラベルなしデータセットの構造を予測して擬似ラベルとし、自己訓練に用いています。

f:id:magattaca:20210904142436p:plain

また、さらにデータの拡張(extra augmentation)として、均一にサンプリングしたMSAを蒸留データセットに追加しています。

以上が、self-distillation datasetの用意です。

2-2. Training data

データセットが準備できたのでトレーニングです(SI 1.2.4 Traning data)

レーニングに用いたデータは「self-distillaion dataset」と「Protein Data Bank」からの構造をハイブリッドしたもので、75%の確率で前者から、25%の確率で後者からとって組み合わせています。

学習時にはこのデータセットを複数回ループして使用していますが、さらに確率的なフィルタリング(SI 1.2.5)、MSA前処理工程(ブロック削除 (SI 1.2.6)、クラスタリング (SI 1.2.7))、残基クロッピング(SI 1.2.8)を加えています。なので、学習のエポックごとに、それぞれ違う標的タンパク質を、違うMSA、異なる領域の切り取り方(cropping)で処理している可能性があります。

ところで、AlphaFold2は「Noisy Student」(文献③)に類似した手順を実行した(SI 1.3)とのことですが、加えたノイズは何でしょうか?

文献③ Figure 1を再度みると、self-distillation datasetに追加したMSAデータ(extra augmentation)や、学習時実行しているドロップアウト(SI 1.11.6)などがノイズに相当しそうです。

以上をふまえて、学習の工程は以下のようになるでしょうか?

f:id:magattaca:20210904142500p:plain

2-3. Confidence metric (KL divergence)

最後に、「SI 1.3 Self-distillation dataset」で出てくる信頼性の指標(confidence metric)もみておきましょう。

これは「カルバック・ライブラー情報量Kullback-Leibler divergence, KL divergence)」によるものです。

公開されたAlphaFold2の指標は著者らがplDDTと呼ぶものですが、自己蒸留の学習の段階ではまだ開発されていなかったらしく、代わりに使っていた指標のようです。(「plDDTをこの段階で使っていたとしても同じような結果が得られただろう」とのこと)

・・・KL divergence。強そうですけどなんなんですかね??

Wikipedia - カルバック・ライブラー情報量」によると、KL divergenceは「2つの確率分布の差異を計る尺度」とのことです。「確率変数(X)に関してデータから得られる情報量の平均値」を表していることから、「情報利得(information gain)」とも呼ばれるそうです。

こんな感じ・・・

f:id:magattaca:20210904142640p:plain
みやすさのため離散確率分布の場合を引用

AlphaFold2の場合は、「予測された残基の各ペア(i, j)」について、「ペアごとの距離の分布(pairwise distance distribution)」と「その残基間隔についての参照分布(reference distribution)」との間でKL divergenceを計算しています(pairwise metric, c_ij)。

このペアごとの値(c_ij)を j について平均をとって、残基 i についての信頼性指標(confidence metirc, c_i)としています。この「残基ごとの信頼性指標大きいことと、予測精度の高さに相関がある」ことがわかったそうです。

f:id:magattaca:20210904150002p:plain

以上のような感じですが、何を計算しているのか?私には難しくてよくわかりません。

特にreference distributionがわからない…

AlphaFold2解体新書の該当箇所の注によると、「参照分布の配列間隔(sequence separation)」は、「アミノ酸配列(1次構造)上でどれだけ離れているか?」のことのようです。確かに「i と jの差分の絶対値」で表されているので、「i - j 間のアミノ酸残基数」と捉えて良さそうです。

ランダムにサンプリングした配列を使って、特定の配列間隔について分布の平均をとっているということは、「(残基の種類関係なく)一定の残基数分、アミノ酸配列上で離れた2残基が、立体的にはどのような距離をとるか?一般的な傾向の分布」を考えている、と解釈できるかもしれません。

一方で、pairwise distance distributionは「(残基の種類を含めた)特定のアミノ酸のペアに着目した場合の立体的な距離」と解釈できそうです。

従って、ここで計算されているKL divergenceは、「配列間隔だけを考えた一般的な距離(reference)」と「アミノ酸の種類も考えたペアの距離(pair-wise)」の2つの分布の差異を計る指標となりそうです。

これなら「指標の値が大きい(分布の差異が大きい)」ほど「特定アミノ酸ペアの残基間相互作用の特徴を学習したモデルによる予測結果」ということになるので、「精度の高さと相関する」というのも納得できるかもしませんね。*12

f:id:magattaca:20210904142756p:plain
PLoS ONE 6(12): e28766より図を一部改変して利用

まとめ

以上、ディープラーニングっぽい用語「Self Distillation」とそのAlphaFold2での利用箇所についてでした。賢い方々の日本語解説記事をよんで、なんとなく雰囲気はわかったものの、具体的にどういう操作をおこなっているか?というところまでは理解できませんでした。まだまだ基礎知識が足りていないです。

とりあえず「精度を上げるための学習の手法で、同じアーキテクチャのモデルで教師役と生徒役の2役をこなす知識蒸留」という雰囲気でした。

また、AlphaFold2では「ラベルありPDBデータ」を使って教師モデルの学習を行い、「ラベルなしの配列のみのデータ」に対する「構造予測(擬似ラベルの生成)」を行うことでSelf-distillationデータセットを用意する、といった利用がなされているようでした。

KL divergenceはよくわかりません!

ところで、バリバリ(?)の深層学習の論文を初めてまともに眺めましたが「ものすごい比喩を多用するなー」という印象をもちました。ひたすら数式やアルゴリズムが書いてあるだけなのかと思っていたので予想外でした。ネーミングも凝ってるしインパクトや親しみやすさも大事な分野なんですかねー???

今回も雰囲気と憶測だけで適当に書いているので間違いが多そうです。ご指摘いただければ幸いです。

ではでは!

f:id:magattaca:20210904160732p:plain
UniClustのキャラクターが可愛い

*1:AIの偉い先生らしい

*2:参考日本語記事① 「Distilling Knowledge in a Neural Network 記事② 「Deep Learningにおける知識の蒸留

*3:文献①では教師・生徒という呼び方はされていませんが、後続の蒸留に関する文献では一般的に使われている呼び方のようです。

*4:図は Alkhulaifi A, Alsahli F, Ahmad I. 2021. Knowledge distillation in deep learning and its applications. PeerJ Comput. Sci. 7:e474を利用して作成 CC-BY 4.0

*5:「dark knowledge」という呼び方をしている文献もあるようです

*6:このあたりは「損失平坦性により過適合を抑える」話として、書籍「[深層学習の原理に迫る](https://www.amazon.co.jp/深層学習の原理に迫る-数学の挑戦-岩波科学ライブラリー-今泉-允聡/dp/4000297031)」今泉允聡 著(第4章)でも取り上げられています。また、元になっている文献は「On Large-Batch Training for Deep Learning: Generalization Gap and Sharp Minimaarxiv:1609.04836v2です。

*7:テストデータでは比較しなくて良いんですかね??

*8:もう一つ「視覚的特徴では分類できない少数の残りのデータは記憶する」という性質も指摘しています。この主張から、単一のモデルの学習では取りこぼす視覚的特徴があるのは、その特徴を学習するのに十分なデータが無いからで、モデルのキャパシティが足りないわけでは無い、ということが示唆されます。学習できた一部の特徴だけでほとんど分類できてしまうので、分類できなかった残り少ない画像だけでは残りの特徴の学習には不十分、というわけです。

*9:初期値が異なることに着目しているのは、「パラメータが初期値からあまり大きくは動かない」という指摘があり、その場合「初期値が異なれば学習後のパラメータもモデル間で異なる」と考えられるからのようです。文献中で比較されている関連(類似?)する話題としてNTK(Neural Tangent Kernel)、レジームというのがあるそう。NTKは線形な特徴選択なのでディープラーニングとは話がかわってくるみたい? NTKの日本語解説記事「Neural Tangent Kernel(NTK)の概要」。よくわかりません。。。

*10:さらにStudentを次のTeacherとして学習を繰り返せば、3つ、4つ、、、と含意されるモデルが多くなることになると思います。

*11:参考:macインフォマティクスUniProtKBデータベースを3つのレベルでクラスタリングしたUniclustデータベース
文献 Milot Mirdita, Lars von den Driesch, Clovis Galiez, Maria J. Martin, Johannes Söding, Martin Steinegger,
Uniclust databases of clustered and deeply annotated protein sequences and alignments,
Nucleic Acids Research, Volume 45, Issue D1, January 2017, Pages D170–D176

*12:雰囲気だけで考えているので、間違っていそう。。。

多重配列アラインメントに含まれる構造情報について ~ Potts modelとついでにTransformer~

AlphaFold2が公開されてひと月ほど経ちましたね。直近でも複合体予測がさらにパワーアップしている等々、どんな発展を遂げていくのか楽しみです。

AlphaFold2に関しては解説の記事や文献等が色々と出ていますが、正直、門外漢にはさっぱりわかりません。

基礎知識が全然足りない!ということで、生物学の専門の方には常識っぽい「多重配列アラインメントから構造情報をとれるよ」ってのはどういう話か調べてみました。 ついでに「AlphaFold2のEvoformerって何を学習しているんですかね?」ってのを妄想します。

f:id:magattaca:20210814145426p:plain

1. 参考文献

以下の文献や記事を参考にしています。

  1. 参考① Protein structure prediction by AlphaFold2: are attention and symmetries all you need? Bouatta, N., Sorger, P. & AlQuraishi, M. (2021). Acta Cryst. D77, 982-991.
  2. 参考② AlphaFold 2 is here: what’s behind the structure prediction miracle
  3. 参考③ Temple University Ron Levy Group Lecture Note 「Potts Models and Protein Covariation」 https://ronlevygroup.cst.temple.edu/courses/2017_spring/chem5412/lectures/Lecture4_SBII_2017.pdf
  4. 参考④ 福永 津嵩. 逆イジング法の生命情報データ解析への応用 JSBi Bioinformatics Review, 1(1), 3-11 (2020)

2. だいたいこんな話?

ものすごい雑にまとめるとこういうことでしょうか???

f:id:magattaca:20210814145502p:plain

・・・怒られそうな図

3. 詳しく

順番にメモしていきます。

3-1. MSAの構造予測への利用

「参考① AlQuraishi先生らの文献」によると、アミノ酸配列からのタンパク質の構造の予測に、配列一つではなく多重配列アラインメント(Multiple sequence alignment、MSA)が用いられるようになったのは、2010年代初めごろからだそうです。背景には、ハイスループットシークエンスが一般化し、関連するタンパク質間における共変異の統計モデルが発達したことなどがあるようです。

AlphaFold2以前の((第一世代)AlphaFoldを含む)構造予測モデルでは、MSAから抽出した残基間の相互作用に関する要約統計量を利用し、残基の共進化(co-evolution)をモデル化していました。ただしこれだけでは色々と無理(物理的に不可能な構造が予測される etc.)があり、物理的手法との組み合わせが必須だったそうです。

対してAlphaFold2では(統計量ではなく)MSAそのものを利用して構造予測をおこなっています。統計物理学に基づくPotts modelを利用するのではなく、生のMSAからTransformerによって学習することで情報を抽出しているそうです。

3-2. MSAってどんな?

MSA、MSAって言うけど実際どんなのなんですかね???ということでPfamから例を貼っておきます。

Pfamは、タンパク質ファミリーのデータベースであり、アノテーション隠れマルコフモデルを用いて生成された多重配列アライメントを含んでいます(Pfam(ウィキペディア))。*1

f:id:magattaca:20210814145543p:plain

図の上半分のように、複数の生物由来の配列を並べたものです。各行がそれぞれの配列で、同じ or 類似の残基が同じ列に並ぶようにギャップが挿入されてます。 3配列以上の多重配列アラインメントは進化的に保存された配列の同定などに用いられるそうです。(シーケンスアラインメント(Wikipedia))

図の下半分のように、各配列の位置で「どのアミノ酸の出現確率が高いか?」みたいなことも分かるようです。・・・なるほど。

3-3. MSAからなぜ立体構造の情報がわかるのか?

さて、MSAの利用が大事なことはわかりましたが、なぜ立体構造情報を抽出することができるのでしょうか?これは共進化の考え方に基づいているようです。

MSAで並べられているのは互いに進化的に関係している配列であり、それぞれのタンパク質のフォールドは類似していると考えられます。 MSAの配列間で異なっている残基は、進化の過程で適当に変化したわけではなく、一定の制約がかかっている考えられます。 例えば、その残基がタンパク質の中で直接他の残基と接触しているのであれば、「残基間の相互作用の親和性を保つように変化する」という制約が考えられます。 逆に言えば、そのような変化をしている残基のペアは「3D構造において直接相互作用する距離関係」にある可能性があります。

f:id:magattaca:20210814145632p:plain

上図で上半分は共進化のイメージで、下にアミノ酸の構造をWikipediaからもってきました。酸性-塩基性、疎水性-疎水性といった相互作用を保つようなペアで変化しているという模式図になっています。

以上のように、MSAの変異のペアには残基間の距離情報が埋め込まれているので、解析すれば逆に「配列から構造を予測する問題」に使えるんじゃない?ってことみたいです。

なお、上図は下記文献を引用していますが、DCA(Direct Coupling Analysis)を利用した立体構造予測を提案した有名な論文だそうです。

3-4. 共進化解析における課題とPotts model

さて、Coevolutionary analysisやContact Predictionといった用語でやりたいことが分かってきました。

「で、実際どうやってやるの?」

「とりあえず残基のペア見つけて相関をスコアにすればいいんじゃない?」

と、なりそうですが、そう単純ではないようです。

参考⑤の文献で指摘されていますが、共進化の解析で「相互作用する残基ペアのみに着目した」方法(localな手法)では、「残基ペアの相互作用が直接のものか、間接なのか見分けるのが難しい」という課題があります。

簡単に言うと、コンタクトする残基のペアA-BB-Cがある時に、ACも一緒に変化するため、実際にはコンタクトしていなくても、A-Cという相互作用があるように見えてしまうということです(偽陽性)。*2

このような直接と間接の区別をつけるため「配列全体を考慮に入れた統計モデル」(globalな手法)が使われ、それがPotts Modelです。

ここでは、対象の配列全体(S)を生成する確率モデルP(S)を考えます。P(S)は「2つの残基間の依存関係の強さを意味するパラメータ」を、「全ての残基ペア分」含むようなモデルです。そしてデータ(MSA)からそのパラメータを学習し、学習されたパラメータを「残基のペア間の依存関係の指標」として捉えます。

ちょっと用語整理。。。

f:id:magattaca:20210814145843p:plain
イジングモデルはボルツマンマシンとも関係しているらしい

タンパク質立体構造では、データがアミノ酸(カテゴリカル、〜20種類)なのでPotts modelになる感じですね。

モデル導出の流れを、参考③をもとに一枚にまとめてみました。最大エントロピーで確率分布を定義した後、ラグランジュの未定乗数法を使って解くことで式を得ています。

f:id:magattaca:20210814145919p:plain
Potts Modelの導出

あとは上記のモデルのパラメーター Jhを求めればOKです。これは最尤推定に基づくMCMCサンプリングや、擬尤度最大化法などなどで推定できるそうです。

モデルができれば以下のようなことが分かるそうです。

f:id:magattaca:20210814151831p:plain
Potts Modelから求められることとその応用

Potts modelから得られる確率的ペアワイズカップリングは、ざっくりしたコンタクトの有無の情報やより細かい距離情報へと変換することができます。 これらはフォールディング予測の幾何的な制約として使うことができるので、Potts modelはタンパク質構造の再構築に有用な情報を持っていることになります。

3-5. もう一度AlphaFold2に戻って

以上、「MSAから構造情報を取得する」というお話の背景をざっとさらってみました。

もう一度、AlphaFold2に話を戻すと、AlphaFold2ではPotts Modelを利用するのではなく、「生のMSA」から直接情報を学習しているとのことでした。 参考①の文献ではTransformerを用いた学習により「事実上独自のPotts modelを導き出しているようなものではないか」と指摘されていました。

これがどういうことなのか?もう少しだけ考えたいと思います。

4. AlphaFold2 のEvoformerは何をやっているのだろうか?

DeepMindの資料とNatureの文献*3をみると、AlphaFold2のアーキテクチャは①Embedding、②Trunk、③Headsの3つに分かれています。

f:id:magattaca:20210814150133p:plain
AlphaFold2全体像

このうちMSAは、直接は前半①と②に関わっています。というわけで、「② Trunk(Evoformer)がPotts model(のようなもの)を学習している」ってのはどういうことか?妄想してみます。

4-1. AlphaFold2の入力とEvoformer

AlphaFold2の前半で何をしているのか、参考②の記事を頼りにざっと眺めます。

入力として配列が渡されると、AlphaFold2は2つの前処理を行います(①Embedding)。 一つは「MSAを作成」することで、もう一つは「類似の構造を持つタンパク質を見つけてtemplateとする」ことです。

templateは「どのアミノ酸が、別のどのアミノ酸とコンタクトしていそうか?」というのをモデル化した「pair representation」という表現に初期化されます。 先に見たContact analysisで作成されている2次元のマップと同じ感じなので、やりたいことがわかりますね。

MSAは進化に基づいてどのアミノ酸とどのアミノ酸がコンタクトしていそうか?という情報を暗に含み、templateは類似(と想定される)構造での明なコンタクト情報を含みます。

EvoformerではMSAについてのTransformerで学習した情報を、template由来のpair representationに渡して学習し、それをまたMSAに戻して改善する、というループを回します。 つまり、Potts Modelを利用する他のモデルであれば「MSAを一度入力として利用してお終い」となりそうなところを、 AlphaFold2では「構造情報(template)をフィードバックしてMSAの精度を改善し、もう一度構造に戻す」という作業を繰り返しながら精度を上げていく、ということが行われているようです。

「独自のPotts Modelを導く」という過程が少しずつ感じられてきました。

4-2. Transformerは全体から学習する

さて、EvoformerのベースになっているTransformer(Attention)ですが、 特徴は「着目している箇所について、それ以外の全体も直接参照した」学習が行われることにあるそうです。

自然言語処理の例では、一つの単語(トークン)について、文全体の各トークンを参照して重みづけ和を計算することになるようです。

f:id:magattaca:20210814150234p:plain
Attention is All You Need

対して、Recurrent Neural Network(RNN)は前から順々に処理するので、手前のトークンの情報に基づいており、また、離れたトークンからの情報が薄まってしまう問題があるそうです。文全体を参照するTransformerの方がより文脈を考慮した学習ができる、という仕組みのようです。

Transformerの「文全体を参照する」仕組みですが、タンパク質のMSAに適応すると「着目する残基について、配列全体の情報を直接参照」することに相当しそうですね。

ところで、MSAのPotts modelは、ペアワイズの解析について直接の相互作用と間接の相互作用とを区別するために、大域的な配列全体を考慮するモデルとして導入されたのでした。配列全体という共通点に着目すると、MSAにTransformerを適用した学習が「独自のPotts Modelを導く」ようなもの、というのも納得できる気がしますね!

・・・しませんか。。。すみません。

4-3. Transformerの用語

ついでにTransformerの用語などまとめておきます。たぶんこういう話。。。*4

f:id:magattaca:20210814150421p:plain
Attentionの用語と作業

Self-Attentionはqueryがkeyとvalueが同じとなっていて、一見、上図の手順とあわないように見えます(queryが1単語だとそもそも次元が一致しない)。これは、文の単語をまとめてバッチでqueryに投げてしまえば、keyやvalueと同じ次元になるからOK、ってな感じだと思います。

もう少し具体的な処理はこんな感じ。。。

f:id:magattaca:20210822202903p:plain
行列形式のAttentionの演算

4-4. EvoformerのMSAに対するAttention

ついでのついでに、EvoformerでMSAに使われているAttentionのアルゴリズムを貼っておきます。

行方向(row-wise)と列方向(column-wise)の2種類があります。どちらもMulti-Head Attentionで、またゲート付き(gated self-attention)です。

ゲートは過去の状態を伝搬させて考慮する仕組みで、RNN(Recurrent Neural Network)における「入力された情報を長時間保持することができない」という課題を解決するために提案された、LSTM(Long-Short Term Memory)というモデルで使われているもののようです。

ではまずcoumn-wiseから。。。

f:id:magattaca:20210815151432p:plain

Nhead=8なので8個のAttentionを使っています。つまり(query, key, value)の組み合わせが8個です。上付きのhでそれぞれのheadを表しています。

(query, key, value)の表記は(q, k, v)です。またこれらの次元(先の説明でd次元になっていたところ)は c (channel,C=32)と表現されています。この表記を踏まえると、入力されたMSA(m)からAttentionを求める流れは先に見たとおりですね。

mをバイアスベクトル無しの線形変換(LinearNoBias)することで(q, k v)を用意しています。q, kとcからsofmaxで重みaを求め、これにvを掛けた重み付け和(Σav)をとっています。

少し違うのは、さらにゲート(g)とのアダマール積(要素積)がとられていることで、これがgated self-attentionとなっている部分だと思います。

また、Multi-Head Attentionは最後に複数のAttentionを一つにまとめる必要があるので、concatで連結して出力しています。

なんとなく操作がわかったぞ!

次に、row-wiseです。こちらはもう少し複雑です。 残基ペアのAttention weightを構築していて、pair representationから得られる情報を追加のバイアスとして取り込んでいます。 これによりペア表現とMSAとの間に一貫性が生まれるとのこと。。。なんのこっちゃ?

f:id:magattaca:20210815151500p:plain

新しく出てきた記号はzで、これは入力のpair representationを表します。 coulmn-wiseとの違いは、Attentionを求める際のsoftmaxで、ここにpair representation由来のバイアス(b)が加えられています。あとはだいたい一緒。

用語をおさえるとなんとなく操作がわかるもんなんですねー。操作の理由はさっぱりわかりませんが。。。

5. まとめ

以上、MSAから構造情報を取得する仕組みと、ついでにAlphaFold2のMSAに関わる部分(Evoformer)について眺めてみました。

AlphaFold2は分野の歴史と発展を踏まえたうえで行きついたもののようですが、門外漢には正直さっぱりわからないです。 賢い方々の解説を色々とつまみ食いしてお勉強してみましたが、背景が膨大すぎて難しいです。

とりあえず配列全体を考えるPotts Modelという方法と、文全体を参照するTransformerにうすーく繋がりが見えた気がするのでヨシ!

今回もいろいろと間違いが多そうなのでご指摘いただければ幸いです。

6. 追記:AlphaFold2の学習で行われるMSAのmaskについて

MSAとTransformerに関する項目として、AlphaFold2の学習(training)で行われているマスクmask)について書き忘れていました。

Transformerをベースとした自然言語処理の非常に性能の良いモデルとして、BERT(Biderectional Encoder Representations from Transformers)というものがあるそうです。

  • Devlin, J., Chang, M.-W., Lee, K. & Toutanova, K. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. arXiv:1810.04805

これもGoogle。。。

このモデルの学習は2段階に分かれており、①大規模な文章コーパスラベル無しデータ)から汎用的な言語のパターンを学習する事前学習と、②個別のタスクのラベル付きデータを用いて、そのタスクに特化させるように学習するファインチューニングからなります。

マスクは事前学習(①)で、ラベル無しデータから学習を行うための工夫として用いられています。ざっくり言うと「隠された単語を当てるゲーム」です。文章の一部をランダムにマスクし、「そのマスクされた単語がなんだったのか?、周りの単語から予測する」というタスクを用いて学習します。*5

f:id:magattaca:20210815151832p:plain

このマスクがAlphaFold2でも使われています。多重配列アラインメント(MSA)の残基の一部分をランダムにマスク(or 変異)し、元々の残基がなんだったのか予測するタスクを学習させています。

関連するMSA入力データの準備の説明について貼っておきます(MSA Clustering (Suppl. 1.2.7))。

入力データを準備する際に、MSAが大きい(sequence数が多すぎる)場合、計算コストを減らすためにランダムに配列を選択して一定の数まで減らします。そのままでは選ばれなかった配列の情報が捨てられてしまうので、選んだ配列をクラスター中心としたグループ化を行うことで残りの配列の情報考慮できるようにしているようです。

f:id:magattaca:20210815151956p:plain

マスクされたMSAに対する予測のタスクは「Supplementary Material: 1.9.9 Masked MSA Prediction」に書かれています。

f:id:magattaca:20210815152026p:plain

このようなタスクを学習することで何が嬉しいか?

AlphaFold2文献本文によると「特定の相関の統計を特徴量にハードコーディングしなくても、ネットワークが系統学的な関係や共変動の関係を学習して解釈できるようになる」ことが期待されているようです。

私は「マスクされた残基の予測」で、Pfamでみたアミノ酸出現確率の図のイメージを思い浮かべました。

f:id:magattaca:20210815152104p:plain
再掲

ネットワークが「配列のどの位置にどの残基が出現しそうか?」を予測できるということは、ある配列全体が生成する確率についての情報を学習できていることになりそうです。

もう一度Potts Modelに戻ると、「対象の配列全体(S)を生成する確率モデルP(S)」として導入されたものでした。

なんだかAlphaFold2が学習の結果、独自のPotts Modelを導いているのに納得感がでてきませんか?・・・やっぱりだめですか。すみません。

おしまい!

*1:参考 TOGO TV 「Pfamを使ってタンパク質のドメインを調べる 2017

*2:因果推論の交絡因子と似た感じのお話でしょうか?

*3:Highly accurate protein structure prediction with AlphaFold

*4:Transformerについては以下を参考にしました。

PROTACの対となるか?〜脱ユビキチン化を促進するキメラ化合物〜

AmgenのKRAS-G12C阻害剤 Sotorasib(AMG510)が承認されましたね!*1 ケミカルバイオロジーの発展に端を発する化合物で、難攻不落と言われていたKRASを標的とする医薬品が誕生したのはとても素晴らしいです。*2

低分子医薬品の世界も、標的となるタンパクが広がったり、PROTACにMolecular glueといった新しいアプローチが実用化に向けてどんどん開発されていて眺めていて楽しいです。

で、コンセプトが明確で、なるほど面白い!そうなるかー!と思ったアプローチを知りました。「Targeted protein stabilization (TPS)」というもので、雑にまとめると「PROTACの逆をいくキメラ化合物」です。

bioRxivに論文が公開されていたのでご紹介。*3

www.biorxiv.org

Targeted protein stabilization

さて、こちらの論文で提案されているアプローチはとても分かりやすいです。

「ユビキチン-プロテアソーム系で分解されるタンパク質に対して、脱ユビキチン化酵素を働かせることで、 分解を抑制して安定化し、機能できる状態を維持しよう」ってな感じです。

変異などにより上手くフォールディングできないタンパク質はユビキチン化を受けプロテアソームで分解されます。 PROTACやMolecular glueはこのシステムを利用して、壊したいタンパク質のユビキチン化を誘導、分解促進するメカニズムでした。

f:id:magattaca:20210606031456p:plain
今更感あふれる図

で、上の図には出てきませんが、 ユビキチン-プロテアソーム系を調節する酵素として脱ユビキチン化酵素(deubiquitinating enzyme、DUB)があります。 その名の通り「タンパク質からユビキチンを切断除去するプロテアーゼ」です。*4

Wikipediaの模式図を貼っておきます。上図の逆向きの反応です。

f:id:magattaca:20210606031533p:plain

さて、著者らのアプローチはこうです。 「PROTACがユビキチンリガーゼを標的タンパクに近接させることで分解を促進するのなら、 脱ユビキチン化酵素を近接させれば分解を抑制・安定化できるだろう。そのような近接化に働くキメラ化合物Deubiquitinase-Targeting Chimera(DUBTAC)と呼ぼう。」

f:id:magattaca:20210606035710p:plain
bioRxiv 2021.04.30.441959 Fig 1aより

比較するとこんな感じの話?

f:id:magattaca:20210606031720p:plain
やっつけの図

なるほど!!でもどうやってそんな化合物作れば良いんですかね???

DUBTACのデザイン戦略

DUBTACによるTPSを達成するためには大きく4つの課題があります。

  1. DUBファミリーのうち、どの脱ユビキチン化酵素を利用するか?
  2. また、その酵素に結合するユニットをどう見つけるか?
  3. 安定化する標的タンパク質として何を選ぶか?
  4. また、その標的タンパク質と結合するユニットをどう見つけるか?

著者らは、前2つの課題に対してはケモプロテオミクスを利用することで、 後2つ課題については嚢胞性繊維症(Csytic Fibrosis, CF)の原因CFTRを標的とすることでコンセプト検証(Proof Of Concept)を行っています。

Nomura Research Groupについて

さて、ここで著者情報です。今回のbioRxiv論文は University of California, BerkeleyDaniel K. Nomura教授らのグループによる報告です(Nomura Research group HP)。 所属としてNovartis-Berkeley Center for Proteomics and Chemistry Technologiesとも記載されています。

こちらはNovartisとUC Berkeleyとが協力し、 「"undruggable"な標的に対する低分子化合物を開発するための次世代の技術」を開発するための組織とのことです*5。新しいモダリティとして深く製薬企業が関わった産学連携の成果なのかもしれません。*6

Nomura研究室の得意とするところは「ケモプロテオミクスを利用した共有結合性リガンドの発見」のようです。

より具体的には「活性ベースタンパク質プロファイリング(Activity-Based Protein Profiling, ABPP)」を利用することで、 「"undruggable"な標的タンパク質に、リガンドが結合可能なユニークなホットスポットがあるか?」、 「ポケットの探索およびそのリガンド(functional covalent ligand)の取得」ができるそうです。*7

今年3月、Acc.Chem.Resに以下のようなレビューも出版されています。

Reimagining Druggability Using Chemoproteomic Platforms Jessica N. Spradlin, Erika Zhang, and Daniel K. Nomura Acc. Chem. Res. 2021, 54, 7, 1801–1813 https://pubs.acs.org/doi/10.1021/acs.accounts.1c00065

では、具体的にどのようにDUBファミリーの探索を行ったのか?、見てみましょう。

DUBファミリーからのOTUB1の選別

さて、DUBファミリーのうち今回の目的に適した脱ユビキチン化酵素の候補をどのように選別すれば良いでしょうか?

著者らは66種の酵素について、ABPPによるケモプロテオミクスのデータを集めて解析を行いました。 このデータは「酵素システイン残基(Cys)について、プローブと反応する残基の位置をラベリングする」というプロファイリングを行ったものです。

目的に合う酵素は以下のような条件を満たすCys残基をもつものです。

  1. アロステリックな位置 … 触媒ドメインのCysでは、TPSで利用したい酵素本来の活性が失われてしまう
  2. 同一酵素内での選択性 … 同じ酵素の複数のCysのうち、特定のCysがよりラベル化されていればホットスポットの可能性が高い
  3. 十分な反応性を示すCys残基 … 複数のケモプロテオミクスデータセットで繰り返しプローブと反応し、ラベリングされているもの

要するに、「触媒活性を邪魔せず、かつ、リガンドがアクセスしやすく反応性が高いCys残基」があればよい、という感じです。

以上を目標に解析を行った結果、OUTB1という脱ユビキチン化酵素が、 触媒C91とは別の箇所C23に適したCys残基を持つことがわかりました。

f:id:magattaca:20210606032603p:plain
bioRxiv 2021.04.30.441959 Fig 1b,c,dより

ところで、シグナル分子としてのユビキチン修飾は多様性があり、ポリユビキチン鎖の連結の仕方によってシグナルの意味合いが異なることが知られています。 プロテアソーム分解のシグナルとなるのは、「ユビキチンの48番目のリジン」を介して連なったポリユビキチン鎖(K48鎖)です。

f:id:magattaca:20210606032724p:plain
ユビキチン修飾の多様性とシグナル

幸いなことに、OUTB1は発現量が高く、かつK48鎖を標的として切断する脱ユビキチン化酵素でした。 まさしく今回の目的にピッタリな酵素だったわけです。

共有結合性リガンドの探索

標的のDUB、OUTB1が定まったので、この酵素に対して適した共有結合性のリガンドを探索しています。ここでもABPPが活躍しています。

探索結果、共有結合部位としてアクリルアミドユニットをもつ化合物EN523を見出しています。 EN523はC23選択的に反応する化合物で、OUTB1の触媒ドメインのCys(C91)との反応はみられませんでした。

EN523の結合によりOUTB1の脱ユビキチン化活性が失われないことは別途in vitro再構成系で確認しています。

f:id:magattaca:20210606032811p:plain
bioRxiv 2021.04.30.441959 Fig. 2a,b,cより

これで、DUBTACデザイン戦略のうち、脱ユビキチン化酵素側の課題2つがクリアできました。 つづいてPOC取得に向けた標的タンパク質とリガンドの選択です。

変異CFTRと低分子医薬品 Lumacaftor

DUBTACを開発するモチベーションとして、疾患の中にはミスフォールドする変異タンパクしか作れなくなり、 ユビキチン-プロテアソーム系で分解されて、機能を発揮できなくなってしまうことでおこるものがある、ということがあげられます。

そのような疾患の例として、CFTRという遺伝子の変異でおこる嚢胞性繊維症(Cystic Fibrosis、CF)があります。 CFでは複数の変異が知られていますが、中でも最も頻度の高いPhe508の欠失したΔF508-CFTRでは、 タンパク質が不安定化し、K48ポリユビキチン化修飾-分解されてしまいます。

ΔF508-CFTRに対してはすでに低分子医薬品Lumacaftorが開発されています。 LumacaftorはCFTRに結合して安定化し、フォールディングを助けるケミカルシャペロンと呼ばれる種類の化合物です。

f:id:magattaca:20210606032910p:plain

しかしながらLumacaftorを用いても完全にCFTRタンパクをレスキューできるわけではなく、多くが分解されてしまうそうです。 つまり、「タンパク質本来の機能を妨げないが結合するリガンド」があり、尚且つまだ効果が不十分というわけです。

このリガンドをDUBTACの標的タンパク側の結合ユニットとして使ってキメラ分子を作成すれば、 さらなるタンパク安定化効果が期待でき、コンセプト検証できるのでは??、ってなわけです。

DUBTACの合成とTPS効果の検証

さて、DUBTACのコンセプト検証に向けた材料が出揃いました。

脱ユビキチン化酵素 OUTB1と結合するリガンド EN523と、標的タンパクCFTRに結合するリガンド Lumacaftorを長さの異なるアルキル鎖(C3、C5)でつなぎ、 2つの化合物NJH-2-056NJH-2-057が合成されました。

f:id:magattaca:20210606032947p:plain
bioRxiv 2021.04.30.441959Fig. 3a,bより

これらの化合物により変異CFTRの安定化は達成されたでしょうか? ΔF508-CFTRを発現するヒト気管支上皮細胞株(CFBE41o-4.7)をもちいてその効果を確かめています。

f:id:magattaca:20210606033142p:plain
bioRxiv 2021.04.30.441959 Fig. 3e,fより

上図の通り、リンカーのアルキル鎖C5としたNJH-2-057において、Lumacaftor単独では観測できなかったCFTRの安定化、増加が確認されました。

興味深いことに、アルキル鎖を短くしたNJH-2-056では効果が確認できず、リンカーの長さの影響を受けることがわかりました。

論文ではこの他に「NJH-2-057が狙い通りのDUBTACとして機能しているか?」メカニズムに関する実験や、「安定化効果が標的タンパク(ΔF508-CFTR)選択的なものか」 といったことに関する実験も行われています。

以上、bioRxivの論文の内容でした。

感想

さて、DUBTACによるTargeted protein stabilization (TPS)、コンセプトが明快で効果検証されており、非常に面白いアプローチだと思いました。

共有結合性のリガンドをうまく活用して、新しい機能を持つハイブリッドな化合物をつくってしまうあたり、 ケミカルバイオロジーの面目躍如という感じでとてもワクワクしますよね!

ただし今回の論文は細胞レベルでの効果検証までなので、実用化にはまだまだハードルはたくさんありそうです。

  • 分子量800近くあるけど、膜透過や薬物動態はどうなの?
  • DUB側が共有結合性だけど、効果の持続とか、DUBの半減期の影響は?
  • ユビキチンを付け外しが繰り返されると小胞体ストレスかかりそうだけど大丈夫?  
  • あれ、そもそもLumacaftorのシャペロン効果が今回の系で確認できてないのは良いんですか?

などなど、私には一読では分からないことがいっぱいです。詳しい方教えてください。

期待

で、いちゃもんつけつつワクワクしているのは、アカデミックな検証だけでおわらないことに期待しているからです。

まず、Novartis社が関わっていそうという点。嚢胞性繊維症の治療薬というと私はVertex社のイメージが強いです。*8 ところが不勉強なので、今更ながらNovartis社がCFTR標的薬 icenticaftor(QWD251)を開発されていることを知りました。

Discovery of Icenticaftor (QBW251), a Cystic Fibrosis Transmembrane Conductance Regulator Potentiator with Clinical Efficacy in Cystic Fibrosis and Chronic Obstructive Pulmonary Disease Journal of Medicinal Chemistry Article ASAP DOI: 10.1021/acs.jmedchem.1c00343 https://pubs.acs.org/doi/abs/10.1021/acs.jmedchem.1c00343

多種多様な研究開発・新規技術を実用化する力のあるメガファーマが治療薬開発に取り組んでいると聞くと、新たなアプローチが可能になるのではないか、と期待してしまいます。

また、逆にベンチャー企業の例として、同様のアプローチを掲げているStablix Therapeuticsという会社があるそうです。 最近できたばかりのようです。*9(シリーズAってなんですか??)

今回の記事で取り上げた論文の著者らとは別のグループの研究者らのようですが、HPで掲げられているアプローチや標的を見る限り類似の戦略をとっているようです。

興味深いのはCo-Founderの一人、Scott A. Kanner博士が以下の論文のファーストオーサーという点です。(課金してないので読んでません。すみません。)

Kanner, S.A., Shuja, Z., Choudhury, P. et al. Targeted deubiquitination rescues distinct trafficking-deficient ion channelopathies. Nat Methods 17, 1245–1253 (2020). https://doi.org/10.1038/s41592-020-00992-6 Targeted deubiquitination rescues distinct trafficking-deficient ion channelopathies | Nature Methods

こちらの論文はAbstractを見る限り人工脱ユビキチン化酵素の開発のようで、キメラ分子のアプローチとは異なりますが、 メカニズムに知見のある研究者が深く関わっているベンチャー企業がどんな化合物を作り出すのか?、楽しみです。*10

まとめ

以上、PROTACの逆をいくかのようなアプローチ、脱ユビキチン化酵素を利用するキメラ分子によるタンパク安定化の文献紹介でした。 私は専門家ではないので、このアプローチがどれくらい実用可能なものか?全く想像がつきません。*11

ですが、非専門家の私でも「なるほど!面白い!」と思える分かりやすい戦略を打ち出していて、それを実際の化合物構造として実現化している、という点で非常に興味深かったです。 また、異なる複数のグループが同じアプローチを目指しているようなのもびっくりしました。私が不勉強なだけでよく知られた手法なのでしょうか?どんどん競争が激しくなりそうです。

いずれにせよ、低分子(中分子?)化合物でできるアプローチもどんどんと進化していて楽しいですね!

毎度毎度、浅い理解で書いているので間違いが多そうです。ご指摘いただけると幸いです。ではでは!

*1:FDA grants accelerated approval to sotorasib for KRAS G12C mutated NSCLC

*2:Fragments in the clinic: AMG 510

*3:CC-BY-NC-ND 4.0 International license

*4:Wikipedia 脱ユビキチン化酵素

*5:Novartis プレスリリースよりNovartis and UC Berkeley collaborate to tackle 'undruggable' disease targets

*6:この辺りの米国の事情は私は詳しくないので間違っていたらすみません。

*7:ABPPについてはChem-Stationさんの記事「活性ベースタンパク質プロファイリング Activity-Based Protein Profiling」の説明と図が分かりやすいです。

*8:本当に格好いいですよね。尊敬しています。

*9: Stablix Therapeutics Launches with $63 Million Series A Financing

*10:なお、上記論文は今回のbioRxivの文献にも引用されています。むしろ私はこちらの文献からbioRxivに辿り着きました。

*11:POCもまだ初期のin vitro細胞系で、そもそもプレプリントの文献です。

RDKit にコントリビュートした

RDKitの2021.03 versionがリリースされましたね!

お気づきの方はいらっしゃるだろうか・・・

f:id:magattaca:20210401233809p:plain

ついに私もRDKit コントリビューターになったよ!

といってもめちゃくちゃ些細な修正をマージして頂いただけです。せっかくなのでご紹介。

最新版のインストール

まずは最新版をインストール。以前と少しコマンドが変わっています。(RDKit documentation installation

# チャンネルがconda-forgeになってる
conda create -c conda-forge -n my-rdkit-env rdkit

作成した環境「my-rdkit-env」にjupyter notebookもインストール。カーネルを選べるように追加しました。*1

# 仮想環境をactivate
conda activate my-rdkit-env

# jupyterをインストール
conda install jupyter

# カーネルを追加
ipython kernel install --user --name=my-rdkit-env --display-name=my-rdkit-env

Jupyter notebookを起動して確認してみます。

from rdkit import rdBase
print('rdkit version: ', rdBase.rdkitVersion)
#  rdkit version:  2021.03.1

インストールできました!*2

修正内容

マージして頂いた修正はMolFromSequenceに関するものです。MolFromSequenceを使うと、アミノ酸一文字表記の配列(シークエンス)からMolオブジェクトを構築することができます。

以前の記事で取り上げていますのでよろしければご参照ください。

magattaca.hatenablog.com

こちらの配列、オプションの引数「flavor=1」とすることでD-アミノ酸も扱えるのですが、なぜかD-セリンD-チロシンだけが含まれておらず利用できませんでした。

無いなら付け加えてもらえばいいじゃない! プルリク!プルリク!

github.com

まずは比較のためL-セリンL-チロシンを確認しておきます。大文字SYでOKです。

from rdkit import Chem
from rdkit.Chem import Draw

L_Ser_L_Tyr = Chem.MolFromSequence('SY')

Draw.MolToImage(L_Ser_L_Tyr)

f:id:magattaca:20210401233841p:plain

さて、問題のD-セリンD-チロシンはどうでしょうか??

flavor=1として小文字syを使えば指定できます。

# D-セリンとD-チロシン
D_Ser_D_Tyr = Chem.MolFromSequence('sy', flavor=1)
Draw.MolToImage(D_Ser_D_Tyr)

f:id:magattaca:20210401233908p:plain

先ほどのL体と立体が逆になった構造が生成されました!

以前のRDKitのバージョンでは、これらの配列は空のオブジェクトとなってしまっていました。きちんと動いてよかった!

おわり

以上、重箱の隅をつつくような些細なお話でした。

しょうもなさすぎて蹴られるかと思いましたが、マージしていただけて嬉しかったので記事にしてしまいました。

RDKitコミュニティー優しくて好き!!

「ベイズ反応最適化」ツールで遊んだ話 ~EVML, auto-QChem, EDBO~

前回、化学合成のための反応条件最適化ツールの文献をご紹介しました。

Shields, B.J., Stevens, J., Li, J. et al. Bayesian reaction optimization as a tool for chemical synthesis. Nature 590, 89–96 (2021).

doi.org

こちらの文献、データ・コードともに公開してくださっています。

せっかくなので今回はこちらで遊んでみたいと思います。

EVMLのゲームを自分で遊んだ後、EDBOに従って遊んだらどうなるか検証してみます。

反応最適化ゲーム(EVML)

概要

まずは、EVML(Expert versus Machine Leaning)(GitHub)です。

これは、開発された反応条件最適化ツール(EDBO)の性能を比較するためのWebゲームアプリです。
化学者に以下のような条件最適化ゲームをしてもらい、その結果をソフトウェアと比較し、ベンチマークとしています。

Reaction optimization game
ボスから司令のあった化合物がパラジウム触媒による直接アリール化で合成できることがわかった。
パラメータ(base, Ligand, solvent, concentration, temperature)を変えて、反応収率を最適化したい。
ボスから与えられた期限は1ヶ月。1日にこなせる実験のキャパシティーは5実験(1 batch)。
つまり最大100実験のうちに最適化しなければならない(全組み合わせの約6%)。
使えるのは、自分の知識と、各実験から得た情報のみ(ネットや本をみたらダメ)。
さあ、やってみよう!*1

f:id:magattaca:20210223215648p:plain
EVMLの反応探索空間

このゲーム、実際の実験データをつかっているのがすごいですよね。

正直、配位子(Ligand)は色々開発されすぎててよくわかりません。・・・トリフェニルホスフィンの安心感。

塩基(Base)も溶媒(Solvent)も研究室ではみなかったものが多いです。製薬企業さんではよく使われるものなんでしょうか???

やってみる

早速あそんでみましょう!

GitHubのファイル(aws_setup.TXT)には、アマゾンウェブサービスAWS)でのインストール方法が書いてありますが、 よくわからないのでMacで試します。

RShinyパッケージでつくられたWebアプリのようなので、Rで実行すれば良さそうです。

以下のRのパッケージが必要そうなのでインストールしました。*2

  1. shiny
  2. dplyr
  3. ggplot2
  4. reshape2
  5. rmarkdown
  6. DT

EVMLのGitHubからフォルダをダウンロードしてきたら、「arylation_appフォルダ」の「app.R」ファイルを実行すればOKです。*3

実行すると、アプリが起動し以下のようなページが開きました!

f:id:magattaca:20210223215805p:plain
Rで実行するとアプリが開く

左側のパネルに、ゲーム利用者の情報をいれたり、反応条件をいれたりします。
プルダウンから条件を選んで「add」し、5つ条件を選んだら「Run Experiments」をクリック!

f:id:magattaca:20210223215905p:plain

選んだ条件の収率(yield)は「Experimet Results」タブに表示されます。 収率を見ながら条件を変えて最適化を目指しましょう!

GitHubの「results」フォルダ「analysis.ipynb」をみると、50人の化学者による結果がまとめられています。

こんな感じです。下図、上が全体で、下3つは属性で分けたサブセットです。

f:id:magattaca:20210223215931p:plain
Expertの結果(EvML GitHubの図表より)

全体では、5バッチ前後で最高収率(100%近く)に到達している方が多い印象でしょうか?

所属の違いをみるとfacultyEngineerは比較的直線的な改善を見せているのに対し、 Med ChemistProcess Chemistは横ばいとなっているところがあるのが面白いですね。 条件検討の進め方にも考え方の違いがありそうです。

ベイズ最適化ソフト(EBDO)で50回行った結果も紹介されています。以下の通りです。

f:id:magattaca:20210223220015p:plain
EDBOの結果(EvML GitHubの図表より)

こちらも5バッチ前後で最高収率に到達しているものが多そうです。ソフトは毎回Maxに到達するまで検討しつづけるのが、ヒトとの違いが見えて面白いですね。

詳細な分析と結果の考察は文献本文をご参照ください。

上から目線でコメントしましたが、私は100%近く(>99%)になるまで7バッチでした。完全にソフトに負けています。
文献で一度結果を読んでいたはずなのに。。。

合成反応の記述子DB(auto-QChem)とベイズ最適化ツール(EDBO)

続いて、ベイズ最適化を実施するためのパッケージauto-QChemEDBOです。 前者は、関連する化学物質を記述子に変換・保管するデータベースで、後者はベイズ最適化パッケージです。

EDBOのGitHubexamplesexperimentsディレクトリのipynbファイルには、 論文中のデータや解析がコードとともに掲載されています。ありがたいですね!

以下の記事はだいたいそこからのコピペ切り貼りです。

この記事の目標は、上記の反応最適化ゲームのデータを題材にして各ツールを使ってみることです。 論文の記述子は量子化学計算(Gaussian)によるものですが、Gaussianは使えないしよくわからないので、 より一般的なケモインフォの記述子(Mordred)を用いたいと思います。

auto-QChem

まずは、合成反応を入力に変換し保管するauto-QChemです。MongoDBタイプのデータベースで、ユーザーガイドやインターフェースについては以下で確認できます。

  1. auto-QChem GitHub
  2. auto-QChem GitHub documentation
  3. auto-QChem Database User Guide
  4. auto-QChem landing page *4

auto-QChemは名前の通り量子化学(Quantum Chemistry)記述子向けで、計算の実行を簡単にする(入力ファイルの作成と出力の解析)というものです。 Mordredを使うのであれば不要ですが、せっかくなのでさわってみます。

インストール

インストールの仕方はGitHubinstall.mdに書いてあります。

condaで仮想環境を作って、関連するパッケージをインストールした後、GitHubから落としてきたauto-QChemのsetup.pyを使えばOKです。

# 新しい仮想環境を作ってactivate
conda create --name autoqchem python=3.7
conda activate autoqchem

# パッケージのインストール
conda install jupyter pandas scipy matplotlib pymongo pyyaml fabric xlrd appdirs openpyxl
conda install -c conda-forge openbabel=2.4.1 # v.3.0.0はよくないらしい
conda install -c rdkit rdkit
python -m pip install imolecule==0.1.13   #3dで描画するのに使う

ついで、git cloneなどでauto-QChemをローカルにもってきたあと、保存したディレクトリに移動し、setup.pyを実行します。

python setup.py install

おしまい!*5

分子の取り扱い(moleculeクラス)

auto-QChemでは分子はmoleculeクラスで扱います。OpenBabelOBMolクラスを継承したクラスになっています。

デフォルトではmoleculeクラスの入力は分子のSMILESで、3次元構造(conformer)を構築します。 コンフォメーションはOpenBabelのOBConformerSearch()で、遺伝的アルゴリズムを使って生成されているようです。

試しにリガンドを一つ読み込んでみましょう!

論文でフィーチャーされてたCgMe-PPhというリガンドです。PubChemではこちらの名前では引っかからずmeCgPPhがシノニムにありました(PubChem CID: 53427790

from autoqchem.molecule import molecule

# Cg-MePPhのSMILESから分子を構築
smiles_str = "CC12CC3(OC(O1)(CC(O2)(P3C4=CC=CC=C4)C)C)C"
mol = molecule(smiles_str)
print(type(mol))
# <class 'autoqchem.molecule.molecule'>

draw()という関数が用意されており、jupyter notebook上でぐりぐり動かせる3D描画が見られます。

# mol.draw()

f:id:magattaca:20210223220249p:plain

かわいい。PubChemと描画スタイル同じですね。PubChemもimolecule使っているんでしょうか?

コンフォマーを複数発生させることもできます。また、loggingでログを表示させられるので、動作確認しておきたいときは表示させても良いかもしれません。

# ログの表示
import logging
logging.basicConfig(level=logging.INFO)

# 最大30のコンフォマー作成
mol_conf = molecule(smiles_str, max_num_conformers=30)

# INFO:autoqchem.molecule:Initializing molecule with canonical smiles: CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C
# INFO:autoqchem.molecule:Creating initial geometry with option 'best'.
# INFO:autoqchem.molecule:Initial geometry created successfully.
# INFO:autoqchem.molecule:Conformer Search generated 6 conformations of CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C molecule

6つのコンフォマーができたみたいです。comformer_numで各コンフォマーにアクセスできます*6

f:id:magattaca:20210223220311p:plain

違いがわからない。。。

その他の使い方や、Gaussianのインプットファイルの作成方法等はFunctional documentationをご参照ください。 また、データベースのインターフェースはFlaskを使っているようですがよくわからないので飛ばします。

Mordredで記述子計算

ベイズ最適化パッケージに進む前に、入力の記述子を用意しておきたいと思います。必要な分子(Ligand、Base、Solvent)のSMILESは以下です。

Ligand_smiles = ["CC(C)C1=CC(C(C)C)=C(C(C(C)C)=C1)C2=C(P(C3CCCCC3)C4CCCCC4)C(OC)=CC=C2OC", "CC(C)(C)P(C1=CC=CC=C1)C(C)(C)C", "CN(C)C1=CC=CC(N(C)C)=C1C2=CC=CC=C2P(C(C)(C)C)C3=CC=CC=C3", "P(C1CCCCC1)(C2CCCCC2)C3CCCCC3", "P(C1=CC=CC=C1)(C2=CC=CC=C2)C3=CC=CC=C3", "CC(C1=C(C2=CC=CC=C2P(C3CCCCC3)C4CCCCC4)C(C(C)C)=CC(C(C)C)=C1)C", "P(C1=CC=CO1)(C2=CC=CO2)C3=CC=CO3", "CP(C1=CC=CC=C1)C2=CC=CC=C2", "CC(OC1=C(P(C2CCCCC2)C3CCCCC3)C(OC(C)C)=CC=C1)C", "FC(F)(F)C1=CC(P(C2=C(C3=C(C(C)C)C=C(C(C)C)C=C3C(C)C)C(OC)=CC=C2OC)C4=CC(C(F)(F)F)=CC(C(F)(F)F)=C4)=CC(C(F)(F)F)=C1", "C[C@]1(O2)O[C@](C[C@]2(C)P3C4=CC=CC=C4)(C)O[C@]3(C)C1", "CP(C)C1=CC=CC=C1"]
Base_smiles = ["O=C([O-])C.[K+]", "O=C([O-])C(C)(C)C.[K+]", "O=C([O-])C.[Cs+]", "O=C([O-])C(C)(C)C.[Cs+]"]
Solvent_smiles = ["CCCCOC(C)=O,CC1=CC=C(C)C=C1", "CCCC#N", "CC(N(C)C)=O"]

それぞれMordredで記述子を計算した後、エラーを取り除きcsvで書き出します。auto-QChemのnotebook Fast featurize with Mordred.ipynbを参考にして進めていきます。

まず、必要なライブラリをインポート。

import pandas as pd
import numpy as np
import mordred
from rdkit import Chem
from mordred import Calculator, descriptors

数が少ないSolventを例に進めます。PandasのDataFrameに記述子を格納したいので、SMILESのlistをDataFrameにしておきます。 ついでにrdkitのrdmolオブジェクトに変換して列を加えておきます。

solvent_df = pd.DataFrame(Solvent_smiles, columns={'smiles'})
solvent_df['rdmol'] = solvent_df['smiles'].map(lambda x: Chem.MolFromSmiles(x))

print(solvent_df.shape)
# (4, 2)

ModredのCalculatorで記述子を計算します。

calc=Calculator(descriptors, ignore_3D=True)
md=calc.pandas(solvent_df['rdmol'])
print(md.shape)
# (4, 1613)

f:id:magattaca:20210223223324p:plain

1613個の記述子が計算されました。エラーが出ている部分をNaNに変換して、dropnaで列を削除します。

# エラーをNaNで置き換え
md=md.applymap(lambda x: np.nan if type(x) in [mordred.error.Missing, mordred.error.Error] else x)

# NaNを含む列の削除
md=md.dropna(axis=1)

print(md.shape)
# (4, 1275)

1275個残りました。あとはcsvとして保存するだけです。

SMILESのDataFrameと記述子のDataFrameをconcatenateで結合します。rdmolオブジェクトの列は不要なので削除しておきます。 また、他の要素と区別するため、列名にprefixsolvent_」をつけておきます。

# 保存用のDataFrameを作成。rdmolは削除しておく。
df_to_save = pd.concat([solvent_df.drop('rdmol', axis=1), md], axis=1)

# prefixをつける
df_to_save = df_to_save.add_prefix('solvent_')

# csvに保存
df_to_save.to_csv("solvent_desc.csv", index=False)

f:id:magattaca:20210223220457p:plain

無事「solvent_desc.csv」ファイルが作成できました!

同様にして「ligand_desc.csv」「base_desc.csv」を作成しました。これで記述子の準備はおしまい。

ベイズ最適化パッケージ(EDBO)

ではいよいよ本命、EDBO(Experimental Design via Basian Optimization)です。

これはベイズ最適化を化学合成に適応するための実用的なパッケージです。

インストール

まずはGitHubREADMEに従ってインストールします。

EDBO自体はpipでインストールできます。必要なライブラリはrdkitMordredPyTorchのようです。

指示に従ってcondaで新しい仮想環境(edbo)を作成し、諸々インストールしました。

# https://github.com/b-shields/edbo
# (0) Create anaconda environment
conda create --name edbo python=3.7.5

# (1) Install rdkit, Mordred, and PyTorch
conda activate edbo
conda install -c rdkit rdkit
conda install -c rdkit -c mordred-descriptor mordred
conda install -c pytorch pytorch=1.3.1

# (2) Install EDBO
pip install edbo

# Running Notebooks
conda install jupyterlab

おしまい。

計算の流れとモデルの中身

EDBOで反応を最適化する流れは以下の通りです。GitHubの例(Bayesian Optimization of a Mitsunobu Reaction)を参考にしています。

  1. 反応をエンコード、入力の記述子を取り込み:edbo.utils
  2. 反応空間(reaction space)のパラメーター設定: components``encodings``descriptor_matrices
  3. ベイズ最適化の設定(ebdoインスタンス化): edbo.bro.BO_express
  4. 計算ループの実施  
    1. 初期化(最初の実験セットを提案): init_sample()
    2. 実験結果の追加 : add_results()
    3. モデルの適合、獲得関数の最適化、次の実験の提案 : run()
    4. 提案実験の出力 : export_proposed()
    5. 提案実験の結果を追加し、収束するまで 2 -> 4を繰り返す

上記手順「3. ベイズ最適化の設定」でedbo.bro.BO_expressというオブジェクトを利用しています。 これは、パラメーターの自動的な読み込みや、全実験空間の組み合わせの計算、必要な前処理の実施といった作業を行ってれるもので、 edboをよりユーザーフレンドリーなプログラムにしてくれています。

BO_expressベイズ最適化のモデルもデフォルトで設定してくれています。デフォルトの設計は以下となっています。(class edbo.bro.BO_express)

ガウス過程はGPyTorchを使って実装されています。*7関連するやつを貼っておきます。*8*9

f:id:magattaca:20210223220719p:plain
カッコイイやつ

では順番に計算していきましょう!

EDBOでEVMLをプレイ

今回は「EDBOの指示に従ってEVMLをプレイする」という体で、計算を実行したいと思います。つまり、EDBOの計算で提案された実験をEVMLで実施し、結果をEDBOにフィードバックして次の実験を再度提案させ、収率の向上を目指します。

ではまず反応のエンコードです。 入力として使う分子の記述子を読み込みます。edbo.utilsを使って、Mordredの計算結果を含むcsvファイルを読み込みます。

計算実行

まず反応のエンコードです。 入力として使う分子の記述子を読み込みます。edbo.utilsを使って、Mordredの計算結果を含むcsvファイルを読み込みます。

# インポート
import pandas as pd
from edbo.utils import Data

# 記述子計算のファイルを読み込み
ligands_desc = Data(pd.read_csv('data/ligand_desc.csv'))
solvents_desc = Data(pd.read_csv('data/solvent_desc.csv'))
bases_desc = Data(pd.read_csv('data/base_desc.csv'))

読み込んだデータはdataで確認でき、PandasのDataFrameとなっています。

print(type(ligands_desc))
print(type(ligands_desc.data))
# <class 'edbo.utils.Data'>
# <class 'pandas.core.frame.DataFrame'>

次に、反応空間(reaction space)のパラメーター設定を行います。今回の合成反応条件で最適化するパラメータ(concentration、temperature)を含めて要素を指定する辞書を作成します。

# 反応空間のパラメーター
components = {'Ligand':'mordred', 'Solvent':'mordred', 'Base':'mordred',
              'Concentration' : [0.057, 0.100, 0.153],
              'Temperature' : [90, 105, 120]}

# Encodingの指定。指定しない場合、自動的にOne-Hot-Encoding (OHE)になる。
encoding = {'Concentration' : 'numeric', 
            'Temperature' : 'numeric'}

# 記述子の設定 (descriptor_matrices) 
mordred = {'Ligand': ligands_desc.data,
           'Solvent' : solvents_desc.data,
           'Base' : bases_desc.data}

descriptor_matricexの設定はEDBO外部で計算した記述子を使う場合に設定する項目で、Gaussianで計算したDFT記述子の場合などに使います。Mordredの計算はEDBOでも行えるので、実際にはここで指定しなくてもよかったのですが、今回は前もって計算したデータを使うことにしました。

次に、EDBOの BO_expressオブジェクトをインスタンス化します。上で指定した要素やモデルの設定など、もろもろを書き込みます。

from edbo.bro import BO_express

# BO object

bo = BO_express(reaction_components = components,
                encoding=encoding,
                descriptor_matrices=mordred,
                acquisition_function='EI',    # 獲得関数はEI
                init_method='rand',            # 初期化はランダム
                batch_size=5,                    # 1ラウンドあたりの実験数は5
                target='yield' )                   #  yieldを最適化

# 事前分布の指定。文献の設定に合わせる

from gpytorch.priors import GammaPrior

bo.lengthscale_prior = [GammaPrior(2.0, 0.2), 5.0]
bo.outputscale_prior = [GammaPrior(5.0, 0.5), 8.0]
bo.noise_prior = [GammaPrior(1.5, 0.5), 1.0]

下準備が完了したので計算ループに進みます。

まずは初期化!最初にする実験をランダムに選択します。

選択された実験はbo.export_proposed()csvファイルとして書き出すことができ、またnotebook上でもbo.get_experiments()で確認できます。

bo.init_sample(seed=0)             # 初期化
bo.export_proposed('data/init.csv')     # CSVファイルに書き出し
bo.get_experiments()               # ノートブック上で確認

f:id:magattaca:20210223220947p:plain
ランダムに選ばれた初期反応条件

5つ実験が選ばれました!この条件で実験を行い、結果を追加して次のループに進むことになります。

実際の環境では、実験に時間がかかるのでPCを付けっ放しにしておくのはもったいないです。そのような場合に設定を保存しておけるよう、EDBOではpickel化するメソッドedbo.bro.BO.save()も用意されています。

bo.save()  # optimizerの保存

作業ディレクトリに「BO.pkl」が作成されました。再度読み込みたいときはedbo.bro.BO.load()を使えばOKです。

from edbo.bro import BO_express

bo = BO_express()
bo.load()

さて、話戻して、出力されたcsvファイルを見ると、以下のようにyield列があり「<Enter Response>」となっています。ここに収率を書き込めば、次のループで使用する入力ファイルとできます。

f:id:magattaca:20210223221039p:plain
出力csvにyieldを書き込む

今回取り上げている反応の結果は、EBDOのGitHub 「experiments/data/direct_arylation/experiment_index.csv」に記載されています。もしくはEVMLで確認することもできます。

今回は「EDBOの指示に従ってEVMLをプレイする」ということにしているので、EVMLに条件を入力し、得られた結果をcsvに追記し、次のラウンドに回しました。

結果を記入したファイルはedbo.bro.BO_express.add_results()メソッドで読み込みます。ついでrun()を実行すれば、ガウス過程モデルのfitと獲得関数(EI)の最適化を行い、次にすべき実験を提案してくれます。

# 結果をboオブジェクトに追加
bo.add_results('data/init_results.csv')

# 計算の実施
bo.run()

# 次にする実験の取り出し
bo.export_proposed('data/round0.csv')  
bo.get_experiments()  

新しい提案の実験条件が手に入りました(round0.csv)。これをEVMLに入れて結果を得て、再度EDBOで計算、次の条件を提案を繰り返します。

・・・結果・・・・

round6の提案実験で、収率100%の反応条件が2つ見つかりました!

最初のinitバッチを含めて計7バッチ35実験)です。だめ押しでround7も実験しましたが、収率は下がったのでここで止めました。

EVML上で確認できる「条件探索の過程」は以下の図のようになりました。

右から左に向かって新しい実験です。赤色EBDOの提案した実験の収率で、緑色ランダムに選んだ場合の実験収率です。

最初は10%にも満たなかった収率が徐々に上がってゆき、最終的に収率100%を達成している様が見て取れます。

f:id:magattaca:20210223221512p:plain

結果の確認

EDBOは結果解析用のツールも提供しています。

例えば、最適化ループを回す過程でoptimizerconvergenceがどのように変化したかプロットすることができます。

bo.plot_convergence()

f:id:magattaca:20210223221542p:plain

順調に収束している様子がわかりますね!

GitHubのnotebookで使われている解析も適用してみましょう!

実施した全40実験(8バッチ)の中で、最も収率の良かった上位5実験について反応条件を提示してみます。

全ての結果を読み込んで統合したDataFrameを作成します。カラムは反応条件(Ligand, Solvent, Base, Concentration, Temperature)と収率(yield)で、yieldで降順に並べ替えます。

results = pd.DataFrame(columns=bo.reaction.index_headers + ['yield'])
for path in ['init', 'round0', 'round1', 'round2', 'round3', 'round4', 'round5', 'round6', 'round7']:
    results = pd.concat([results, pd.read_csv('data/' + path + '_results.csv', index_col=0)], sort=False)

results = results.sort_values('yield', ascending=False)
results.head()

f:id:magattaca:20210223221834p:plain

できました!

上位5つのリガンド、塩基はどれも同じですね。温度はより高温、濃度は濃い方が良さそうですが、溶媒の影響もありそうです。

DMAc溶媒では薄め(0.057)でも収率96.6%なのが興味深いですね。

棒グラフも書いて見ます。

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, len(results.columns.values[:-1]), figsize=(30,5))

for i, feature  in enumerate(results.columns.values[:-1]):
    results[feature].iloc[:5].value_counts().plot(kind="bar", ax=ax[i]).set_title(feature)
plt.show()

f:id:magattaca:20210223221902p:plain

また、EDBOは構造式描画のためのChemDrawというモジュールも用意されています。

上位5反応に出てきた分子の構造式を描画しましょう。

from edbo.chem_utils import ChemDraw

for col in results.iloc[:,:3].columns.values:
    print('\nComponent:', col, '\n')
    cdx = ChemDraw(results[col].iloc[:5].drop_duplicates())
    cdx.show()

Component: ligand_smiles_index

f:id:magattaca:20210223221925p:plain

Component: solvent_smiles_index

f:id:magattaca:20210223221954p:plain

Component: base_smiles_index

f:id:magattaca:20210223222015p:plain

なるほど便利!!

やっぱり構造式でみるのがが一番しっくりくる!!

感想

以上、今回は前回の記事で取り上げたNatureの文献のツールで遊んで見ました。データとツールが公開されているというのはとてもありがたいですね!

特にベイズ最適化ツールEDBOは変更が必要な設定項目も少なく、コピペするだけで簡単に使えてとても便利な印象でした。

EDBOの実践として、記述子にMordredを使い、EVMLをプレイして見ました。最初は低い収率で類似の実験を検討しており、「大丈夫かなー、本当にこれでいいのかなー」と半信半疑でしたが、収率上昇の兆候が見え始めると一気に条件検討が進んでいき最終的に収率100%の条件を見つけることができました!

私は自分でEVMLをプレイした際、99%までしか行かなかったので、100%の結果を見てびっくりしました。・・・完敗です。

ところで、文献中では最終的に記述子として量子化学計算によるものが用いられていました。今回はMordredでも良い結果に到達できることが確認できたので良かったです。Gaussianなんて一生触れそうそうにないですからね!

溶媒の沸点を考慮して反応温度に拘束をかけるとか、追加の設定を加えたらより探索効率上がるかもなー、とか思いましたがその辺りどうなんでしょう???

いずれにせよ、面白いツールがどんどん出てきて公開されるのは見ていて楽しいですね!

ほとんどGitHubのコードそのままコピペした記事ですみませんでした。。。ではでは!

*1:READMEをそれっぽっくかえてみました

*2:参考「aws_setup.TXT」の「5. Install required packages」

*3:私はRStudioで実行しました。Rはほとんど触ったことがなかったのでRStudioの使い方に手間取りました・・・

*4:ランディングページ(LP)というのは、ユーザーがWebで最初にアクセスするページのことらしいです。知らなかった。

*5: 仮想環境でjupyter notebookの立ち上げにエラー(Bad config ~~)が出るときは、「pip install environment-kernels」をすれば動くようになるかもしれません。参考

*6:コンフォマーの数以上の数字を与えると、最後のコンフォマーが選ばれるそうです

*7:GPyTorchについてはHELLO CYBERNETICSさんの記事「GPyTorchでガウス過程を見てみよう」がとても参考になります。

*8:参考 講談社ガウス過程と機械学習

*9:参考ガウス過程回帰の基礎 赤穂昭太郎 システム/制御/情報 2018(62)390

Macで独習!量子化学計算入門!! 〜インストールでつまずいた話〜

電子書籍「独習 量子化学計算」を購入し、一通りポチポチ遊んでみました。

ゼロからわかる!!」の謳い文句の通り、ど素人でもソフトウェアを動かすことができ、とてもオススメです!

WindowsでもMacでも実践できるように解説してくださっているのですが、 私のMacではソフトウェアのインストールなど、環境構築でいくつかつまずきました。

備忘録としてメモを残します。

尚、本記事は以下の環境での事例です。

書籍の紹介とオススメポイント

まず書籍について、初学者からみた感想とオススメポイントをご紹介します。

ゼロからわかる!!独習 量子化学計算」は、めちゃくちゃわかりやすくて素晴らしいサイト「PC CHEM BSICS.COM」さんが出されている電子書籍です。

こちらのサイト、量子化学の用語を検索すると何度もトップにヒットしました。いずれの記事も簡潔でわかりやすく、初学者の私にとって非常に参考になりました。

また、理論・用語の解説だけでなく、「無料で楽しむ量子化学計算サイト」とのことで、 自宅のPCで試すことができるソフトウェアとその利用方法についてもたくさん紹介してくださっています。

気軽に試せるのは初心者にとってとてもありがたいですよね!

で、この素晴らしいサイトの電子書籍独習 量子化学計算」ですが、 サブタイトルは「理論からはじめない新しい量子化学計算の本」!!!

「難しいことよくわかんないけど、とりあえずやってみたい!」という、私にとってピッタリすぎる謳い文句です!

早速購入して一通りやってみましたが、「初学者が独学でつまずくポイントがよくわかっていらしゃる・・・」という感じでとても分かりやすかったです。

具体的には、以下のような構成となっているのが良かったです。

  1. 実際に計算を経験することを通して学習できるように構成
  2. ソフトウェアのインストールから使い方結果の解析まで解説
  3. 一つ一つの章が簡潔で、初学者でも力尽きずに取り掛かかれる
  4. ステップバイステップでより複雑な計算の設定がわかるようになる

OSはWindowsMac両方に対応しており、一方でしか使えないソフトについては、その旨と別の方法の解説がのっています。(使われているソフトはちょっとWindowsより?)

また、操作方法は実際の画面の図で解説されているのでわかりやすくなっています。

個人的にとても助かったのは、操作方法を繰り返し説明してくださっている点です。

私は記憶力が貧弱なので、書籍をよんでいるとよく「基本忘れた!また、前の章に戻って確認しなきゃ。面倒くさいからもういいや。」と投げ出してしまうのですが、今回はそうならずにすみました!*1

ちょっと気になった点は、操作画面の写真の解像度が少し低く、 拡大してもボタンや設定の文字がよくみえないということがありました。*2

あと、「文字検索できたら便利なのになー」と思いました。*3

また、ソフトウェアのインストールや様々なサイトの紹介に際してリンクが記載されていますが、 Kindleから直接Webページに飛んだり、文字列のコピーをしたり、ということができませんでした。ただし、この点は「書籍のサポートページ」にリンク集を用意してくださっているので、そちらで解決しました。

さらにリンクだけでなく、サポートページには「サンプルファイル」も用意してくださっています。「自分で作った初期構造だと計算が収束しない!振動数めっちゃ負!もう嫌だ!!」となった時には、サンプルファイルを利用して計算を進めていくことができます。

「できの悪い人間の諦めポイント」がよくわかっていらっしゃる!ありがとうございます!!

というわけで、「量子化学計算興味はあるけど、どこからスタートすれば良いかわからない」という方は この書籍から始めるのがとてもオススメだと思います!

Macで環境構築する際につまずいた点

では本題、Macで環境構築する際につまずいた点とその回避策の備忘録です。

Avogadroの英語化

まず、Avogadroについてです。 Avogadroはさまざまな量子化学計算ソフトウェアの入力ファイルを作成することができ、 本書でも繰り返し利用方法が解説されています。

MacでAvogadroをダウンロードすると、メニューが日本語化されているのですが、 ちょっと微妙で、英語・日本語の重複があったりして混乱します。

こんな感じ・・・

f:id:magattaca:20210204232001p:plain
Avogadroの日本語がちょっと微妙

書籍の表記と合わせるためにも英語化した状態で利用するのがオススメです。

Avogadroのメニューを英語にする方法は、PC CHEM BASICS.COMさんの記事 「Avogadroを使ってみよう」 でも解説されていますが、私の環境ではうまくいきませんでした。

サイト記載の方法は以下の通り。

  1. Avogadroのアプリケーションをcontrolボタンを押した状態でクリック
  2. 「パッケージ内容を表示」を選択
  3. Contents⁩>Resources⁩内の「ja.lproj」をフォルダごと削除

この通りにしても日本語のままでした。

そこで、Avogadroのメールフォーラム(?)(こちら)の記述を参考にして、さらに「xxxxx.qm ファイル」を削除しました。

f:id:magattaca:20210204232510p:plain

こちらのファイルもAvogadroを右クリックして「パッケージの中身を表示」すれば確認できます。

このqmファイルが翻訳に関与しているそうなので、削除すれば自動的にデフォルトの英語となります。

こんな感じ・・・

f:id:magattaca:20210204232138p:plain
英語にしてしまった方がわかりやすい

メニューの重複も無くなってわかりやすくなりました。

Fireflyのテスト失敗(Wineの設定、prefix)

次にFirefly(PC GAMESS)(version 8.2.0)の計算を実行しようとしてつまずいた点です。*4

Fireflyインストール後、以下のコマンドでテストを実行しましたが失敗しました。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT-DualCore-Only.sh

ターミナルで実行すると、以下のようなメッセージがでて計算が走りませんでした。

wine: 'ホームディレクトリ/.wine' is a 64-bit prefix, it cannot be used with 32-bit Wine.

「テストファイルが問題?」と、別に用意した「水(H2O)の構造最適化計算」を試しましたが、 こちらもアウトプットがファイルが作成されず、Fireflyがすぐに閉じてしまいました。

エラーメッセージのWINEですが、「Windows用のソフトウェアをMacで利用するための環境をつくる」ためのものだそうです。

Mac用のFireflyをインストールすると、内部に「WINE」というディレクトリもあり、Fireflyパッケージの中に含まれているようです。

で、ソフトウェア起動時に、WINEの設定がホームディレクトリ下の「.wine」ディレクトリにつくられて読み込まれるようなのですが、 どうもこのprefix(?)のバージョンの違いが問題となったようです。

そこでこの「.wine ディレクトリを一度全て削除」してしまったあとに、再度Fireflyのテストを試したところ無事実行することができました!

Firefly実行後に再度確認すると、同じ場所に新しい「.wine」ディレクトリができていました。

注) この方法だとWINEに依存する既存の他のパッケージに影響があるかもしれません。もっと良い方法をご存知の方はご教示いただければ幸いです。

私の場合、以前に別のWindows用ソフト(WinMOPAC)を利用しようとして、Wineを個別にインストールしたことがあり、 その時に作られた設定が今回のFireflyと合わなかったのが原因だったようです。FireflyのWINEはwine 1.1.33 ソースに基づいているようですが、10年くらい前のもののようなので、以前のWINEが64-bitで、Fireflyが32-bitだったのかなぁと推測しています。

尚、Firefly付随のインストールマニュアルに書かれているテスト用のコマンドは、現在のバージョンとずれているので修正が必要です。

マニュアルでは、Dual coreの場合ターミナルで以下を実行するように記載されています。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT.sh

実際には、配布されているファイルは以下でした。

/Applications/Firefly/BENCHMARKS-THREADED-TESTJOBS-RUNSCRIPT-DualCore-Only.sh

テストされる際にはご注意を!

Fireflyでインストールが推奨されているアプリ

Firefly 8.0.0(for Mac)で推奨されているソフトウェアとして以下のものがあります。

  1. TextWrangler
  2. wxMacMolPlt

Firefly 7.1.G update 以降、これらのソフトウェアへの依存性(dependency)はなくなっているようなので、 別に入れなくても問題はないようです。

ですがFireflyのインストールマニュアルに「便利だし無料だからマジでおすすめ!(“strongly recommended”)」と書かれていたので導入しました。 デフォルトの「/Applications」ディレクトリに入れておくことがオススメされています。

TextWranglerは無料のテキストエディタです。 Firefly計算実行した際に出力の「.outファイル」がTextWranglerで自動的に開かれ、 逐次内容が更新されるので計算が走っているかどうかの確認に便利です。・・・が、邪魔といえば邪魔です。。。

無料のテキストエディタは色々あるので、別にTextWranglerに拘らなくても良いのかなーという感じです。 オススメな理由をご存知な方がいらっしゃったら教えていただければ幸いです。

もう一つのwxMacMolPltは、計算結果の解析に使うことができ、「独習 量子化学計算」でも何度もでてきます。こっちはインストールしておきましょう!

GamessQ立ち上げ時のエラー

2021年現在、「GAMESS for MacOS X」にはGamessQが含まれていません。別途インストールする必要があります。

私のMacOS環境では「GamessQ (v1.2.2)」を使おうとすると以下のエラーが出ました。バックエンドの立ち上げに失敗しているみたいです。

f:id:magattaca:20210204232717p:plain
GamessQのエラー

「Add」などのボタンがグレーになっていて使うことができません。・・・困った。

回避策(?)はこちらのGitHubページに書かれていました。(→Error launching backend on fresh install under OS X #19

  1. GamessQのディレクトリにあるコマンドラインバイナリ「gamessqd」をダブルクリックして実行
  2. アプリ「GamessQ.app」を立ち上げ

黒い四角のアイコン gamessqdをクリックするとターミナルが立ち上がります。 この状態でGamessQのアプリを立ち上げれば先のエラーが出ず、正常(?)に起動できました!

f:id:magattaca:20210204235848p:plain
「GamessQ.app」の前に「gamessqdd」を立ち上げる

「Add」や「Clean up」といったボタンが使えるようになっています(アクティブ?)。

f:id:magattaca:20210204232915p:plain

アプリを閉じるとターミナルにもlogoutと表示されてプロセスが終了します。

Gamess計算実行時のエラー(HOST_NOT_FOUND)

GamessQが使えるようにできたので、計算を実行したところ結果が得られずログに以下のエラーが出ました。

 Error: Gethostbyname(XYZXYZ) returned HOST_NOT_FOUND.
 ddikick.x: Fatal error detected.

Gethostbyname(XYZXYZ)のXYZXYZの箇所には、実際には私のMacBookに設定した名前が入っています。*5

上記のエラーはGamessQを使わずに、Gamess単独で以下のテストを実行しても再現しました。(GAMESS(US)は少し古い2018年のビルドです)*6

  1. ターミナルでGamess実行ファイルを格納したディレクトリに移動
  2. ./gms tests/standard/exam01
  3. output file name? [tests/standard/exam01.log]と聞かれるのでReturnキーを押す
  4. 出力されたexam01.logを確認

Gethostbyname()は指定したホスト名から、ホストに関する情報を得るための命令のようです。
つまり「PC(MacBook)に設定した名前のホストが見つからないよ!!(HOST_NOT_FOUND)」みたいなエラー???

色々調べた結果、macOShostsファイルの設定を追加してやればよいのではないか、となりました。

hostsは「./etc」(あるいは「./private/etc」)にあります。テキストエディタで開くとこんな感じでした。

##
# Host Database
#
# localhost is used to configure the loopback interface
# when the system is booting.  Do not change this entry.
##
127.0.0.1 localhost
255.255.255.255 broadcasthost
::1             localhost

「よくわからないけどローカルだからローカルホストでしょ!」ということで、 7行目の後ろにPCの名前(XYZXYZ)を追加して以下のように書き換えました。(書き換えには権限が必要です)

127.0.0.1 localhost XYZXYZ(←MacBookにつけた名前)

hostsファイル変更後に、Gamessのテストを再度実行すると以下のようなポップアップが出てきました。

f:id:magattaca:20210204233240p:plain

ddikick? OK!OK!・・・知らんけど。

「ネットワーク受信接続」を許可すると「exam01.dat」ファイルが作成され、無事計算実行できました!

「exam01.log」の中身も 「EXECUTION OF GAMESS TERMINATED NORMALLY」となっていました。

こちらのページによると、 「127.0.0.1ローカル・ループバック・アドレスと呼ばれ、自分自身を指す特別なIPアドレスである。「localhost」という名前でも参照できる。自分自身の上で動作しているサービスへ接続する場合は、このIPアドレスを利用できる。」だそうです。

GamessはローカルPCで計算を実行する場合でも、内部でネットワーク接続処理(TCP/IP?)を行う構成になっているのでしょうか???

localhostに追加するのが正しいのか全くわかりませんがとりあえず動いたのでヨシ!!

IPアドレスをもう少し検証

せっかくなのでもう少し実験してみました。

先ほどのページには以下の記載がありました。*7

一般的には「127.0.0.1」というIPアドレスIPv4の場合)が利用されるが、実際にはIPアドレスの最上位のバイト(最上位の8bit)の内容が「127」でありさえすればよいので、「127.0.0.1~127.255.255.254」の範囲内ならばどのIPアドレスでも利用できる(127.0.0.0と127.255.255.255の2つはブロードキャスト・アドレスのため除外される)。例えば127.0.0.2でもよいし、127.1.2.3でもよいが、一般的な用途では127.0.0.1だけが利用される。

別のアドレスでもOKなら、localhostの後ろに加えずに「127.0.0.2 XYZXYZ」としてもよいのでは???

ということでhostsファイルを以下のようにして再テストしてみました。

127.0.0.1 localhostt
127.0.0.2 XYZXYZ(←MacBookにつけた名前)

127.0.0.2を使っても計算は終了し、「exam01.log」は EXECUTION OF GAMESS TERMINATED NORMALLYとなりました。

ですが、データファイル「exam01.dat」(PUNCH出力)が作成されず、代わりに「exam01.F05」というファイルができました。 また、ターミナルもコマンド実行中のままで次に移行しませんでした。

logには計算結果が書かれていたものの、「.datファイル」(PUNCHファイル)への出力、計算プロセスの終了がうまくいっていない感じでしょうか?

正常な動作のためにはlocalhostと同じIPアドレスにしておいた方が良さそうです。

ソフトウェアの通信やコマンドラインの処理、何が行われているのかさっぱり意味がわからないです。。。。

wxMacMolPltによる結果の可視化の違い(FireflyとGamessQ)

「独習 量子化学計算 Chapter 8 フロンティア軌道で反応を予測しよう」では、MacMolPltにより分子軌道を可視化する方法が解説されています。

この操作で、書籍では「Fireflyで得た計算結果(.outファイル)」を解析していましたが、私の環境では再現できませんでした。

一方で、「GamessQで同じinputファイルを使って得た計算結果(.logファイル)」では、可視化操作がうまくいきました。

行いたかった作業は以下のように「Select Orbital Set : Molecular Eigenvectors」で、 「Select Orb : Energy」から目的の分子軌道を選択し描画することです。

f:id:magattaca:20210204233550p:plain
やりたかった作業

上図は、GamessQで得たlogファイルを用いたものです。対して、Fireflyのoutファイルでは以下となりました。

f:id:magattaca:20210204233625p:plain
Fireflyの出力ファイルだとできない

上図のように、そもそも選択できる描画「Surfce Type」が少なく、「Atomic Orbitals」しか選択できませんでした。

Fireflyの計算結果(.outファイル)の中身を確認すると、以下のように「MOLECULAR ORBOTALS」の情報が記載されていたので、 おそらく描画ソフト側(wxMacMolPlt)の問題ではないかと推測しています。

f:id:magattaca:20210206210437p:plain
Fireflyでも分子軌道は出力されている

私のPC環境(or MacOS)特有の問題かもしれませんが、同じ問題に出くわした方のためにご参考までに紹介しておきます。

おまけ Gatekeeperへのアプリ許可の追加

GamessQのbackend エラー回避で参照したGitHubGateKeeperパーミッションについても言及されていました。 「administratorの権限がないと、先の回避策が使えないのでは???」って感じの話みたいです。

こちらのページによると、「Gatekeeperとは、Mac App Storeからインストールするかユーザーの右クリック操作により許可されたアプリだけを起動できるようにする、セキュリティ強化の仕組み」だそうです。

「開発元本当に信頼できる???」って聞いてくるあれですね!

この許可をターミナルのコマンド「spctl -add パス」を使って処理することもできるとのことです。

例えば、GamessQの場合なら以下のようなコマンドとなります。(パスは適宜修正してください)

spctl --add "/Applications/GamessQ/GamessQ.app"
spctl --add "/Applications/GamessQ/gamessqd"

許可リストから外す場合は「spctl -remove パス」とすれば良いそうです。

spctl --remove "/Applications/GamessQ/GamessQ.app"
spctl --remove "/Applications/GamessQ/gamessqd"

インストールの際にsudoと組み合わせておけば、administrator権限がない場合でも問題を回避できるようになるみたいです。(GitHubでは実習で生徒が使うPCの設定が議論されてました)

まとめ

以上、「独習 量子化学計算」の感想と、Mac量子化学計算の環境構築につまずいた点のご紹介でした!

繰り返しになりますが、とりあえず計算を始めてみたい初学者にとてもオススメな書籍という感想です。*8

ハートリー・フォックやDFTという用語をなんとなく知ったものの、実際にはどんな感じか全くわからなかったので非常に勉強になりました。

動かしてみると、「SCF計算ってわかったつもりになってたけど、具体的にどの部分の何を計算してるんだっけ???」とか、 「DFTで計算したら、HFでできたアセトンの構造最適化失敗するやん!負の振動とは一体???」などなど、新たな疑問がたくさん出てきました。

すこし初期構造を変えるだけで計算がうまくいかなくなったりして、計算化学の実験科学的な側面をちょっと体感できたのでよかったです。たぶん論文に載っているような結果は、膨大な試行錯誤の結果なのでしょうね。格好いい図くらいに思っててすみませんでした。

あと、途中で北浦-諸熊解析という用語が出てきて、おおっよくわからないけど格好いい!!!これが噂のONIOM??FMO??とテンション上がりましたが全然関係なかったです。

水分子1つの計算でつまずいているのに、タンパク質の計算を理解できる日はくるのでしょうか。。。。

ゆれるアセトン

f:id:magattaca:20210204234525g:plain

*1:レベル低い感想ですみません

*2:外付けモニタと私の視力のせいかも・・・

*3:私のKindleアプリが古いせいかも

*4:Fireflyのライセンス認証のメールが中々来ないと思っていたら、Gmailの迷惑フォルダに入っていました。受信日、翌日だったのでレンスポンスめちゃくちゃ早くてびっくりです。

*5:恥ずかしいので隠した

*6:テスト方法参考 量子化学計算ソフト GAMESS のエラー

*7:ローカル・ループバック・アドレス(127.0.0.1)とは?

*8:回し者ではないですよ!ノーアフィリエイトですよ!

PyMOLをJupyter Notebookと連携させる方法について

前回の記事で、PyMOLのスクリプト機能を利用して3D PSAを求めました。『PyMOL=生物学向けの描画ソフト』と思い込んでいたのですが、まだまだ便利な機能がたくさんあって、化学方面にも応用できそうです。

素敵なPyMOLですが、プログラミングができない私にとってコマンドラインでの操作はとっつきづらいです。コマンドと一緒にコメントや出力結果を残しやすいJupyter notebookで使えた方が後でたどりやすいのになー、と、いうことでJupyterから使う方法を調べました。

PyMOL WikiJupyterという、そのものな項目があったので記載の方法を試していきたいと思います。ほぼそのままです。

検証する操作の確認

Jupyter notebookを試す前に、検証で行う操作を念のためPyMolアプリのGUIで確認します。

  1. ラニンを読み込む(fragmentcmd.fragment)
  2. 描画を画像(PNG)で保存(pngcmd.png)

アミノ酸残基はBuildメニューのResidueから選択して作成することができます。コマンドラインfragment alaとするか、cmd.fragment('ala')としてもOKです。

f:id:magattaca:20210109154316p:plain

画像の保存はFileメニューのExport Image Asあるいは右上のDraw/Rayパネルからできます。コマンドラインで行うには、png testとすると作業ディレクトリの中に「test.png」として保存できます。cmd.png('test')でもOKです。

f:id:magattaca:20210109154344p:plain

詳しい使い方はPyMol Bookが日本語でわかりやすいです。なお、cmdPyMOL Python APIを使うためのモジュールで、PyMol上のコマンドラインを使うのであればわざわざ使わなくても良いと思いますが、以降の内容との関係するので併記しました。

JupyterからPyMOL

では本題、Jupyter notebookからPyMOLを利用する方法です。PyMol Wikiには大きく2つ、直接PyMOL APIを使う方法と、サードパーティーのライブラリを使う方法が書かれています。それぞれ制限があるようです。

概要はこんな感じです。

直接PyMOL APIを使う

同じスレッド

まず、同じスレッドで使う方法です。こちらはGUIの画面をみることはできません(headless)。

操作は簡単。Jupyter notebookを起動し、PyMOLのモジュールを読み込むだけです。

# cmd モジュールのインポート
from pymol import cmd 

# アラニンを構築し、拡大した後に画像を保存
cmd.fragment('ala')
cmd.orient()
cmd.png('./test.png', ray=1)

PyMOL not running, entering library mode (experimental)」とかいう出力がでてきますが問題なく使えます。PyMOL WikiLaunching From a Scriptによると、PyMol 2.1から、cmdを使うと自動的にメインスレッドでバックエンドプロセスをGUI無しで起動するらしいので、その関係だと思います。*1

さて、上記で無事「test.png」が作成されたので中身を確認してみます。

from IPython.display import Image
Image(filename='./test.png')

f:id:magattaca:20210109154444p:plain

うまくいきました!

非同期スレッド

次に、非同期スレッド(asynchronous thread)で使う方法です。こちらはGUIが使えます。

が、残念ながらmacOSではこの操作を行えません(GitHub Issue)。

コマンドだけコピペしておきます。

# PyMOLウィンドウを開く
import sys
import pymol
_stdouterr = sys.stdout, sys.stderr
pymol.finish_launching(['pymol', '-q'])
sys.stdout, sys.stderr = _stdouterr

# PyMOLウィンドウにアラニンを読み込む
from pymol import cmd
cmd.fragment('ala')

Macでこれを実行するとカーネルが死にます。。。

コマンドラインモード

PyMOL Wikiから外れますが、コマンドラインモードであればMacでもJupyterからPyMOLを立ち上げることができます。起動のオプション-cを使うと、GUIなしのコマンドラインのみのモードとなり、バッチプロセスに使うことができます(PyMOLWiki Command Line Options)。

import pymol

# PyMolをコマンドラインモードで立ち上げる
pymol.pymol_argv = ['pymol','-c']
pymol.finish_launching()
# Detected 4 CPU cores.  Enabled multithreaded rendering.

ターミナルにPyMOL起動時のバージョン情報が表示され、起動がうまくいったことがわかります。オプション-qをつければこの表示が出ないようにすることもできます(quietモード)。

注意点として、このPyMOLの立ち上げ方はCPUを100%使うらしいです*2Macが悲鳴をあげるので不要であればcmd.quit()で終了しておいた方が良さそうです。

# コマンドラインモードのPyMOLを終了
cmd.quit()

RPCコミュニケーション

つづいて、Jupyter notebookとは別に立ち上げたPyMol(スタンドアローン)をRPCコミュニケーションで操作する方法です。RPC(Remote procedure call、遠隔手続き呼び出し)は「プログラムから別のアドレス空間にあるサブルーチンや手続きを実行することを可能にする技術」で、「遠隔相互作用の詳細を明示的にコーディングする必要がない」そうです(Wikipedia 遠隔手続き呼び出し)。・・・よくわかりません。

まず、pymol -Rで立ち上げておきます。-RオプションはRPCサーバーを立ち上げるためのオプションです。

pymol -R
# xml-rpc server running on host localhost, port 9123

PyMOL GUIアプリが立ち上がりますが、上記のように表示され、xml-rpc serverlocalhostのポート番号9123で起動していることがわかります。

このRPCサーバーにJupyter notebookからアクセスするにはPythonxmlrpc.clientモジュールを使えば良いそうです。

# Jupyter notebookからPyMOL にRPCコミュニケーション
import xmlrpc.client as xmlrpclib
cmd = xmlrpclib.ServerProxy('http://localhost:9123')
cmd.fragment('ala')

以上を実行すると、PyMOLのGUIウィンドウにアラニンが読み込まれます。

xmlrpc.client.ServerProxy()でRPCコミュニケーションを管理するためのオブジェクトを生成し、これにPyMOL Python APIを扱うためのモジュールcmdと同じ名前をつけておくことで、APIを扱う感覚で外部のPyMolを操作ができているようです。・・たぶん。

以上が、直接PyMol APIを使う方法です。

サードパーティーライブラリ

つづいてサードパーティーのライブラリを使った方法についてです。ライブラリとしてIPyMOLRDKitが言及されていますがPyMOL Wikiには詳細は書かれていません。

ライブラリの説明やGitHubをたよりに試してみましょう。

IPyMOL

まずIPyMOLです。IPythonからPyMolのセッションを制御するためのライブラリで、こちらもRPCコミュニケーションを利用します(MITライセンス)。

「Jupyter Notebookできちんと見せたい時や、PyMOLスクリプトのプロトタイプを簡便に作成したい時に理想的なツール」とのことなので、今回の目的にまさにぴったりなライブラリです。

pipでインストール可能(ipymol - PyPI)ですが、古いバージョンらしく不具合があるらしいです(GitHub issue #21)。GitHubバージョンは修正されているらしいので、直接インストールしてみます。

pip install git+https://github.com/cxhernandez/ipymol.git

IPyMOLでは、さきと同様に、あらかじめpymol -Rで立ち上げておいたPyMOLに接続できますが、IPyMolからstart()PyMOL RPCサーバーを起動することもできます。閉じる時はstop()です。

from ipymol import viewer as pymol
pymol.start()

PyMOLのGUIが立ち上がりました!今回は「xml-rpc server running on host localhost, port 9124」と出たので、ポート番号が先とは異なるようです。

ラニンを読み込むには以下のようにします。

pymol.fragment('ala')

display()を使うことでPyMOL GUIのディスプレイをJupyter Notebook上に描画することができます。

pymol.display()

f:id:magattaca:20210109154735p:plain

IPyMolで利用可能なAPIは以下を打ち込むと確認できます。出力は多いので省略します。

print(pymol._server.system.listMethods())

上記のリストに含まれていないコマンドを使いたい時はpymol.do("実行したいコマンド")とすればOKです。

使い方の実例はこちらのノートブックが参考になります。

終わったらstop()でウィンドウを閉じられます。

# IPyMOLで立ち上げたPyMOLの終了
pymol.start()

RDKit PyMolモジュール

続いてRDKitからPyMolに RPCコミュニケーションする方法です。Chem.PyMolモジュールを利用すればOKで、化学の新しいカタチさんが「RDKitからPyMOLを利用する」で詳しい使い方をご紹介してくださっています。

参考までにこれまでと同じ処理を確認します。

まず、pymol -RXML-RPC サーバーを立ち上げておきます。

RDKitで必要なモジュールをインポートし、MolViewerオブジェクトを作成します。

from rdkit import Chem
from rdkit.Chem import PyMol

v = PyMol.MolViewer()

MolViewerに含まれていない操作を行うには、sever.do("実行したいコマンド")とすれば良いそうです。アラニンフラグメントを読み込み、GetPNG()でJupyter notebook上に描画します。

# アラニンの読み込み
v.server.do("fragment ala")

# ズーム
v.Zoom('ala')

# PyMOL GUIの画面をPNGで取得
v.GetPNG()

f:id:magattaca:20210109154938p:plain

SaveFile()を使うことで画像として保存することもできます。

折角ですので、ついでに前回の記事で取り上げた3D PSAのスクリプトを実行してみましょう!

まずDeleteAll()で一度全て消した後に、LoadFile()でファイルを読み込みます。「conformer-00.xyz」というxzyファイルを「conf1」という名前のオブジェクトとしてロードしています。

# viewewrのオブジェクトを全て消す
v.DeleteAll()

# ファイルの読み込み
v.LoadFile('./conformer-00.xyz', 'conf1')

別途作成したスクリプトpsa3d.py)をつかえるように読み込み(run)、計算を実行します。

# スクリプトの実行
v.server.do('run ./psa3d.py')

# 計算の実行
v.server.do('set solvent_radius=1.4; set dot_density=4; set dot_solvent=1; psa3d')

無事計算が実行され、PyMOL上にSA 3D PSAの値「201.973」が表示されました。

以上、サードパーティライブラリの動作も確認できました。

まとめ

以上、PyMOLをJupyter notebookで利用するための方法の紹介でした。ほとんどPyMOL Wikiそのままですが。。。

PyMOLのスクリプトは色々な方が作成、公開してくださっているようなので、Jupyter notebookとうまく連携させればプログラミング苦手な私でももっと活用できるかもしれません。

Jupyter notebookで操作して出力された結果を、PyMOLから受け取ってJupyter Notebookに表示することもしたかったのですが、やり方がよくわかりませんでした。サブプロセス(?)を使えば良いのかもしれませんが私のキャパシティでは無理でした。。。

ではでは。

*1:どういうことかはさっぱりわかりません。。。

*2:Launching From a Script