機械学習メモ

化学を機械学習で何か

五日目　iQSPR

Xenonpy

iQSPR

今回の検討で分かったことがある。

「SMILESの微小改変に追従するぐらい予測モデルの値が変化しないとXenonPyによる構造改変が機能しない。」

少ない記述子で十分に予想できるモデルを作成したとしても、iQSPRを実行するとターゲット領域の構造を探索しきれないことに繋がってしまう。これは記述子だけではなく、学習モデル自体に対しても要求されるはずだ。

今回は比較的手軽な機械学習モデルを用いているが、Fingerprint・機械学習モデルの組み合わせによっては、iQSPRの探索範囲が限定される結果をもたらす可能性がある。そのため、化学構造探索数・精度の観点で、Fingerprintと機械学習モデルの最適化を行うのが、実用上望まれるであろう。

16から256にn_Bitsを変更することで、標的領域に化学構造が現れるようになった。

まず最初に、作成したN-gramを読み込む。

prd_mdls=like_mdl

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# load a pre-trained n-gram from the pickle file
n_gram = joblib.load('ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'rb') as f:
#     n_gram = pk.load(f)
    
# load a pre-trained n-gram from the pickle file
n_gram2 = joblib.load('ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'rb') as f:
#     n_gram2 = pk.load(f)

_ = n_gram.merge_table(n_gram2)

次に、初期の探索点を選定する。

# set up initial molecules for iQSPR
np.random.seed(201906) # fix the random seed
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi)) for i, smi in enumerate(data['SMILES'])
       if (data['HOMO-LUMO gap'].iloc[i] > 3.5)]
init_samples = np.random.choice(cans, 25)
print(init_samples)

また、アニーリング（書籍参照のこと）パラメーターを設定する。

# set up annealing schedule in iQSPR
beta = np.hstack([np.linspace(0.01,0.2,20),np.linspace(0.21,0.4,10),np.linspace(0.4,1,10),np.linspace(1,1,10)])
print('Number of steps: %i' % len(beta))
print(beta)

ここから、N-gramとLogLikelihood関数をIQSPRクラスに設定して、探索を行う。

# library for running iQSPR in XenonPy-iQSPR
from xenonpy.inverse.iqspr import IQSPR

# update NGram parameters for this exampleHOMO-LUMO gap
n_gram.set_params(del_range=[1,20],max_len=500, reorder_prob=0.5, sample_order=(1,20))

# set up likelihood and n-gram models in iQSPR
iqspr_reorder = IQSPR(estimator=prd_mdls, modifier=n_gram)
    
# main loop of iQSPR
iqspr_samples1, iqspr_loglike1, iqspr_prob1, iqspr_freq1 = [], [], [], []
for s, ll, p, freq in iqspr_reorder(init_samples, beta, yield_lpf=True):
    iqspr_samples1.append(s)
    iqspr_loglike1.append(ll)
    iqspr_prob1.append(p)
    iqspr_freq1.append(freq)
# record all outputs
iqspr_results_reorder = {
    "samples": iqspr_samples1,
    "loglike": iqspr_loglike1,
    "prob": iqspr_prob1,
    "freq": iqspr_freq1,
    "beta": beta
}
# save results
with open('iQSPR_results_reorder.obj', 'wb') as f:
    pk.dump(iqspr_results_reorder, f)

探索する際に、likelihood関数がどのような構造が得られたのか、履歴を確認する。Target領域近傍の点がサンプルできた。

# plot the likelihood evolution

# set up the min and max boundary for the plots
tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]
flat_list = np.asarray([item for sublist in tmp_list for item in sublist])
y_max, y_min = max(flat_list), min(flat_list)

_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Log-likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), ll ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_loglike_reorder.png',dpi = 500)
#plt.close()

y_max, y_min = np.exp(y_max), np.exp(y_min)
_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), np.exp(ll) ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_like_reorder.png',dpi = 500)
#plt.close()


探索された点に対して、改めて予測モデルを適用してどの程度の予測ブレがあるのか調べる。結構大きい結果となったので、境界にいる構造は外れている可能性もある。

# re-calculate the property values for the proposed molecules
x_mean, x_std, y_mean, y_std = [], [], [], []
r_std = []
FPs_samples = []
for i, smis in enumerate(iqspr_results_reorder["samples"]):
    tmp_fps = FP_E.transform(smis)
    FPs_samples.append(tmp_fps)
    
    tmp1, tmp2 = custom_mdls1["E"].predict(tmp_fps)
    x_mean.append(tmp1)
    x_std.append(tmp2)
    
    tmp_fps1 = FP_H.transform(smis)
    tmp1, tmp2 = custom_mdls2["HOMO-LUMO gap"].predict(tmp_fps1)
    y_mean.append(tmp1)
    y_std.append(tmp2)
    
    r_std.append([np.sqrt(x_std[-1][i]**2 + y_std[-1][i]**2) for i in range(len(x_std[-1]))])

# flatten the list for max/min calculation
flat_list = [item for sublist in r_std for item in sublist]
print('Range of std. dev.: (%.4f,%.4f)' % (min(flat_list),max(flat_list)))

Range of std. dev.: (0.7748,0.8449)

探索履歴を、二つの物性値のグラフ上で動かした履歴を確認できるようにする。

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_prd/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

flat_list = np.asarray([item for sublist in r_std for item in sublist])
s_max, s_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["E"],
    np.asarray([item for sublist in x_mean for item in sublist])))
x_max, x_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["HOMO-LUMO gap"],
    np.asarray([item for sublist in y_mean for item in sublist])))
y_max, y_min = max(flat_list), min(flat_list)
tmp_beta = iqspr_results_reorder["beta"]

for i in range(len(r_std)):
    dot_size = 45*((np.asarray(r_std[i])-s_min)/(s_max-s_min)) + 5
    
    _ = plt.figure(figsize=(5,5))
    rectangle = plt.Rectangle((0,0),2.5,3,fc='y',alpha=0.1)
    _ = plt.gca().add_patch(rectangle)
    _ = plt.scatter(data["E"], data["HOMO-LUMO gap"],s=3, c='b', alpha=0.2)
    _ = plt.scatter(x_mean[i], y_mean[i],s=dot_size, c='r', alpha=0.5)
    _ = plt.title('Step: %i (beta = %.3f)' % (i,tmp_beta[i]))
    _ = plt.xlim(x_min,x_max)
    _ = plt.ylim(y_min,y_max)
    _ = plt.xlabel('Internal energy')
    _ = plt.ylabel('HOMO-LUMO gap')
    #plt.show()
    _ = plt.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

化学構造の変遷を見てみる。

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_smiles/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]    

n_S = 25
for i, smis in enumerate(iqspr_results_reorder['samples']):
    tmp_smis = iqspr_results_reorder['samples'][i][
        np.argsort(tmp_list[i])[::-1]]
    fig, ax = plt.subplots(5, 5)
    _ = fig.set_size_inches(20, 20)
    _ = fig.set_tight_layout(True)
    for j in range(n_S):
        xaxis = j // 5
        yaxis = j % 5
        try:
            img = Draw.MolToImage(Chem.MolFromSmiles(tmp_smis[j]))
            _ = ax[xaxis, yaxis].clear()
            _ = ax[xaxis, yaxis].set_frame_on(False)
            _ = ax[xaxis, yaxis].imshow(img)
        except:
            pass
        _ = ax[xaxis, yaxis].set_axis_off()
    _ = fig.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

ターゲット領域の設定からして、分子量が小さくて、電子共鳴しやすい構造であるなというのが想像つくのだが、おおよそそうなっていたので、うまく行ったということだろう。

まとめ

iQSPRについては、ほとんど寄り道がなかったので、まとめは実施しない。コード自体はチュートリアルの域を出なかったが、パラメーターを色々変えてみた場合にどうなるか、実施してみることで、次の応用への理解が深まった。

github.com


"Thank you for using XenonPy-iQSPR. We would appreciate any feedback and code contribution to this open-source project.The sample code can be downloaded at: https://github.com/yoshida-lab/XenonPy/blob/master/samples/iQSPR.ipynb"

結局のところ、フィードバックを与えるような試みもしてはいなかった。

likelihoood関数とか記述子くらいか。一次元探索の問題にしたいときは、ダミーの物性予測（RDKitで算出できる分子量とか）をもう一次元に設定して、予測モデルもRDKitの関数にするなりして、範囲も(0.0,np.inf)とかにすれば良さそう。

興味を持った読者は、自身で実行しながら、githubサイト、書籍を読んではいかがだろうか。サンプルコードも寄り道をしながらわかりやすいように説明してくれているので、XenonPyで用意されたフレームワークが柔軟に設計されていることは大変ありがたい。記述子周りも、すっきりさせることができるので、表現方法探索をするのにもこのコードを応用してやれば、泥臭いコードでエラーが頻発するなんてことも提言できそうなので、仕事にも役立つかもしれない。

書籍版のiQSPR用のコードでは、学習モデルのハイパーパラメーター探索もしていたようなので、後で確認してみようと思う。

最後にコード全体がどうなっていたのかを再表示してみる。

run tools.ipynb

from xenonpy.descriptor import Fingerprints
# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
# make histogram
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()

# set target range
target_range = {'E': (0,2.5), 'HOMO-LUMO gap': (-np.inf, 3)}

# prepare indices for cross-validation data sets
from xenonpy.datatools import Splitter
sp = Splitter(data.shape[0], test_size=0, k_fold=10)

# set property list
prop=list(data.columns)
prop.remove("SMILES")

from copy import deepcopy
from sklearn.metrics import r2_score

# initialize model dictionary
mdls = {key: [] for key in prop}

def CrossValidation(sp,propid,Desc,model):
    y_trues_test=[]
    y_preds_test=[]
    y_trues_train=[]
    y_preds_train=[]
    out_mdls=[]
    # cross-validation test
    for iTr, iTe in sp.cv():

        x_train = data['SMILES'].iloc[iTr]
        x_test = data['SMILES'].iloc[iTe]
    
        fps_train = Desc.transform(x_train)
        fps_test = Desc.transform(x_test)
    
        y_train = data[prop].iloc[iTr]
        y_test = data[prop].iloc[iTe]

        mdl = deepcopy( model )
        
        mdl.fit(fps_train, y_train.iloc[:,propid])
        prd_train = mdl.predict(fps_train)
        prd_test = mdl.predict(fps_test)

        y_trues_test.append(y_test.iloc[:,propid].values)
        y_trues_train.append(y_train.iloc[:,propid].values)
    
        y_preds_test.append(prd_test)
        y_preds_train.append(prd_train)
        out_mdls.append(mdl)
    
    y_true_test = np.concatenate(y_trues_test)
    y_pred_test = np.concatenate(y_preds_test)
    y_true_train = np.concatenate(y_trues_train)
    y_pred_train = np.concatenate(y_preds_train)
    xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
    xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
    xy_diff = xy_max - xy_min

    r2=r2_score(y_true_test,y_pred_test)

    _ = plt.figure(figsize=(5,5))
    _ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
    _ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
    _ = plt.title('Property: ' + Target)
    _ = plt.xlim(xy_min,xy_max)
    _ = plt.ylim(xy_min,xy_max)
    _ = plt.legend(loc='upper left')
    _ = plt.xlabel('Prediction')
    _ = plt.ylabel('Observation')
    _ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
    _ = plt.show()
    return out_mdls

from sklearn.ensemble import GradientBoostingRegressor
Target="E"
FP_E = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=256,counting=True, input_type='smiles')
model=GradientBoostingRegressor(loss='squared_error')
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_E,model)

from sklearn.ensemble import RandomForestRegressor
Target="HOMO-LUMO gap"
FP_H= Fingerprints(featurizers=['PatternFP',"ECFP","MACCS",'LayeredFP'],radius=10,n_bits=256,counting=True, input_type='smiles')
model=RandomForestRegressor(100, n_jobs=-1, max_features=None)
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_H,model)


# Forward model template in XenonPy-iQSPR 
from xenonpy.inverse.iqspr import GaussianLogLikelihood

class bootstrap_fcn():
    def __init__(self, ms, var):
        self.Ms = ms
        self.Var = var
    def predict(self, x):
        val = np.array([m.predict(x) for m in self.Ms])
        return np.mean(val, axis = 0), np.sqrt(np.var(val, axis = 0) + self.Var)
    
# include basic measurement noise
c_var = {'E': 0.3, 'HOMO-LUMO gap': 0.3}

# generate a dictionary for the final models used to create the likelihood
custom_mdls1 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="E"}
custom_mdls2 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="HOMO-LUMO gap"}

# import descriptor calculator and forward model to iQSPR
prd_mdls1 = GaussianLogLikelihood(descriptor=FP_E, targets={'E': target_range["E"]}, **custom_mdls1)
prd_mdls2 = GaussianLogLikelihood(descriptor=FP_H, targets={'HOMO-LUMO gap': target_range['HOMO-LUMO gap']}, **custom_mdls2)


from xenonpy.inverse.base import BaseLogLikelihoodSet
from sklearn.linear_model import BayesianRidge

custom_mdls = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() }
desc={"E":FP_E,'HOMO-LUMO gap':FP_H}

class MyLogLikelihood(BaseLogLikelihoodSet):
    def __init__(self):
        super().__init__()
        self.loglike = prd_mdls1
        self.loglike = prd_mdls2
        self._mdl=custom_mdls
        self.descriptor=desc
        self.targets=target_range
    def predict(self, smiles, **kwargs):
        tmp = {}
        for k, v in self._mdl.items():
            fps = self.descriptor[k].transform(smiles, return_type='df')
            fps_ = fps.dropna()
            tmp[k + ': mean'], tmp[k + ': std'] = v.predict(fps_)
        tmp = pd.DataFrame(data=tmp, index=fps_.index)
        return pd.DataFrame(data=tmp, index=fps.index)

like_mdl = MyLogLikelihood()

pred=like_mdl.predict(data['SMILES'])

# calculate log-likelihood for a given target property region
tmp_ll = like_mdl(data['SMILES'])
tmp = tmp_ll.sum(axis = 1, skipna = True)

_ = plt.figure(figsize=(8,5))
_ = plt.hist(tmp, bins=50)
_ = plt.title('ln(likelihood) for SMILES in the data subset')
_ = plt.xlabel('ln(likelihood)')
_ = plt.ylabel('Counts')
_ = plt.show()

# plot histogram of likelihood values
_ = plt.figure(figsize=(8,5))
_ = plt.hist(np.exp(tmp), bins=50)
_ = plt.title('Likelihood for SMILES in the data subset')
_ = plt.xlabel('Likelihood')
_ = plt.ylabel('Counts')
_ = plt.show()

# check the predicted likelihood
dot_scale = 0.5
l_std = np.sqrt(pred['E: std']**2+pred['HOMO-LUMO gap: std']**2)

_ = plt.figure(figsize=(6,5))
rectangle = plt.Rectangle((0,0),2.5,3,fc='y',alpha=0.1)
_ = plt.gca().add_patch(rectangle)
im = plt.scatter(pred['E: mean'], pred['HOMO-LUMO gap: mean'], s=l_std*dot_scale, c=np.exp(tmp),alpha = 0.2,cmap=plt.get_cmap('cool'))
_ = plt.title('Pubchem data')
_ = plt.xlabel('E')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.colorbar(im)
_ = plt.show()

prd_mdls=like_mdl

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# load a pre-trained n-gram from the pickle file
n_gram = joblib.load('ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'rb') as f:
#     n_gram = pk.load(f)
    
# load a pre-trained n-gram from the pickle file
n_gram2 = joblib.load('ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'rb') as f:
#     n_gram2 = pk.load(f)

_ = n_gram.merge_table(n_gram2)

# set up initial molecules for iQSPR
np.random.seed(201906) # fix the random seed
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi)) for i, smi in enumerate(data['SMILES'])
       if (data['HOMO-LUMO gap'].iloc[i] > 3.5)]
init_samples = np.random.choice(cans, 25)
print(init_samples)

# set up annealing schedule in iQSPR
beta = np.hstack([np.linspace(0.01,0.2,20),np.linspace(0.21,0.4,10),np.linspace(0.4,1,10),np.linspace(1,1,10)])
print('Number of steps: %i' % len(beta))
print(beta)

# library for running iQSPR in XenonPy-iQSPR
from xenonpy.inverse.iqspr import IQSPR

# update NGram parameters for this exampleHOMO-LUMO gap
n_gram.set_params(del_range=[1,20],max_len=500, reorder_prob=0.5, sample_order=(1,20))

# set up likelihood and n-gram models in iQSPR
iqspr_reorder = IQSPR(estimator=prd_mdls, modifier=n_gram)
    
# main loop of iQSPR
iqspr_samples1, iqspr_loglike1, iqspr_prob1, iqspr_freq1 = [], [], [], []
for s, ll, p, freq in iqspr_reorder(init_samples, beta, yield_lpf=True):
    iqspr_samples1.append(s)
    iqspr_loglike1.append(ll)
    iqspr_prob1.append(p)
    iqspr_freq1.append(freq)
# record all outputs
iqspr_results_reorder = {
    "samples": iqspr_samples1,
    "loglike": iqspr_loglike1,
    "prob": iqspr_prob1,
    "freq": iqspr_freq1,
    "beta": beta
}
# save results
with open('iQSPR_results_reorder.obj', 'wb') as f:
    pk.dump(iqspr_results_reorder, f)


with open('iQSPR_results_reorder.obj', 'rb') as f:
    iqspr_results_reorder = pk.load(f)

# plot the likelihood evolution

# set up the min and max boundary for the plots
tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]
flat_list = np.asarray([item for sublist in tmp_list for item in sublist])
y_max, y_min = max(flat_list), min(flat_list)

_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Log-likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), ll ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_loglike_reorder.png',dpi = 500)
#plt.close()

y_max, y_min = np.exp(y_max), np.exp(y_min)
_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), np.exp(ll) ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_like_reorder.png',dpi = 500)
#plt.close()



# re-calculate the property values for the proposed molecules
x_mean, x_std, y_mean, y_std = [], [], [], []
r_std = []
FPs_samples = []
for i, smis in enumerate(iqspr_results_reorder["samples"]):
    tmp_fps = FP_E.transform(smis)
    FPs_samples.append(tmp_fps)
    
    tmp1, tmp2 = custom_mdls1["E"].predict(tmp_fps)
    x_mean.append(tmp1)
    x_std.append(tmp2)
    
    tmp_fps1 = FP_H.transform(smis)
    tmp1, tmp2 = custom_mdls2["HOMO-LUMO gap"].predict(tmp_fps1)
    y_mean.append(tmp1)
    y_std.append(tmp2)
    
    r_std.append([np.sqrt(x_std[-1][i]**2 + y_std[-1][i]**2) for i in range(len(x_std[-1]))])

# flatten the list for max/min calculation
flat_list = [item for sublist in r_std for item in sublist]
print('Range of std. dev.: (%.4f,%.4f)' % (min(flat_list),max(flat_list)))


import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_prd/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

flat_list = np.asarray([item for sublist in r_std for item in sublist])
s_max, s_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["E"],
    np.asarray([item for sublist in x_mean for item in sublist])))
x_max, x_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["HOMO-LUMO gap"],
    np.asarray([item for sublist in y_mean for item in sublist])))
y_max, y_min = max(flat_list), min(flat_list)
tmp_beta = iqspr_results_reorder["beta"]

for i in range(len(r_std)):
    dot_size = 45*((np.asarray(r_std[i])-s_min)/(s_max-s_min)) + 5
    
    _ = plt.figure(figsize=(5,5))
    rectangle = plt.Rectangle((0,0),2.5,3,fc='y',alpha=0.1)
    _ = plt.gca().add_patch(rectangle)
    _ = plt.scatter(data["E"], data["HOMO-LUMO gap"],s=3, c='b', alpha=0.2)
    _ = plt.scatter(x_mean[i], y_mean[i],s=dot_size, c='r', alpha=0.5)
    _ = plt.title('Step: %i (beta = %.3f)' % (i,tmp_beta[i]))
    _ = plt.xlim(x_min,x_max)
    _ = plt.ylim(y_min,y_max)
    _ = plt.xlabel('Internal energy')
    _ = plt.ylabel('HOMO-LUMO gap')
    #plt.show()
    _ = plt.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_smiles/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]    

n_S = 25
for i, smis in enumerate(iqspr_results_reorder['samples']):
    tmp_smis = iqspr_results_reorder['samples'][i][
        np.argsort(tmp_list[i])[::-1]]
    fig, ax = plt.subplots(5, 5)
    _ = fig.set_size_inches(20, 20)
    _ = fig.set_tight_layout(True)
    for j in range(n_S):
        xaxis = j // 5
        yaxis = j % 5
        try:
            img = Draw.MolToImage(Chem.MolFromSmiles(tmp_smis[j]))
            _ = ax[xaxis, yaxis].clear()
            _ = ax[xaxis, yaxis].set_frame_on(False)
            _ = ax[xaxis, yaxis].imshow(img)
        except:
            pass
        _ = ax[xaxis, yaxis].set_axis_off()
    _ = fig.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

次はiQSMDとnn_modelで学習モデルを作って、transfer_learningを試してみるかな。

四日目　XenonPyを使ったモデル作成

Xenonpy

モデル作成

XenonPyの記述子クラスと、これまで作成した交差検証関数をつかって、モデルを作成し、HOMO-LUMOギャップに対する予測精度を上げる方法を探索した。

目指す目的は、短い学習時間で最良の結果を得ることである。

人間の頭でパラメーターを変更して行った結果、誤差・決定係数ともに改善させることに成功した。学習モデルはRadon Forestを用いている。(はてなの仕様だと思うけど、表を横にスクロールしてください。R2とMAE載せてますよ)


FP radius n_bits n_estimators max_features R2 MAE
PatternFP',"MACCS" 2048 500 None 0.92 0.24
PatternFP',"MACCS" 128 500 None 0.91 0.27
PatternFP',"MACCS" 16 500 sqrt 0.90 0.29
PatternFP',"ECFP","MACCS",'LayeredFP' 10 16 100 None 0.90 0.30
PatternFP',"MACCS" 16 500 sqrt 0.89 0.33
PatternFP',"ECFP","MACCS",'LayeredFP' 10 16 10 None 0.89 0.31
PatternFP',"MACCS" 16 10 sqrt 0.87 0.36
LayeredFP',"MACCS" 16 10 sqrt 0.87 0.36
AtomPairFP',"MACCS" 16 10 Sqrt 0.86 0.38
TopologicalTorsionFP',"MACCS" 16 10 sqrt 0.86 0.37
ECFP',"MACCS" 3 16 10 Sqrt 0.85 0.39
FCFP',"MACCS" 3 16 10 sqrt 0.85 0.38
MHFP',"MACCS" 16 10 sqrt 0.85 0.38
RDKitFP',"MACCS" 16 10 sqrt 0.84 0.40


記述子の要素数を2048にし、RadomForestのn_estimatorsの値を500にし、max_featuresをNoneにすると、かなり学習に時間がかかってしまう。

２分程度で終わる条件は下記となった。

PatternFP',"ECFP","MACCS",'LayeredFP' 10 16 100 None 0.90 0.30

つまり、学習モデルを単純化させた場合、記述子側をある程度複雑にする必要がある。これは下記のURLにて最後に書かれている言葉と同じであるが、どれだけ複雑にする必要があるのか、必要条件を得ることは現時点ではできない。

chemintelligence.com

RDKitやMordredなどに頼らずに、別の方法で記述子を作成する方法は、量子化学計算を毎回実行し、記述子を得ること、こういう方法もあるだろう。しかし、逆解析や多数の未知化合物を予測する場合、計算時間が足枷となってしまう問題がある。

別の方法も存在しつつあり、DScribeなども探索範疇に入れてもいいかもしれないし、もっと良いものがあるかもしれない。

https://www.sciencedirect.com/science/article/pii/S0010465519303042

DScribe is a Python package for transforming atomic structures into fixed-size numerical fingerprints. These fingerprints are often called “descriptors” and they can be used in various tasks, including machine learning, visualization, similarity analysis, etc. To get started you can check the basic tutorial.

singroup.github.io

トレンドとしてはGraph convolutional neural networkがよく使われることが多いだろうし、転移学習もあるだろう。

RandomForestではなく、Gradientboostを使うと高エネルギー側で記述が格段に悪くなる。

lightgbmをXenonPyで使うために、wrapperを準備してやれば、利用しやすい。

import lightgbm as lgb
Target="HOMO-LUMO gap"
FP_H= Fingerprints(featurizers=['PatternFP',"ECFP","MACCS",'LayeredFP'],radius=10,n_bits=32,counting=True, input_type='smiles')

class lightgbm_class():
    def __init__(self):
        self.params= {"objective": "regression","metrics": "mae",}  
        return
    def fit(self,train,teacher):
        lgb_train = lgb.Dataset(train.to_numpy(), teacher)
        self.model=lgb.train(self.params,lgb_train)
    def predict(self,data):
        return self.model.predict(data)
      

model=lightgbm_class()
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_H,model)

xgboostも同じようなノリで記載はできるはずだ。lightgbmでやった結果、やはり、高エネルギー側で記述が悪い結果となってしまった。

ChatGptで

「Xenonpyを使ったiQSPRコードを生成して」と聞いて生成されたコード。多分違う。

import pandas as pd
from xenonpy.descriptor import Compositions
from xenonpy.model import XGBRegressor
from xenonpy.utils import settings
from xenonpy.datatools import preset

# データの読み込み
data = pd.read_csv('solubility_data.csv')

# 分子記述子の生成
comp = Compositions(n_jobs=-1)
X = comp.transform(data['smiles'])

# モデルの構築
settings.general.seed = 1
model = XGBRegressor()
model.fit(X, data['logS'])

# モデルの評価
preset.regression(model, X, data['logS'])

これではダメなので、
「Xenonpyを使って、ElasticNetでモデリングして、LogLikelihood関数を作成して、N-gramを生成して、iQSPRで逆解析して、化学構造を得るコードを書いて」
と適当なことを質問してみた。

import pandas as pd
from xenonpy.descriptor import Fingerprints
from xenonpy.model import ElasticNet
from xenonpy.utils import settings
from xenonpy.datatools import preset

# データの読み込み
data = pd.read_csv('solubility_data.csv')

# Fingerprintsの生成
fp = Fingerprints(n_jobs=-1, type='morgan', n_bits=512, length=2)
X = fp.transform(data['smiles'])

# ElasticNetのモデルの構築
settings.general.seed = 1
model = ElasticNet()
model.fit(X, data['logS'])

# モデルの評価
preset.regression(model, X, data['logS'])

# LogLikelihood関数の作成
from xenonpy.inverse.iqsp import LogLikelihood
llh = LogLikelihood(model, fp, data['logS'], scale='standard')
n_gram = llh.get_ngram(5, 0.1)

# iQSPRで逆解析して化学構造を得る
from xenonpy.inverse.iqsp import iQSPR
iqspr = iQSPR(model, fp, scale='standard')
smiles = iqspr.inverse(n_gram)
print(smiles)

ChatGpy 3.0くらいなのかな、トライアル版だと。

Materials informaticsって

学習モデルについては、かなりな勢いで進化を遂げていて、深層学習は任意の写像を表現できる強力な手法でもある。

 

ニューラルネットワークと深層学習 (nnadl-ja.github.io)

「まず、この表現はニューラルネットワークは任意の関数を完全に計算できる事を意味するのではありません。 そうではなく、好きなだけ近い近似関数を得られるという意味です。 隠れ層のニューロンを増やすことで、近似の精度を上げることができます。」

 

他方で、記述子選びが大変で、一筋縄で行かない場合がほとんど。

 

Fingerprintを利用する問題点も書籍で指摘されていた通り。分子自体の記述子の導入という手も確かにあるだろう。それでも、一筋縄で行かない場合の方策があれば良いのだが、まだ出会ったことはない。

 

グラフ折りたたみニューラルネットワークも、良いのだけれど、最良の方法でもない場合がある。

【機械学習+分子設計】分子生成モデルの主要トレンド - ころがる狸 (dajiro.com)

 

現状考えうる、分子表現は以下のようにレビューが多数ある。

 

Molecular representations for machine learning applications in chemistry - Raghunathan - 2022 - International Journal of Quantum Chemistry - Wiley Online Library

 

A review of molecular representation in the age of machine learning - Wigh - 2022 - WIREs Computational Molecular Science - Wiley Online Library

 

Molecular representations in AI-driven drug discovery: a review and practical guide | Journal of Cheminformatics | Full Text (biomedcentral.com)

 

次はMOLGANも見てみるかな。

 

magic3007/MolGAN-pytorch: 🦑 Pytorch implementation of MolGAN: An implicit generative model for small molecular graphs. (github.com)

 

 

 

三日目　N-gramを作る。

Xenonpy

前回まで

dango0dan5.hatenablog.com


N-gram

ようやく、書籍を手にすることができた。iQSPRの基本的な概念が日本語で記載されているので、コードの理解にも役に立つ。Xenonpyはいろいろと便利なクラスが用意されているし、ドキュメントもコード内にも細かく書かれている。新しいことをしようとすると、ちゃんと調べてやらないと、かなり難しい。

www.kyoritsu-pub.co.jp

次に、N-gramを作成する。

まずN-gramの概念については、下記を参照のこと。

mojitoba.com

qiita.com

SMILES文字列の場合は、単語ではなく、文字で区切りを設けている。
推奨文字数として5文字以上が経験的に推奨されている。最大で10文字が許されるが、以下のようにすると、最大が5文字として設定される。

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# initialize a new n-gram
n_gram = NGram()

# train the n-gram with SMILES of available molecules
n_gram.fit(data['SMILES'],train_order=5)

NGram(ngram_table=[[[ C ( O [N+] ! & N Cl =C =O ... \
['C'] 26254 11168 2899 58 3496 13361 2666 237 1251 215 ...
[')'] 11471 2204 5707 488 0 0 3750 1333 0 4 ...
['O'] 4430 2 40 22 3926 0 39 1 0 1 ...
['[N+]'] 0 589 0 0 0 14 0 0 0 0 ...
['N'] 1937 1537 61 6 1609 940 135 1 247 116 ...
... ... ... ... ... ... ... ... ... ... ... ...
['=[Se]'] 0 1 0 0 0 0 0 0 0 1 ...
['[Ti]'] 0 1 0 0 0 0 0 0 0 0 ...
['=[Si]'] 0 0 0 0 0 0 0 0 0 1 ...
['[Cu]'] 0 0 0 0 0 0 0 0 0 0 ...
['[O+]'] 0 0 0 0 0 0 0 0 0...
['C', '(', '=C', 'C', 2] 1 0 0 0 0 0
['(', '=C', 'C', 2, '=C'] 0 0 0 1 0 0
['=C', 'C', 2, '=C', 'C'] 0 0 0 0 0 1
['C', '=C', 'C', '=C', 0] 0 0 0 1 0 0
['=C', 'C', '=C', 0, 'C'] 0 1 0 0 0 0,
=C C 2 0
['&', '=C', '&', 'C', '('] 1 0 0 0
['=C', '&', 'C', '(', '=C'] 0 1 0 0
['&', 'C', '(', '=C', 'C'] 0 0 1 0
[1, 'C', '&', '=C', '('] 0 1 0 0
['C', '&', '=C', '(', 'C'] 1 0 0 0
['&', '=C', '(', 'C', '=C'] 0 1 0 0
['=C', '(', 'C', '=C', 'C'] 1 0 0 0
['(', 'C', '=C', 'C', '=C'] 0 0 0 1]]],
sample_order=(1, 5))

データの最初の５個目までの構造で、N-gramを使って構造発生すると、以下のようなものが得られる。

# perform pure iQSPR molecule generation starting with 5 initial molecules
n_loop = 5
tmp = data['SMILES'][:5]
for i in range(n_loop):
    tmp = n_gram.proposal(tmp)
    print('Round %i' % i,tmp)

Round 0 ['CC(=O)OC(CC(=O)[O-])C', 'CC(=O)OC(CC(=O)O)C[N+](C)(C)N', 'C1=CC(C(C(=C1)C)(C(=O)O)CBr)O', 'C1COCCN1S(=O)(=O)OCCCl', 'C(C(=O)COP(=O)(O)O)O']
Round 1 ['CC(=O)OC(CC(=O)[O-])C1=CC=C(C=C1)CC=C', 'CC(=O)OC(CC(=O)O)C[N+](C)(C)C', 'C1=CC(C(C(=C1)C)(C(=O)O)O)C(=O)C=CC(=O)C', 'C1COCCN1S(=O)N1CCC2=CC=CC=C21', 'C(C(=O)COP(=O)(O)OP(=O)(O)OC1(C(N(C2=CC3=C(CCC4=CC(=O)CCC4(C3(C(CCC2C1)O)O)C)O)C1=CC=CC=C1)O)C)NC1=C(C=CC(=C1)O)Cl']
Round 2 ['CC(=O)OC(CC(=O)[O-])C1=CC=C(C=C1)[N+](=O)[O-]', 'CC(=O)OC(CC(=O)O)C(C)C1=CC=CC=C1O', 'C1=CC(C(C(=C1)C)(C(=O)O)O)C(=O)O', 'C1COCCN1S(=O)N1CCC2=CC=CC=C21', 'C(C(=O)COP(=O)(O)OP(=O)(O)OC1(C(N(C2=CC3=C(CCC4=CC(=O)CCC4(C3(C(CCC2C1)O)O)C)O)C1=CC=CC=C1)O)C)NC1=C(C=CC(=C1)OP(=O)(O)O)OC(=O)C=C(C#N)N=CC1=CC=C(C=C1)N']
Round 3 ['CC(=O)OC(CC(=O)[O-])C1=CC=C(C=C1)[N+](=O)[O-]', 'CC(=O)OC(CC(=O)O)C(C)C1=CC=CC=C1OC', 'C1=CC(C(C(=C1)C)(C(=O)O)O)C(=O)OC', 'C1COCCN1S(=O)N1CCC2=CC=CC=C2S1', 'C(C(=O)COP(=O)(O)OP(=O)(O)OC1(C(N(C2=CC3=C(CCC4=CC(=O)CCC4(C3(C(CCC2C1)O)O)C)O)C1=CC=CC=C1)O)C)NC1=C(C=CC(=C1)OP(=O)(O)O)OC(=O)C=C(C#N)N=CC1=CC=C(C=C1)C(=O)NC1=CC=CS1']
Round 4 ['CC(=O)OC(CC(=O)[O-])C1=CC=C(C=C1)OCC=CC1=CC2=C(C1=O)C(=C2)C', 'CC(=O)OC(CC(=O)O)C(C)C1=CC=C(C=C1)Br', 'C1=CC(C(C(=C1)C)(C(=O)O)O)O', 'C1COCCN1S(=O)N1CCC2=C(C1=O)C=CC(=C2O)Br', 'C(C(=O)COP(=O)(O)OP(=O)(O)OC1(C(N(C2=CC3=C(CCC4=CC(=O)CCC4(C3(C(CCC2C1)O)O)C)O)C1=CC=CC=C1)O)C)NC1=C(C=CC(=C1)OP(=O)(O)O)OC(=O)C=C(C#N)N=CC1=CC=C(C=C1)C(=O)OCC(=O)N(C1=CC=C(C=C1)Cl)C1=CC=C(C=C1)C1=CC=CC2=CC3=CC=CC=C3C=C12']

N-gramの挙動を調べるために、パラメーターを与えて違いを見てみている。最大長さを500文字として、末尾から原子を削除する個数（ランダムに可変する）の範囲を1-10とし、並び替えが行われる確率を0.5としている。同じことをしているのにも関わらず、違う構造が得られているのがわかる。

# change internal parameters of n_gram generator
_ = n_gram.set_params(del_range=[1,10], max_len=500, reorder_prob=0.5)

# perform pure iQSPR molecule generation starting with 10 PG molecules
n_loop = 5
tmp = data['SMILES'][:5]
for i in range(n_loop):
    tmp = n_gram.proposal(tmp)
    print('Round %i' % i,tmp)

Round 0 ['C[N+](C)(C)CC(CC(=O)[O-])OC(C1=CC=CC=C1)CC', 'O(C(C)=O)C(CC(=O)O)C[N+](C)(C)O', 'C1=CC=C(C(=O)NC(=O)NC1=O)[N+](=O)[O-]', 'C1=NC2=C(C1=O)C(=O)C=CC2=O', 'C(C(=O)COP(=O)(O)OC1CCN(CC1)CN1CCOCC1)C(=O)N']
Round 1 ['C[N+](C)(C)CC(CC(=O)[O-])OC(C1=CC=CC=C1)COP(=O)(O)OP(=O)(OC)O', 'O=C(C)OC(CC(=O)O)CNCCCN', 'O=C1NC(=O)NC(=O)C=CC=C1C(=O)OCC', 'C1=NC2=C(C1=O)C(=O)C=CC2=O', 'N1(CN2CCOCC2)CCC(OP(=O)(O)OCC(=O)CC(O1)C1=CC=CC=C1)(C(=O)O)C1=CC=CC=C1O']
Round 2 ['C[N+](C)(C)CC(CC(=O)[O-])OC(C1=CC=CC=C1)COP(=O)(O)OP(=O)(O)OC', 'C(CCN)NCC(CC(=O)O)OC(CC(CCCO)O)C(C)N', 'N1C(=O)C=CC=C(C(=O)OCC)C(=O)NC(=O)N1C1=CC=CC=C1', 'C1=NC2=C(C1=O)C(=O)C=CC2=O', 'N1(CN2CCOCC2)CCC(OP(=O)(O)OCC(=O)CC(O1)C1=CC=CC=C1)(C(=O)O)C1=CC=CC=C1O']
Round 3 ['c1ccc(C(COP(=O)(O)OP(=O)(O)OC)OC(CC(=O)[O-])C[N+](C)(C)C)NC(=O)O1', 'OC(=O)CC(CNCCCN)OC(CC(O)CCCO)C(=O)C', 'c1(-n2[nH]c(=O)cccc(C(=O)OCC)c(=O)[nH]c2=O)N(C1=O)C', 'C1=NC2=C(C1=O)C(=O)C=CC2=O', 'C1COCCN1CN1CCC(C(=O)O)(c2ccccc2O)OP(=O)(O)OCC(=O)CC(c2ccccc2)C=CC=C1']
Round 4 ['O=P(O)(OC)OP(=O)(O)OCC(OC(CC(=O)[O-])C[N+](C)(C)C)C1=CC(=C(C=C1)C)C', 'C(NCCCN)C(CC(=O)O)OC(CC(O)(C)O)(C)C', 'c1(-n2[nH]c(=O)cccc(C(=O)OCC)c(=O)[nH]c2=O)N(C1=O)C', 'C1=NC2=C(C1=O)C13CCCCC1C2CC2=C3C=C(C=C2Cl)C=CC2C3C(CCC1C(C2CC(C3)O)O)CCOC1(C)C', 'C1COCCN1CN1CCC(C(=O)O)(c2ccccc2O)OP(=O)(O)OCC(=O)CC(c2ccccc2)C=CC=C1']

並び替えを行う方が、多彩な化学構造を得ることができるので、逆解析の際に有利となる。次の例では文字列の並び替えを行わない例を示す。デフォルトのtrain_order=10なので、10文字が現れている。

# Method 1: use canonical SMILES in RDKit with no reordering
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi)) for smi in data['SMILES']]
n_gram_cans = NGram(reorder_prob=0)
n_gram_cans.fit(cans)

# save results
with open('ngram_cans.obj', 'wb') as f:
    pk.dump(n_gram_cans, f)


NGram
NGram(ngram_table=[[[ C ( O [N+] ! =C & N P =[N+] ... \
['C'] 24476 10218 3356 49 2126 1755 2752 2379 35 1 ...
[')'] 6594 1984 2590 68 0 284 0 1429 22 0 ...
['O'] 3240 0 34 21 2355 1782 0 64 302 277 ...
['[N+]'] 0 271 0 0 0 0 7 0 0 0 ...
['=C'] 1292 1760 45 6 46 8 550 35 1 0 ...
... ... ... ... ... ... ... ... ... ... ... ...
['[GeH4]'] 0 0 0 0 1 0 0 0 0 0 ...
['[Ti]'] 0 1 0 0 0 0 0 0 0 0 ...
['[SeH]'] 0 0 0 0 1 0 0 0 0 0 ...
['[Cu]'] 0 0 0 0 0 0 0 0 0 0 ...
['=[O+]'] 0 0 0 0 0 0 0 0 0 0 ....
['c', '(', '=O', ')', 'c', '&', 'c', 'c', 'c', ... 1 0 0 0
['(', '=O', ')', 'c', '&', 'c', 'c', 'c', 'c', ... 0 1 0 0
['c', 'c', 'c', '&', 'c', 'c', 'c', '(', 'c', '&'] 1 0 0 0
['c', 'c', '&', 'c', 'c', 'c', '(', 'c', '&', 'c'] 1 0 0 0
['c', '&', 'c', 'c', 'c', '(', 'c', '&', 'c', 'c'] 1 0 0 0
['&', 'c', 'c', 'c', '(', 'c', '&', 'c', 'c', 'c'] 0 0 1 0
['c', 'c', 'c', '(', 'c', '&', 'c', 'c', 'c', '('] 1 0 0 0
['c', 'c', '(', 'c', '&', 'c', 'c', 'c', '(', 'c'] 0 0 0 1]]])

以上は肩慣らしで、N-gram作成はここからが本番となる。まず、CanonicalなSMILESに変換しなおし、それをmolオブジェクトに変換する。SMILESへの逆変換の際に、SMILESの最初の原子を順繰りに変えた文字列を作成する。重複をなくしたリストを作成し、これをsmi_reorder_allに入れる。重複をなくしたリストの個数がn_reorder以下であれば、ランダムに並び替えて、sample_reorder
入れる。そして、生成したSMILESの文字列数に応じて、ヒストグラムを表示する。

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# Method 2: expand n-gram training set with randomly reordered SMILES
# (we show one of the many possible ways of doing it)
n_reorder = 15 # pick a fixed number of re-ordering

# convert the SMILES to canonical SMILES in RDKit (not necessary in general)
cans = []
for smi in data['SMILES']:
    # remove some molecules in the full SMILES list that may lead to error
    try:
        cans.append(Chem.MolToSmiles(Chem.MolFromSmiles(smi)))
    except:
        print(smi)
        pass

mols = [Chem.MolFromSmiles(smi) for smi in cans]
smi_reorder_all = []
smi_reorder = []
for mol in mols:
    idx = list(range(mol.GetNumAtoms()))
    tmp = [Chem.MolToSmiles(mol,rootedAtAtom=x) for x in range(len(idx))]
    smi_reorder_all.append(np.array(list(set(tmp))))
    smi_reorder.append(np.random.choice(smi_reorder_all[-1], n_reorder,
                                        replace=(len(smi_reorder_all[-1]) < n_reorder)))

n_uni = [len(x) for x in smi_reorder_all]
_ = plt.hist(n_uni, bins='auto')  # arguments are passed to np.histogram
_ = plt.title("Histogram for number of SMILES variants")
_ = plt.show()

次に、ここで改めて、新しいNGramクラスを定義する。ここでは、train_orderは10文字、並び替え確率は0.5としている。その後、20文字に拡張してモデルを作っている。このモデル作成にはかなり時間がかかる。M1proで数時間くらい。

# flatten out the list and train the N-gram
flat_list = [item for sublist in smi_reorder for item in sublist]

# first train up to order 10
n_gram_reorder1 = NGram(reorder_prob=0.5)
_ = n_gram_reorder1.fit(flat_list, train_order=10)
# save results
_ = joblib.dump(n_gram_reorder1, 'ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'wb') as f:
#     pk.dump(n_gram_reorder1, f)
    
# Then, train up from order 11 to 20
n_gram_reorder2 = NGram(reorder_prob=0.5, sample_order=(11, 20))
_ = n_gram_reorder2.fit(flat_list, train_order=(11, 20))
# save results
_ = joblib.dump(n_gram_reorder2, 'ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'wb') as f:
#     pk.dump(n_gram_reorder2, f)

データをロードするときには、以下のようにする。二つのN-gramをマージすることで一つにして、次のiQSPR解析に用いる。

# load a pre-trained n-gram from the pickle file
n_gram = joblib.load('ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'rb') as f:
#     n_gram = pk.load(f)
    
# load a pre-trained n-gram from the pickle file
n_gram2 = joblib.load('ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'rb') as f:
#     n_gram2 = pk.load(f)

_ = n_gram.merge_table(n_gram2)

まとめ

この作業自体はかなり長くかかるため、別解析として考えた方が良い。セーブするところまでを一つの区切りとする。

前述の書籍によるとそのほかの逆解析方法として遺伝的アルゴリズム、VAE、GANがある。Xenonpyが終わったら、VAEか、GANはここで、試してみたい。

さて、次のiQSPR部分は、表示が面倒なだけで、コーディングの面では、かなり楽だと思っている。

一番の難所は、
Descriptorの定義
学習モデルの準備
LogLikelihood関数の定義
である。

そのため、先にN-gram生成を行なっておくと、実際の研究には良いだろう。

run tools.ipynb

# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# Method 2: expand n-gram training set with randomly reordered SMILES
# (we show one of the many possible ways of doing it)
n_reorder = 15 # pick a fixed number of re-ordering

# convert the SMILES to canonical SMILES in RDKit (not necessary in general)
cans = []
for smi in data['SMILES']:
    # remove some molecules in the full SMILES list that may lead to error
    try:
        cans.append(Chem.MolToSmiles(Chem.MolFromSmiles(smi)))
    except:
        print(smi)
        pass

mols = [Chem.MolFromSmiles(smi) for smi in cans]
smi_reorder_all = []
smi_reorder = []
for mol in mols:
    idx = list(range(mol.GetNumAtoms()))
    tmp = [Chem.MolToSmiles(mol,rootedAtAtom=x) for x in range(len(idx))]
    smi_reorder_all.append(np.array(list(set(tmp))))
    smi_reorder.append(np.random.choice(smi_reorder_all[-1], n_reorder,
                                        replace=(len(smi_reorder_all[-1]) < n_reorder)))

n_uni = [len(x) for x in smi_reorder_all]
_ = plt.hist(n_uni, bins='auto')  # arguments are passed to np.histogram
_ = plt.title("Histogram for number of SMILES variants")
_ = plt.show()


# flatten out the list and train the N-gram
flat_list = [item for sublist in smi_reorder for item in sublist]

# first train up to order 10
n_gram_reorder1 = NGram(reorder_prob=0.5)
_ = n_gram_reorder1.fit(flat_list, train_order=10)
# save results
_ = joblib.dump(n_gram_reorder1, 'ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'wb') as f:
#     pk.dump(n_gram_reorder1, f)
    
# Then, train up from order 11 to 20
n_gram_reorder2 = NGram(reorder_prob=0.5, sample_order=(11, 20))
_ = n_gram_reorder2.fit(flat_list, train_order=(11, 20))
# save results
_ = joblib.dump(n_gram_reorder2, 'ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'wb') as f:
#     pk.dump(n_gram_reorder2, f)

二日目　学習モデルとLogLikelihood関数について

Xenonpy

前回のおさらい

dango0dan5.hatenablog.com


XenonPyを使って記述子を設定する方法を学んだ。記述子にECFP6を選んだ。

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'])

学習モデルの作成

二日目の最終ゴールは、異なる学習方法で作成したモデルから、LogLikelihood関数を作成することである。

まずは、物性名を辞書として使うために、リストを先に生成しておく必要がある。

# property name will be used as a reference for calling models
prop = ['E','HOMO-LUMO gap']

コード使い回し時に、めんどくさければ、dataから"SMILES"を取り除くようにすることでも設定できる。

print(data.columns)
prop=list(data.columns)
prop.remove("SMILES")
print(prop)

Index(['SMILES', 'E', 'HOMO-LUMO gap'], dtype='object')
['E', 'HOMO-LUMO gap']

XenonPyの交差検証用にデータを分割する方法を利用して以下ではモデリングをしていく。一つ一つの物性値に対してまずは実施して、やり方を確かめてみる。k_foldの値は分割数を表す。

from xenonpy.datatools import Splitter
# prepare indices for cross-validation data sets
sp = Splitter(data.shape[0], test_size=0, k_fold=10)

次に、学習データ(train)に対応する教師データ・予測データ、テストデータ(test)に対する教師データ・予測データを格納する配列を定義する。
また、mdlsは辞書ファイルで、mdl["物性名"]とするとモデルオブジェクトが返されるものと定義する。

# initialize output variables
y_trues_test, y_preds_test = [[] for i in range(len(prop))], [[] for i in range(len(prop))]
y_trues_train, y_preds_train = [[] for i in range(len(prop))], [[] for i in range(len(prop))]
mdls = {key: [] for key in prop}

Sampleコードを参照しながら、prop["E"]に対するモデリングを行う。とりあえず、ECFP6_FPsとElasticNetの組み合わせをここでは使ってみた。

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import ElasticNet

Target="E"
Desc=ECFP6_FPs
# cross-validation test
for iTr, iTe in sp.cv():

    x_train = data['SMILES'].iloc[iTr]
    x_test = data['SMILES'].iloc[iTe]
    
    fps_train = Desc.transform(x_train)
    fps_test = Desc.transform(x_test)
    
    y_train = data[prop].iloc[iTr]
    y_test = data[prop].iloc[iTe]

    propid=prop.index(Target)

    #mdl = RandomForestRegressor(500, n_jobs=-1, max_features='sqrt')
    #mdl = ElasticNetCV(cv=5)
    #mdl = GradientBoostingRegressor(loss='absolute_error')
    mdl = ElasticNet(alpha=0.5)
        
    mdl.fit(fps_train, y_train.iloc[:,propid])
    prd_train = mdl.predict(fps_train)
    prd_test = mdl.predict(fps_test)

    y_trues_test[propid].append(y_test.iloc[:,propid].values)
    y_trues_train[propid].append(y_train.iloc[:,propid].values)
    
    y_preds_test[propid].append(prd_test)
    y_preds_train[propid].append(prd_train)
    
    mdls[Target].append(mdl)

from sklearn.metrics import r2_score

y_true_test = np.concatenate(y_trues_test[propid])
y_pred_test = np.concatenate(y_preds_test[propid])
y_true_train = np.concatenate(y_trues_train[propid])
y_pred_train = np.concatenate(y_preds_train[propid])
xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
xy_diff = xy_max - xy_min
r2=r2_score(y_true_test,y_pred_test)
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
_ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
_ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
_ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
_ = plt.title('Property: ' + Target)
_ = plt.xlim(xy_min,xy_max)
_ = plt.ylim(xy_min,xy_max)
_ = plt.legend(loc='upper left')
_ = plt.xlabel('Prediction')
_ = plt.ylabel('Observation')
_ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
_ = plt.show()

sampleの予測結果とあまり変わらないように見えている。ここで、mdlsには、各検証でFittingした予測モデルが入っている。この予測モデルを調べてみる。一つのSMILESに対して、記述子を計算し、入れてみると、それぞれのモデルで値が若干異なっていることがわかる。これは、LogLikelihood関数を作成するときに大事なことになる。

tmpFP = ECFP6_FPs.transform([data["SMILES"][1]])
for x in mdls["E"]:
    print(x.predict(tmpFP))

[1.95511605]
[1.95277628]
[1.95631381]
[1.94501778]
[1.95411992]
[1.94175693]
[1.94580894]
[1.94818088]
[1.95047658]
[1.95249834]


一個一個個別に作って行ってもいいのだが、結局は同じ操作をすることになります。それぞれの物性値の予測モデルの、学習モデルと記述子を個別に与えたいために、以下の関数を導入します。deepcopyを使わなければならないのに注意。

from copy import deepcopy

# initialize model dictionary
mdls = {key: [] for key in prop}

def CrossValidation(sp,propid,Desc,model):
    y_trues_test=[]
    y_preds_test=[]
    y_trues_train=[]
    y_preds_train=[]
    out_mdls=[]
    # cross-validation test
    for iTr, iTe in sp.cv():

        x_train = data['SMILES'].iloc[iTr]
        x_test = data['SMILES'].iloc[iTe]
    
        fps_train = Desc.transform(x_train)
        fps_test = Desc.transform(x_test)
    
        y_train = data[prop].iloc[iTr]
        y_test = data[prop].iloc[iTe]

        mdl = deepcopy( model )
        
        mdl.fit(fps_train, y_train.iloc[:,propid])
        prd_train = mdl.predict(fps_train)
        prd_test = mdl.predict(fps_test)

        y_trues_test.append(y_test.iloc[:,propid].values)
        y_trues_train.append(y_train.iloc[:,propid].values)
    
        y_preds_test.append(prd_test)
        y_preds_train.append(prd_train)
        out_mdls.append(mdl)
    
    y_true_test = np.concatenate(y_trues_test)
    y_pred_test = np.concatenate(y_preds_test)
    y_true_train = np.concatenate(y_trues_train)
    y_pred_train = np.concatenate(y_preds_train)
    xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
    xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
    xy_diff = xy_max - xy_min

    r2=r2_score(y_true_test,y_pred_test)

    _ = plt.figure(figsize=(5,5))
    _ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
    _ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
    _ = plt.title('Property: ' + Target)
    _ = plt.xlim(xy_min,xy_max)
    _ = plt.ylim(xy_min,xy_max)
    _ = plt.legend(loc='upper left')
    _ = plt.xlabel('Prediction')
    _ = plt.ylabel('Observation')
    _ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
    _ = plt.show()
    return out_mdls

n_bitsを小さくして、予測モデル作成を早く終わるようにしてみます。16の値を入れたことはないのだが、やってみました。

Target="E"
ECFP6_FPs = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
Desc=ECFP6_FPs
model=GradientBoostingRegressor(loss='squared_error')
mdls[Target]=CrossValidation(sp,prop.index(Target),Desc,model)

結果として、良いR2の値と、MAEの値となっていましたね。次に、予測モデルを使って、一個のSMILESから予想される値を調べると、バラついていることがわかります。

tmpFP = ECFP6_FPs.transform([data["SMILES"][1]])
for x in mdls["E"]:
    print(x.predict(tmpFP))

[1.55619993]
[1.4899335]
[1.57128927]
[1.58779223]
[1.56506169]
[1.51193849]
[1.61289896]
[1.54215519]
[1.62158862]
[1.52752974]


次に、バンドギャップについても作成してみますが、違う記述子と違うモデルを作成してみました。

from sklearn.ensemble import RandomForestRegressor

Target="HOMO-LUMO gap"
Desc= Fingerprints(featurizers=['PatternFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=RandomForestRegressor(10, n_jobs=-1, max_features='sqrt')
mdls[Target]=CrossValidation(sp,prop.index(Target),Desc,model)

これもまぁまぁなモデルが作れています。Sampleと違って、10eV以上の点で改善されているように見えます。RadomForestがベストなものではないので、いろいろ変えてみると、さらに良くなるかもしれません。

tmpFP = Desc.transform([data["SMILES"][1]])
for x in mdls["HOMO-LUMO gap"]:
    print(x.predict(tmpFP))

LogLikelihood関数

LogLikelihood関数を定義します。ここでは、少し特殊なことをしています。custom_mdlsを定義するところまでは、sampleと同じです。今回は別々の記述子を使うので、それぞれの予測モデルを組み込んだcustom_mdls1とcustom_mdls2を定義しています。それぞれの記述子を指定するために、GaussianLogLikelihood関数をつかって、prd_mdls1とprd_mdls2を定義しています。

# Forward model template in XenonPy-iQSPR 
from xenonpy.inverse.iqspr import GaussianLogLikelihood

class bootstrap_fcn():
    def __init__(self, ms, var):
        self.Ms = ms
        self.Var = var
    def predict(self, x):
        val = np.array([m.predict(x) for m in self.Ms])
        return np.mean(val, axis = 0), np.sqrt(np.var(val, axis = 0) + self.Var)
    
# include basic measurement noise
c_var = {'E': 0.1, 'HOMO-LUMO gap': 0.1}

# generate a dictionary for the final models used to create the likelihood
custom_mdls = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() }
custom_mdls1 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="E"}
custom_mdls2 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="HOMO-LUMO gap"}

# import descriptor calculator and forward model to iQSPR
prd_mdls1 = GaussianLogLikelihood(descriptor=ECFP6_FPs, targets={'E': (0, 2)}, **custom_mdls1)
prd_mdls2 = GaussianLogLikelihood(descriptor=Desc, targets={'HOMO-LUMO gap': (-np.inf, 3)}, **custom_mdls2)

このprd_mdls1が予測モデルとして使えるかどうか、調べてみましょう。

prd_mdls1.predict(data["SMILES"])

次からは難問です。オリジナルのLogLikelihoodを定義するのですが、predictorの設定や、それに伴ってsampleにあったコードを書き換える必要があります。基本的には、GaussianLogLikelihood関数のソースコードを参考にしながら、改変することでうまくいきました。
github.com

from xenonpy.inverse.base import BaseLogLikelihoodSet
from sklearn.linear_model import BayesianRidge

target_range = {'E': (0,2), 'HOMO-LUMO gap': (-np.inf, 3)}
desc={"E":ECFP6_FPs,'HOMO-LUMO gap':Desc}

class MyLogLikelihood(BaseLogLikelihoodSet):
    def __init__(self):
        super().__init__()
        self.loglike = prd_mdls1
        self.loglike = prd_mdls2
        self._mdl=custom_mdls
        self.descriptor=desc
        self.targets=target_range
    def predict(self, smiles, **kwargs):
        tmp = {}
        for k, v in self._mdl.items():
            fps = self.descriptor[k].transform(smiles, return_type='df')
            fps_ = fps.dropna()
            tmp[k + ': mean'], tmp[k + ': std'] = v.predict(fps_)
        tmp = pd.DataFrame(data=tmp, index=fps_.index)
        return pd.DataFrame(data=tmp, index=fps.index)

like_mdl = MyLogLikelihood()

これでpredictorが動くかどうか調べてみましたが、うまく行っていますね。

pred=like_mdl.predict(data['SMILES'])
pred

では、最後にいろいろとプロットさせていきます。

# calculate log-likelihood for a given target property region
tmp_ll = like_mdl(data['SMILES'])
tmp = tmp_ll.sum(axis = 1, skipna = True)

_ = plt.figure(figsize=(8,5))
_ = plt.hist(tmp, bins=50)
_ = plt.title('ln(likelihood) for SMILES in the data subset')
_ = plt.xlabel('ln(likelihood)')
_ = plt.ylabel('Counts')
_ = plt.show()

# plot histogram of likelihood values
_ = plt.figure(figsize=(8,5))
_ = plt.hist(np.exp(tmp), bins=50)
_ = plt.title('Likelihood for SMILES in the data subset')
_ = plt.xlabel('Likelihood')
_ = plt.ylabel('Counts')
_ = plt.show()



# check the predicted likelihood
dot_scale = 0.5
l_std = np.sqrt(pred['E: std']**2+pred['HOMO-LUMO gap: std']**2)

_ = plt.figure(figsize=(6,5))
rectangle = plt.Rectangle((0,0),2,3,fc='y',alpha=0.1)
_ = plt.gca().add_patch(rectangle)
im = plt.scatter(pred['E: mean'], pred['HOMO-LUMO gap: mean'], s=l_std*dot_scale, c=np.exp(tmp),alpha = 0.2,cmap=plt.get_cmap('cool'))
_ = plt.title('Pubchem data')
_ = plt.xlabel('E')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.colorbar(im)
_ = plt.show()

ここで作成したモデルは、高いバンドギャップの領域(10eV以上）での予想の傾向が向上していたので、そちらの探索にも使えるかもしれませんね。次回以降は、領域指定の値を変えて実施してもいいかなと思います。

二日目のまとめ

これらの結果をまとめると、下のようなコードになります。jupyter-notebookでやるときには適当なところで、セルを分割してください。

run tools.ipynb

from xenonpy.descriptor import Fingerprints

# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
# make histogram
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()

# set target range
target_range = {'E': (0,2), 'HOMO-LUMO gap': (-np.inf, 3)}

# prepare indices for cross-validation data sets
from xenonpy.datatools import Splitter
sp = Splitter(data.shape[0], test_size=0, k_fold=10)

# set property list
prop=list(data.columns)
prop.remove("SMILES")

from copy import deepcopy
from sklearn.metrics import r2_score

# initialize model dictionary
mdls = {key: [] for key in prop}

def CrossValidation(sp,propid,Desc,model):
    y_trues_test=[]
    y_preds_test=[]
    y_trues_train=[]
    y_preds_train=[]
    out_mdls=[]
    # cross-validation test
    for iTr, iTe in sp.cv():

        x_train = data['SMILES'].iloc[iTr]
        x_test = data['SMILES'].iloc[iTe]
    
        fps_train = Desc.transform(x_train)
        fps_test = Desc.transform(x_test)
    
        y_train = data[prop].iloc[iTr]
        y_test = data[prop].iloc[iTe]

        mdl = deepcopy( model )
        
        mdl.fit(fps_train, y_train.iloc[:,propid])
        prd_train = mdl.predict(fps_train)
        prd_test = mdl.predict(fps_test)

        y_trues_test.append(y_test.iloc[:,propid].values)
        y_trues_train.append(y_train.iloc[:,propid].values)
    
        y_preds_test.append(prd_test)
        y_preds_train.append(prd_train)
        out_mdls.append(mdl)
    
    y_true_test = np.concatenate(y_trues_test)
    y_pred_test = np.concatenate(y_preds_test)
    y_true_train = np.concatenate(y_trues_train)
    y_pred_train = np.concatenate(y_preds_train)
    xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
    xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
    xy_diff = xy_max - xy_min

    r2=r2_score(y_true_test,y_pred_test)

    _ = plt.figure(figsize=(5,5))
    _ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
    _ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
    _ = plt.title('Property: ' + Target)
    _ = plt.xlim(xy_min,xy_max)
    _ = plt.ylim(xy_min,xy_max)
    _ = plt.legend(loc='upper left')
    _ = plt.xlabel('Prediction')
    _ = plt.ylabel('Observation')
    _ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
    _ = plt.show()
    return out_mdls

from sklearn.ensemble import GradientBoostingRegressor
Target="E"
FP_E = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=GradientBoostingRegressor(loss='squared_error')
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_E,model)

from sklearn.ensemble import RandomForestRegressor
Target="HOMO-LUMO gap"
FP_H= Fingerprints(featurizers=['PatternFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=RandomForestRegressor(10, n_jobs=-1, max_features='sqrt')
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_H,model)


# Forward model template in XenonPy-iQSPR 
from xenonpy.inverse.iqspr import GaussianLogLikelihood

class bootstrap_fcn():
    def __init__(self, ms, var):
        self.Ms = ms
        self.Var = var
    def predict(self, x):
        val = np.array([m.predict(x) for m in self.Ms])
        return np.mean(val, axis = 0), np.sqrt(np.var(val, axis = 0) + self.Var)
    
# include basic measurement noise
c_var = {'E': 0.1, 'HOMO-LUMO gap': 0.1}

# generate a dictionary for the final models used to create the likelihood
custom_mdls1 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="E"}
custom_mdls2 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="HOMO-LUMO gap"}

# import descriptor calculator and forward model to iQSPR
prd_mdls1 = GaussianLogLikelihood(descriptor=FP_E, targets={'E': target_range["E"]}, **custom_mdls1)
prd_mdls2 = GaussianLogLikelihood(descriptor=FP_H, targets={'HOMO-LUMO gap': target_range['HOMO-LUMO gap']}, **custom_mdls2)


from xenonpy.inverse.base import BaseLogLikelihoodSet
from sklearn.linear_model import BayesianRidge

custom_mdls = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() }
desc={"E":FP_E,'HOMO-LUMO gap':FP_H}

class MyLogLikelihood(BaseLogLikelihoodSet):
    def __init__(self):
        super().__init__()
        self.loglike = prd_mdls1
        self.loglike = prd_mdls2
        self._mdl=custom_mdls
        self.descriptor=desc
        self.targets=target_range
    def predict(self, smiles, **kwargs):
        tmp = {}
        for k, v in self._mdl.items():
            fps = self.descriptor[k].transform(smiles, return_type='df')
            fps_ = fps.dropna()
            tmp[k + ': mean'], tmp[k + ': std'] = v.predict(fps_)
        tmp = pd.DataFrame(data=tmp, index=fps_.index)
        return pd.DataFrame(data=tmp, index=fps.index)

like_mdl = MyLogLikelihood()

pred=like_mdl.predict(data['SMILES'])

# calculate log-likelihood for a given target property region
tmp_ll = like_mdl(data['SMILES'])
tmp = tmp_ll.sum(axis = 1, skipna = True)

_ = plt.figure(figsize=(8,5))
_ = plt.hist(tmp, bins=50)
_ = plt.title('ln(likelihood) for SMILES in the data subset')
_ = plt.xlabel('ln(likelihood)')
_ = plt.ylabel('Counts')
_ = plt.show()

# plot histogram of likelihood values
_ = plt.figure(figsize=(8,5))
_ = plt.hist(np.exp(tmp), bins=50)
_ = plt.title('Likelihood for SMILES in the data subset')
_ = plt.xlabel('Likelihood')
_ = plt.ylabel('Counts')
_ = plt.show()

# check the predicted likelihood
dot_scale = 0.5
l_std = np.sqrt(pred['E: std']**2+pred['HOMO-LUMO gap: std']**2)

_ = plt.figure(figsize=(6,5))
rectangle = plt.Rectangle((0,0),2,3,fc='y',alpha=0.1)
_ = plt.gca().add_patch(rectangle)
im = plt.scatter(pred['E: mean'], pred['HOMO-LUMO gap: mean'], s=l_std*dot_scale, c=np.exp(tmp),alpha = 0.2,cmap=plt.get_cmap('cool'))
_ = plt.title('Pubchem data')
_ = plt.xlabel('E')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.colorbar(im)
_ = plt.show()

一日目　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'])