機械学習メモ

化学を機械学習で何か

一日目　Xenonpy内の記述子関連クラスについて

Xenonpy

動作確認

Jupyter notebookをvisual studio codeの中で利用すると、わざわざコマンドで実行しなくても良いので便利です。

atmarkit.itmedia.co.jp

最低限の機能があるので、利用には必要十分だと思います。早速、Xenonpyをインストールしたので、github上のサンプルに従って実行してみます。最初はimportして、versionを確認する、RDkitが動くかどうかも確認してみます。

import xenonpy
from rdkit import Chem
from rdkit.Chem import Draw
xenonpy.__version__

'0.6.7'

mol=Chem.MolFromSmiles("C1CCOC1")
Draw.MolToImage(mol)

しっかりと実行できたと思います。最後に構造も描画できていましたね。

iQSPR.ipynb:記述子クラス関連

次は、Xenonpyのgithubサイト(https://github.com/yoshida-lab/XenonPy)のsampleフォルダにあるiQSPR.ipynbを使って実行し、コードを理解した上で、将来的には改変していこうと思います。

github.com

最初に、tools.ipynbを実行します。こちらには、実行に必要なライブラリのimport文やmatplotlibによる定型プロッターが準備されております。実行するために、sampleフォルダ内のtools.ipynbをこれから作るコードと同じ場所に格納してください。

run tools.ipynb

次に、トレーニング用のサンプルデータベースをダウンロードします。

https://github.com/yoshida-lab/XenonPy/releases/download/v0.4.1/iQSPR_sample_data.csv

このデータには、有機化合物の化学構造を表現したSMILES文字列と、内部エネルギー(internal energy E(kJ/mol))と、バンドギャップ（HOMO-LUMO gap (eV, or eV/particle))が格納されています。これは量子化学計算パッケージであるGAMESSを用いて計算されたもので、用いた理論はB3LYP/6-31G(d)です。データ数として16674個あるので、これを用いてANN/GcNNなどの学習を行って転移学習用途にも使えるかもしれません。

このブログでは、エネルギーの単位が違うので、eVに統一したいと思います。

また、未取得データがある場合にNaNという文字列を与えれば、Xenonpy上で無視してくれたりする機能があるみたいです。

ダウンロードしたcsvはここで述べるコードと同じフォルダに格納してあるものとします。

# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)

# take a look at the data
data.columns
data.head()

1 KJ/mol = 0.0104 eV/particleの関係があるので、次のようにします。

data["E"]=data["E"]*0.0104
data

これで、物性値として使っている値の単位があったので、数字の大きさが似通ってきました。次に、データベース内の内容を理解するために、matplotlibのhistコマンドを用いて、SMILESの文字数に対する分布グラフを描いてみましょう。ピークがおおよそ20-30くらいなので、iQSPR実行時の文字列として20文字くらいあれば十分かなという目算を行うことができます。ここで、ピークの位置を求めてみました。そうすると、おおよそ23.2くらいにあることがわかります。binの数を変えてもほとんど値は変わりません。(一部誤りがあり、2023/03/19修正)

count = [len(x) for x in data['SMILES']]
y,x,_=plt.hist(count, bins=20)
y=list(y)
x=list(x)
print("Peak position of SMILES length histogram"
      ,x[y.index(max(y))+1])

二つの物性値が格納されているので、どのようにデータが分布しているか調べめてみましょう。

# check target properties: E & HOMO-LUMO gap
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(data['E'],data['HOMO-LUMO gap'],s=15,c='b',alpha = 0.1)
_ = plt.title('iQSPR sample data')
_ = plt.xlabel('Internal energy')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.show()

内部エネルギーは0-6eV、バンドギャップは0-18eVくらいに分布していることがわかります。薄い色のところにデータがない領域があることに着目しましょう。sampleフォルダのコードでは、内部エネルギーが低く、バンドギャップも低いエネルギーを返す、SMILES文字列を探しています。

記述子クラス

XenonpyにはRDKit由来のmolecular finger printを利用するためのラッパーがあります。チュートリアルに従い、ECFPとMACCSを連結したFinger printを用いて記述子クラスを定義しています。サンプルとして下記のようなものが与えられています。

from xenonpy.descriptor import Fingerprints
RDKit_FPs = Fingerprints(featurizers=['ECFP', 'MACCS'], input_type='smiles')

finger printに関する脱線

上記までの簡素な方法で記述子を得ることができます。この記述子を自分の好きなように変えたいとなると、これだけでは情報が足りません。

Xenonpyのソースコードを確認し、Fingerprintsとして何が定義されているのかみてみましょう。

github.com


コードの冒頭部分に記載されているように、RDkit由来で利用できるFingerprintは以下のものです。

  • 'RDKitFP'
  • 'AtomPairFP'
  • 'TopologicalTorsionFP'
  • 'MACCS'
  • 'FCFP'
  • 'ECFP'
  • 'PatternFP'
  • 'LayeredFP'
  • 'MHFP'

Fingerprintsクラスは750行目くらいから記載されています。

class Fingerprints(BaseDescriptor):
    """
    Calculate fingerprints or descriptors of organic molecules.
    Note that MHFP currently does not support parallel computing, so n_jobs is fixed to 1.
    """

利用できる主要な引数として、以下のものがあります（詳しくはソースコードの英語を読んでください）。

  • n_jobs:計算するスレッド数
  • radius:CircularタイプのFinger printであるMorgan finger print向けの半径
  • n_bits:平たくいうとFingerprintの要素数
  • counting:0/1の値を要素数に入れる場合はFalse
  • featurizers:デフォルトは全てのfingerprintを考慮
  • input_type:SMILES記法を入力とするのか、RDKitのmolオブジェクト（rdkit.Chem.rdchem.Mol）をそのまま入力とするのか。

これがわかったら、次のようにfinger printを設定することもできます。

ECFP6_FPs = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=128*4,counting=True, input_type='smiles')

このクラスを理解するにはBaseDescriptorから理解する必要があります。

https://github.com/yoshida-lab/XenonPy/blob/master/xenonpy/descriptor/base.py

from sklearn.base import TransformerMixin, BaseEstimator
class BaseDescriptor(BaseEstimator, TransformerMixin, metaclass=TimedMetaClass):

このBaseDescriptorクラスを利用することで、各featurizerをタスクリストのように格納して、Xenonpy内で利用できるようです。このクラスのtransform関数を用いて、SMILESなどから、finger printを算出することができるようになります。

ECFPの例の場合、

https://github.com/yoshida-lab/XenonPy/blob/master/xenonpy/descriptor/fingerprint.py

class ECFP(BaseFeaturizer):

    def __init__(self, n_jobs=-1, *, radius=3, n_bits=2048, counting=False,
                 input_type='mol', on_errors='raise', return_type='any', target_col=None):
        """
        Morgan (Circular) fingerprints (ECFP)
        The algorithm used is described in the paper Rogers, D. & Hahn, M. Extended-Connectivity Fingerprints.
        JCIM 50:742-54 (2010)
        Parameters

が呼ばれます。

さらに、クラス継承元のBaseFeaturizerを見ると、

https://github.com/yoshida-lab/XenonPy/blob/master/xenonpy/descriptor/base.py

このクラスのtransform内で変換が行われており、結果としてECFPのfeaturize関数を呼び出している部分があります。ここまでわかったら、SMILESから、ある数値リストを出力するような関数を作り、Xenonpyの中に組み込むことは簡単です。次に、DummyFuncと、DummyClassを定義して使うこともできます。この仕組みを利用することで、ANN/GcNNによる学習モデルから転移学習させる際にDummyFuncに組み込んでやれば、Xenonpy向けの記述子クラスを定義できるようになります。以下に実装例を記載します。適当に作ったので、使用する際にはしっかりと作り込むことが必要です。

from xenonpy.descriptor.base import BaseDescriptor, BaseFeaturizer
def DummyFunc(x, nBits=1):
    # mol オブジェクトからcanonical smilesに直して、n_bitsの数だけの長さの数列を返す無意味な関数。
    dummy=Chem.MolToSmiles(x)
    return np.ones(nBits)*len(dummy)

class DummyClass(BaseFeaturizer):
    def __init__(self, n_jobs=-1, *, radius=3, n_bits=2048, counting=False,
                 input_type='mol', on_errors='raise', return_type='any', target_col=None):
        """
        """
        super().__init__(n_jobs=n_jobs, on_errors=on_errors, return_type=return_type, target_col=target_col)
    def featurize(self, x):
        if self.input_type == 'smiles':
            x_ = x
            x = Chem.MolFromSmiles(x)
            if x is None:
                raise ValueError('cannot convert Mol from SMILES %s' % x_)
        if self.input_type == 'any':
            if not isinstance(x, Chem.rdchem.Mol):
                x_ = x
                x = Chem.MolFromSmiles(x)
                if x is None:
                    raise ValueError('cannot convert Mol from SMILES %s' % x_)
        if self.counting:
            return count_fp(dummyFunc(x,nBits=self.n_bits) )
        else:
            return list(dummyFunc(x,nBits=self.n_bits) )

    @property
    def feature_labels(self):
        if self.counting:
            return [f'Dummy_c:' + str(i) for i in range(self.n_bits)]
        else:
            return [f'Dummy:' + str(i) for i in range(self.n_bits)]

# prepare descriptor function from XenonPy (possible to combine multiple descriptors)
class DummyDesc(BaseDescriptor):
    def __init__(self, n_jobs=-1, on_errors='nan'):
        super().__init__()
        self.n_jobs = n_jobs
        self.rdkit_fp = DummyClass(n_jobs, on_errors=on_errors, input_type='smiles')

Desc = DummyDesc()

記述子の適用

csvファイル内に"FBr(F)(F)(F)F"があるので、このデータは抜く必要があります。最初の記述子に話を戻して、SMILES文字列から、記述子を算出します。

data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
tmpFP = ECFP6_FPs.transform(data['SMILES'])
tmpFP


まとめ

記述子関連のやれたことを短くまとめると、下記になります。

run tools.ipynb

from xenonpy.descriptor import Fingerprints
# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
count = [len(x) for x in data['SMILES']]
y,x,_=plt.hist(count, bins=20)
y=list(y)
x=list(x)
print("Peak position of SMILES length histogram"
      ,x[y.index(max(y))+1])

# Convert energy unit
data["E"]=data["E"]*0.0104

# check target properties: E & HOMO-LUMO gap
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(data['E'],data['HOMO-LUMO gap'],s=15,c='b',alpha = 0.1)
_ = plt.title('iQSPR sample data')
_ = plt.xlabel('Internal energy')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.show()

# Descriptor
ECFP6_FPs = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=128*4,counting=True, input_type='smiles')
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
tmpFP = ECFP6_FPs.transform(data['SMILES'])