1. openbabelによるMOPACのインプットファイルの作成
複数分子
smilesのリストの複数分子を変換するやり方
複数分子の記載があるdataframeから分子を作成する
!wget https://raw.githubusercontent.com/chemoinformatics-lecture/lecture-beginner/main/lesson02_jupyter/data/compounds-structures.csv
dfに読み込みます。
import pandas as pd
df = pd.read_csv("compounds-structures.csv")
お試しなので、3個で行ってみたいと思います。先頭の3行だけを実行します。
df = df.head(3)
df
3つの分子が表示されます。
smilesの列があるので、これで分子を読み込んでみます。また、名前はidにしてみます。
前のページで行ったことを合わせて行ってみます。
(forループの中で複数のリストの要素を同時に取得して使いたいので、zip()関数を使っています)
for id, smiles in zip(df['id'], df["smiles"]):
print(id)
print(smiles)
m = pybel.readstring("smiles", smiles)
m.make3D(forcefield='mmff94', steps=1000)
mopac_format = m.write("mop")
mopac_format = '\n'.join(mopac_format.splitlines()[3:])
memo1 = m.write("smiles") #任意の文字_メモ書き 。今回は分子のsmilesを書き込んだ。
memo2 = m.write("inchi") #任意の文字_メモ書き。今回は分子のinchiを書き込んだ。
mopac_code = mopac_keyword + memo1 + memo2 + mopac_format
path = str(id) + ".dat"
with open(path, mode='w') as f:
f.write(mopac_code)
これを表示させると、左側のフォルダに3.dat, 5.dat, 6.datというファイルが作成されているのがわかると思います。また、このファイルをクリックすると、ファイルの中身が確認できて、MOPAC形式のファイルが出力されているのが確認できます。
sdfファイルから複数分子を変換するやり方
複数分子の計算方法の仕方の2番めとしてsdfから読み込む方法を行います。また、作成したファイルをgoogle driveに保存することを行ってみます。(たくさんのファイルを行うときにはgoogle driveに保存しておくほうが良いと思います)
google driveに保存しますので、接続します。
from google.colab import drive
drive.mount('/content/drive')
sdfファイルは、以下のものを使用します。
!wget -P "/content/drive/MyDrive/data" https://raw.githubusercontent.com/chemoinformatics-lecture/lecture-beginner/main/lesson05_MOPAC/PubChem_compound_list_antioxidant.sdf
pybel.readfile
pybel.readfileでsdfを指定すると、sdfのまま読み込むことができます。
molecules = pybel.readfile("sdf","/content/drive/MyDrive/data/PubChem_compound_list_antioxidant.sdf")
# どれくらいの分子が入っているかを確認
len(list(molecules))
今回は、122個の分子のsdfファイルになります。
試しに読み込んだ分子構造をsmilesで表示させてみます。
for molecule in pybel.readfile("sdf","/content/drive/MyDrive/data/PubChem_compound_list_antioxidant.sdf"):
print(molecule.write())
これまでのまとめとして、複数分子のmopac用の計算ファイルを出力します。
今回は、molecule.titleにinchikeyをいれて、inchikeyをファイル名にしました。
smilesやInChIは分子によっては、長くなりすぎたり特定の記号が使われてソフトによって使われている文字によりエラーになる(たとえば/
, @
, []
などの文字)のでinchikeyをファイル名にしました。(inchikeyは、25文字の固定長でデジタル表現なので人間には読むことができません。ただし、InChIそのものとは異なり、InChIKeyは一意ではなく、非常に稀ではあるが重複が発生するのでmopacのメモの欄にsmilesやinchiを書き込んでおきます。)
# 今まで行ってきたことのまとめで複数分子を同時に出力する。
for molecule in pybel.readfile("sdf","/content/drive/MyDrive/data/PubChem_compound_list_antioxidant.sdf"):
molecule.make3D(forcefield='mmff94', steps=1000)
mopac_format = molecule.write("mop")
mopac_format = '\n'.join(mopac_format.splitlines()[3:])
mopac_cartesian = mopac_format.replace('PUT KEYWORDS HERE\n\n\n', '')
mopac_keyword = "PM7 EF\n" #ここにキーワードを書く
memo1 = molecule.write("smiles") #任意の文字_メモ書き 。今回は分子のsmilesを書き込んだ。
memo2 = molecule.write("inchi") #任意の文字_メモ書き。今回は分子のinchiを書き込んだ。
mopac_code = mopac_keyword + memo1 + memo2 + mopac_cartesian
molecule.title = molecule.write("inchikey") #ファイル名をinchikeyにするために、分子のタイトルをつける。
path = "/content/drive/MyDrive/data/mopac/"+molecule.title+".dat"
with open(path, mode='w') as f:
f.write(mopac_code)
ファイルを実行すると、MyDrive > data > mopac の中にdatファイルが作成できてくる。
合計122個を出力するので、少し時間はかかります。
他の方法(グラフィカルで入力ファイルを作成する)
- windowsUserの人は、winmostarという日本の会社が作成しているソフトがおすすめです。
無料でも使用できる部分があります。
- Chem3Dのライセンスを持っている人は、Chem3Dでも作成できます。
MOPACのキーワード
MOPACのキーワードは色々なものが設定できます。詳しくは、以下のサイトなどを参考にしてください。
英語
日本語