magattacaのブログ

日付以外誤報

SDFを読み込んでPythonのジェネレータを理解しようとした話

RDKitにはSDFを読み込む方法としてSDMolSupplierとFowardSDMolSupplierの大きく二つの方法があるそうです。

後者を使えば圧縮されたSDFを直接読み込めるということで、大きなデータを読み込むため使用してきましたが、今一両者の違いについて理解できない点がありました。問題の箇所はこの部分です。

毎度引用させていただいて恐縮ですがいますが、「化学の新しいカタチ」さんの参考記事1 から

SDMolSupplierとForwardSDMolSupplierのもう１つの違いは，前者はMolオブジェクトのリストを作成するのに対し，後者は異なるという点が挙げられます．そのため上の例ではsuppl[0]とすると１つめのMolオブジェクトにアクセス可能ですが，fsuppl[0]ではエラーを返します．

ForwardSDMolSupplierはリストではない・・・では一体なんなのか？

回答

「ジェネレータ 」

twitterで@iwatobipen先生に教えていただきました。「メモリの消費が少ないのがメリット」だそうです。 ありがとうございました！

・・・・・・格好つけてジェネレータとか呟いてしまいましたが、私、オブジェクトよく分からない。

ということでもうちょっと調べてみました。*1

入門 Python3の記述

まずは手元にあった書籍「入門 Python3 」から関係のありそうな記述を引用します。

ジェネレータは一度しか実行できない。リスト、集合、文字列、辞書はメモリ内にあるが、ジェネレータは一度に一つずつその場で値を作り、イテレータに渡して行ってしまうので、作った値を覚えていない。そのため、ジェネレータをもう一度使ったり、バックアップしたりすることはできない。入門 Python 3 オライリージャパン　Bill Lubanovic 著 斎藤　康毅　監訳　長尾　高弘　訳(p109)

ジェネレーターは、Pythonのシーケンスを作成するオブジェクトである。ジェネレータがあれば、シーケンス全体を作ってメモリに格納しなくても、（巨大になることがある）シーケンスを反復処理できる。・・・（省略）・・・ジェネレータは、反復処理のたびに、最後に呼び出されたときにどこにいたかを管理し、次の値を返す。これは、以前の呼び出しについて何も覚えておらずいつも同じ状態で１行目を実行する通常の関数とは異なる　同上 (p125)

上記と@momijiameさんのこちらの記事（参考記事2） Python: ジェネレータをイテレータから理解する を合わせると以下のようになりそうです。

コンテナオブジェクト ジェネレータオブジェクト
⑴ リスト、集合、辞書
⑵ ランダムアクセスできる できない
⑶ 値を取り出した後も値を覚えている 取り出すと値を忘れる
⑷ 何度でも反復して使用できる 一回きりしか使えない
⑸ 処理後、最初の場所に戻る 処理ごとに順番に一つずつ進んでいく
⑹ メモリ内に全ての要素を格納したまま 忘れるのでメモリを節約できる

以上の性質をふまえると、SDMolSupplierで読み込むとコンテナ、ForwardSDMolSupplierで読み込むとジェネレータ、となっていそうです。一つずつ違いを確認したいと思いますが、その前に両者の違いをまとめます。

どちらも分子のリストを作ることができますが、それぞれ

SDMolSupplier FowardSDMolSupplier
① 圧縮ファイルからは読み込めない 読み込める（file-like オブジェクトも使える）
② ランダムアクセスできる できない
③ 分子を取り出しても分子を覚えている 取り出すと忘れる
④ くり返し同じ分子のリストを作ることができる 一回きりしか使えない
⑤ 処理後、最初の場所に戻る 処理ごとに順番に一つずつ進んでいく（状態を覚えている）

実際のデータで検証してみる

検証用のデータとして73個の構造を含むSDF（'test.sdf'）を使用しました。 上記「参考記事1」 で例として用いられている、東京化成から購入可能なサリチル酸メチル誘導体のリストです。

from rdkit import rdBase, Chem
print(rdBase.rdkitVersion) 
#2018.09.1

①圧縮ファイルの読み込み

まずは圧縮されていないSDFで確認
#SDMolSupplierの場合
SDMSup = Chem.SDMolSupplier('./test.sdf')
SDMols = [x for x in SDMSup if x is not None]
len(SDMols)

「73」と出力されました。無事読み込めたようです。

#ForwardSDMolSupplierの場合
FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
FSDMols = [x for x in FSDMSup if x is not None]
len(FSDMols)

こちらも「73」と出力されました。

圧縮されたSDFの場合
import gzip
test_gz = gzip.open('./test.sdf.gz')
#SDMolSupplierの場合
SDMSup = Chem.SDMolSupplier(test_gz)

#エラーが返ってきた
---------------------------------------------------------------------------
ArgumentError                             Traceback (most recent call last)
<ipython-input-30-eb17d9fedebe> in <module>()
      1 #SDMolSupplierの場合
----> 2 SDMSup = Chem.SDMolSupplier(test_gz)

ArgumentError: Python argument types in
    SDMolSupplier.__init__(SDMolSupplier, GzipFile)
did not match C++ signature:
    __init__(_object*, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > fileName, bool sanitize=True, bool removeHs=True, bool strictParsing=True)
    __init__(_object*)

SDMolSupplierでは上記のエラーがでました。

#ForwardSDMolSupplierの場合
FSDMSup = Chem.ForwardSDMolSupplier(test_gz)
FSDMols = [x for x in FSDMSup if x is not None]
len(FSDMols)

「73」と出力されました。gzip.open()関数を使って圧縮ファイルを開いて得られたファイルライクオブジェクトは、SDMolSupllierでは読み込めず、FowardSDMolSupplierでは読み込むことができました。*2

②ランダムアクセス

次に、インデックスを使って好きな要素に自由にアクセスできるか検証します。

#SDMolSupplierの場合
SDMSup = Chem.SDMolSupplier('./test.sdf')
m10 = SDMSup[10]

forループを回してリストを作成しなくてもアクセスできました。

味気ないので構造を出してみます。

from rdkit.Chem import Draw
Draw.MolToImage(m10)

f:id:magattaca:20190119164523p:plain

構造もきちんと認識されています。次にForwardSDMolSupplierを試します。

#ForwardSDMolSupplierの場合
FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
fm10 = FSDMSup[10]

#エラーが返ってきた
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-45-bb94ea24d000> in <module>()
      1 #ForwardSDMolSupplierの場合
      2 FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
----> 3 fm10 = FSDMSup[10]

TypeError: 'ForwardSDMolSupplier' object does not support indexing

怒られました。indexによる呼び出しはサポートされていないそうです。

ひょっとしてSDMolSupllierはそのままリストととして扱えるのでしょうか？

読み込まれた分子の総数を確認するため、長さを取得できるか検証します。

len(SDMSup)
#73

おお！　forループを回さなくても長さを確認できました。

複数の分子を一度に取り出すことはできるでしょうか？

SDMSup[1:5]

#エラーが返ってきた
---------------------------------------------------------------------------
ArgumentError                             Traceback (most recent call last)
<ipython-input-56-66d817f86e32> in <module>()
----> 1 SDMSup[1:5]

ArgumentError: Python argument types in
    SDMolSupplier.__getitem__(SDMolSupplier, slice)
did not match C++ signature:
    __getitem__(RDKit::SDMolSupplier*, int)

怒られました。スライスでの取り出しには対応していないみたいです。SDMolSupplierそのものでリスト型と同じというわけではないみたいです。

結局、どういった操作ができるのでしょうか？ APIに果敢にアクセスしましたがさっぱりわかりません。多分こんな感じ。

SDMolSupplier FowardSDMolSupplier
Python API SDMolSupplier FowardSDMolSupplier
C++ API SDMolSupplier FowardSDMolSupplier
next できる できる
length できる できない
[ ] operator ("idx") できる できない

以上のようです。

SDMolSupplierならindexによる任意の場所のアクセスと、長さの取得はできますが、ForwardSDMolSupplierではできません。

イテレータ？

上記の表にnextという関数があります。こちらは一つずつ順番に取り出していくようです。

SDMSup = Chem.SDMolSupplier('./test.sdf')
m0 = next(SDMSup)
m1 = next(SDMSup)
m2 = next(SDMSup)

Draw.MolsToImage([m0, m1, m2])

以下の構造が取り出されたようです。

f:id:magattaca:20190119164836p:plain

indexで取り出したものと比較します。

Draw.MolsToImage([SDMSup[0], SDMSup[1], SDMSup[2]])

f:id:magattaca:20190119164836p:plain

indexでとりだしたものも、nextで取り出したものと同じ結果を返しました。

こちらの組み込み関数 next() はPythonのイテレータと深く関連しているそうです。

上記の「参考記事２」の説明によるとイテレータは「要素を一つずつ取り出ことのできるオブジェクト」であり「for 文こそ最も身近なイテレータ」とのことだそうです。

「イテレータ？そんな意味不明なもの一生使わんやろう。」と思っていましたが、forループを使ってSDMolSupplierをMolオブジェクトのリストに変換している時点でとてもお世話になっていたみたいです。

知らなかった・・・

いずれにせよSDMolSupplierのままでもいくつかの操作は行えるみたいですが、リストに変換しておいた方が後々便利そうです。

コンテナ？　ジェネレータ？

forループはイテレータということでしたが、「イテレータオブジェクトをつくるもと」はなんなのか？　というとコンテナオブジェクトとジェネレータオブジェクトということになるようです。

・・・やっと冒頭のテーブルに話が戻ってきました。

それではコンテナとジェネレータの性質の違いを順番に見ていきたいと思います。

記憶力（③）と反復利用（④）の可否

まずはSDMolSupplierで読み込んだ場合

#一回め
SDMSup = Chem.SDMolSupplier('./test.sdf')
SDMols = [x for x in SDMSup if x is not None]
len(SDMols)
#73

イテレータの作成に使用したオブジェクト、SDMSupはもう一度使用できるでしょうか？

#２回め
SDMols_2 = [x for x in SDMSup if x is not None]
len(SDMols_2)
#73

２回めも同数の結果が返ってきました。SDMolSupplierによるオブジェクトは元の状態を覚えており反復利用できました。（コンテナ）

では、ForwardSDMolSupplierではどうでしょうか？

#１回め
FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
FSDMols = [x for x in FSDMSup if x is not None]
len(FSDMols)
# 73

#2回め
FSDMols_2 = [x for x in FSDMSup if x is not None]
len(FSDMols_2)
# 0

0個！　

一度イテレータに渡して全ての分子を読み込むと、FSDMSup は空っぽになってしまいました。ForwardSDMolSupplierによるオブジェクトは、使用の度に忘れ、同じ内容は一度しかつかえないようです。（ジェネレータ）

⑤処理後の状態

先の処理ではイテレータに一度にSDFの73個全ての分子を渡してしまいました。もし、一度途中で渡すのをやめたらどうなるでしょうか？

まずはSDMolSupplierの場合です。

#念のため全部再読み込み
SDMSup = Chem.SDMolSupplier('./test.sdf')
#最初の30個を読み込み
SDMol_list_1 = []
for i, j in zip(range(30), SDMSup):
    SDMol_list_1.append(j)
len(SDMol_list_1)
#30

さらに10個取り出して見ます。

#続けて10個取り出し新しいリストに格納
SDMol_list_2 = []
for i, j in zip(range(10), SDMSup):
    SDMol_list_2.append(j)
len(SDMol_list_2)
#10

それぞれ最初の分子の構造を描画します。

Draw.MolsToImage([SDMSup[0], SDMol_list_1[0], SDMol_list_2[0]], 
                 legends=['SDMSup[0]', 'SDMol_list_1[0]', 'SDMol_list_2[0]'])

f:id:magattaca:20190119165436p:plain

すべて同じ構造となりました。一度構造を取り出したあと、次に取り出すときも最初に戻って取り出しているようです。

次に同様の処理をFowardSDMolSupplierで行います。

#念のため全部再読み込み
FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
#最初の30個を読み込み
FSDMol_list_1 = []
for i, j in zip(range(30), FSDMSup):
    FSDMol_list_1.append(j)
len(FSDMol_list_1)
#30

#続けて10個取り出し新しいリストに格納
FSDMol_list_2 = []
for i, j in zip(range(10), FSDMSup):
    FSDMol_list_2.append(j)
len(FSDMol_list_2)
#10

#それぞれ最初の分子の構造を描画
Draw.MolsToImage([FSDMol_list_1[0], FSDMol_list_2[0], SDMSup[30]], 
                 legends=['FSDMol_list_1[0]', 'FSDMol_list_2[0]', 'SDMSup[30]'])

f:id:magattaca:20190119165551p:plain

２度めに呼び出した10個の構造のリスト（FSDMol_list_2）の最初の分子は、初めの30個のリストの（FSDMol_list_1）の最初の分子と異なりました。元々のSDFの31番目の分子（SDMSup[30]）と構造が一致していることから、30個呼び出したあとで、２度めは31番めからの呼び出しが始まっていることがわかります。

SDMolSupplierは処理後、最初の場所に戻り、FSDMolSupplierでは、処理ごとに順番に一つずつ進み、処理後の状態を覚えているということが確認できました。

メモリの消費

以上から、SDMolSupplierはコンテナ（シーケンス(?)）の特徴をもっており、FSDMolSupplierはジェネレータの特徴を持っていることがなんとなくわかってきました。

では、後者を使うメリットの「メモリの消費」というのはどういうことなのでしょうか？

参考記事2を再度引用すると、コンテナオブジェクトではいつでも取り出せる状態にしておくため、「あらかじめ全ての要素をメモリに格納しなければならない」、一方、「そのつど生成した値を使い終わったら後は覚えておく必要はない」場合、「変数から参照されなくなったら要素はガーベジコレクションの対象となるためメモリの節約につながる」とのことだそうです。

ガベージコレクション（garbage collection、GC）というのは、参照されなくなったと判断されたメモリ領域を自動的に解放する仕組み、だそうです。C言語などでは、プログラムの作成者がメモリのマネジメントも行う必要があるそうですが、Pythonではその必要がなく、メモリの利用の宣言や、解放といったことを明示的にプログラミングする必要がない、とのことです。 （参考書籍2 (p328)）

FowardSDMolSupplierを使う時の注意点

以上見てきた違いから、非常に多数の分子を取り扱う必要がある場合、

  • 圧縮ファイルを読み込める
  • メモリの消費が少ない

という点で、FowardSDMolSupplierを使う方が良さそうです。ただし、使用にあたって気をつけるべきこともあります。

FowardSDMolSupplierは読み込んだ順番に分子を忘れていきますが、これは処理に失敗した場合でも同じです。

どういうことか、わざと失敗して確かめて見ます。

FSDMSup = Chem.ForwardSDMolSupplier('./test.sdf')
#存在しないリストを指定してエラーを出す
for x in FSDMSup:
    dummy.append(x)

---------------------------------------------------------------------------
NameError: name 'dummy' is not defined

#間違いに気づき空のリストを用意する
dummy = []
for x in FSDMSup:
    dummy.append(x)
len(dummy)

#72

73個の分子を含むはずのSDFなのに、72個しかリストに格納されていません。

一番最初の分子を確認します。

Draw.MolsToImage([dummy[0], SDMSup[1]],
                 legends=['dummy[0]', 'SDMSup[1]'])

f:id:magattaca:20190119165837p:plain

dummyリストの最初の分子は、元のSDFの2番めの分子と構造が一致します。

元のSDFの最初の分子は、エラーを出した際に忘れ去られてしまい、あとから正しい処理を書き直しても、もうその処理には送られません。なので、再度もともとのSDFをFowardSDMolSupplierで読み込むところからやり直す必要が生じます。

私はこのことに気づかず、誤った処理を書いて直すたびに分子の総数が変化し、混乱するということに陥っていました。

まとめ

以上、今回RDKitからSDFを読み込むための二つの方法、SDMolSupplierとFowardSDMolSupplierの違いを検証することでPythonのコンテナ、イテレータ、ジェネレータといった用語の雰囲気をつかむことを試みました。Pythonの入門書を読んで用語は目にしたことがあっても結局なにがしたいのかよくわかっていなかったので、実際に色々触ってみるとなんとなくやりたいことがわかってきました。

私はRDKitもPythonも初心者なので、用語の使い方や理解など間違っていることも多いと思います。ご指摘いただければ幸いです。

~~~~~~~~~~~~~~~~~~~~~~~~~

本記事の作成にあたって参考にさせていただいページ、書籍

参考記事1: 「化学の新しいカタチ」 RDKitでケモインフォマティクスに入門

参考記事２: 「CUBE SUGAR CONTAINER」 Python: ジェネレータをイテレータから理解する - CUBE SUGAR CONTAINER

参考書籍１: 入門 Python 3 オライリージャパン　Bill Lubanovic 著 斎藤　康毅　監訳　長尾　高弘　訳

参考書籍2: 科学技術計算のためのPython入門 技術評論社　中久喜　健司 著

*1:* ジェネレータと書くとジェネレータ関数かジェネレータオブジェクトかわからないと、ネットで言われていましたが、両者の違いを判別できるまで理解できていないので、以下ジェネレータとのみ表記します。

*2:gzipの使い方にしてはこちらの記事を参考にしました。

失われたIDを探せ！

先の２つの記事でrdkitを使った前処理を試みてきましたが、sdf出力の際の設定が間違っており、各エントリのidnumberが新しいsdfには出力できていませんでした。

SDWriterのSetPropsの設定を修正することで、idnumberをsdfに含めることはできたのですが、そもそもMolオブジェクトからidnumberがなくなってしまっているものがありました。どうやら、複数のフラグメントをもつオブジェクトについて、脱塩処理（一番大きいサイズのフラグメントのみ残す）をした場合に、idnumberが失われている、というように見えました。

構造を頼りに化合物の数を絞りこんだとしても、もともとの化合物のidnumberがわからないと出所がわからず困ったことになりそうです。そこで、再度idnumberを保持した形で前処理を実施したいと思います。折角なのでついでに分子量（Moleculer weight）も計算して、情報を付け足して見たいと思います。

idnumberを保持したSDFの生成

懲りずにGoogle Colaboratoryでやりますが、処理までの操作は変わらないので割愛します。

from rdkit import Chem
import gzip

#分子の前処理のためのMolStandardize、
#および分子量計算のためのDescriptorsをimport
from rdkit.Chem import MolStandardize, Descriptors

EPc_gz = gzip.open('Enamine_Premium_collection.sdf.gz')
EPc_suppl = Chem.ForwardSDMolSupplier(EPc_gz)  
EPc_mols_pro = []

for x in EPc_suppl:
  if x is not None:
    mol = x
  
  #処理の前にidnumberを取り出しておく(前回との違い)
  ID = mol.GetProp('idnumber') 
  
  #構造の標準化
  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)
  
  #以下は前回と異なる追加の処理
  
  #処理後のMolオブジェクトに、取り出しておいた元々のidnumberをoriginal_idとして付与
  mol_neu.SetProp('original_id', ID)
  
  #分子量を計算
  MW_value = Chem.Descriptors.MolWt(mol)
  
  #分子量をプロパティとして持たせる。
  #小数点を含むfloat型のプロパティなので SetDoubleProp を使って情報を付加する
  mol_neu.SetDoubleProp('MW', MW_value)
  
  #新しいリストに追加
  EPc_mols_pro.append(mol_neu)

かかった時間は + 1分くらいです。 　CPU times: user 3min 47s, sys: 13.9 s, total: 4min 1s 　Wall time: 4min 1s

念のため最初の分子でそれぞれのプロパティが格納されているか確認したいと思います。

test = EPc_mols_pro[0]
print('idnumber:', test.GetProp('idnumber'))
print('original_id:', test.GetProp('original_id'))
print('MW:', test.GetDoubleProp('MW'))

idnumber: Z1498649509
original_id: Z1498649509
MW: 210.28099999999995

と出力されました。うまくいってそうです。

Molオブジェクトのプロパティの設定で、idnumberと分子量では異なるものを使いました。 ケモインフォマティクスのオンライン入門書によると、

関数 データの型
SetProp('名前',値) str型
SetIntProp('名前',値) int型
SetUnsignedProp('名前',値) int型（非負）
SetDoubleProp('名前',値) float型
SetBoolProp('名前',値) bool型

と、なっているそうです。値を取り出すときはSetの部分をGetにすれば良さそうです。

sdfで出力し、自分のPCで確認します。

#構造とidnumber、original_id、MWをもつsdfファイルを作成
writer = Chem.SDWriter('Enamine_Premium_collection_pro_id_MW.sdf')

#プロパティとして持たせたいものをリストにしておく
prop_names = ['idnumber', 'original_id', 'MW']

writer.SetProps(prop_names)
for mol in EPc_mols_pro:
  writer.write(mol)
writer.close()

#大きいので圧縮
!gzip -c Enamine_Premium_collection_pro_id_MW.sdf > Enamine_Premium_collection_pro_id_MW.sdf.gz

#Googleドライブへ出力
upload_file = drive.CreateFile()
upload_file.SetContentFile('Enamine_Premium_collection_pro_id_MW.sdf.gz')
upload_file.Upload()

MarvinViewで開くとこのような見た目です。

f:id:magattaca:20190114233416p:plain

今度こそ上手くいったのではないでしょうか。

idnumberがなくなっていたものについて原因を検証

該当の構造の抜き出し

idnumberが失われてしまった構造は元々どのようなものだったのか、気になります。indexとして抜き出してみます。

lost_id_index = []

for i in range(len(EPc_mols_pro)):
    mol = EPc_mols_pro[i]
    
    #Molオブジェクトの持つプロパティの一覧を取得する
    props = mol.GetPropNames()
    
    #idnumberを持たない場合はそのindex（化合物idではない）をリストに加える
    if 'idnumber' not in props:
        lost_id_index.append(i)

idnumberが失われた構造の総数を確認してみます。

len(lost_id_index)

「2486」個と出ました。

これらの化合物について処理前と後で構造の変化を眺めてみたいと思います。

PandasToolsを用いた処理前後の構造の比較

比較のため再度、処理前のsdfからMolオブジェクトのリストを作成します。

EPc_gz2 = gzip.open('Enamine_Premium_collection.sdf.gz')
EPc_mols2 = [mol for mol in Chem.ForwardSDMolSupplier(EPc_gz2) if mol is not None]

構造の比較のためPandsToolsを利用してみたいと思います。

RDKitのPandasToolsについては「化学の新しいカタチ」さんのこちらの記事などを参考にしてください。

from rdkit.Chem import PandasTools
import pandas as pd

#処理前の構造から比較に使いたい構造を取り出したリストを作成する。
before_list = []

for i in lost_id_index:
    before_list.append(EPc_mols2[i])
    
#リストからSDFを作成したのち、PandasToolsのLoadSDFで読み込みDataFrameとする。
writer_b = Chem.SDWriter('before.sdf')
writer_b.SetProps(['idnumber'])
for mol in before_list:
  writer_b.write(mol)
writer_b.close()

df_before = PandasTools.LoadSDF('./before.sdf')

#処理後の構造についても同じことをする。
after_list =[]

writer_a = Chem.SDWriter('after.sdf')
writer_a.SetProps(['original_id'])
for mol in after_list:
  writer_a.write(mol)
writer_a.close()

df_after = PandasTools.LoadSDF('./after.sdf')

#2つのDataFrameを連結する。
df_compare = pd.concat([df_before, df_after], axis=1)

#全て取り込めているかを確認
df_compare.shape

(2486, 6) と返ってきました。無事2486個処理されたみたいです。

それでは早速構造の比較を見てみましょう

df_compare.head()

f:id:magattaca:20190114233756p:plain

塩酸塩となっていました。最初の考え通り、フラグメントを取り除く処理をしたものについてidnumberのプロパティが失われていた可能性が高そうです。

ついでに、最後の方の番号も見てみましょう。

df_compare.tail()

f:id:magattaca:20190114233904p:plain

後ろの方はカルボキシル基がイオン化してナトリウム塩となっているエントリーでした。 処理後の構造ではカルボキシル基は水素が付加されており、中和処理も問題なく行われているようです。

以上、idnumberが失われる場合の検証を行ってみました。予想外に処理が行われているものについて、実際に構造を比較し、処理の結果を検証することまでできました。欠損値（?）っていうのも大事なんですね。

たぶん、本当は処理の段階で、処理が必要になったものについては別のリストを作るとかそういう対応をするんだと思います。

すごい遠回りしてる感ある・・・

メモリの節約でColaboratoryの壁を越えようとする話

先の記事でGoogle colaboratoryの制限により、大きなファイルの処理ができなかったという結果になりました。

「先にMolオブジェクトをまとめて作るのではなくてforの中で一つずつ作ればメモリの節約にはなる」と教えていただいたので早速試してみたいと思います。 @yamasasaKit_先生ありがとうございました！

処理としては先の記事のままで、sdfファイルの読み込み以降の操作が異なります。

Google colaboratoryで諸々準備をした後に以下を実行しました。

#必要なモジュールの読み込み
from rdkit import Chem
import gzip
from rdkit.Chem import MolStandardize

#Googleドライブから取ってきた圧縮ファイルの読み込み
EPc_gz = gzip.open('Enamine_Premium_collection.sdf.gz')

今回もFowardSDMolSupplierで読み込みますが、Molオブジェクトをforの中でつくるため、リスト内包表記にはいれません。（←前回との違い）

#ForwardSDMolSupplierで読み込む。まだMolオブジェクトのリストにしない。
EPc_suppl = Chem.ForwardSDMolSupplier(EPc_gz)

それでは脱塩処理してMolオブジェクトのリストを作成します。 時間も測ってみます。

%%time

#molオブジェクトのリストを作る段階で前処理を実行してメモリを節約する。
#空のリストを作成

EPc_mols_pro = []

#ループ!!!
for x in EPc_suppl:
  if x is not None:
    mol = x
  
  #構造の標準化
  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)
  
  #新しいリストに追加
  EPc_mols_pro.append(mol_neu)

かかった時間・・・

CPU times: user 3min 37s, sys: 12.5 s, total: 3min 50s Wall time: 3min 50s

ここまでの「使用したRAM 1.85GB」でした。

念のため処理できた数を確認して見ます。

len(EPc_mols_pro)

「128816」とでました。前回も「128816」だったので、同じ数だけできているようです。

一応出力して、ローカルで確認します。

#構造とidnumberのみを残したsdfファイルを作成
writer = Chem.SDWriter('Enamine_Premium_collection_processed_2.sdf')
writer.SetProps(['idnumber'])        #(01/14)追記 SetPropsの引数をリストに変更
for mol in EPc_mols_pro:
  writer.write(mol)
writer.close()

#圧縮して出力
!gzip -c Enamine_Premium_collection_processed_2.sdf > Enamine_Premium_collection_processed_2.sdf.gz

upload_file = drive.CreateFile()
upload_file.SetContentFile('Enamine_Premium_collection_processed_2.sdf.gz')
upload_file.Upload()

出力まで行って、「使用したRAM 1.92GB」。前回は「3.91GB」だったので、約半分です。 Molオブジェクトのリストを２回作っていたのが、１回になったからでしょうか？

まさかこんなに簡単にメモリの節約ができるとは・・・・。

ちなみに自分のPCでsdfを確認すると前回と同じような構造が出ていました。

f:id:magattaca:20190114122649p:plain

注 ) 上の図は修正前のものです。見た目は変わってしまいますが、こちらの方がわかりやすいので残しておきます。

より大きなファイルで同じ処理をした結果は以下の通り。
ファイルごとにセッションを新しくして実行しています。

Enamine_Advanced_collectionの場合

  • 処理の時間 CPU times: user 13min 27s, sys: 54.8 s, total: 14min 22s Wall time: 14min 22s
  • 化合物総数 486322
  • 全体のメモリ 6.21GB

UOS_HTS_collectionの場合

  • 処理の時間 CPU times: user 15min 17s, sys: 54.5 s, total: 16min 11s Wall time: 16min 12s
  • 化合物総数 516664
  • 全体のメモリ 7.34GB

Enamine_HTS_collectionの場合

残念ながら最後まで至らず・・・ファイルサイズがEnamine_Advanced_collectionの５倍くらいあるので仕方ないかもしれません。

私のパソコンでやると総数「1921120」でした。

創薬ちゃんさんによると、ライブラリの構成は以下のようになっているそうです。

About 2.5 million・・・・ここからどう絞り込めば・・・

RDKitとcolabolatoryで前処理しようとして挫折した話

RDKitの使い方を調べたりインフォマティクスに関する本を読んでみたりしましたが、結局あまり実際には動かせていません。
おそらく、写経（？）したりして身につける必要があるのでしょうが、いまいちそれだと楽しくない（ダメな入門者の典型）。

何か面白いこと・・・ということで、調子に乗って創薬ちゃんさんの創薬レイドバトルに登録してしまいました。

どうしよう・・・

とりあえずどんな化合物から絞り込めば良いのか化合物ライブラリを眺めてみることにしました。

化合物ライブラリをダウンロード、解凍してみると以下の４つのsdfファイルが入っていました。

f:id:magattaca:20190113215243p:plain

試しにsdfを開いてみたらMarvinViewが起動し手軽に見ることができました。
いつのまにこんなソフトが？と思ったら、ChemDrawの代わりに使えるものをということでインストールしたMarvinSketchと一緒にダウンロードされていたみたい・・・ChemAxon社様様です。

f:id:magattaca:20190113215514p:plain

中身は上記のような感じ・・・。構造に加えてidnumber、heavy_atoms、LogSなどなどのカラムがあるようです。

思ったよりシンプルな化合物です。もっとドラッグライクな化合物はあるのでしょうが、なにぶんデータが多すぎてよくわからない。 これは前処理（？）というやつで、ある程度数を減らさないといけなさそうです。

しかし何をすれば良いかわからない・・・。

まずは専門家の真似をしよう！ということで、@yamasakit_ 先生の記事（ raziのDocker-composeで創薬レイドバトル2018用のJupyter Notebookからアクセスできる化合物データベースを作った話 ）とケモインフォマティクス若手の会の ハンズオン資料 を参考に前処理を行ってみました。*1

方針としては

  • Google Colaboratoryを使ってみる（サイズ大きいから）
  • 処理としてRDKitを使って
    1. 構造の標準化
    2. 脱塩（一番大きい左図のフラグメントだけを残す）
    3. 電荷の中和を行う
  • idnumberと構造だけ残してsdfファイルにする

という感じです。

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上に表示されます。

f:id:magattaca:20190113221309p:plain

それでは処理を行います。

# 分子の標準化を行うためのモジュールを読み込む
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
f:id:magattaca:20190113221438p:plain f:id:magattaca:20190113221503p:plain f:id:magattaca:20190113221520p:plain

無事、塩酸塩が脱塩されました！ 標準化と電荷の中和が必要かいまいち理解できませんが、とりあえず処理に含めておきます。

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で見るとこんな感じ・・・

f:id:magattaca:20190113222354p:plain

注 ) 上の図は修正前のものです。見た目は変わってしまいますが、こちらの方がわかりやすいので残しておきます。

どれくらいメモリをつかったのか？

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:

お正月読書感想文

１月も１２日目ですが、あけましておめでとうございます。

インシリコ創薬に興味を持ち勉強を始めてみましたが、そもそもインフォマティクス（？）に関する知識が全然無い・・・ということで、年末年始の休暇を利用して数冊本を流し読みしてみました。

 

コードも数式も理解できないので読み飛ばしてしまいましたが、折角なので、備忘録を兼ねて感想を書きます。

 

１冊目：よくわかるバイオインフォマティクス 入門 

よくわかるバイオインフォマティクス入門 (KS生命科学専門書)

よくわかるバイオインフォマティクス入門 (KS生命科学専門書)

  • 作者: 岩部直之,川端猛,浜田道昭,門田幸二,須山幹太,光山統泰,黒川顕,森宙史,東光一,吉沢明康,片山俊明,藤博幸
  • 出版社/メーカー: 講談社
  • 発売日: 2018/11/19
  • メディア: 単行本（ソフトカバー）
  • この商品を含むブログを見る
 

 

ケモインフォマティクス についての入門書はどれを読めば良いか分からなかったので、まずはバイオインフォマティクス の本をということでTwitterでオススメされていたこちらの本を購入しました。*1

 

目次をご覧いただくのが一番早いと思うので講談社サイエンティフィク さんのページより拝借・・・

 　1章　配列解析

 　2章　分子進化

 　3章　タンパク質の立体構造解析

 　4章　ncRNA解析

 　5章　NGSデータ概論

 　6章　ゲノム解析

 　7章　トランスクリプトーム解析

 　8章　エピゲノム解析

 　9章　メタゲノム解析

  10章　プロテオーム解析

  11章　データベース

  12章　バイオのための機械学習概論

 

バイオインフォマティクスに関し全くの門外漢だったのですが、分子生物学を中心として、遺伝子（上流、１次元）からタンパク質の構造・機能解析（下流、３次元）まで、さまざまな段階においてデータ解析を活用した研究を行うといった印象でした。

各トピックに関して「そもそもそれをなぜ研究したいのか」という分野の意義の説明から、その分野において「インフォマティクスの技術はどう使えるか」、また「実際に使うにはどうすれば良いか」という具体的なデータベースのアクセス・解析の原理と方法についてカラーの図で解説してあり非常に面白かったです。

 

　ところで、分子生物学というと細胞内部や、個体内部といったミクロな視点がメインという印象を勝手にいだいていたのですが、こちらの本ではよりマクロな視点も取り扱われていました。具体的には、「１章 配列解析」における、生物種間でのタンパク質の差異の検証に基づき進化の分岐点を探る、といった歴史を辿るような試みや、「９章 メタゲノム解析」における、ある場所に存在する多種多様な微生物の集まりをひとまとめにして、遺伝子プールという観点からその集団の持つ特性を評価する、という群としての取り扱いなども取り上げられています。

　私は薬学部出身なのですが、大学に進学した際に「高校の生物学ではエンドウ豆のシワの話やら、淡水魚・海水魚と浸透圧、みたいな話をしていてわかりやすかったのに、なぜ大学はタンパク質の話しかしないんだ？？？　もっとマクロな話はどうなった？？」と、「生体内分子の機能を明らかにすることで生物を理解する」という趣旨を理解するのにとても時間がかかったダメ学生でした・・・農学部や理学部に行ったら違ったのかも。なので、この本を読んで「マクロな生物学の解析にこんな風に進んでいたのか！」と非常に面白かったです。*2

 

　個人的に最もオススメしたいのは「６章 ncRNA解析」です。セントラルドグマのせいかRNAに対して「DNAとタンパク質の間をつなぐ影の存在」的な印象を抱いていたので、タンパク質には翻訳されずRNA自身が生体内で機能する（noncoding RNA: ncRNA）の研究が近年急速に進歩していること、機能と結びつくRNAの２次・高次構造解析、RNA-RNA / RNA-タンパク質 相互作用解析 といった研究が広がっていることを知りませんでした。20ページもありませんが、その中でncRNAの基礎から最新の研究への流れと、今後の展望について一気に語られており、「急速に知識を詰め込まれていく感じとても気持ち良い！」*3という感じです。

 

　総じて、各章「そもそもそれをなぜ研究したいのか」という分野の意義の説明から、その分野において「インフォマティクスの技術はどう使えるか」、また「実際に使うにはどうすれば良いか」という具体的なデータベースのアクセス・解析の原理と方法についてカラーの図で解説してあり非常に面白かったです。カラーなのでお値段高めですが、理工書としてはお手頃ですしこれから入門したいという方にはオススメ！

 

２冊目：Pythonデータサイエンスハンドブック  

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習

Pythonデータサイエンスハンドブック ―Jupyter、NumPy、pandas、Matplotlib、scikit-learnを使ったデータ分析、機械学習

 

 

　過去の記事でPandasのdataframeを使ってみたりしましたが、「そもそもPandasって何だ？ハードル高いエクセルか？もうエクセルでよくない？？？」ぐらいの低レベルなので購入してみました。　エクセルもsumしか使えないんですけどね・・・

原書は無料で公開されているそうですが、英語だと多分１年かかるし紙の本じゃないと140字以上読めない・・・。

*こちらブログの書評が素晴らしいので有用な情報が欲しい方はこちらをお読みください。

  →  『Pythonデータサイエンスハンドブック』は良書（NumPy, pandasほか） | note.nkmk.me

 

目次はこんな感じ

　１章 IPython : Pythonより優れたPython

　２章 NumPyの基礎  

　３章 pandasを使ったデータ操作

　４章 Matplotlibによる可視化

　５章 機械学習

 

　Pythonを使ってデータ解析を行う際に基礎となるライブラリが紹介されているそうです。データ解析をやったことないのでよくわかりませんが、とりあえず素人理解では、

・IPython → インタラクティブに検証可能なので解析が捗る

・NumPy → 科学計算でPythonを扱うのに有用な数値配列のパッケージ

・Pandas → NumPyの上に構築されたデータ操作に有用なパッケージ

・Matplotlib → 綺麗なグラフ描ける

・機械学習 → scikit-learnを使えば原理をよくわからなくても試すことができる

ということみたいです。

章を読み進めるに従って、

・NumPyがPythonの配列よりもなぜ効率がよく効果的か

・NumPyの制限をPandasがいかに解決したか

・Pandasで扱うデータをMatplotlibを使えばどう効果的にビジュアル化できるか

順番に理解していくことができるため、初心者にとって非常にわかりやすい構成になっていました。　（読み飛ばしましたが）コピペして使えそうなコードが多数あり、また「４章 Matplotlib」では描きたグラフを描くための細かな調整の方法が詳しく書いてあるので、手元に置いておいて困った時に参照しよう、という内容になっています。

 

　また、先の「よくわかるバイオインフォマティクス 入門」の「12章　バイオのための機械学習概論」は名前の通り概論で機械学習の具体的な内容まではよくわからなかったのですが、こちらの「Pythonデータサイエンスハンドブック」ではデータサイエンスを謳っているだけあり、機械学習の事例・実行するためのコードについて多数紹介されていて先の補完の意味でもとても勉強になりました。

 

　とりあえず大きなデータを扱ってやりたい分析をできるようにするには、この本で紹介されているような内容を把握なければいけなさそう・・・ということで出発点として非常に良い本なのかな、という感想です。

 

３冊目：PythonとKerasによるディープラーニング

PythonとKerasによるディープラーニング

PythonとKerasによるディープラーニング

 

 

　やはりディープな話題にも触れてみなければならないのではないのか、ということで購入。ネットで購入して届いたら予想以上に分厚くて読む前から挫折しそうになったのですが、非常に読みやすくなんとか読み通すことができました。

まあ、コードを読まずに、図を追いかけただけなんですけどね。何となく雰囲気はつかめたので良し。

 

目次をマイナビブックスさんから・・・（サポートサイトへのリンクもあります）

　Part 1　ディープラーニングの基礎

　　1章　ディープラーニングとは何か

　　2章　予習：ニューラルネットワークの数学的要素

　　3章　入門：ニューラルネットワーク

　　4章　機械学習の基礎

   Part 2　ディープラーニングの実践

      5章　コンピュータビジョンのためのディープラーニング

      6章　テキストとシーケンスのためのディープラーニング

      7章　高度なディープラーニングのベストプラクティス

      8章　ジェネレーティブディープラーニング

      9章　本書のまとめ

 

　ディプラーニングの基礎的な話から、具体的に行うにはどうすれば良いのかKerasを使って説明されてます。原書の著者はKerasの作成者（！）だそうで、「LEGOブロックを組み立てる」ように簡単な「ユーザーフレンドリなライブラリ」により「ディープラーニングが大衆化」し、新規参入者が増えたことが分野を後押ししている、とのことです。もちろん作成者だけに現状での限界や、できないことについても記述してあり、その点が非常に良いと思いました。

 

　ケモインフォマティクス の分野から勉強を始めたため、あまり認識していなかったのですがディープラーニングの強い分野はやはり画像認識、ということで様々な興味深い事例が紹介されていました。悪夢的な画像を生成するDeepDreamや、任意の画像をゴッホが描いたかのように変換するニューラルスタイル変換など、やはり画像に関する技術はビジュアル的なインパクトが強く、門外漢にとっては読んでいてワクワクする純粋に楽しめる内容でした。

　また、ケモインフォ の記事をみていて見かける化合物生成に関連する技術、変分オートエンコーダ（VAE）や敵対生成ネットワーク（GAN）についてもとりあげられており、元々はこういう技術だったのか！と出所がわかり面白かったです。

 

とりあえず、ネットワークの損失関数をオプティマイズすれば良いってさ　（←よくわからなかった）

 

４冊目：データ解析のための統計モデリング入門

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

 

 

　流行りのディープラーニングについて浅すぎる理解ができたので、次は機械学習に関連して各所でオススメされている本を・・・ということで緑本（で、あってます？）を購入しました。ハードカバーと難しそうな単語が並ぶタイトルからしてハードルが高かったのですが、驚くほど分かり易く読み易い語り口で一気に読み通すことができました。いろいろなブログでオススメされているのも納得という良書・・・もちろんコードは追えていない

 

著者の久保拓弥先生のページで出版のもとになった講義のーと が公開されており、章立てと全体の流れの解説図があるため目次は割愛します。

 

　さて、先の「PythonとKerasによるディープラーニング」の中で、ディープラーニングと従来のシャローラーニング（表層学習）との違いとして、ニューラルネットワークは生のデータから有益な特徴量を自動的に抽出できるため、特徴エンジニアリング（特徴量をより単純な方法で表現することで、問題を容易にすること）を必要とせずに行うことができる、ということが挙げられていました（４章 p.106）。だからと言って特徴エンジニアリングを疎かにすべきではないし、重要性として

　・リソースの消費を抑えた上で問題を的確に解決できる

　・より少ないデータで解決できる

があげられていましたが、やはりディープラーニングの大きな長所として特徴抽出に関して何度も述べられていました。

　他方、「データ解析のための統計モデリング入門」はまさに上記と対になるような内容という印象でした（間違ってたらすみません）。与えらえれたデータをもとにどのようにしてモデルを構築すれば良いか？　最小二乗法によるシンプルな線形モデルから、より高度な一般化線形モデル、階層ベイズモデルへと順をおってより複雑なモデルが説明されています。データのバックグラウンドを考慮しながら、複雑な自然のデータを解釈するためのモデルを作り上げていく流れは、”ドメイン知識に基づく特徴の抽出”といった感じで、研究者の思考をたどるような面白さがありました。

　残念ながらベイズ統計（というより統計全体）を理解していないため、後半、話の流れを追えなかった部分もあるのですが、とりあえず初心者でも一度は目を通した方が良い本だ、と思いました。

 

５冊目：これならわかる最適化数学

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

 

 

　統計モデリング入門は大変分かりやすかったのですが、やはり数学的素養があまりにもなさすぎてところどころ理解できない部分がありました。そこで「化学の新しいカタチ」さんのこちらの記事 でDr. Tomさんがオススメされていた、こちらの本を購入しました。

目次はこのような感じ・・・

　第1章 数学的準備

　第2章 関数の極値

　第3章 関数の最適化

　第4章 最小二乗法

　第5章 統計的最適化

　第6章 線形計画法

　第7章 非線形計画法

　第8章 動的計画法

 

　やはりDr. Tomさんがオススメされているだけあって非常に良い本で、読み進めるための数学的基礎の準備から、いかに関数を処理し最適化を行えば良いか、非常に分かりやすく解説されています。特にたまに見かけるけどもよくわからないまま放置していた「ラグランジュの未定乗数法」についても取り扱われていたので、個人的にはそれだけで読む価値がありました（・・・普通の理系なら当たり前の知識かも）。

　私の浅い理解では、機械学習は「モデルの作成 → 最適化問題を解く」という感じなので、これから勉強する上で色々と役に立ちそう・・・と思っています。また、本書のまえがきで述べられているところによると、最適化数学（数理計画）は「経営学やオペレーションズリサーチ（OR）の分野の中心テーマ」ということのようなので、「ある制約の下で、できるだけ最大限の利益を得る（＝ 損失を最小にする）ための最適な解を求める」という考え方は、機械学習に限らず様々な局面で役立ちそうです。

 

・・・っていうか大学で教えて欲しかった

６冊目：プログラミングのための線形代数 

プログラミングのための線形代数

プログラミングのための線形代数

 

 

　さて、「これなら分かる最適化数学」は大変良い本だったのですが、致命的に数学ができない私はこれでもわからない・・・。「今、何で行列のカッコ取れた？」「転置どこいった？」という調子で、とにかく式変形が追えない。これはダメだ・・・ということで線形代数に再入門すべくこちらの本を購入しました。

　amazonの評価が高いだけあって、またしても良書でした！　本書の「はじめに」で述べられているように、

　・「数学のプロではない」読者を対象に

　・「数学のための数学ではなく」

　・「何の役にたつのか」の説明に気を配った

との言葉通り、多数の図を使って線形代数の様々な操作がどういう意味をもつのか説明されていて、私でも「線形代数が何をしたいのか」わかった気になりました。

 

目次は割愛し、重要なメッセージを！

・「行列は写像だ！」

・「行列は写像だ！」

・「行列は写像だ！」

 そして、

・「空間で発想だ！」

・「たいていのものはズームアップすればほとんどまっすぐだ！」

 

以上！ 

まとめ

　コツコツと読み進めるといった計画的な勉強ができないダメ人間なので、長期休暇を活用すべく思い切って６冊も本を購入してしまいました。しかも、どうしても電子書籍に馴染めないという旧世代の人間なので全て紙版・・・。正直、金銭的にかなり辛かったのですが、各所でオススメされている本を中心に購入しただけあって、どの本も大変素晴らしく読んでよかった！と、思える本ばかりでした。

やはり「3日で分かる・・・」とかそんなタイトルにつられて何冊も購入するよりも、評価の高い本を腰を据えて読むのが大事だなあ、と今更ながらに実感した次第です。

 

　当初は「大学で統計や線形代数は一度勉強してるしー」みたいな軽い気持ちだったのですが、思った以上に何も理解していなかっです。大学院まで行ってこの低レベルの私は一体・・・。というよりジョルダン標準形やらラグランジュの未定乗数法やらやった記憶が一切ないんですが、私は本当に理系なのか？？？

・・・

というわけで、本年も雰囲気理系として、ぼちぼちやっていきたいと思いますのでどうぞよろしくお願いいたします。

 

 

 

 

*1:ケモインフォ はこちらみたいな本をがあったら是非読みたい！

*2:腸内細菌叢と疾患の関わりに関する研究も熱いのでマクロは薬学と無関係という意味ではないです・・・念のため

*3:気持ちの悪い発言

RDKit直訳 総括

RDKit直訳

勢いだけで始めた RDKit直訳 Advent Calendar 2018 - Adventar ですが、なんとか日程を全て終えることができました。ご支援いただいた皆様ありがとうございました。

 

ご覧いただいているという嬉しいコメントもいただいたので、調子に乗って試訳をまとめたテストサイトを作成しました。よろしければご参照ください。こんな感じです↓ 

f:id:magattaca:20181230210225p:plain

 

以下、備忘録を兼ねてRDKit直訳アドベントカレンダーの総括を記載しておきます。

 

 1. 振り返り

　当初の予定では、

①「 The RDKit Documentation 」のうち「Getting Started with the RDKit in Python」を対象とし、

② 専門用語についてはtwitterでハッシュタグを使いながら修正していこう

と、考えていました。

 　最終的に、

①については「The RDKit Documentation」のうち、「Python API Reference」を除く主要な項目の試訳を作成することができました。

②については、試訳の作成だけで手いっぱいとなり、SNS活用能力の低さも相まって最初の１週間もしないうちに諦めてしまいました。

　結果として、素人が専門用語もわからずにとりあえず日本語化した、というなんとも歪なものとなってしまいましたが、たたき台にでもしていただければ幸いです。

 

2. 日本語版サイトの追加項目について

　さて、試訳をまとめたサイトを作成するにあたって、日本語版の追加の項目として「Getting Started with the RDKit in Python」の補遺と、翻訳時に訳につまづいた単語の訳語対応表を追加しました。

1) 日本語版補遺を追加した意図

　オンラインドキュメンテーションを訳している間、具体的にどのような操作が行われているのかイマイチわかりませんでした。理由としては

・ 私がPython初心者でコードが読めない

ということに加えて

・SMILES、SMARTSといった表記法に馴染みがないため、操作の意味が想像できない

ということがあると思いました。

そこで、各コードで使われている線形表記を構造式におこしたFigureを作成し、日本語版の補遺として追加することとしました(*1)。

2) 訳語対応表を追加した意図

　用いた訳語の日英対応を掲載した理由は、

・単純に不適切な訳語を使用している可能性が高いため、後で修正する際の参考としたい

ということに加えて

・モジュール、関数の名称が専門用語を反映しているように見えたので、英語の言い回しについて辿れるようにしておけば、RDKitを用いたコードを理解するうえで役に立つのではないか、

と考えたためです。

各日程で訳出に困った単語を集めて並べただけなのですが、まとめてみると

・「非専門家の私が訳出に困る単語」 ＝ 「ケモインフォマティクス の専門用語の可能性が高い」

というわけで、「専門用語集的な役割も果たせるのでは？？」という気もしてきました。　（自分で褒めていくスタイル！！）

 

3. テクニカルな話

　翻訳の記法として、はてなブログが対応しているからという理由でMarkdown記法を採用しJupyter notebookでぼちぼち作成していたのですが、日程最終段階でもっといい方法があるということがわかりました。念のため備忘録。

① エディタは大事

　ふとしたきっかけでAtomというエディタというのを使って見たのですが、こちらのプレビュー機能というのがとても便利でした。こんな感じ↓

f:id:magattaca:20181230193800p:plain

Atomの見た目（左半分：テキスト、右半分：対応箇所のプレビュー）

　左半分のMarkdownを編集すると、右側のプレビューが随時変化していくため、見た目を確認しながら作業を進めることができます。また、コードや段落、リンクといった箇所をそれぞれ色分けして表示してくれるため、記載のおかしい部分を随時把握することができます(*2)。

② reStructuredText 

　試訳をまとめたサイトを作成するにあたってMarkdownからreStructuredTextに記法を変更しました (*3)。

 Markdownよりも優れていると感じた点は、

・脚注の記載方法（リンクが相互なので、行ったり来たりしやすい）

・文書間のリンクを作成しやすい

と、いった点です。一方で

・テーブルの表記方法は字数を合わせないとすぐに崩れてしまうため、全角半角混合なものを書こうとすると幅を合わせるのが難しい

といった点は少し使いづらいと感じました(*4)。

必要に応じて table generator のtext table作成などを活用するのが良いかと思います。

 

　MarkdownからreStructuredTextへの変換にあたってはこちら の記事を参考にPandocを使用しました。とても簡単につかえて便利なのでオススメです。また、AtomでreStructuredTextを使うには拡張（reSTへの対応 と reSTのプレビュー機能 の追加）が必要なのでご注意ください(*5)。

③ Sphinx

　公開するページの作成にあたってはこちらの記事などを参考にさせていただきつつ、Sphinxを使って見ました。Sphinxの日本ユーザ会 のページはとても内容が豊富でわかりやすいのでオススメです。

　またサイトの見た目（theme）にはこちら を参考にRead the Docs「sphinx_rtd_theme」を使用しました。themeを変えるためのconf.pyの設定については岐阜大学藤沢研究室のこちらのページを参考にさせていただきました(*6)。

 

以上、備忘録を兼ねてざっと便利だと思ったソフトウェアや注意点を振り返りました。

 

4. まとめ

　思いつきで始めた翻訳で思った以上に大変だったのですが、RDKitだけでなく翻訳や技術文書に関する周辺についても色々と知ることができとても勉強になりました。また、必要に迫られないと調べないモチベーションの低い人間なので、アドベントカレンダーを利用するというのは、逃げ道を断つというのと、強制的にスケジュール管理せざるを得ないという点で、とても有効だと思いました。

・・・残念ながら翻訳に必死で肝心のRDKitのPythonのコードを読む余裕はなかった 

 

　ふとしたきっかけでこのブログを始めてから、今年はとてもたくさんの方にお世話になりました。この場を借りてお礼申し上げます。

　今更ながらにとても驚いたのがツイッターにはすごい方々がたくさんいらっしゃって、情報をたくさん発信してくださっている、というのが今年の１番の発見でした。田舎だと情報手に入らないなーとぼんやり過ごしていたのですが、まさかこんな手段があったとは！ 来年も気が向いたらぼちぼちやっていこうと思いますので、皆様どうぞよろしくお願い申し上げます。

 

良いお年を！

 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

*1)  構造式の作成にMarvinSketchを使って、RDKitを使っていない時点で「１ヶ月やって何も学んでない！」と怒られてしまいそうですが、MarvinSketch自由度高くて使い勝手良さすぎた・・・許してください

 *2)  初日に知りたかった！！Jupyter notebookで「Shift+ Enter」「セルの追加」を連打し続けた日々を返して欲しい・・・

*3) RDKitオリジナルサイトはMarkdownのものとreSTのものがどちらもありましたが、「技術書、翻訳」で検索したらreSTをオススメする記事がたくさん出てきたのでこちらに統一しました。

*4) RDKitのオリジナルサイトを流用したため表記が難しかったですが、

csv形式など幅を気にしなくてもいい表記方法もあるみたいです。

*5) reST表記を理解していなかったため変換ミスと思って削除した「\」が、実はエスケープ文字で必須だったと後で気付いた時は発狂しそうになりました・・・思い込みで作業するのは良くない

*6) Sphinxの使い方以外にも生物学的対称性の説明など面白そうな記事がいっぱいありました。

      研究室のHP →藤澤研究室

 

 

Backwards incompatible changes 〜RDKit 直訳 Day25〜

RDKit直訳

こちらはRDKit直訳 Advent Calendar 2018 - Adventar 25日目の記事です。基本的な進め方は1日目の記事をご覧ください。

本日の訳出に困った用語
backward compatibility: 後方互換性
break compatibility: 互換性を破る
break down: 分類する

以下、訳

後方互換性の無い変更（Backwards incompatible changes）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#backwards-incompatible-changes

リリース間の後方互換性を維持しようと本当に努めていますが、互換性を破る必要がある事例がまれにあります。このドキュメントの目的は、RDKitで行われた後方互換性の無い変更についての情報を補足することです。リリースサイクルによって分類されており、だいたいにおいてバグの修正で生じた結果による変更は含んでいません; バグの修正はリリースノートで注意喚起するようにしています。

Release 2018.03

[LInk] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#release-2018-03

デフォルトではMolToSmiles()はisomeric SMILESを生成します（MolToSmiles() generates isomeric SMILES by default）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#moltosmiles-generates-isomeric-smiles-by-default

以前のリリースでは、出力されるSMILESに立体化学や同位体ラベルについての情報を含めたい場合、オプションのisomericSmiles引数をtrueに設定する必要がありました。現在、このデフォルトの値はtrueとなっています。昔の動作に戻し、立体化学の情報を持たないSMILESを取得したい場合は、isomericSmilesをfalseに設定するだけです。

MolToMolBlock()はincludeStereoフラグがセットされると2Dコンフォメーションを生成する（MolToMolBlock() generates a 2D conformation when the includeStereo flag is set）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#moltomolblock-generates-a-2d-conformation-when-the-includestereo-flag-is-set

Mol blockの立体化学を取得したい場合、出力が座標の情報を有する必要があります。以前のRDKitのバージョンでは、座標の生成を忘れずに自分で行う必要がありました。現在では、includeStereoがセットされている場合、コンフォメーションを持たない分子に対してはデフォルトで座標が生成されます。

現在のデフォルトではコンフォメーション生成コードはETKDGを使う（The conformation generation code now uses ETKDG by default）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#the-conformation-generation-code-now-uses-etkdg-by-default

以前のRDKitのリリースでは、デフォルトでは、標準的なデジスタンスジオメトリー法を使ってコンフォメーションを生成していました。新しいデフォルトの設定では、Sereina RinikerによるETKDGアルゴリズムを使います。少し遅いですが、ずっと良い結果を出すことが示されています。

Release 2018.09

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#release-2018-09

デフォルトではGetAtomSmile()はisomeric SMILESを生成する（GetAtomSmiles() generates isomeric SMILES by default）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#getatomsmiles-generates-isomeric-smiles-by-default

以前のリリースでは、出力されるSMILESに立体化学や同位体ラベルについての情報を含めたい場合、オプションのisomericSmiles引数をtrueに設定する必要がありました。現在、このデフォルトの値はtrueとなっています。昔の動作に戻し、立体化学の情報を持たないSMILESを取得したい場合は、isomericSmilesをfalseに設定するだけです。

MCSコードのringMatchesRingOnlyオプションの変更点（Changes to ringMatchesRingOnly option in the MCS code）

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#changes-to-ringmatchesringonly-option-in-the-mcs-code

FindMCS()関数のringMatchesRingOnlyオプションは結合と結合のマッチングと同様に、原子と原子のマッチングにも適用されます。

現在のデフォルトではPythonから呼び出すとコンフォメーション生成コードはETKDGを使う（The conformation generation code now uses ETKDG by default when called from Python）

[Link]

現在デフォルトでは、Python関数のEmbedMolecule()とEmbedMultipleConfs()は標準的なデジスタンスジオメトリーの代わりにETKDGアルゴリズムを使用します。

License

[Link] https://www.rdkit.org/docs/BackwardsIncompatibleChanges.html#license

この文書の著作権は copyright (C) 2013-2018 by Greg Landrum に所属しています。

この文書はCreative Commons Attribution-ShareAlike 4.0 Licenseのもとでライセンスされています。このライセンスを見るためにはhttp://creativecommons.org/licenses/by-sa/4.0/ にアクセスするか、Creative Commons, 543 Howard Street, 5th Floor, San Francisco, California, 94105, USA.　に手紙を送ってください。

このライセンスの意図はRDKitそのものの意図と似ています。簡単に言えば "これを使ってなんでもやりたいことをやっていいけど、私たちの功績にも言及してね”

12/25/2018

このブログ記事のライセンス

以上はRDKit Documentationの試訳です。
ライセンスはCC BY-SA 4.0で、Greg Landrum氏の功績の元に作成しています。
Creative Commons — Attribution-ShareAlike 4.0 International — CC BY-SA 4.0