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

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

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

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

計算するのに化合物が多すぎたので
１．分子量300以下は削除
２．オルト位に置換基をもつビフェニルに限定
３．それでも多かったのでざーっとドッキングして-8でカット
４．ファーマコフォアをビフェニルに指定して試行回数増やして再計算

— mayudqx (@mayudqx) 2019年1月15日

@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')

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

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

PandasTools.AddMoleculeColumnToFrame(df, 'smiles')

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

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

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

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

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

分子量の分布がどのようになっているか知りたいので、まずは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)

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

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

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

分子量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)

分子量が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()

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

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

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