RDKitで化合物の芳香族度合を示す分子記述子を計算する

分子記述子とは分子の性質を決める指標のことです．特にケモインフォマティクスの分野では分子の構造から導くことが可能なものを記述子と呼ぶことが多いです．

一般的には，複数の記述子の組み合わせによって生理活性などの実験データを理解・予測するモデル（QSARモデル）の作成が記述子を用いる目的になります．

さまざまな記述子が報告されており，分子記述子に特化した分厚い専門書が売られている程度には，多数の記述子が既に存在しています．すべての性質に万能な記述子は今のところ存在しないために，目的とする性質に最適な記述子を目指して，「化学者が記述子を作り込んできた」というのが実情でしょう．

これまでに本ブログでも

• 「RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ」
• 「RDKitで立体的な分子の形を記述する」
• 「RDKitで合成難易度を評価して化合物をスクリーニング」
• 「RDKitで天然物らしさをスコア化して化合物ライブラリーの指標にする」

などの記事でいろいろな分子記述子を扱ってきました．本記事では分子の性質のうち，特に「化合物の芳香族度合」に着目した記述子についてみていきたいと思います．

分子への芳香環の導入は，クロスカップリング反応の発展に伴って格段に容易になりました．一方で，複数の芳香環を有する化合物の物性はあまり好ましいものではないことも知られています．

そのため化合物の芳香族の程度を数値的に把握することは創薬化学研究における有用な指針と言えます．

実際さまざまな記述子が提案されていますので，今回は

• 芳香族度合の指標となる記述子の定義や性質
• RDKitを用いた計算の仕方

について概説していきたいと思います．

RDKitの分子Molオブジェクトを扱う

「RDKitでケモインフォマティクスに入門」という記事では，RDKitの分子であるMolオブジェクトの作成方法を中心に学習しました．今回はそのMolオブジェクトに対してどのような操作が行えるかを学んでいきます． 化学は分子を中心に扱う学問ですので，Molオブジェクトの扱いに習熟す...

future-chem.com

2021.01.09

芳香族指数（Aromatic Indicator）：Ar/HA, AP

芳香族指数とは「化合物の芳香族度合」を表す値として導入され，

と定義されています．しばしば「Ar/HA」や「AP」などと表記されています（HAはHeavy Atomですね）．

RDKitを用いたAr/HAの計算

RDKitのMOLオブジェクトには

• 全重原子の数
• 芳香族原子の数

を求めるメソッドが用意されています．それらを用いることで容易にAr/HAを計算することが可能です．

芳香族比率（Aromatic Ratio）：ARR

Ar/HAとほぼ同等の記述子が芳香族比率（ARR）になります．この指標は

と定義されます．Ar/HAのような原子数ではなく，結合数に着目した記述子になります．

分子記述子の計算によく用いられる商用ソフトウェアDragonで計算可能なことから様々なQSARモデルで用いられる傾向にあります．

RDKitを用いたARRの計算

分子中の全結合数はGetNumBondsで取得可能です．また結合が芳香族か否かはBondオブジェクトのGetIsAromaticで調べることが可能です．

これらを組み合わせることでARRを計算することが可能になります．

def Calc_ARR(mh):
    m = Chem.RemoveHs(mh)
    num_bonds = m.GetNumBonds()
    num_aromatic_bonds = 0
    for bond in m.GetBonds():
        if bond.GetIsAromatic():
            num_aromatic_bonds += 1
    ARR = num_aromatic_bonds/num_bonds
    return ARR

芳香環カウント（Aromatic Ring Count）：NAR, AROM

芳香環の数も，芳香族具合を示す際によく使われる分子記述子になります．注意点としては

• ベンゼンは環１つ
• ナフタレンは２つ
• アントラセンは３つ

のように，縮環構造をそれぞれ環の数として数える点があげられます．

芳香環の数が増えると一般的に脂溶性も向上する傾向にあります．しかしCLogPを一定にした化合物群において芳香環数が化合物の水溶性に悪い影響を与えるといった報告もされていることから，芳香環の数が脂溶性とは無関係な何らかの効果を持っていると考えてもよさそうです．

RDKitを用いて芳香環の数をカウント

Descriptorsモジュールには芳香環の数を返すメソッドがありますので，そちらを用いて計算可能です．

ここではRingInfoオブジェクトを用いて自分で関数を書いてみましょう．

MOLオブジェクトのGetRingInfoで得られるRingInfoオブジェクトを利用することで分子中の環構造に関する情報を取得できます．

AtomRingsでは環を形成する原子IDの，BondRingsでは環を形成する結合IDのタプルのタプルが得られます．これらの原子や結合が芳香族か否かを調べることで芳香環の数を数え上げていきます．

以下のコードでは個々の環について環内の全ての原子が「Aromatic」である場合に環を芳香族として数え上げています．

def Calc_AROM(mh):
    m = Chem.RemoveHs(mh)
    ring_info = m.GetRingInfo()
    atoms_in_rings = ring_info.AtomRings()
    num_aromatic_ring = 0
    for ring in atoms_in_rings:
        aromatic_atom_in_ring = 0
        for atom_id in ring:
            atom = m.GetAtomWithIdx(atom_id)
            if atom.GetIsAromatic():
                aromatic_atom_in_ring += 1
        if aromatic_atom_in_ring == len(ring):
            num_aromatic_ring += 1
    return num_aromatic_ring

炭素芳香族とへテロ芳香族数

上記のAROMでは全ての芳香環を一律にカウントしていきました．ただし物性への影響などを考えた際に，「ベンゼン」と「ピリジン」を同じように芳香環１つと数えることに感覚的に違和感を覚える人もいると思います．

そのような場合には芳香環を

• 炭素芳香環
• へテロ芳香環

に区分けして数え上げるとより適切に分子の性質を記述できることがあります．

RDKitを用いて炭素芳香族とへテロ芳香族を数え上げ

Descriptorsモジュールには炭素芳香環，へテロ芳香環の数を返すメソッドもありますのでそれらを使うことですぐに計算可能です．

ここでは先ほどと同様に自分で関数を定義してみます．AROMの計算に用いたコードを改変することでへテロ芳香族の識別が可能です．

具体的には環内の構成原子に炭素以外の元素が含まれているかどうかをチェックしていきます．

def Calc_Carbo_Hetero_Aromatic(mh):
    m = Chem.RemoveHs(mh)
    ring_info = m.GetRingInfo()
    atoms_in_rings = ring_info.AtomRings()
    num_Caromatic_ring = 0
    num_Hetaromatic_ring = 0
    for ring in atoms_in_rings:
        aromatic_atom_in_ring = 0
        heteroatom_in_ring = 0
        for atom_id in ring:
            atom = m.GetAtomWithIdx(atom_id)
            if atom.GetIsAromatic():
                aromatic_atom_in_ring += 1
            if atom.GetSymbol() != 'C': ### 環内の原子が炭素かどうかをチェック
                heteroatom_in_ring += 1
        if aromatic_atom_in_ring == len(ring):
            if heteroatom_in_ring == 0:
                num_Caromatic_ring += 1
            else:
                num_Hetaromatic_ring += 1
    return (num_Caromatic_ring, num_Hetaromatic_ring)

脂肪族指数：Fsp3

芳香族と脂肪族化合物は表と裏の関係ですから，芳香族指数が高い化合物は，必然的に脂肪族性が低い化合物ということになります．

と定義されるFsp3は化合物の脂肪属性の度合を示す指標ですので，小さいFsp3の値が高い芳香族指数となります．

Fsp3は化合物の立体的な特徴を記述する記述子として用いられることが多いです．「RDKitで立体的な分子の形を記述する」という記事ではFsp3をはじめとした，分子の立体的な性質を示す記述子を取り上げています．

RDKitで立体的な分子の形を記述する

「RDKitのおける記述子の扱い方をリピンスキーのルール・オブ・ファイブを通して学ぶ」という記事では，分子の性質を表現する指標として記述子（デスクリプター）について説明し，RDKitを用いて実際に計算を行ってみました．その際に扱った記述子は分子量や水素結合ドナーの数など，2次元の...

future-chem.com

2021.01.17

RDKitにおけるFsp3の計算

RDKitではFsp3を計算するメソッドがChem.Descriptorsに実装されています．

物性予測指数（Property Forecast Index）：PFI

という式で定義されるPFIは，logDやAROMを単体で用いた時よりも様々な化合物の性質をよく記述できることが報告されています．つまりAROMの項でも述べたように，芳香環の数が与える影響が

• 脂溶性と関連する影響
• 脂溶性とは無関係な影響

に切り分けられることを示唆しています．

PFIでは各々を１：１の重み付けで和を取っていますが，これはある値のlogDとAROMを持った化合物を調べた際に，各々の記述子が１変わったさいに化合物の溶解性への影響がほぼ同じであったことに起因しています．

溶解性以外の性質へと展開する場合には，特に根拠はないわけですが計算しやすいことから１：１で和を取っていると考えてよさそうです．

残念ながら2018年12月現在のところ，logDの計算はRDKitには実装されていませんのでPFIは求められません．

芳香族・脂肪族バランス：Ar-sp3

化合物の芳香族性と脂肪族性は表裏一体という説明をしましたが，この２つのバランスを表す指標として

で定義される値が使われることがあります．この指標は芳香族原子の数やsp3炭素の数がゼロの場合にもパラメータ化可能な点で，Ar/HAやFsp3と比べて利点があります．

RDKitを用いて芳香族・脂肪族バランスを計算する

先ほど述べたようにMolオブジェクトのGetAromaticAtomsを用いて芳香族原子の数が求められます．sp3炭素の数はAtomオブジェクトのGetHybridizationとGetSymbolを組み合わせることで数え上げています．

def Calc_Ar_Alk_balance(mh):
    m = Chem.RemoveHs(mh)
    num_aromatic_carbon = len(m.GetAromaticAtoms())
    num_sp3_carbon = 0
    for atom in m.GetAtoms():
        if str(atom.GetHybridization()) == 'SP3' and atom.GetSymbol() == 'C':
            num_sp3_carbon += 1
    ar_alk_balance = num_aromatic_carbon - num_sp3_carbon
    return ar_alk_balance

ベンゾイド指数：BI

先ほども述べたように炭素芳香環とへテロ芳香環が物性に与える影響が異なると考えられるため，

としてベンゾイド指数を定義します．

大きいBI値は高い炭素芳香環の割合を示し，開発成功率が低下する傾向があります．

RDKitを用いてベンゾイド指数を計算する

ベンゾイド指数はDescriptorsモジュールの

• NumAromaticCarbocycles
• NumAromaticRings

を使うことで容易に計算可能です．

芳香族系分子記述子の相関関係について

今回扱ってきた記述子は全て「芳香族度合」を表すように設計されたものです．そのため各記述子間ではある程度の相関があることが予想されます．

同じような情報を含んだ記述子を複数用いてモデル構築を行うべきではありませんので，ここでは今回取り上げた記述子の間の相関係数を眺めてみることとします．

分子としてはZINC15データベースから取得した，酵素をターゲットとしてin vivoで活性のある化合物を抽出したSDFを使ってみます．このSDFには653化合物が含まれていまして，これらをpandasのデータフレームとして読み込みます．

以下のコードでは

1. 分子をpandasのデータフレームとして読み込み
2. pandasのmap関数を用いて各種記述子の一括計算
3. データの標準化
4. 各記述子間の相関係数を計算

という作業を行っています．

相関係数を計算するだけなら標準化（データの平均を0，標準偏差を１に揃える作業）は必要ありませんが，今後のモデル構築を考えると標準化を行う癖をつけておいた方がよいと思います．

### ライブラリのimport
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, PandasTools, Descriptors
import pandas as pd
print(rdBase.rdkitVersion) ### 2018.09.1
### 分子の準備
df = PandasTools.LoadSDF('enzyme-in-vivo.sdf')
len(df) ### 653
### 記述子の計算
df['Ar/HA'] = df.ROMol.map(lambda x: len(x.GetAromaticAtoms())/x.GetNumHeavyAtoms())
df['ARR'] = df.ROMol.map(Calc_ARR)
df['AROM'] = df.ROMol.map(Descriptors.NumAromaticRings)
df['Fsp3'] = df.ROMol.map(Descriptors.FractionCSP3)
df['Ar-sp3'] = df.ROMol.map(Calc_Ar_Alk_balance)
### データの標準化
df['Ar/HA_STD'] = (df['Ar/HA'] - df['Ar/HA'].mean()) / df['Ar/HA'].std()
df['ARR_STD'] = (df['ARR'] - df['ARR'].mean()) / df['ARR'].std()
df['AROM_STD'] = (df['AROM'] - df['AROM'].mean()) / df['AROM'].std()
df['Fsp3_STD'] = (df['Fsp3'] - df['Fsp3'].mean()) / df['Fsp3'].std()
df['Ar-sp3_STD'] = (df['Ar-sp3'] - df['Ar-sp3'].mean()) / df['Ar-sp3'].std()

df[['Ar/HA_STD', 'ARR_STD', 'AROM_STD', 'Fsp3_STD', 'Ar-sp3_STD']].corr().round(2)

計算された相関係数は以下の通りでした．Ar/HAとARRは芳香族の原子か結合に注目するかの違いだけですので，ほぼ同じ情報を含んでいることがわかります．

終わりに

今回は「RDKitで化合物の芳香族度合を示す分子記述子を計算する」という話題について

• 分子の芳香族具合を表すさまざまな記述子の定義と説明
• それら記述子をRDKitを用いてどのように計算するか
• 計算した記述子はお互いにどのような相関があるのか

について説明してきました．

いろいろと手の込んだ記述子が発表されていますが，記述子は基本的には分子の部分的な特徴を表現しているに過ぎません．一方で分子全体の構造的特徴を網羅するような表現方法として分子の「フィンガープリント」があります．

フィンガープリントはビット配列の羅列ですので，人間には理解しにくいですがコンピューターにはとても扱いやすい分子の表現方法の１つといえます．

＞＞次の記事：「RDKitでフィンガープリントを使った分子類似性の判定」