RDKitとcolabolatoryで前処理しようとして挫折した話
RDKitの使い方を調べたりインフォマティクスに関する本を読んでみたりしましたが、結局あまり実際には動かせていません。
おそらく、写経（？）したりして身につける必要があるのでしょうが、いまいちそれだと楽しくない（ダメな入門者の典型）。
何か面白いこと・・・ということで、調子に乗って創薬ちゃんさんの創薬レイドバトルに登録してしまいました。
どうしよう・・・
とりあえずどんな化合物から絞り込めば良いのか化合物ライブラリを眺めてみることにしました。
化合物ライブラリをダウンロード、解凍してみると以下の４つのsdfファイルが入っていました。
試しにsdfを開いてみたらMarvinViewが起動し手軽に見ることができました。
いつのまにこんなソフトが？と思ったら、ChemDrawの代わりに使えるものをということでインストールしたMarvinSketchと一緒にダウンロードされていたみたい・・・ChemAxon社様様です。
中身は上記のような感じ・・・。構造に加えてidnumber、heavy_atoms、LogSなどなどのカラムがあるようです。
思ったよりシンプルな化合物です。もっとドラッグライクな化合物はあるのでしょうが、なにぶんデータが多すぎてよくわからない。 これは前処理（？）というやつで、ある程度数を減らさないといけなさそうです。
しかし何をすれば良いかわからない・・・。
まずは専門家の真似をしよう！ということで、@yamasakit_ 先生の記事（ raziのDocker-composeで創薬レイドバトル2018用のJupyter Notebookからアクセスできる化合物データベースを作った話 ）とケモインフォマティクス若手の会の ハンズオン資料 を参考に前処理を行ってみました。*1
方針としては
- Google Colaboratoryを使ってみる（サイズ大きいから）
- 処理としてRDKitを使って
- 構造の標準化
- 脱塩（一番大きい左図のフラグメントだけを残す）
- 電荷の中和を行う
- idnumberと構造だけ残してsdfファイルにする
という感じです。
- Google colaboratoryで使うためのデータセットの準備
- Google colaboratoryでの作業
- どれくらいメモリをつかったのか？
- Google ドライブを使いたくない場合？
- 残念なお知らせ
- まとめ
Google colaboratoryで使うためのデータセットの準備
Google colaboratoryからデータセットを使うには、Googleドライブへのアップロードが必要そうでした。
まずは、サイズが大きいのでデータセットを圧縮します。
gzip -c Enamine_Advanced_collection.sdf > Enamine_Advanced_collection.sdf.gz gzip -c Enamine_HTS_collection.sdf > Enamine_HTS_collection.sdf.gz gzip -c Enamine_Premium_collection.sdf > Enamine_Premium_collection.sdf.gz gzip -c UOS_HTS.sdf > UOS_HTS.sdf.gz
Googleドライブにログインして上記のファイルをアップロードしましたが、圧縮してあっても時間がかかるので、この間にGoogle colaboratoryの準備をしました。
Google colaboratoryでの作業
1. Google colaboで新しいnotebookを作成
Google Colaboにアクセスして、「ファイル」から「Python 3の新しいノートブック」をクリックするだけで使えます。すごい・・・ *2
2. rdkitをインストール(ハンズオン資料そのまま)
!curl -Lo rdkit_installer.py https://git.io/fxiPZ
import rdkit_installer
%time rdkit_installer.install()
3. Google ドライブのファイルにアクセスするための準備
Google colaboratoryからGoogleドライブのファイルを使用する方法は脚注の記事（→*3 )を参考にしました。
!pip install -U -q PyDrive from pydrive.auth import GoogleAuth from pydrive.drive import GoogleDrive from google.colab import auth from oauth2client.client import GoogleCredentials auth.authenticate_user() gauth = GoogleAuth() gauth.credentials = GoogleCredentials.get_application_default() drive = GoogleDrive(gauth)
Google Cloud SDKの認証リンクが出てくるので、認証を行います。
4. GoogleドライブからファイルをColaboratoryのローカルにファイルを持ってくる
先の脚注の参考記事に従ってGoogle ドライブから共有リンク('id='以降)を取得します。 まずは一番ファイルサイズの小さいEnamine_Premium_collection から進めます。
id = 'Googleドライブのファイルのidをここに書く' downloaded = drive.CreateFile({'id': id}) downloaded.GetContentFile('Enamine_Premium_collection.sdf.gz')
これでColaboratoryからsdfファイルを使う準備ができました。
5. RDKitを使ってSDFファイルを読み込む
今回は圧縮したファイルなので、SDMolSupplierではなくForwardSDMolSupplierを使用します。
読み込み方法の詳細は脚注のページを参考にしました。（→ *4 ）
from rdkit import Chem import gzip EPc_gz = gzip.open('Enamine_Premium_collection.sdf.gz') EPc_mols = [mol for mol in Chem.ForwardSDMolSupplier(EPc_gz) if mol is not None]
読み込むことができた分子の総数を確認します。
len(EPc_mols)
「128816」と出力されました。 戦闘力１３万・・・強い・・・
6. プロパティの確認
SDFから読み込んだMolオブジェクトに、元々のSDFの情報がどのように紐づいているか確認するため、構造以外の情報を確認してみます。 Molオブジェクトのリストから最初の一つを取り出して、プロパティを取得しました。
詳細は「化学の新しいカタチ」さんの記事 *5 などを参考にしてください。
EPc_test = EPc_mols[0] names=list(EPc_test.GetPropNames()) print(names)
プロパティのリストが出力されました。
['idnumber', 'heavy_atoms', 'LogS', 'LogP', 'rotating_bonds', 'PSA', 'hb_acceptors', 'hb_donors']
MarvinViewでみたときとプロパティが違う・・・（linkというプロパティがない）。
先ほどMarvinViewでみたのはEnamine_HTS_collectionだったので、どうやらsdfファイルによって含むプロパティは違うみたいです。
7. 前処理の試し打ち
いよいよ本題です。
ハンズオン資料によると、MolStandardizeというモジュールを使うことで前処理ができるみたいです。
まずは試しに一つやってみます。 塩酸塩は無事、脱塩処理されるでしょうか？
まずは元の構造の確認
#構造式をnotebook上で表示させるための設定 from rdkit.Chem.Draw import IPythonConsole from IPython.display import display #塩酸塩のエントリー例 EPc_mols[2]
こんな構造がnotebook上に表示されます。
それでは処理を行います。
# 分子の標準化を行うためのモジュールを読み込む from rdkit.Chem import MolStandardize #標準化 normalizer =MolStandardize.normalize.Normalizer() test_mol_norm = normalizer.normalize(test_mol) #一番大きいサイズのフラグメントのみ残す（ここで脱塩されるみたい） lfc = MolStandardize.fragment.LargestFragmentChooser() test_mol_desalt = lfc.choose(test_mol_norm) #電荷の中和 uc = MolStandardize.charge.Uncharger() test_mol_neu = uc.uncharge(test_mol_desalt)
処理後の構造をそれぞれnotebook上に描画した結果は以下の通り・・・
| test_mol_norm | test_mol_desalt | test_mol_neu |
|---|---|---|
 |
 |
 |
無事、塩酸塩が脱塩されました！ 標準化と電荷の中和が必要かいまいち理解できませんが、とりあえず処理に含めておきます。
8. 前処理本番
時間がかかりそうなので、時間も測定してみます。
#セルの処理の時間測定 %%time #前処理を実行して新しいMOlオブジェクトのリストを作る #空のリストを作成 processed_EPc_mols = [] #ループを回せ!!! for i in range(len(EPc_mols)): mol = EPc_mols[i] #構造の標準化 normalizer =MolStandardize.normalize.Normalizer() mol_norm = normalizer.normalize(mol) #一番大きいサイズのフラグメントのみ残す lfc = MolStandardize.fragment.LargestFragmentChooser() mol_desalt = lfc.choose(mol_norm) #電荷の中和 uc = MolStandardize.charge.Uncharger() mol_neu = uc.uncharge(mol_desalt) #新しいリストに追加 processed_EPc_mols.append(mol_neu)
かかった時間・・・
CPU times: user 3min 6s, sys: 12.4 s, total: 3min 19s Wall time: 3min 19s
これならカップ麺作ってる間に前処理できる！！
念のため処理した後に残った分子の数を確認
len(processed_EPc_mols)
「128816」と出力されました。全部うまくいったみたいです。
9. 出力用のSDFを作成
SDFとして出力したいと思いますが、処理をおこなったので元々のプロパティのうち、
logSやlogPなどは意味のない値になってしまってそうです。
そこでidnumberと構造だけを含むSDFとしたいと思います。
SDFの出力にはSDWriterを使い、SetPropで紐づけたいプロパティを指定します。
#構造とidnumberのみを残したsdfファイルを作成 #SDWriterを使用する writer = Chem.SDWriter('Enamine_Premium_collection_processed.sdf') #プロパティの設定 #右だとうまくいきません。プロパティはリストで渡す必要があります。→ writer.SetProps('idnumber') writer.SetProps(['idnumber']) #ループを回せ！！ for mol in processed_EPc_mols: writer.write(mol) #そっ閉じ writer.close()
これでSDFができました。あとは自分のPCにもってくるだけ・・・ とりあえずGoogleドライブに出力します。
10. Googleドライブに出力
Google Colaboratoryローカルに入力する際に参考にした記事に、出力方法も載っていました。
念のため圧縮して出力します。
#大きいので圧縮 #colaboratoryでは"!"を先につけるらしい !gzip -c Enamine_Premium_collection_processed.sdf > Enamine_Premium_collection_processed.sdf.gz #Googleドライブへ出力 upload_file = drive.CreateFile() upload_file.SetContentFile('Enamine_Premium_collection_processed.sdf.gz') upload_file.Upload()
Googleドライブで無事出力されていることが確認できました！
ローカルPCにもってきてMarvinViewで見るとこんな感じ・・・
注 ) 上の図は修正前のものです。見た目は変わってしまいますが、こちらの方がわかりやすいので残しておきます。
どれくらいメモリをつかったのか？
Google colaboratoryは一度のセッションで12GBまでしか使えないそうです。
上記の処理では、最もファイルの大きさの小さいEnamine_Premium_collectionを使いましたが、
出力まで終えて「使用したRAM 3.91GB」でした。
Google ドライブを使いたくない場合？
今回、「Google colaboratoryのローカルにGoogleドライブからファイルをもってくる」、 という入力を行いましたが、後から考えたら「colabolatory上でファイルをダウンロードしてもよかったかな」、と思いました。
一応確認・・・
!wget https://xxxxxx #創薬レイドバトルのページにある化合物ライブラリのリンクをはる #library.tar.gz というのがcolaboratoryローカルにダウンロードされる #解凍する !tar -xvf library.tar.gz
解凍されて以下ができました。
- souyakuchan_library/
- souyakuchan_library/Enamine_Advanced_collection.sdf
- souyakuchan_library/Enamine_HTS_collection.sdf
- souyakuchan_library/Enamine_Premium_collection.sdf
- souyakuchan_library/UOS_HTS.sdf
これをさらにcolaboratory上でダウンロードしても使えるみたいなのですが、
試しにやってみたらとても遅いので、先述のGoogleドライブを経由する方法の方が良さそうでした。
残念なお知らせ
よし！残りのファイルも同じ感じでやるぞ！と思ったのですが
Enamine_HTS_collectionで同じ操作をおこなったところ、
メモリの使いすぎ（？）のためかエラーが出てしまいました。
ランタイムをリセットしてやり直しても12GBにおさまらず、前処理の途中で力尽きてしまいました。
自分のPCでは遅いからこっちでやろうと思ったのですが・・・残念
まとめ
以上、RDKit と Google colaboratory を使ってみた話でした。
お手元にRDKitをインストールしていない場合でもGoogle colaboratory上で遊べるっていうのはとても良いですね。
SDFをいい感じに分割したら、他のファイルも同じ処理ができるのではないかと思うのですが、私の能力では無理でした・・・
わかりやすくするために省きましたが、実際には例で検証したりせず、いきなり全構造処理したり、その他様々なトラブルで見事に丸一日を無駄にしましたよ！！！ sdf一つしか処理できなかった・・・
追記（01/14） SDWriterの設定が間違っていたので修正しました。 SetPropsはプロパティをリストかタプルで渡す必要があるので、 writer.SetProps('idnumber') を writer.SetProps(['idnumber']) のように修正する必要があります。
また、この記事の通りにやった場合、脱塩処理された化合物（フラグメントの大きい方が取り出されたもの）については idnumberがなくなってしまいます。
こちらは修正方法を検討して別の記事で書きます（いつか・・・）
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*1:
ハンズオン資料はtwitterをさまよってたらありました。とても素晴らしい資料ですよ！
@yamasakit_ 先生の記事の通りにやならかったのは、Dockerの使い方がわからなかったからです。。。
*2: 環境構築不要でPython入門！Google Colaboratoryの使い方を分かりやすく説明
*3:
*4:
「化学の新しいカタチ」さんの記事 RDKitでケモインフォマティクスに入門
-
*5:
