ライブラリの指標計算結果まとめ

前回の記事で、活性化合物群の情報をもとにライブラリから候補化合物を絞り込むため、指標の閾値（フィルター）を決めようと試みました。とりあえずの結果として、以下の値を求めました。

	指標		small_molecules min	small_molecules max	フィルタリング目安
①	分子量	MW	215	835	> 300
②	cLogP	MolLogP	-2.5	9.5	-3 ~ 10
③	水素結合供与体数	NumHDonoes	1	7	< 10
④	水素結合受容体数	NumHAcceptors	3	9	< 10
⑤	回転可能結合数	NumRotatableBonds	2	15	< 15
⑥	極性表面積	TPSA	48	250	< 250

こちらのフィルタリング目安を、前処理済みのライブラリEnamine_Premium_collection（分子量300より小さいものは除去済み）に適用すべく、各化合物の指標を計算、その統計値を求めました。

結果は以下の通り・・・

f:id:magattaca:20190126163131p:plain

なんと、欲張って幅を広くとりすぎたためか、分子量以外の指標では化合物の数をうまく減らせそうにありません。

絞り込み対象のライブラリにおける指標の値の分布を検証せずに、手抜きして活性化合物群のみから閾値を決めようとしたことに原因があるかもしれません。

まずは、相手を知る必要がありそうです。

そこで各指標値の値を残りのライブラリ３つ（Enamine_Advanced_collection、Enamine_HTS_collection、UOS_HTS）に対しても求め分布を把握してみることにしました。これらのライブラリに対しては分子量 300以下を除く処理はしていません。

まずは結果から・・・

① Enamine_Advanced_collection

f:id:magattaca:20190126163312p:plain

② Enamine_HTS_collection

f:id:magattaca:20190126163415p:plain

③ UOS_HTS

f:id:magattaca:20190126163429p:plain

分子量300以下を除く処理をしていないためか、分子量の最小値の小ささにすこしびっくりしました。

このままではよくわからないので、値のばらつきを見るため箱ひげ図を用いてばらつきを可視化しました。以下は、Enamine_Advanced_collectionの場合です。

f:id:magattaca:20190126163614p:plain

こちらを見ると先日設定したフィルタリング目安は、最大値よりも大きく、外れ値とも言えそうな化合物も全て残してしまう、ということがわかります。

先にこっちを見ておくべきだった！！

以下に、上記を求めた際につまづいた点などを記載しておきます。

MolLogPの計算
前処理のコード全体
処理に失敗した分子の確認
SDFの持つ情報をとりあえず確認する方法
統計量の計算
箱ひげ図を用いたばらつきの可視化

MolLogPの計算

Enamine_Premium_collectionを計算する際には問題とならなかったのですが、Enamine_Advanced_collectionを処理する際に、MolLogPの計算でエラーが生じました。指標の要約統計量を確認するとEnamine_Advanced_collectionのMolLogPのみ、count(値の総数:486321)が分子の総数(486322)よりも1小さくなっています。

エラーの内容は以下の通り。

---> 42   MolLogP_value = Chem.Descriptors.MolLogP(mol_neu)
     43   mol_neu.SetDoubleProp('MolLogP', MolLogP_value)
     44   MolLogP_list.append(MolLogP_value)

/root/miniconda/lib/python3.6/site-packages/rdkit/Chem/Crippen.py in <lambda>(*x, **y)
    168 _pyMolMR.version = "1.1.0"
    169 
--> 170 MolLogP = lambda *x, **y: rdMolDescriptors.CalcCrippenDescriptors(*x, **y)[0]
    171 MolLogP.version = rdMolDescriptors._CalcCrippenDescriptors_version
    172 MolLogP.__doc__ = """ Wildman-Crippen LogP value

ValueError: Sanitization error: Explicit valence for atom # 6 B, 5, is greater than permitted

Value Error と書かれたエラー内容から推測すると、この分子はホウ素（B : 原子番号 6）を含み、その原子に対する結合の数が、通常よりも多いため分子としてうまく認識できなくなっているようです。構造は以下の通り。

f:id:magattaca:20190126163852p:plain

上記のようにホウ素に3つのフッ素原子とピロールが結合した分子のカリウム塩となっています。４本の結合＋負電荷で「valence # 6 B, 5」となりエラーとなったようです。どうやらこのライブラリにはクロスカップリング反応の試薬（有機トリフルオロボレート塩 (Chem-Station さんより)）となるような化合物も含まれているようです。

前処理のコード全体

上記エラーのためforループが途中でとまってしまい、前処理ができなくなってしまいました。そこで、MolLogPの計算については例外処理（try…except…）を追加しました。前処理のコード全体としては下記となります。

from rdkit import Chem
import gzip
from rdkit.Chem import Descriptors, MolStandardize

EAc_gz = gzip.open('Enamine_Advanced_collection.sdf.gz')
EAc_suppl = Chem.ForwardSDMolSupplier(EAc_gz) 

#Molオブジェクト、および各指標を計算し格納するためのからのリストを作成
EAc_mols = []
MW_list = []
MolLogP_list = []
NumHDonors_list = []
NumHAcceptors_list = []
NumRotatableBonds_list = []
TPSA_list = []

for x in EAc_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_neu)

  #分子量をプロパティとして持たせる。(小数点を含むfloat型)
  mol_neu.SetDoubleProp('MW', MW_value)
  
  #値のみのリストにも追加する
  MW_list.append(MW_value)
  
  #MolLogPを計算(MolLogPは対応しない元素を含むとエラーを返すので例外処理として記入)
  try: 
    MolLogP_value = Chem.Descriptors.MolLogP(mol_neu)
    mol_neu.SetDoubleProp('MolLogP', MolLogP_value)
    MolLogP_list.append(MolLogP_value)
  except:
    MolLogP_list.append(None)
   
  #NumHDonorsでも同じことをする。
  NumHDonors_value = Chem.Descriptors.NumHDonors(mol_neu)

  #非負の整数値（int型）のプロパティなので SetUnsignedIntProp でも良いが面倒なのでSetDoublePropにする。
  #たぶんメモリの無駄遣い。
  mol_neu.SetDoubleProp('NumHDonors', NumHDonors_value)
  NumHDonors_list.append(NumHDonors_value)  

  #NumHAcceptors
  NumHAcceptors_value = Chem.Descriptors.NumHAcceptors(mol_neu)
  mol_neu.SetDoubleProp('NumHAcceptors', NumHAcceptors_value)
  NumHAcceptors_list.append(NumHAcceptors_value)  

  #NumRotatableBonds
  NumRotatableBonds_value = Chem.Descriptors.NumRotatableBonds(mol_neu)
  mol_neu.SetDoubleProp('NumRotatableBonds', NumRotatableBonds_value)
  NumRotatableBonds_list.append(NumRotatableBonds_value)    

  #TPSA
  TPSA_value = Chem.Descriptors.TPSA(mol_neu)
  mol_neu.SetDoubleProp('TPSA', TPSA_value)
  TPSA_list.append(TPSA_value)   
 
  #Molオブジェクトのリストも忘れずに
  EAc_mols.append(mol_neu)

以後の処理をできるだけ簡単にするため、指標の数値のみを含むcsvファイルを作成しました。

import pandas as pd
property_names = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA'] 
property_values = [MW_list, MolLogP_list, NumHDonors_list, NumHAcceptors_list, NumRotatableBonds_list, TPSA_list]

df = pd.DataFrame(index=[], columns=property_names)

for i, j in zip(property_names, property_values):
  df[i] = j
  
df.to_csv('Enamine_Advanced_collection_desc.csv', sep=',')

また別途、構造情報と各指標の値を紐づけた上でsdfを作成しました。

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

writer.SetProps(property_names)   #01/29追記 このままだとidが抜けるので適当に追加してください。またやってしまった・・・
for mol in EAc_mols:
  writer.write(mol)
writer.close()

以上が、各ライブラリファイルに対して実行した前処理の全体の流れとなります。

処理に失敗した分子の確認

先の処理でMolLogPの計算ができなかったものについてはNoneを格納するように例外処理を指定しました。このリストをpandasで読み込むとNoneは自動的にNaNとしてDataFrameに取り込まれるそうです。

どのような化合物が処理に失敗したかを確認するため、

上記で出力したcsvファイルをPandasのDataFrameとして読み込む
MolLogPの列で値がNaNとなっている行のindexを取得する
取得したindexに相当する分子の構造を確認するという手順を実施しました。

（上述のValue Errorとなった構造の描画のために行った処理です）

# csvをDataFrameとする
df_EAc = pd.read_csv('./Enamine_Advanced_collection_desc.csv')

# MolLogPの列がNaNとなっている(query)行番号(index)を取り出し、リスト化する(list)
df_EAc_MolLogP_None = list(df_EAc.query('MolLogP == "NaN"').index)

print(df_EAc_MolLogP_None)
# [8135] と出力された

index 8135の分子一つのみが、処理に失敗したもののようです。コードは省略しますが、この分子をsdfから取り出し描画することで、先の構造（ピロールのトリフルオロボレート塩）を確認しました。

SDFの持つ情報をとりあえず確認する方法

上記、前処理ですが Enamine~~ と名前のついた3つのsdfではうまくいったものの、UOS_HTSではうまくいきませんでした。原因はidnumber というプロパティが存在しなかったことです。どうやらこちらのライブラリのみ、化合物のIDが別の名称の属性として格納されているようです。

TextエディタやMarvinViewなどで開いてSDFの持つ構造情報以外の属性を確認する、という手段もありますが、なにぶんファイルサイズが大きく開くだけでも大変です。ここは　RDKitをうまく使いたい・・・ということで、とりあえずsdfから最初の１分子のみを読み込んでみることとしました。

以前こちらの記事で検証したように、ForwardSDMolSupplierで読み込んだSDFからはnext()を使うことで順番に一つずつ取り出すことができます。

UH_suppl = Chem.ForwardSDMolSupplier('./UOS_HTS.sdf') 

#組み込み関数next()を使って最初の一つの分子を取り出す
test_mol = next(UH_suppl)

#どんなPropertyが含まれているか取り出し、リスト化、出力
prop_names = list(test_mol.GetPropNames())
print(prop_names)
#['ID']と出力された

UOS_HTSについては、化合物のIDはIDという属性に格納されているようです。これでidnumberの代わりにIDとして、他のライブラリに行ったのと同じ処理を行えば良いことがわかりました。

なんどもクドイようですが、上記処理に続けて前処理のコードを実行すると失敗します。

理由は属性の確認のため、ForwardSDMolSupplierのオブジェクトから最初の一つをとりだしてしまったため、再度sdfを読み込むところから始めないと、最初の一つを捨ててしまったことになるからです。（私は以前記事を書いておきながらもこの失敗をまた繰り返しました・・・学習能力低すぎ・・・）

統計量の計算

pandasのdescribeを用いると項目数が増えるので、総数（count）、最小値(min)、平均値(mean)、最大値(max)を別個に計算しました。

df_count = df_EAc[descriptors].count() 
df_min = df_EAc[descriptors].min()
df_mean = df_EAc[descriptors].mean()
df_max = df_EAc[descriptors].max()

df_desc_view_EAc =  pd.DataFrame([df_count, df_min, df_mean, df_max], index=['count', 'min', 'mean', 'max'])

df_desc_view_EAc.round(2)

これで冒頭の統計量を出しました。

箱ひげ図を用いたばらつきの可視化

要約統計量のみでは値の広がりがわからないので箱ひげ図を用いてばらつきを可視化しました。

import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

一つのプロットに全ての指標を乗せてしまうと、値の範囲の違い（分子量と水素結合受容体の数では値の幅が全然違う）のため、潰れてしまったのでsubplotを指定しました。

descriptors = ['MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA']

df_EAc[descriptors].plot(kind = 'box', subplots=True, figsize=(20, 6))

これで冒頭の図を得ました。

以上です。