失われたIDを探せ！ - magattacaのブログ

先の２つの記事で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