Python手遊び(化合物のSDFをRDBへInsert)

#この記事、何?
機械学習をあれこれやりたいので、化合物に関する情報をRDB(postgres)に入れたい。
ひとまずうまくいったので、備忘録かねてそのデータ・コマンド・などの記録。

#どういう人向け?
「ケモインフォマティクス」を始めようとする方向け。
多分、コピペして実行すれば動くはず・・・というつもり。
まあ・・・この後、scikit-learnでSVMとかRandomForestとかを試してみたい。

#ざっくりの方針

化合物ファイルを入手し、RDKitで処理してPostgreSQLに格納。

今までの記事のいくつかを参照。

##参考URL
■化合物を300万個目指していっぱい集めてみた
https://qiita.com/siinai/items/7b4196f4271448093a1f

■化合物の記述子：RDKitを計算してみた
https://qiita.com/siinai/items/1803142350df41185e73

■PostgreSQLを入れてみた
https://qiita.com/siinai/items/457977b61b6c4cb15ac4

・・・全部自分の投稿ですけど。。。

#できあがり

##データ

今後の検証用もかねて、化合物を20個格納したSDF作ってた。
ちょっと長いけどデータも・・・と思ったけどさすがに平文だと長すぎるので割愛。
PubChemから以下のファイルを取得

■PubChem
ftp://ftp.ncbi.nlm.nih.gov/pubchem/Compound/CURRENT-Full/SDF/
Compound_000000001_000025000.sdf

##テーブル定義
テーブル定義はこれ。
200記述子分を連番で作っておいた。あとSMILESとInChIも。

```sql
--22-02.sql
create table t_sample (
sampleid integer
,samplename text
,molblock text
,smiles text
,inchi text
,c001 double precision
,c002 double precision
,c003 double precision
----さすがに長いので割愛---
,c199 double precision
,c200 double precision
,primary key (sampleid)
)
;

```

##コード

```python
#22-01.py
import rdkit
import pandas as pd
from rdkit import Chem
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors

import psycopg2

print('--- import sdf has 20 compounds ---')
mols = Chem.SDMolSupplier('Compound_000000001_000025000.sdf')

print('--- descriptors ------------------------')
descs = [desc_name[0] for desc_name in Descriptors._descList]
descs_cnt = len(descs)
if descs_cnt > 200:
    descs_cnt = 200
print(descs)

print('--- calculate ---')
desc_calc = MoleculeDescriptors.MolecularDescriptorCalculator(descs)

print('--- define connection ---')
conn = psycopg2.connect('postgresql://postgres:manager@localhost:5432/postgres')
cur = conn.cursor()

sql = 'delete from t_sample'
cur.execute(sql)

for mol in [x for x in mols if x is not None]:
    sampleid = mol.GetProp('PUBCHEM_COMPOUND_CID')
    samplename = mol.GetProp('_Name')
    molblock = Chem.MolToMolBlock(mol)
    smiles = Chem.MolToSmiles(mol)
    inchi = Chem.MolToInchi(mol)
    values = desc_calc.CalcDescriptors(mol)

    sql = 'insert into t_sample (sampleid, samplename, molblock, smiles, inchi'
    for i in range(descs_cnt):
        sql += ',c' + str(i+1).zfill(3)
    sql += ')'
    sql += ' values ('
    sql += "%s, %s, %s, %s, %s"
    for i in range(descs_cnt):
        sql += ', %s'
    sql += ')'

    paras = []
    paras.append(sampleid)
    paras.append(samplename)
    paras.append(molblock)
    paras.append(smiles)
    paras.append(inchi)
    for i in range(descs_cnt):
        paras.append(values[i])
    
    cur.execute(sql, paras)

conn.commit()
cur.close()
conn.close()

```

##データ
```sql
select sampleid, smiles, inchi, c001, c002, c003 from t_sample order by sampleid limit 5;
```

|sampleid|smiles|inchi|c001|c002|c003|
|:-:|:-:|:-:|:-:|:-:|:-:|
|1|CC(=O)OC(CC(=O)[O-])C[N+](C)(C)C|InChI=1S/C9H17NO4/c1-7(11)14-8(5-9(12)13)6-10(2,3)4/h8H,5-6H2,1-4H3|10.67899282|-1.203981481|10.67899282|
|2|CC(=O)OC(CC(=O)O)C[N+](C)(C)C|InChI=1S/C9H17NO4/c1-7(11)14-8(5-9(12)13)6-10(2,3)4/h8H,5-6H2,1-4H3/p+1|10.69940098|-0.953981481|10.69940098|
|3|O=C(O)C1=CC=CC(O)C1O|InChI=1S/C7H8O4/c8-5-3-1-2-4(6(5)9)7(10)11/h1-3,5-6,8-9H,(H,10,11)|10.33703704|-1.30787037|10.33703704|
|4|CC(O)CN|InChI=1S/C3H9NO/c1-3(5)2-4/h3,5H,2,4H2,1H3|8.236111111|-0.337962963|8.236111111|
|5|NCC(=O)COP(=O)(O)O|InChI=1S/C3H8NO5P/c4-1-3(5)2-9-10(6,7)8/h1-2,4H2,(H2,6,7,8)|10.26513889|-4.509158951|10.26513889|

```sql
select count(*) from t_sample;
```
※23255件でした。



#感想
いい意味であまり面白くなかった。
これまでにあれこれ用意してたので結構すんなり。
あと、23000件も処理してもトランザクションもてるのね・・・
多分落ちるだろうから、それから何とかしようと思ってたので最後にコミットしてそのまま。


#今後
まずはこの後、処理時間計測用のスクリプトを作ろうかな、と。
デコレートでいい感じにしようかと思ったけど、処理側に一切干渉せずに測れるようにしたいな。。。
と、もちろん、冒頭に書いた機械学習あれこれを。