magattacaのブログ

日付以外誤報

活性化合物群を分析して指標の閾値を決めようとする話

前回、約13万個の化合物を含むEnamine_Premium_collectionから、分子量300を基準とした分子の絞り込みにより約11万個まで化合物数を減らすことができました。

まだまだ多い・・・と思っていたら@iwatobipen先生に分子を絞り込む指標を考える上でとても参考になる記事 SAR News Np.19 をご紹介いただきました。

SAR newsは日本薬学会 構造活性相関部会のHPの発行しているニュースレターで、だいたい年に２本のペースで発行されており、部会に所属していなくても過去のニュースレターを含めて自由に閲覧可能となっています（バックナンバーはこちらから閲覧可能です）。

SAR Newsで用いられている指標

さて、ご紹介いただいたニュースレターの 株式会社理論創薬研究所 吉森篤史 氏による記事では「標的タンパク質指向型フラグメントライブラリの構築」を適用例として、RDKitを用いたフィルタリングの方法が実際のコード含めて解説されています。

用いられているフィルタリング基準は、フラグメント指標「Rule of Three」で、具体的には以下となっています。

指標 基準 記事中のコード 現在のコード？
① 分子量 ≦300 AvailDescriptors.descDict['MolWt'](mol) Descriptors.MolWt(mol)
② cLogP ≦3 AvailDescriptors.descDict['MolLogP'](mol) Descriptors.MolLogP(mol)
③ 水素結合供与体数 ≦3 AvailDescriptors.descDict['NumHDonors'](mol) Descriptors.NumHDonors(mol)
④ 水素結合受容体数 ≦3 AvailDescriptors.descDict['NumHAcceptors'](mol) Descriptors.NumHAcceptors(mol)
⑤ 回転可能結合数 ≦3 AvailDescriptors.descDict['NumRotatableBonds'](mol) Descriptors.NumRotatableBonds(mol)
⑥ 極性表面積 ≦60 AvailDescriptors.descDict['TPSA'](mol) Descriptors.TPSA(mol)

記事が2010年のもののためか、現在のRDKitのモジュールを眺めてもAvailDescriptorsというものが見つけられませんでした。「化学の新しいカタチ」さんの記事「 RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ」 に同じ目的と思われるモジュール（Descriptors）が紹介されていましたので、こちらも上記の表に追加しました。

上記、指標の基準はフラグメント（"活性が弱いながらも適切な結合を形成できる小分子"）を対象としているためか、かなり小さな値となっているように見えます。今回の創薬レイドバトルの目標は、ある程度活性が強いものを取得することなので、この基準をそのまま適用するのは良くなさそうです。

また、他の指標としてリピンスキーの法則を用いる、というのも良さそうですが、前回活性化合物群の分子量のデータを眺めた印象では、分子量 500以上の大きなものも多く、すでにRule of Five を外れてしまっています。

一般的な指標をそのまま適用するだけでは、既知活性化合物群ですら外れてしまいそうなので、まずは活性化合物群の値がどのようになっているかを検証したいと思います。

PD/PD-L1阻害活性　化合物群の指標の分析

指標の計算

まずは必要なものをimportします。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools
import pandas as pd
 
print('rdkit version: ',rdBase.rdkitVersion)
print('pandas version: ', pd.__version__)
#rdkit version:  2018.09.1
#pandas version:  0.23.4

前回同様、創薬レイドバトルの共闘用シェアデータを使用します。

#csv形式のデータをPandasのDataFrameで読み込む
df = pd.read_csv('../pd1_inhibitor_dataset-master/PD1_inhibitor_dataset.csv')

#smilesカラムのSMILES情報を使ってをMolオブジェクトを追加する
PandasTools.AddMoleculeColumnToFrame(df, 'smiles')

#ROMolカラムに追加されているMolオブジェクトを使って記述子を計算し、値を追加する。
#①分子量
df['MW'] = df.ROMol.map(Descriptors.MolWt)
#②LogP
df['MolLogP'] = df.ROMol.map(Descriptors.MolLogP)
#③水素結合供与体数
df['NumHDonors'] = df.ROMol.map(Descriptors.NumHDonors)
#④水素結合受容体数
df['NumHAcceptors'] = df.ROMol.map(Descriptors.NumHAcceptors)
#⑤回転可能結合数
df['NumRotatableBonds'] = df.ROMol.map(Descriptors.NumRotatableBonds)
#⑥極性表面積
df['TPSA'] = df.ROMol.map(Descriptors.TPSA)

計算した指標の要約統計量を求めてみます。

#計算した指標のリスト
Descriptors_list = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA']

#要約統計量を計算。小数点以下の桁数が多いので丸める
df[Descriptors_list].describe().round(2)

f:id:magattaca:20190120180844p:plain

相関係数は以下のようになっています。

#相関係数の計算
df[Descriptors_list].corr().round(2)

f:id:magattaca:20190120180932p:plain

ざっと見る限りMolLogPのみが他の指標との相関係数が負かつ絶対値が小さくなっています。

指標をプロット

数字を見てもイマイチ感触がつかめません。前回分子量を眺めた際に描画すると値の分布がわかりやすかったので今回も描画を行います。一つ一つみるのは面倒なのでSeabornを使ってペアプロットを作成します。

import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline

#Seabornを用いたペアプロットの作成
sns.pairplot(df[Descriptors_list])

f:id:magattaca:20190120181042p:plain

確かに、MolLogPだけが他の指標とのプロットした際に相関が悪そうです。

見やすさのためMolLogPを外して再度プロットして見ます。

#計算した指標のリストから MolLogPを除いたもの
Descriptors_list_2 = ['MW', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA']
sns.pairplot(df[Descriptors_list_2])

f:id:magattaca:20190120181154p:plain

残りの指標間は直線上に乗っているようにみえわかりやすくなりました。

外れ値を取り除く

各プロットを見ると、最大値となっている分子が１つだけ他の分子から遠く、外れ値のようになっています。特にNumRotatableBondsで顕著にみえるので、こちらの最大値がどのような分子なのか眺めて見たいと思います。

#NumRotatableBondsが最大値となる化合物を眺める。
#まずは最大値となる化合物の行番号を取得する
print(df['NumRotatableBonds'].idxmax())
#5と出力

#Molオブジェクトの格納されているcolumnの番号を確認
print(df.columns.get_loc('ROMol'))
#6と出力

DatafFrameの(5, 6) がみたい構造のMOlオブジェクトがある場所のようです。これを描画します。描画の手順は省きます。

f:id:magattaca:20190120181634p:plain

環状ペプチドにさらに長い鎖状構造がつながっています。他の構造と比べてかなり特異な構造なので、こちらは外れ値として省いてしまっても良さそうです。

新しいDataFrameを作成し、再度ペアプロットを作成します。

#行番号5を除去する
df_new = df.drop(5, axis=0)
sns.pairplot(df_new[Descriptors_list_2])

f:id:magattaca:20190120181824p:plain

かなり分布がわかりやすくなってきました。いずれの指標においても大きく２つのグループに分子を分けることができそうです。

これは・・・クラスタリングというやつが使えるのではないか？？？

クラスタリングによる分類

Morganフィンガープリントとタニモト係数による類似性評価

「化学の新しいカタチ」さんのこちらの記事（RDKitとk平均法による化合物の非階層型クラスタリング ）を参考に化合物を2つのクラスターに分けたいと思います。以下の操作は上記記事の処理手順をほぼそのまま使用させていただいています。

各分子のMorganフィンガープリントを計算し、分子間の類似性をタニモト係数によって比較、クラスタリングに用いるための距離行列を作成します。

from rdkit import DataStructs
from rdkit.Chem import AllChem
import numpy as np

#計算に使用するためMolオブジェクトのリストを作成
mol_list = list(df_new['ROMol'])
print(len(mol_list))
#52個の分子を含む

#Morganフィンガープリントの計算
morgan_fp = [AllChem.GetMorganFingerprintAsBitVect(x,2,2048) for x in mol_list]

#総数52の分子間で距離行列の計算
dis_matrix = [DataStructs.BulkTanimotoSimilarity(morgan_fp[i], morgan_fp[:52],returnDistance=True) for i in range(52)]

### numpy.ndarrayへの変換
dis_array = np.array(dis_matrix)
dis_array.shape

52 x 52 の距離行列が得られました。

k平均法によるクラスタリング

k平均法で２つのクラスターに分類します。

from sklearn.cluster import KMeans
import pandas as pd
kmeans = KMeans(n_clusters=2)
kmeans.fit(dis_array)
pd.value_counts(kmeans.labels_)

f:id:magattaca:20190120182815p:plain

10個の分子を含むクラス0と、42個の分子を含むクラス1に分類できました。 このクラス分類もDataFrameに追加しておきます。

class_list = []
for i in range(len(kmeans.labels_)):
    if kmeans.labels_[i] == 0:
        class_list.append('class_0')
    elif kmeans.labels_[i] == 1:
        class_list.append('class_1')
        
df_new['kmeans_class'] = class_list

ペアプロットへのクラス分類の反映

先に作成したpairplotをkmeans法で分類したクラスに応じて色分けするとどうなるか、見て見たいと思います。

もともと他の指標との相関が悪そうだったMolLogPについても、色分けしたらどうなるか気になるので、今回はpairplotに加えたいと思います。

#計算した指標のリスト（MolLogPを含む）にkmeansで形成したクラスも加える
Descriptors_class = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA', 'kmeans_class']
sns.pairplot(df_new[Descriptors_class], hue='kmeans_class')

f:id:magattaca:20190120183026p:plain

Morganフィンガープリントをつかって２つのクラスターに分けた結果で色分けすると、指標の分布でも比較的良好に２つのクラスターに分かれているように見えます。

とくに興味深いのが、MolLogPで、ヒストグラムでは２つのクラスターが重なっていますが、ペアプロット上では割と良好に分布が分かれているように見えます。

Morganフィンガープリントがどのように計算しているのか、まだ理解していないのですが指標の分布をうまく説明するようなクラスタリングができているのは面白い結果です。

各クラスターに含まれる構造のインタラクティブな描画

プロット上での描画

各点の構造がどのようなものか、分布の広がりがありつつ２つのクラスターが綺麗に分かれているMolLogPとTPSAの２軸によるプロットを用いて可視化して見たいと思います。

@iwatobipen先生の記事 Plot Chemical space with d3js based libraryを参考にさせていただきました。

import mpld3  #別途インストールしておかないとimportできない
from mpld3 import plugins
mpld3.enable_notebook()

#描画を作成するためのモジュール
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

#svgを格納するための空のリスト
svgs_list = []


for i in range(len(df_new)):
    view = rdMolDraw2D.MolDraw2DSVG(300,300)
    
    #Molオブジェクトを含むカラムの番号
    ROMol_col = df_new.columns.get_loc('ROMol')
    
    #i行目のMolオブジェクトを取得
    mol = df_new.iloc[i, ROMol_col]
    
    #描画する分子の前処理
    prepro_mol = rdMolDraw2D.PrepareMolForDrawing(mol)
    
    #分子をコンテナに追加
    view.DrawMolecule(prepro_mol)
    
    #コンテナをファイナライズ
    view.FinishDrawing()
    
    #コンテナに書き込んだデータを取り出す
    svg = view.GetDrawingText()
    
    #データをリストに追加
    svgs_list.append(svg.replace('svg:,',''))

以上で、見たい化合物の構造情報（svg形式）を含むリストが作成できました。あとは、プロットと組み合わせて、プロット上の各点がどのような構造に相当するのかを眺めます。

fig, ax = plt.subplots()
ax.set_xlabel('MolLogP')
ax.set_ylabel('TPSA')
points = ax.scatter(df_new['MolLogP'], df_new['TPSA'])

tooltip = plugins.PointHTMLTooltip(points, svgs_list)
plugins.connect(fig, tooltip)

f:id:magattaca:20190120183609g:plain

y軸（TPSA）値が 500よりも上に位置するのがclass_1、500よりも下に位置するのがclass_0となっています。

カーソルをうまくGIFに写すことができずわかりづらくなってしまいましたが、プロット上を右上から反時計回りに順にカーソルを移動させ、構造を表示させています。 class_1は全て環状ペプチド、class_0はほとんどが環状ペプチドではない低分子ですが、左下に位置する３つは環状ペプチド構造となっています。

スライドバー形式での描画

少しわかりづらかったので、もうひとつのインタラクティブな描画方法を試します。こちらはRDKitのBlogで紹介されていたものです。

まずは　class_0 のみを要素とするDataFrameを作成し、さらにMolオブジェクトのリストとします。

df_class_0 = df_new[df_new['kmeans_class']== 'class_0']
class_0_list = list(df_class_0['ROMol'])

以下の方法ではipywidgetsを使って可視化します。

from ipywidgets import interact,fixed,IntSlider
import ipywidgets

#class_0 に分類された構造を順番に眺めていく
def class_0_viewer(idx):
    mol = class_0_list[idx]
    return(display(Draw.MolToImage(mol)))
    
interact(class_0_viewer, idx=ipywidgets.IntSlider(min=0,max=len(class_0_list)-1, step=1));

f:id:magattaca:20190120184006g:plain

class_0 のうち idx = [29, 30, 32] は環状ペプチド構造となっており、その他の構造とは少し違った形となっています。

ついでにclass_1もやってみます。同様の操作なのでコードは省略します。

f:id:magattaca:20190120184131g:plain

順番に構造をおくるとかなりわかりやすくなりました。

環状ペプチドと低分子とで２つに再度クラス分け

環状ペプチドと、それ以外で分けた方がより良いクラス分けとなりそうです。 そこでclass_0 に含まれる環状構造 をclass_1に移動し、新しいデータフレームを２つ作成します。

df_peptides = pd.concat([df_class_1, df_class_0.iloc[[29,30,32]]])
df_small_molecules = df_class_0.drop([40,41,43])  # idx=[29, 30, 32]番目の行のindexは[40,41,43]

print(len(df_peptides))  # 13と出力
print(len(df_small_molecules)) # 39と出力

合計52（13 + 39）個の２つのデータフレームに分けることができました。　処理の最初の段階で、外れ値とも思える大きな分子を一つ取り除いているので、もともとのデータセット（53個）よりも一つ少なくなっています。

活性化合物群に基づき指標の閾値を考察

さて、いろいろと操作してきたため、本来の目的を見失ってきました。

もともとの目的は、化合物ライブラリから化合物の数を絞り込むための指標の閾値を、活性があると報告されている化合物群のデータをもとに決めよう、というものでした。

活性化合物群を可視化した結果、２つのクラスターに分類することができ、これをもとに環状ペプチドと、それ以外の低分子との２つに分けることができました。このうち、「PD/PD-L1の相互作用を阻害する低分子を見つける」という目標として参照するのに適しているのは、後者となりそうです。

以下のテーブルでは、低分子化合物クラスター(small_molecules)と環状ペプチドクラスター（peptides）における各指標の最小値と最大値、およびLipinski の Rule of Five の値を並べています。

化合物を絞り込む際に用いる指標の閾値としては、値の範囲よりも少し広めにとって最右列のフィルタリング目安とするのが手かもしれません。

指標 small_molecules
min		 small_molecules
max		 peptides
min		 peptides
max		 Lipinski
Rule of 5 		フィルタリング目安
① 分子量 MW 215 835 491 2096 <= 500 > 300
（実施済）
② cLogP MolLogP -2.5 9.5 -6.35 3.20 <= 5 -3 ~ 10
③ 水素結合供与体数 NumHDonoes 1 7 9 22 <= 5 < 10
④ 水素結合受容体数 NumHAcceptors 3 9 10 26 <= 10 < 10
⑤ 回転可能結合数 NumRotatableBonds 2 15 3 39 < 15
⑥ 極性表面積 TPSA 48 250 237 713 < 250

念のため今回クラス分類した２つのデータフレームを出力しておきます。

properties = df_small_molecules.columns.values
Chem.PandasTools.WriteSDF(df_small_molecules, 'small_molecules.sdf', molColName='ROMol', idName='compound_id', properties=properties)
Chem.PandasTools.WriteSDF(df_peptides, 'peptides.sdf', molColName='ROMol', idName='compound_id', properties=properties)

まとめ

以上、本記事では分子量以外の化合物絞り込みの基準を探すべく、活性化合物群の情報を分析して見ました。

Morganフィンガープリントとk平均法を組み合わせた分類では、環状ペプチドのような構造もいくつかは鎖状低分子と同一のクラスにふくまれるなど、想定外のクラス分けがされて興味深い結果になりました。 また、ある程度構造が限られた化合物セットであれば順番に表示させるような機能を利用することで、化合物セットの雰囲気をつかむことができることがわかりました。

本記事を作成するにあたって、さまざまなページを参照させていただきました。また、多くのコードを拝借させていただいております。誤り等ございましたらご指摘いただければ幸いです。

追記

前回の記事で分子量で絞り込んだあとの化合物群（約11万個）の指標を計算してみたら、上記のテーブルのフィルタリング目安ではほとんど絞り込みにならなさそうです。 もっとちゃんと考えないとダメそうです・・・私の休日・・・

f:id:magattaca:20190121000021p:plain

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

参考記事

① SAR News No.19
Chemoinformatics Toolkits を用いた創薬システム開発におけるラピッドラタイピング　吉森　篤史
http://bukai.pharm.or.jp/bukai_kozo/SARNews/SARNews_19.pdf

② 化学の新しいカタチ
RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ

③ 化学の新しいカタチ
RDKitとk平均法による化合物の非階層型クラスタリング

④ IS LIFE WORTH LIVING? Plot Chemical space with d3js based library #RDKit #Chemoinformatics | Is life worth living?

⑤ The RDKit Blog
RDKit: Using the new fingerprint bit rendering code

⑥ ipywidgetsのインストール方法
ipywidgetsの使い方

⑦ ipywidgets の User Guide
ipywidgets — Jupyter Widgets 7.4.2 documentation

分子量で絞り込もうとした話

なんとか化合物を読み込むということはできましたが、総数は200万以上・・・

このままでは全くどうして良いかわかりません。

どうしたものかと思っていたらTwitterで処理の方法を紹介してくださっている方がいらっしゃいました。

@mayudqxさんのされていた処理の順序は、

処理 意図
① 分子量の範囲を設定 粗い絞り込み
② 部分構造で絞り込み 活性化合物を参照したLBVS（？）
③ 粗いドッキング タンパク質との相互作用を予測するSBVS（？）
④ ファーマコフォアを指定したより詳細なドッキング SBVS ２回目（？）

という感じです。

このうち①と②なら、RDKitを使って行うことができそうです。とりあえずプロの真似をしよう！

今回は、まずは①分子量での絞り込みに挑戦したいと思います。

絞りこみの対象となる化合物群の構造式をいくつか眺めた感想として、確かにとても小さな分子（試薬？？）もライブラリに含まれていることが意外でした。今回の目標は、生理活性のある分子を予測しよう、というものなので、@mayudqxさんがされているように、あまりにも小さく活性の期待できなさそうなものを省くというのは、最初に行う絞り込みとして良さそうです。

ところで、絞り込みの閾値を分子量300とされていましたが、こういう値をどうやって設定して良いものかわかりません。

何か参考になるデータはないか・・・と思っていたら創薬レイドバトルのgithubページ に「共闘用シェアデータ」という、活性化合物群のデータを掲載してくださっていました！

様々な特許から「低分子PD-1阻害剤と思われる 低分子化合物の構造を抽出し」たもので、「各特許の記述に基づいて最も高活性とされる化合物を1例ずつ抽出し」たものだそうです。

こんな大変そうなデータセットを公開してくださるとは・・・。ありがとうございます！

活性のある化合物は閾値を決める根拠としてとても良さそうです。これを使わない手はない！ということでこちらを参照して見たいと思います。

前回に引き続き、RDKitのPandasToolsを使って解析を行って見たいと思います。「化学の新しいカタチ」さんのこちらの記事 を参考にさせていただきました。

まずは色々必要なものをimportします。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools
import pandas as pd
 
print('rdkit version: ',rdBase.rdkitVersion)
print('pandas version: ', pd.__version__)

#rdkit version:  2018.09.1
#pandas version:  0.23.4

データはcsv形式で提供されているので、まずはpandasで読み込みます。

githubからダウンロードしてくるとpd1_inhibitor_dataset-masterというフォルダの中にPD1_inhibitor_dataset.csvという名前でデータが保存されていました。

df = pd.read_csv('../pd1_inhibitor_dataset-master/PD1_inhibitor_dataset.csv')
print(df.shape)
#(53, 6)

print(df.columns)
#Index(['compound_id', 'patent_no', 'example_no', 'schembl_id', 'applicant', 'smiles'], dtype='object')

f:id:magattaca:20190119223825p:plain

全部で53個の化合物が含まれており、構造情報はSMILESとして格納されています。

SMILESのままではよくわからないので、構造式を眺めるためにMolオブジェクトに変換します。

PandasTools.AddMoleculeColumnToFrame(df, 'smiles')

f:id:magattaca:20190119223956p:plain

ROMolという列にMolオブジェクトが保存されており、小さくて見づらいですが構造式として表示されます。 最初の8個を構造式だけを取り出して並べて見ます。

PandasTools.FrameToGridImage(df[:8], column='ROMol', legendsCol='compound_id', molsPerRow=4, subImgSize=(300,300))

f:id:magattaca:20190119224146p:plain

BMS-##というcompoun_idのついた化合物はマクロサイクル（環状ペプチド）のような分子も含まれているようです。このあたりは分子量が大きそうです。

今回知りたいのは、これらの活性化合物の分子量がどのような値となっているか？、なので分子量を計算しDataframeに加えます。

df['MW'] = df.ROMol.map(Descriptors.MolWt)

f:id:magattaca:20190119224433p:plain

分子量の分布がどのようになっているか知りたいので、まずはpandasのdescribe()で要約統計量を取得します。

df['MW'].describe()

'''
count      53.000000
mean      819.726075
std       626.128326
min       215.213000
25%       419.525000
50%       574.899000
75%       753.340000
max      2699.165000
Name: MW, dtype: float64
'''

平均値 820 で最大値 2699 ！　予想以上に分子量が大きいです。分子量の大きそうな巨大なマクロサイクルで、分子量の統計量が大きな値へと引っ張られていそうです。

要約統計量だけから判断すると、大きな値に引きづられてしまいそうなので、念のためヒストグラムを作成し、分子量の値の広がりを確認したいと思います。

import matplotlib.pyplot as plt

#notebookにグラフを描くためのマジックコマンド
%matplotlib inline

#MWのMaxが約2700なので 0 から 2800 までの幅50のヒストグラムを作成
edges = range(0,2800,50)
n, bins, patches = plt.hist(df['MW'], bins=edges)

f:id:magattaca:20190119225412p:plain

分子量500前後と分子量2000を超える分子の大きく２つに山が分かれています。 今回の創薬レイドバトルの目標は低分子のPD/PD-L1の阻害剤を見つけるということですから、分子量500前後の山に相当するような分子を見つけ出すことが目的になりそうです。（そもそも分子量2000を超える分子がライブラリに含まれているか疑問・・・）

分子量が小さい山から、分子量による絞り込みの閾値を決めたいと思います。より具体的な値をみるため各ビンに含まれる数を取り出して見ます。

hist_df= pd.DataFrame(n, bins[1:]).T

f:id:magattaca:20190119230452p:plain

分子量300以下は一つしかないため、300を足切りにして絞り込むというのは非常に妥当な値のように見えます。閾値を一瞬で見定めるとは・・・プロすごい・・・

それでは、スクリーニング対象の化合物を検証したいと思います。

データサイズが大きいのでまたしてもGoogle colaboratoryを利用したいと思います。

ちょうど以前の記事で前処理の際に分子量を計算し、圧縮ファイル（Enamine_Premium_collection_pro_id_MW_sdf.gz）としてGoogleドライブに保存していたので、こちらを利用します。

もろもろの設定は省略し、SDFを読み込むところから始めます。

from rdkit import Chem
import gzip
from rdkit.Chem import Descriptors

EPc_gz = gzip.open('Enamine_Premium_collection_pro_id_MW_sdf.gz')

#ForwardSDMolSupplierで読み込む。
EPc_suppl = Chem.ForwardSDMolSupplier(EPc_gz)  

#molオブジェクトのリストを作る段階で、分子量を取得して分子量のみのリストを作成する。
#空のリストを作成
EPc_MW_list = []
EPc_mols = []

#ループを回せ!!!
for x in EPc_suppl:
  if x is not None:
    mol = x
  
  #MolオブジェクトにMWという名前で格納している分子量の値（プロパティ）を取得し、リストに追加する。
  #小数点を含むfloat型のプロパティなので GetDoubleProp を使う
  mw = mol.GetDoubleProp('MW')
  
  EPc_MW_list.append(mw)
 
  #Molオブジェクトのリストも作成
  EPc_mols.append(mol)
  
#総数を確認
len(EPc_MW_list)

# 128816と出力された

分子量のリスト（EPc_MW_list）が作成できたので、pandasで処理します。

import pandas as pd
df_MW = pd.DataFrame(EPc_MW_list, columns=['MW'])
df_MW['MW'].describe()

"""
count    128816.000000
mean        332.064740
std          32.450643
min         204.277000
25%         312.377000
50%         332.426000
75%         352.398000
max         527.150000
Name: MW, dtype: float64
"""

平均 332 で 204 から 527　の範囲の分子量となっているようです。

ヒストグラムを描いてみます。

import matplotlib.pyplot as plt
%matplotlib inline

#分子量 100 から 600 までの幅50のヒストグラムを作成
edges = range(100,600,50)
n, bins, patches = plt.hist(df_MW['MW'], bins=edges)

f:id:magattaca:20190119230828p:plain

分子量が300を超えるものについてのみのMolオブジェクトのリストを作成します。

EPc_mols_300 = []

for mol in EPc_mols:
    MW = mol.GetDoubleProp('MW')
    
    if MW > 300:
        EPc_mols_300.append(mol) 

len(EPc_mols_300)

#109602

残った分子の総数は約11万個。元々のEnamine_Premium_collectionには約13万の分子が含まれていたので、とりあえず２万個程度の分子を減らすことができたようです！

粗い絞り込みでまだまだ残った分子の数は多いですが、最初のフィルタリングとしてはそれなりにうまくいったのではないでしょうか。

sdfとして出力したいと思いますが、前回気づかぬ間にMolオブジェクトのプロパティが失われているというミスを犯したので、まずはプロパティを確認します。

#今あるプロパティのリスト
prop_names = [p for p in EPc_mols_300[0].GetPropNames()]
print(prop_names)
# ['idnumber', 'original_id', 'MW']

きちんとプロパティが保持されていそうです。安心してSDFを作成できます。

#構造と各種プロパティを残したsdfファイルを作成
writer = Chem.SDWriter('Enamine_Premium_collection_MW_filter.sdf')

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

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

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

以上、最初のステップとして分子量による絞り込みでした。

まだまだ分子数多い・・・しかもこれより多いデータセットがまだ３つある・・・

締め切りに間に合うのでしょうか？？？っていうか締め切りいつ？？？

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:気持ちの悪い発言