Pythonで分子の性質を調べる　ー飲み薬の成分候補の探索ー（情報Ⅰプログラミング×化学②）

1.　分子の性質（分子量、水素結合ドナー・アクセプター数、LogP）

今回はRDKitやPubChemPyのライブラリを使って、分子の性質を調べていきます。

今回は新しい医薬品を作る際に、その候補成分として着目されることがある、以下の4つの性質をPythonで情報を集めたいと思います💡

・分子量: 分子に含まれる全ての原子の原子量の総和です。この原子量とは、¹²C（質量数12の炭素原子）の質量を12と定め、これを基準として表した相対的な質量のことです。
・水素結合ドナー数: 分子の中で、他の分子に水素原子を提供して水素結合を作ることができる部分の数です。主に-OH基や-NH基などが該当します。
・水素結合アクセプター数: 分子の中で、他の分子から水素原子を受け取って水素結合を作ることができる部分の数です。主に酸素原子(O)や窒素原子(N)が該当します。
・LogP (分配係数): オクタノール/水 分配係数の常用対数値のことです。これは分子が水と油のどちらに馴染みやすいかの度合いを示す数値指標で、LogPの値が大きいほど油に溶けやすく、小さいほど水に溶けやすい傾向があります。

2.　Pythonで分子の性質を計算する準備

Pythonを使って分子の性質を計算する準備を始めます。
Google Colaboratory (Colab) などのPython実行環境を使います。

◆Step 1：必要なライブラリをインストール＆インポート
今回使うRDKitとPubChemPy（PubChemから情報を取得するためのライブラリ）をインストールします。Colabでは、以下のコードをセルに入力して実行してください。（もし既にインストール済みであれば、このステップは不要です。）

!pip install rdkit-pypi pubchempy

Pythonでこれらのライブラリを使えるようにします。今回は、性質計算に使うRDKitのモジュールも一緒にインポートします。

from rdkit import Chem
from rdkit.Chem import Draw # 分子構造の描画用
from rdkit.Chem import Descriptors # 分子量やLogPなどの計算用
from rdkit.Chem import Lipinski # 水素結合ドナー/アクセプター数などを計算用
import pubchempy as pcp

#Jupyter NotebookやGoogle Colabで構造をきれいに表示
from rdkit.Chem.Draw import IPythonConsole
IPythonConsole.ipython_useSVG=True # SVG形式で表示（推奨）

◆Step 2：調べたい化合物名を入力しよう！
次に、性質を調べたい化合物の名前（英語名）を指定します。

今回は、カフェイン（コーヒーにも含まれ、眠気防止や痛み止めにも使用される成分）、アスピリン（痛み止めとして使用される成分）、オレイン酸（食用油として使用される）の3つについて調べていきます。

# ここに調べたい化合物の英語名を入力してください
compound_name = "caffeine"

上記のコードの “caffeine” の部分を、調べたい化合物の英語名に書き換えると、以降のコードで、調べたい成分の情報を得ることができます。
“aspirin”（アスピリン）、”oleic acid”(オレイン酸)、”ethanol”（エタノール）、”water”（水）など

◆Step 3：PubChemから分子情報を取得し、RDKitで扱える形に！
指定した化合物名を使って、PubChemデータベースからその分子のSMILESを取得します。そのSMILESを元にRDKitが扱える分子オブジェクトを作成します。この分子オブジェクトで、様々な性質を計算していきます。

mol = None # RDKitの分子オブジェクトを格納する変数 (最初は空)
smiles = None # SMILES文字列を格納する変数 (最初は空)

try:
    # PubChemで化合物名を検索
    compound_data_list = pcp.get_compounds(compound_name, 'name')
    
    if compound_data_list: # もし化合物が見つかったら
        pubchem_compound = compound_data_list[0] # 最初の検索結果を使う
        # SMILES文字列を取得 (isomeric_smilesがあればそれを、なければcanonical_smiles)
        smiles = pubchem_compound.isomeric_smiles if pubchem_compound.isomeric_smiles else pubchem_compound.canonical_smiles
        
        if smiles: # SMILESが取得できたら
            print(f"取得したSMILES: {smiles}")
            # SMILESからRDKitの分子オブジェクトを作成
            mol = Chem.MolFromSmiles(smiles)
            
            if mol: # 分子オブジェクトがちゃんとできたら
                print(f"'{compound_name}' の分子オブジェクト作成に成功しました！")
                # 分子構造も表示してみましょう
                display(Draw.MolToImage(mol, size=(300,300))) # Colab/Jupyter環境で表示
            else:
                print(f"エラー: SMILES '{smiles}' から分子オブジェクトを作成できませんでした。")
        else:
            print(f"エラー: '{compound_name}' のSMILES情報がPubChemで見つかりませんでした。")
    else:
        # PubChemで指定された化合物名が見つからなかった場合
        print(f"エラー: '{compound_name}' がPubChemデータベースで見つかりませんでした。英語の名称か、スペルを確認してください。")

except Exception as e: 
    print(f"処理中に予期せぬエラーが発生しました: {e}")

# 出力結果
取得したSMILES: CN1C=NC2=C1C(=O)N(C(=O)N2C)C
'caffeine' の分子オブジェクト作成に成功しました！

このコードを実行すると、指定した化合物のSMILESと、その構造式（画像）が表示されます。

もしエラーが出た場合は、エラーメッセージをよく読んで、化合物名が正しいか、インターネットに繋がっているかなどを確認してみてください。

・compound_name =aspirinの場合
取得したSMILES: CC(=O)OC1=CC=CC=C1C(=O)O
‘aspirin’ の分子オブジェクト作成に成功しました！
・oleic acidの場合
取得したSMILES: CCCCCCCC/C=C\CCCCCCCC(=O)O
‘oleic acid’ の分子オブジェクト作成に成功しました！

3.　分子量の計算

作成した分子オブジェクト (mol) を使って、RDKitのDescriptors.MolWt() で分子量を計算します。

# caffeineの場合
Descriptors.MolWt(mol) 

# 出力結果
194.194

・aspirinの場合: 180.15899999999996
・oleic acidの場合: 282.4679999999999

num_h_donors = Lipinski.NumHDonors(mol)
num_h_acceptors = Lipinski.NumHAcceptors(mol)
print(f"水素結合ドナーの数 ({compound_name}): {num_h_donors}")
print(f"水素結合アクセプターの数 ({compound_name}): {num_h_acceptors}")

# 出力結果
水素結合ドナーの数 (caffeine): 0 
水素結合アクセプターの数 (caffeine): 6

log_p_value = Descriptors.MolLogP(mol)
print(f"計算されたLogP ({compound_name}): {log_p_value:.2f}")

# 出力結果
計算されたLogP (caffeine): -1.03

pubchem_xlogp = "N/A"

try:
    # 化合物名でPubChemを検索し、Compoundオブジェクトを取得
    compounds = pcp.get_compounds(compound_name, 'name')
    if compounds:
        target_compound_obj = compounds[0]

        # Compoundオブジェクトの xlogp 属性を取得
        if hasattr(target_compound_obj, 'xlogp') and target_compound_obj.xlogp is not None:
            pubchem_xlogp = f"{float(target_compound_obj.xlogp):.2f}" # 小数点以下2桁で表示
        else:
            pubchem_xlogp = "xlogp属性なし"

        print(f"PubChemの計算値 (XLogP3相当): {pubchem_xlogp}")
    else:
        print(f"PubChemで {compound_name} が見つからず、比較できませんでした。")
except Exception as e:
    print(f"PubChemからの情報取得中にエラーが発生しました: {e}")

# 出力結果
PubChemの計算値 (XLogP3相当): -0.10