In [39]:
plot_the_pred_distribution_for_four_systems(reg, features_df, bar_ene)
0.3580137787084095
import numpy as np
import pandas as pd
# load the preprocessing package(mainly sklearn)
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.linear_model import LinearRegression
# load ploting packages
import matplotlib.pyplot as plt
import seaborn as sns
# load features
# d1 Phe72 H𝜁 (Tyr72 H𝜂 ) – Glu166 O𝜀2
# d2 Lys73 H𝜁2 – Glu166 O𝜀2
# d3 Water H1 – Glu166 O𝜀2
# d4 Water H1 – Glu166 O𝜀1
# d5 Asn170 H𝛿2 – Glu166 O𝜀1
# d6 Water H2 – Asn170 O𝛿
# d7 Lys73 H𝜁1 – Water O
# d8 Lys73 H𝜁2 – Asn132 O𝛿
# d9 Water O – IPM C7
# d10 Lys73 H𝜁1 – Ser70 O𝛾
# d11 IPM 6𝛼OH – Water O
# d12 IPM 6𝛼OH – Asn132 O𝛿
# d13 IPM 6𝛼OH – Glu166 O𝜀1
# d14 IPM 6𝛼OH – Glu166 O𝜀2
# d15 Lys73 H𝜁1 – Ser130 O𝛾
# d16 Ser130 H𝛾 – IPM N4
features_name = ['d1', 'd2', 'd3', 'd4', 'd5', 'd6', 'd7', 'd8', 'd9', 'd10', 'd11','d12', 'd13', 'd14', 'd15', 'd16']
wt_d1_16_distances = np.load("/users/chaoy/scratch/0.proj_kpc/2.data_anal/2.disanal/3.new_features_16_distances/kpc_wt_d1_16_features.npy")
wt_d2_16_distances = np.load("/users/chaoy/scratch/0.proj_kpc/2.data_anal/2.disanal/3.new_features_16_distances/kpc_wt_d2_16_features.npy")
y72_d1_16_distances = np.load("/users/chaoy/scratch/0.proj_kpc/2.data_anal/2.disanal/3.new_features_16_distances/kpc_y72_d1_16_features.npy")
y72_d2_16_distances = np.load("/users/chaoy/scratch/0.proj_kpc/2.data_anal/2.disanal/3.new_features_16_distances/kpc_y72_d2_16_features.npy")
# concatenate four systems together
features = np.concatenate((wt_d1_16_distances, wt_d2_16_distances, y72_d1_16_distances, y72_d2_16_distances), axis = 1).T
# build a dataframe for features
features_df = pd.DataFrame(features, columns = ['d1', 'd2', 'd3', 'd4', 'd5', 'd6', 'd7', 'd8', 'd9', 'd10', 'd11','d12', 'd13', 'd14', 'd15', 'd16'])
# load energy
system_lables = ["wt_d1"] * 200 + ["wt_d2"] * 200 + ["y72_d1"] * 200 + ["y72_d2"] * 200
# load the barrier energy for every pathway
wt_d1_bar_ene = np.load("/users/chaoy/scratch/0.proj_kpc/1.sample/1.kpc_wt_d1/7.sp/0.b3lyp_d3/charge_ene/enes_all_paths.npz")['ene_barrier']
wt_d2_bar_ene = np.load("/users/chaoy/scratch/0.proj_kpc/1.sample/0.kpc_wt_d2/7.sp/0.b3lyp_d3/charge_ene/enes_all_paths.npz")['ene_barrier']
y72_d1_bar_ene = np.load("/users/chaoy/scratch/0.proj_kpc/1.sample/3.kpc_y72_d1/7.sp/0.b3lyp_d3/charge_ene/enes_all_paths.npz")['ene_barrier']
y72_d2_bar_ene = np.load("/users/chaoy/scratch/0.proj_kpc/1.sample/2.kpc_y72_d2/7.sp/0.b3lyp_d3/charge_ene/enes_all_paths.npz")['ene_barrier']
# merge barrier energy together
bar_ene = np.concatenate((wt_d1_bar_ene, wt_d2_bar_ene, y72_d1_bar_ene, y72_d2_bar_ene))
np.mean(wt_d1_16_distances, axis = 1)
array([3.3374522, 2.1554048, 1.6497445, 2.5110903, 1.9213964, 1.7693936,
2.8241355, 1.8590629, 2.8657162, 2.116153 , 3.468609 , 3.612612 ,
4.1744184, 4.5902443, 3.0784116, 1.8957262], dtype=float32)
np.mean(wt_d2_16_distances, axis = 1)
array([3.1083362, 2.022509 , 1.652809 , 2.4155698, 1.8823148, 1.7635456,
3.1083996, 1.9095639, 2.7960742, 2.1741593, 3.2451575, 3.6568122,
4.087621 , 4.6142855, 3.1971235, 2.5870457], dtype=float32)
np.mean(y72_d1_16_distances, axis = 1)
array([1.6614423, 2.3930242, 1.6546925, 2.4763231, 1.9471045, 1.7859906,
2.916861 , 1.8526888, 2.948042 , 2.0655007, 2.8675365, 3.3877716,
3.6053724, 4.165021 , 2.8541772, 1.9278173], dtype=float32)
np.mean(y72_d2_16_distances, axis = 1)
array([1.6480212, 2.1525745, 1.7016983, 2.3271143, 1.8654382, 1.7868141,
3.157306 , 1.89794 , 2.8289635, 2.0891602, 2.7356691, 3.3523095,
3.1997576, 4.021411 , 3.0606377, 2.5454714], dtype=float32)
X_train, X_test, y_train, y_test = train_test_split(features_df, bar_ene, test_size= 0.1, random_state=42, stratify = system_lables)
reg = LinearRegression().fit(X_train, y_train)
def plot_the_pred_distribution_for_four_systems(model, bonds_new_df, bar_ene):
# features for four systems
data_wt_d1 = bonds_new_df[0:200]
data_wt_d2 = bonds_new_df[200:400]
data_y72_d1 = bonds_new_df[400:600]
data_y72_d2 = bonds_new_df[600:800]
# energy for four systems
y_wt_d1 = bar_ene[0:200]
y_wt_d2 = bar_ene[200:400]
y_y72_d1 = bar_ene[400:600]
y_y72_d2 = bar_ene[600:800]
y_wt_d1_pred = model.predict(data_wt_d1).flatten()
y_wt_d2_pred = model.predict(data_wt_d2).flatten()
y_y72_d1_pred = model.predict(data_y72_d1).flatten()
y_y72_d2_pred = model.predict(data_y72_d2).flatten()
wt_d1_error = y_wt_d1_pred - y_wt_d1
wt_d2_error = y_wt_d2_pred - y_wt_d2
y72_d1_error = y_y72_d1_pred - y_y72_d1
y72_d2_error = y_y72_d2_pred - y_y72_d2
# print("MAX error in the training set:", max(train_error))
# print("MIN error in the training set:", min(train_error))
# print("MAE for the training set:", np.mean(np.abs(train_error)))
rstr_r2 = r'$R^2$'
rstr_mae = r'$MAE$'
rstr_rmse = r'$RMSE$'
rstr_nsamp = r'$n_{samples}$'
print(np.mean(wt_d1_error))
# WT_D1
r2_wt_d1 = round(r2_score(y_wt_d1, y_wt_d1_pred),4)
mae_wt_d1 = np.mean(np.abs(y_wt_d1 - y_wt_d1_pred))
rmse_wt_d1 = (np.sum((y_wt_d1 - y_wt_d1_pred)**2)/y_wt_d1.shape[0]) ** 0.5
# n_train = y_train.shape[0]
# WT_D2
r2_wt_d2 = round(r2_score(y_wt_d2, y_wt_d2_pred),4)
mae_wt_d2 = np.mean(np.abs(y_wt_d2 - y_wt_d2_pred))
rmse_wt_d2 = (np.sum((y_wt_d2 - y_wt_d2_pred)**2)/y_wt_d2.shape[0]) ** 0.5
# Y72_D1
r2_y72_d1 = round(r2_score(y_y72_d1, y_y72_d1_pred),4)
mae_y72_d1 = np.mean(np.abs(y_y72_d1 - y_y72_d1_pred))
rmse_y72_d1 = (np.sum((y_y72_d1 - y_y72_d1_pred)**2)/y_y72_d1.shape[0]) ** 0.5
# Y72_D2
r2_y72_d2 = round(r2_score(y_y72_d2, y_y72_d2_pred),4)
mae_y72_d2 = np.mean(np.abs(y_y72_d2 - y_y72_d2_pred))
rmse_y72_d2 = (np.sum((y_y72_d2 - y_y72_d2_pred)**2)/y_y72_d2.shape[0]) ** 0.5
# plot a figure for the predition and true value
sns.set_style(style='white')
sns.set_style("ticks")
fig, axs = plt.subplots(2,2,figsize = (7,7), dpi =300, facecolor = 'w', edgecolor = 'k')
# KPC-WT/IPM-D1
plt.subplots_adjust(left=.3/7, right=6.7/7, top=6.9/7, bottom=.1/7, wspace = 0.3, hspace = 0.3)
axs[0,0].tick_params(which='major', direction='in', length=6, width=1)
axs[0,0].scatter(y_wt_d1, y_wt_d1_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
axs[0,0].axis([-5, 85, -5, 85])
xpoints = ypoints = axs[0,0].set_xlim()
axs[0,0].plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
axs[0,0].set_title(r'KPC-WT/IPM-$\Delta$1')
axs[0,0].set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
axs[0,0].set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
axs[0,0].xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[0,0].yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[0,0].text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_wt_d1: <8.4f}\n{rstr_mae: >10} = {mae_wt_d1: <8.4f}\n{rstr_rmse: >10} = {rmse_wt_d1: <8.4f}\n {rstr_nsamp: >10} = {200: <8}',
ha='right', va='bottom', transform=axs[0,0].transAxes, fontsize=10)
#. KPC-WT/IPM-D2
axs[0,1].tick_params(which='major', direction='in', length=6, width=1)
axs[0,1].scatter(y_wt_d2, y_wt_d2_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
axs[0,1].axis([-5, 85, -5, 85])
xpoints = ypoints = axs[0,1].set_xlim()
axs[0,1].plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
axs[0,1].set_title(r'KPC-WT/IPM-$\Delta$2')
axs[0,1].set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
axs[0,1].set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
axs[0,1].xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[0,1].yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[0,1].text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_wt_d2: <8.4f}\n{rstr_mae: >10} = {mae_wt_d2: <8.4f}\n{rstr_rmse: >10} = {rmse_wt_d2: <8.4f}\n {rstr_nsamp: >10} = {200: <8}',
ha='right', va='bottom', transform=axs[0,1].transAxes, fontsize=10)
# KPC-F72Y/IPM-D1
axs[1,0].tick_params(which='major', direction='in', length=6, width=1)
axs[1,0].scatter(y_y72_d1, y_y72_d1_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
axs[1,0].axis([-5, 85, -5, 85])
xpoints = ypoints = axs[1,0].set_xlim()
axs[1,0].plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
axs[1,0].set_title(r'KPC-F72Y/IPM-$\Delta$1')
axs[1,0].set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
axs[1,0].set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
axs[1,0].xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[1,0].yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[1,0].text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_y72_d1: <8.4f}\n{rstr_mae: >10} = {mae_y72_d1: <8.4f}\n{rstr_rmse: >10} = {rmse_y72_d1: <8.4f}\n {rstr_nsamp: >10} = {200: <8}',
ha='right', va='bottom', transform=axs[1,0].transAxes, fontsize=10)
# KPC-F72Y/IPM-D2
axs[1,1].tick_params(which='major', direction='in', length=6, width=1)
axs[1,1].scatter(y_y72_d2, y_y72_d2_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
axs[1,1].axis([-5, 85, -5, 85])
xpoints = ypoints = axs[1,1].set_xlim()
axs[1,1].plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
axs[1,1].set_title(r'KPC-F72Y/IPM-$\Delta$2')
axs[1,1].set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
axs[1,1].set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
axs[1,1].xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[1,1].yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
axs[1,1].text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_y72_d2: <8.4f}\n{rstr_mae: >10} = {mae_y72_d2: <8.4f}\n{rstr_rmse: >10} = {rmse_y72_d2: <8.4f}\n {rstr_nsamp: >10} = {200: <8}' ,
ha='right', va='bottom', transform=axs[1,1].transAxes, fontsize=10)
# fig.suptitle(f"Linear Regression Model", fontsize = 12)
plt.tight_layout()
fig.savefig(f'./SI/0.Linear_Regression_four_systems_results.png')
plot_the_pred_distribution_for_four_systems(reg, features_df, bar_ene)
0.3580137787084095
def plot_the_pred_distribution(model, X_train, X_test, y_train, y_test):
# Training set
y_train_pred = model.predict(X_train).flatten()
y_test_pred = model.predict(X_test).flatten()
train_error = y_train_pred - y_train
# print("MAX error in the training set:", max(train_error))
# print("MIN error in the training set:", min(train_error))
# print("MAE for the training set:", np.mean(np.abs(train_error)))
rstr_r2 = r'$R^2$'
rstr_mae = r'$MAE$'
rstr_rmse = r'$RMSE$'
rstr_nsamp = r'$n_{samples}$'
r2_train = round(r2_score(y_train, y_train_pred),4)
mae_train = np.mean(np.abs(y_train - y_train_pred))
rmse_train = (np.sum((y_train - y_train_pred)**2)/y_train.shape[0]) ** 0.5
n_train = y_train.shape[0]
r2_test = round(r2_score(y_test, y_test_pred),4)
mae_test = np.mean(np.abs(y_test - y_test_pred))
rmse_test = (np.sum((y_test - y_test_pred)**2)/y_test.shape[0]) ** 0.5
n_test = y_test.shape[0]
print(r2_train, mae_train, rmse_train, n_train)
# Testing set
y_test_pred = model.predict(X_test).flatten()
test_error = y_test_pred - y_test
# print("MAX error in the testing set:", max(test_error))
# print("MIN error in the testing set:", min(test_error))
# print("MAE for the testing set:", np.mean(np.abs(test_error)))
r2_test = round(r2_score(y_test, y_test_pred),4)
# plot a figure for the predition and true value
sns.set_style(style='white')
sns.set_style("ticks")
fig, (ax1,ax2) = plt.subplots(1,2,figsize = (7,3.5), dpi =300, facecolor = 'w', edgecolor = 'k')
# plt.setp(ax1.spines.values(), linewidth=2)
# plt.subplot(121)
plt.subplots_adjust(left=.3/7, right=6.7/7, top=3.4/3.5, bottom=.1/3.5, wspace = 2/7)
ax1.tick_params(which='major', direction='in', length=6, width=1)
ax1.scatter(y_train, y_train_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
ax1.axis([-5, 85, -5, 85])
xpoints = ypoints = ax1.set_xlim()
ax1.plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
ax1.set_title('Training Set')
ax1.set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
ax1.set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
ax1.xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
ax1.yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
# ax1.text(0.06, 0.94, syslbl, ha='left', va='top', transform=ax.transAxes, fontsize=8)
ax1.text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_train: <8.4f}\n{rstr_mae: >10} = {mae_train: <8.4f}\n{rstr_rmse: >10} = {rmse_train: <8.4f}\n{rstr_nsamp: >10} = {n_train: <8} ',
ha='right', va='bottom', transform=ax1.transAxes, fontsize=10)
ax2.tick_params(which='major', direction='in', length=6, width=1)
ax2.scatter(y_test, y_test_pred , color=sns.color_palette('colorblind')[0], zorder=2, s = 3)
ax2.axis([-5, 85, -5, 85])
xpoints = ypoints = ax2.set_xlim()
ax2.plot(xpoints, ypoints, linestyle='-', color='r', lw=1, scalex=False, scaley=False)
ax2.set_title('Validation Set')
ax2.set_xlabel(r"$\Delta E_{QM/MM} (kcal$ $ mol^{-1})$")
ax2.set_ylabel(r"$\Delta E_{ML} (kcal$ $ mol^{-1})$")
ax2.xaxis.set_major_locator(plt.MultipleLocator(base=20.0))
ax2.yaxis.set_major_locator(plt.MultipleLocator(base=20.0))
# ax1.text(0.06, 0.94, syslbl, ha='left', va='top', transform=ax.transAxes, fontsize=8)
ax2.text(0.94, 0.06,
f'{rstr_r2: >10} = {r2_test: <8.4f}\n{rstr_mae: >10} = {mae_test: <8.4f}\n{rstr_rmse: >10} = {rmse_test: <8.4f}\n{rstr_nsamp: >10} = {n_test: <8} ',
ha='right', va='bottom', transform=ax2.transAxes, fontsize=12)
# fig.suptitle(f"Linear Regression Model", fontsize = 12)
plt.tight_layout()
# ax1.set_aspect(1.)
fig.savefig(f'./SI/0.Linear_Regression_training_and_validation_set_results.png')
plot_the_pred_distribution(reg, X_train, X_test, y_train, y_test)
0.8375 3.171696872300365 4.167423910835283 720
