Python による RDKit 練習ノート

インフォメーション

オンラインドキュメント

RDKit: Open-Source Cheminformatics Software 内に見つかる，下記ドキュメントで勉強しました．

目次（ページ内リンク）

最初の例

サンプルコード 1

「始めよう」の冒頭部分から抜粋・改変して単純なサンプルコードを作成してみました．
このコードは，ベンゼンの骨格を SMILES 形式の文字列から作成し，MDL MOL 形式で出力します．
「始めよう」には，「基本的な分子を取り扱う機能の多くは rdkit.Chem モジュールに含まれています」と書いてあります．
サンプルコード 1 でも rdkit.Chem モジュールの関数を使っています．

from rdkit import Chem

tmp_mol = Chem.MolFromSmiles(SMILES='c1ccccc1')  #rdkit.Chem.rdchem.Mol オブジェクトを生成．キーワード引数で引数を与えている
tmp_str = Chem.MolToMolBlock(mol=tmp_mol)  #MDL MOL 形式の文字列を生成．キーワード引数で引数を与えている
print(tmp_str)

Chem.MolFromSmiles() 関数

RDKit API Reference を調べると，rdkit.Chem.rdmolfiles module 内に Chem.MolFromSmiles() が見つかります．
SMILES 形式の文字列から rdkit.Chem.rdchem.Mol オブジェクトを生成する関数で，引数を 2 個取る関数と 3 個取る関数が記載されています．

1. rdkit.Chem.rdmolfiles.MolFromSmiles((AtomPairsParameters)SMILES, (SmilesParserParams)params) → Mol :
2. MolFromSmiles( (AtomPairsParameters)SMILES [, (bool)sanitize=True [, (dict)replacements={}]]) -> Mol :

後者の関数は，第二引数以降にデフォルト引数が設定されているので省略可能です．サンプルコード 1 では引数が 1 個なので，後者の関数が呼び出されるようです．
Python と RDKit の初心者（== 筆者）には分かりにくいので，後者の関数を [] に注目して分解してみます．

• MolFromSmiles( (AtomPairsParameters)SMILES
• [, (bool)sanitize=True
• [, (dict)replacements={}]]
• -> Mol :

サンプルコード 1 では，第一引数 (AtomPairsParameters)SMILES の部分に 'c1ccccc1' を充てて SMILES='c1ccccc1' としています．
第二引数と第三引数は省略しているので，デフォルト引数が充てられます．

関数の戻り値は，"→ Mol" で示されています．Mol オブジェクトの作成に失敗した場合は，None を返すとのことです．
この Mol クラスに関するページは，RDKit API Reference では rdkit.Chem.rdchem.Mo に記されています．
Mol クラスには，分子，原子，結合の状態を取得したり，分子を編集したりする関数が見つかります．

Chem.MolToMolBlock() 関数

上と同じく，RDKit API Reference では，rdkit.Chem.rdmolfiles.MolToMolBlock((Mol)mol[, (bool)includeStereo=True[, (int)confId=-1[, (bool)kekulize=True[, (bool)forceV3000=False]]]]) → str :．
分解してみると，

• MolToMolBlock((Mol)mol
• [, (bool)includeStereo=True
• [, (int)confId=-1
• [, (bool)kekulize=True
• [, (bool)forceV3000=False]]]])
• → str :

サンプルコード 1 では，第一引数 (Mol)mol の部分にコード内で作成した tmp_mol を充てて mol=tmp_mol としています．
第二引数以降はデフォルト引数が与えられているので，省略しています．
戻り値は"→ str" すなわち，文字列です．もう少し詳しく言うと，MDL MOL 形式で分子が表現された文字列です．

実行

サンプルコード 1 を sample1.py という名前で保存し，そのディレクトリで，
~$ python3 sample1.py
としたら，下の内容が表示されました．

（ここは空行）
     RDKit          2D
（ここは空行）
  6  6  0  0  0  0  0  0  0  0999 V2000
    1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0
  2  3  1  0
  3  4  2  0
  4  5  1  0
  5  6  2  0
  6  1  1  0
M  END

最初の例に引数を設定してみる

サンプルコード 2

上のサンプルコード 1 を改変して，引数に値を与える練習をしました．
サンプルコード 2 では，デフォルト値を与えているので，出力はサンプル 1 と同じです．
コードの前半では，引数のひとつ params がどのようなオブジェクトなのかを調べています．

from rdkit import Chem

#MolFromSmiles() 関数の第二引数となるオブジェクトを調べる
#まず，prm オブジェクトを作成
prm = Chem.SmilesParserParams()

#prm のデフォルト値を出力
print("prm のデフォルト値")
print(prm.allowCXSMILES, prm.debugParse, prm.parseName, prm.removeHs, prm.sanitize, prm.strictCXSMILES)

#代入可能なことを，デフォルト値を代入して確認
prm.allowCXSMILES = True
prm.debugParse = False
prm.parseName = True
prm.removeHs = True
prm.sanitize = True
prm.strictCXSMILES = True


#上で調べた prm を第二引数にした
mol0 = Chem.MolFromSmiles(SMILES='c1ccccc1', params=prm)

#第二引数以降は，RDKit API Reference に記されているデフォルト引数を与えた
mol1 = Chem.rdmolfiles.MolToMolBlock(mol=mol0, includeStereo=True, confId=-1, kekulize=True, forceV3000=False)

print("mol1")
print(mol1)

MolFromSmiles() 関数の第二引数 params について

第二引数の SmilesParserParams を直訳すると，SMILES 形式の解析を制御するためのパラメータ，ということです．
RDKit API Reference を検索すると，SmilesParserParams へのリンクが見つかります．
class rdkit.Chem.rdmolfiles.SmilesParserParams((object)arg1) というクラスなのでした．
このクラスは，次の property を保有しています．
RDKit API Reference にはデフォルト値が記されていないので，サンプルコード 2 ではデフォルト値を調べるコードを入れてあります．

• property allowCXSMILES：controls whether or not the CXSMILES extensions are parsed
• property debugParse：controls the amount of debugging information produced
• property parseName：controls whether or not the molecule name is also parsed
• property removeHs：controls whether or not Hs are removed before the molecule is returned
• property sanitize：controls whether or not the molecule is sanitized before being returned
• property strictCXSMILES：controls whether or not problems in CXSMILES parsing causes molecule parsing to fail

さらに RDKit API Reference をページ内検索すると，以下のセクションへのリンクを見つけることができます．執筆時点で，まだ理解が不充分です．

• CXSMILES/CXSMARTS extensions：SMILES/SMARTS の ChemAxon 拡張
• Reading molecule names：RDKit Mol オブジェクトの Name の構文解析拡張
• Molecular Sanitization：RDKit Mol オブジェクトのデータ構造を矛盾のないものに調整する

実行

サンプルコード 2 を sample2.py という名前で保存し，そのディレクトリで，
~$ python3 sample2.py
としたら，下の内容が表示されました．

prm のデフォルト値
True 0 True True True True
mol1
（ここは空行）
     RDKit          2D
（ここは空行）
  6  6  0  0  0  0  0  0  0  0999 V2000
    1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0
  2  3  1  0
  3  4  2  0
  4  5  1  0
  5  6  2  0
  6  1  1  0
M  END

デフォルト値以外の引数でどのように動作が変わるのか，調べる準備が整いました．

ファイルを開く

「始めよう」の冒頭で，MDL MOL 形式のファイルを開く関数が紹介されています．そのまま引用すると，
m = Chem.MolFromMolFile('data/input.mol')
この関数を RDKit API Reference で探すと，次の項目が見つかります．
rdkit.Chem.rdmolfiles.MolFromMolFile((str)molFileName[, (bool)sanitize=True[, (bool)removeHs=True[, (bool)strictParsing=True]]]) → Mol :

この関数は，名前が示すように MDL MOL 形式のファイルを読み込みます．
戻り値は，"→ Mol" で示されているように，rdkit.Chem.rdchem.Mol オブジェクトです．
上と同じように，[] に注目して引数を書き下すと下のようになります．第二引数以降は省略可能で，デフォルト引数はすべて True です．

1. MolFromMolFile((str)molFileName
2. [, (bool)sanitize=True
3. [, (bool)removeHs=True
4. [, (bool)strictParsing=True]]])
5. → Mol :

サンプルコード 3

この関数はファイルを開く関数なので，カレントディレクトリに input0.mol という MDL MOL 形式のファイルを作成しておきます．すなわち，

（ここは空行）
     RDKit          3D
（ここは空行）
  2  1  0  0  0  0  0  0  0  0999 V2000
   -0.6290   -0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
M  END

inut0.mol を読み取るサンプルコードを作成しました．
このコードは，input0.mol を読み込んで，その内容を端末に出力します．

from rdkit import Chem

tmp_mol = Chem.MolFromMolFile(molFileName='input0.mol')  #rdkit.Chem.rdchem.Mol オブジェクトを生成
tmp_str = Chem.MolToMolBlock(mol=tmp_mol)  #MDL MOL 形式の文字列を生成
print(tmp_str)

実行

サンプルコード 3 を sample3.py という名前で，カレントディレクトリに保存したものとします．
そのディレクトリで，
~$ python3 sample3.py
としたら，input0.mol の内容が端末に表示されました．

ファイルに保存

サンプルコード 4

「始めよう」の冒頭部分から，ファイル出力している文をピックアップして，サンプルコードにまとめました．

from rdkit import Chem

#SMILES 形式から rdkit.Chem.rdchem.Mol オブジェクトを生成．サンプルコード 1 の同じ方法
mol0 = Chem.MolFromSmiles(SMILES='C')  #メタン
mol1 = Chem.MolFromSmiles(SMILES='CC')  #エタン
mols = {mol0, mol1}  #mol0 と mol1 からなるリスト

#MOL MOL 形式で mol0 をファイル出力
print(Chem.MolToMolBlock(mol0),file=open('mol0.mol','w'))

#SDF 形式で mols（mol0 と mol1）をファイル出力
write = Chem.SDWriter('mols1.sdf')
for m in mols:
  write.write(m)

#SMILES 形式で mols（mol0 と mol1）をファイル出力
write = Chem.SmilesWriter('mols2.sml')
for m in mols:
  write.write(m)

これまでに RDKit API Reference を読む練習をしたので，RDKit API Reference へのリンクを作成すれば充分かと思います．
MolToMolBlock() は関数，SDWriter と SmilesWriter はクラスです．

• rdkit.Chem.rdmolfiles.MolToMolBlock((Mol)mol[, (bool)includeStereo=True[, (int)confId=-1[, (bool)kekulize=True[, (bool)forceV3000=False]]]]) → str :
• class rdkit.Chem.rdmolfiles.SDWriter((AtomPairsParameters)arg1, (AtomPairsParameters)arg2)
• class rdkit.Chem.rdmolfiles.SmilesWriter((AtomPairsParameters)arg1, (AtomPairsParameters)fileObj[, (str)delimiter=' '[, (str)nameHeader='Name'[, (bool)includeHeader=True[, (bool)isomericSmiles=True[, (bool)kekuleSmiles=False]]]]])

実行

サンプルコード 4 を sample1.py という名前で保存し，そのディレクトリで，
~$ python3 sample4.py
としたら，mol0.mol，mols1.sdf，および mols2.sml というファイルが生成しました．

mol0.mol

（ここは空行）
     RDKit          2D

  1  0  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
M  END

mols1.sdf

（ここは空行）
     RDKit          2D

  2  1  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.2990    0.7500    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
M  END
$$$$
（ここは空行）
     RDKit          2D

  1  0  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
M  END
$$$$

mols2.sml

SMILES Name 
CC 0
C 1

補遺 関数のキーワード引数と位置引数

引数の使い方で混乱したので，この場を利用して整理しておきます．
このセクションで引用しているページは，Python 3.12.2 ドキュメントの下にあります．

Python における実引数

Python 用語集を見ると，実引数は 2 種類あります．キーワード引数と位置引数です．ほぼそのまま引用します．

キーワード引数: 関数呼び出しの際に引数の前に識別子がついたものや， ** に続けた辞書の中の値として渡された引数．
例えば，次の complex() の呼び出しでは， 3 と 5 がキーワード引数です:

• complex(real=3, imag=5)
• complex(**{'real': 3, 'imag': 5})

位置引数: キーワード引数以外の引数．
位置引数は引数リストの先頭に書くことができ，また * に続けた iterable の要素として渡すことができます．
例えば，次の例では 3 と 5 は両方共位置引数です:

• complex(3, 5)
• complex(*(3, 5))

引数の書式

もう少し具体的な説明が，Python チュートリアルにあります．4.8. 関数定義についてもう少し．
いくつか文言を引用します．

• キーワード引数は位置引数の後でなければならない
• キーワード引数は全て，関数で受け付けられる引数のいずれかに対応していなければならない
• キーワード引数の順序は重要ではない
• いかなる引数も値を複数回は受け取れない

サンプルコード

これらを踏まえて関数定義の練習をしたときのコードを，多少整理して載せておきます．

def tmp_func(num:int=123, string:str='qwerty'):
  if type(num) != (int):
    print("num 型エラー")
    return
  if type(string) != (str):
    print("string 型エラー")
    return

  print(type(num), type(string))
  print(num, string)


# キーワード引数を使うと，引数の順番を入れ替えてもエラーにならないし，任意の引数を省略できる

print("Case 0：キーワード引数")
tmp_func(num=456, string="Hello")

print('\n',"Case 1：キーワード引数")
tmp_func(string="Hello", num=456)

print('\n',"Case 2：キーワード引数")
tmp_func(string="Hello")

print('\n',"Case 3：キーワード引数")
tmp_func(num=456)

print('\n',"Case 4：位置引数")
tmp_func(456, "Hello")

print('\n',"Case 5：位置引数")
tmp_func(456)


# キーワード引数と位置引数の混在
print('\n',"Case 6：キーワード引数と位置引数の混在")
tmp_func(456, string="Hello")

# これはエラーになる
#print('\n',"Case 7：キーワード引数と位置引数の混在")
#tmp_func(string="Hello", 456)


# 位置引数の型エラー
# def tmp_func() 関数の if 文をコメントアウトするとエラーにならず，端末に出力される
# もちろん，引数で何らかの処理をするとエラーになる可能性が高い
print('\n',"Case 8：位置引数")
tmp_func("Hello", 456)