空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

上机作业1: 不对称催化中的非共价作用（学生版）副本副本

机器学习及其在化学中的应用

Void

发布于 2024-03-29

推荐镜像 :Basic Image:bohrium-notebook:2023-04-07

推荐机型 :c2_m4_cpu

数据集

nci_birman(v1)

背景概述

对映选择性(enantioselectivity)是有机化学重要的研究主题之一. 在有机反应中, 手性催化剂的对映选择性往往与反应过渡态中一些重要的非共价相互作用(non-covalent interactions, NCIs)有关. 以Birman等报道的酰基转移反应为例:

alt image.png

在Birman等假设的机理中, 手性催化剂的对映选择性由 $π$ - $π$ 堆积或 $π$ -阳离子堆积作用介导, 特定手性的过渡态因位阻作用导致两个芳环不平行, 不利于堆积作用的形成.

alt image.png

本次上机作业中, 我们将以线性回归模型定量地探索所用手性催化剂的性质对对映选择性的影响, 试图揭示: (1) 哪些性质参数会影响对映选择性? (2) 如何优化性质参数可以提升对映选择性 (例如减少还是增加某参数)? 为此, 我们将考察如下图所示的模型体系:

alt image.png

其中, 底物的芳香环(A, a, b)按照特定的几何约束条件与催化剂(D, c, d)排列成复合物, 而芳环与底物的其余分子骨架(包含醇羟基和烷基)在图中箭头所示位置相连.

参考文献: J. Am. Chem. Soc. 2017, 139, 20, 6803–6806

代码

文本

数据探索与清洗

代码

文本

数据的初步探索

在数据集nci_birman.csv中, 包含如下三组(用作特征的)参数:

$π$ - $π$ 堆积作用参数: d_pi_d, d_pi_D, e_pi_d, e_pi_D. 根据特定的电子结构计算方法(B97-D/def2TZVP水平), 对复合物进行几何优化, 给出几何结构信息与互作能(复合物能量减去各组分能量和). 其中:
- 每组复合物中, 底物芳环由于具体朝向的不同, 存在2种可能构象(比如下图的例子), 分别以后缀_d和_D作区分;
- 在计算中, 对底物与催化剂分子的两个芳环施加约束条件, 使二者完全“正对彼此”(环平面平行, 且两环中心的连线与环平面垂直). 此时, 我们找到使电子能量(前缀e_, 单位kcal/mol)最低的、也就是使复合物最稳定的环间距(前缀d_, 单位Å)
烷基几何参数: L_Alk, B1_Alk, B5_Alk, 分别代表长度、最大宽度、最小宽度.
芳环几何参数: L_Ar, B1_Ar, B5_Ar, 分别代表长度、最大宽度、最小宽度. 后两组参数称为Sterimol参数, 其具体含义可从下图中做简单的理解:

alt image.png

我们的线性回归模型将以上述参数作为特征(在化学信息学中, 这些特征也常称为描述符, descriptor), 试图预测特定催化剂针对特定底物的对映选择性. 在数据集中, 选择性以对映比 $er$ 描述, 这是一个实验可测的量. 通过热力学公式, 我们可以将其与两种对映体的活化自由能之差 $ΔΔ G^{‡}$ 之间建立定量联系:

$ΔΔ G^{‡} = - RT ln (er) .$

代码

文本

[2]

import pandas as pd

nci_birman = pd.read_csv("/bohr/MLInChem-licv/v1/nci_birman.csv", sep="\t")

nci_birman.iloc[:5, :]

	d_pi_d	d_pi_D	e_pi_d	e_pi_D	L_Alk	B1_Alk	B5_Alk	L_Ar	B1_Ar	B5_Ar	er(%)
0	3.54	3.54	-10.24	-10.24	4.36	2.92	3.35	6.38	1.77	3.15	1.71
1	3.54	3.54	-10.24	-10.24	4.35	2.09	3.34	6.38	1.77	3.15	3.36
2	3.54	3.54	-10.24	-10.24	4.38	1.73	3.33	6.38	1.77	3.15	3.78
3	3.54	3.54	-10.24	-10.24	3.08	1.70	2.20	6.38	1.77	3.15	5.03
4	3.49	3.49	-11.88	-11.88	3.08	1.70	2.20	6.38	1.88	4.52	5.03

代码

文本

数据的清洗与加工

对于堆积相互作用参数, 我们不妨对其进行降维: 取两个可能构象_d和_D的参数加权和: $D π_{w} E π_{w} := c_{d} D π_{d} + c_{D} D π_{D}, := c_{d} E π_{d} + c_{D} E π_{D},$ 其中, 权重 $c_{X} (X = d, D)$ 取为(按电子能量计算的)Boltzmann因子: $c_{X} = \frac{exp ( - \frac{E π _{X}}{RT} )}{exp ( - \frac{E π _{d}}{RT} ) + exp ( - \frac{E π _{D}}{RT} ) .}$

代码

文本

任务1: 对数据nci_birman进行进一步加工处理.

在函数prepare_data()中, 输入:
- 直接读取自csv的原始数据集in_df;
- 用于进行加权和与 $ΔΔ G^{‡}$ 计算的温度temperature.
返回: 处理好的特征和标签.
特征为pd.DataFrame格式, 包括:
- 降维后的(加权)参数 $D π_{w}, E π_{w}$ , 以及二者的交叉项(乘积) $D π_{w} E π_{w}$ , 分别命名为d_pi_w, e_pi_w与de_pi_w;
- 所有原有的Sterimol几何参数, 名称不变.
标签为pd.Series格式, 为自由能差 $ΔΔ G^{‡} = - RT ln (er)$ , 单位取为kcal/mol, 命名为delta_delta_G. 这里给出单位换算因子: 1 kcal/mol = 4.184 kJ/mol.

代码

文本

[3]

from typing import Tuple

import numpy as np

R = 8.314 # ideal gas constant

kcal_to_kj = 4.184 # unit conversion

def prepare_data(in_df: pd.DataFrame, T: float=298.0) -> Tuple[pd.DataFrame, pd.Series]:

### BEGIN YOUR SOLUTION ###

D_pi_d = in_df["d_pi_d"]

D_pi_D = in_df["d_pi_D"]

E_pi_d = in_df["e_pi_d"]

E_pi_D =in_df["e_pi_D"]

er = in_df["er(%)"]

I=np.exp((-1*E_pi_d)/R*T)+np.exp((-1*E_pi_D)/R*T)

c_d = np.exp((-1*E_pi_d)/R*T)/I

c_D = np.exp((-1*E_pi_D)/R*T)/I

D_pi_w = c_d*D_pi_d+c_D*D_pi_D

E_pi_w = c_d*E_pi_d+c_D*E_pi_D

DE_pi_w = D_pi_w*E_pi_w

delta_delta_G = -4.184*R*T*(er)

df1=in_df.drop(['d_pi_d','d_pi_D','e_pi_d','e_pi_D','er(%)'],axis=1)

df2=pd.DataFrame({'d_pi_w':D_pi_w,'e_pi_w':E_pi_w,'de_pi_w':DE_pi_w})

df3 = pd.concat([df2,df1],axis=1)

Turple = (df3,delta_delta_G)

return(Turple)

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 请务必运行下面的代码块先做初步检查.

代码

文本

[4]

X, y_true = prepare_data(nci_birman)

X.iloc[:5, :], y_true[:5]

(   d_pi_w  e_pi_w  de_pi_w  L_Alk  B1_Alk  B5_Alk  L_Ar  B1_Ar  B5_Ar
 0    3.54  -10.24 -36.2496   4.36    2.92    3.35  6.38   1.77   3.15
 1    3.54  -10.24 -36.2496   4.35    2.09    3.34  6.38   1.77   3.15
 2    3.54  -10.24 -36.2496   4.38    1.73    3.33  6.38   1.77   3.15
 3    3.54  -10.24 -36.2496   3.08    1.70    2.20  6.38   1.77   3.15
 4    3.49  -11.88 -41.4612   3.08    1.70    2.20  6.38   1.88   4.52,
 0   -17726.135734
 1   -34830.301793
 2   -39184.089517
 3   -52141.791077
 4   -52141.791077
 Name: er(%), dtype: float64)

代码

文本

完成初步检查后, 可以运行下面的代码块对该函数进行测试. 测试通过情况将关系到该任务的得分. 请勿修改该代码块中的任何内容.

代码

文本

[5]

### CAUTION: DO NOT MODIFY THIS CELL. ###

test_in_df = pd.concat([

pd.DataFrame({

"d_pi_d": [4.89, 4.91, 6.48, 7.00],

"d_pi_D": [7.00, 7.21, 7.22, 7.20],

"e_pi_d": [-3.48, -2.98, -1.80, -1.45],

"e_pi_D": [-1.48, -2.00, -1.96, -2.26]

}),

nci_birman.iloc[:4, 4:]

], axis=1)

test_labels = np.array([2.409283, 2.009314, 1.939568, 1.770392])

test_weighted_features = np.array([

[4.959643, -3.413988, -16.932160],

[5.279007, -2.822771, -14.901427],

[6.899685, -1.890743, -13.045529],

[7.159408, -2.095601, -15.003264]

])

def test_prepare_data_names(in_df):

X_, y_true_ = prepare_data(in_df)

assert set(X.columns) == {"d_pi_w", "e_pi_w", "de_pi_w", "L_Alk", "B1_Alk", "B5_Alk", "L_Ar", "B1_Ar", "B5_Ar"}, "wrong column names."

assert y_true_.name == "delta_delta_G", "wrong label name."

def test_prepare_data_labels(in_df, labels):

_, y_true_ = prepare_data(in_df)

np.testing.assert_allclose(

y_true_, labels, rtol=1.0e-3,

err_msg="wrong label values."

)

def test_prepare_data_weighted_features(in_df, weighted_features):

X_, _ = prepare_data(in_df)

np.testing.assert_allclose(

X_.loc[:, ["d_pi_w", "e_pi_w", "de_pi_w"]], weighted_features, rtol=1.0e-3,

err_msg="wrong weighted feature values."

)

try:

test_prepare_data_names(test_in_df)

print("names: PASSED")

except Exception as err:

print(err)

try:

test_prepare_data_labels(test_in_df, test_labels)

print("labels: PASSED")

except Exception as err:

print(err)

try:

test_prepare_data_weighted_features(test_in_df, test_weighted_features)

print("weighted features: PASSED")

except Exception as err:

print(err)

wrong label name.

Not equal to tolerance rtol=0.001, atol=0
wrong label values.
Mismatched elements: 4 / 4 (100%)
Max absolute difference: 52143.56146944
Max relative difference: 29453.11629822
 x: array([-17726.135734, -34830.301793, -39184.089517, -52141.791077])
 y: array([2.409283, 2.009314, 1.939568, 1.770392])

Not equal to tolerance rtol=0.001, atol=0
wrong weighted feature values.
Mismatched elements: 12 / 12 (100%)
Max absolute difference: 1.268736
Max relative difference: 0.084564
 x: array([[  4.89    ,  -3.48    , -17.0172  ],
       [  4.91    ,  -2.98    , -14.6318  ],
       [  7.217617,  -1.959485, -14.142809],
       [  7.2     ,  -2.26    , -16.272   ]])
 y: array([[  4.959643,  -3.413988, -16.93216 ],
       [  5.279007,  -2.822771, -14.901427],
       [  6.899685,  -1.890743, -13.045529],
       [  7.159408,  -2.095601, -15.003264]])

代码

文本

拆分与归一化

原数据没有直接分出训练集和测试集, 我们可以用sklearn.model_selection.train_test_split对数据集进行手动拆分.

任务2: 将28个样本先按照23:5的份额划分训练集与测试集, 再在各个数据集上分别对特征作均值-方差归一化:

$x_{n}^{'} := \frac{x _{n} - Mean _{n = 1}^{N} ( x )}{Std _{n = 1}^{N} ( x )}, n = 1, \dots, N,$

返回(归一化后的)训练集与测试集. 注意: 测试集上归一化时, 正确的做法是使用训练集上算出的均值和方差, 而非在测试集上重新计算.

任务2.1: 完成函数split_and_normalize()的编写, 输入:
- 原始数据集的特征X与标签y_true;
- 测试集大小test_size(默认为5).
返回拆分好的X_train, X_test, y_true_train, y_true_test, 其中, X_train与X_test进行了归一化处理. 为了和StandardScaler保持统一, 请确保返回的四个值均为np.array格式.
任务2.2: 在以下代码块的注释区回答问题: 为什么均值-方差归一化需要在训练集和测试集上分开进行? 如果先归一化再拆分, 会导致什么后果?

代码

文本

提示

归一化(在sklearn中叫做标准化, standardization)既可以像上机实习1演示的那样手动计算, 也可以通过函数sklearn.preprocessing.StandardScaler实现.
- 注意: 如果自己用pandas库计算mean和std, 请设置计算std时的参数ddof=0. 这是StandardScaler方法的计算原则.
警惕机器学习模型在数据集上“作弊”.

代码

文本

[87]

from sklearn.model_selection import train_test_split

from sklearn.preprocessing import StandardScaler

import numpy

def split_and_normalize(

X: pd.DataFrame, y_true: pd.Series, test_size: int=5

) -> Tuple[np.array, np.array, np.array, np.array]:

### BEGIN YOUR SOLUTION ###

X_train,X_test = train_test_split(X,train_size = 5)

y_true_train,y_true_test=train_test_split(y_true,train_size=5)

X_train1 = StandardScaler().fit_transform(X_train)

X_test1 = StandardScaler().fit_transform(X_test)

y_true_train1 = y_true_train

y_true_test1 = y_true_test

Tuple1 = (X_train1,X_test1,y_true_train1,y_true_test1)

return(Tuple1)

### END YOUR SOLUTION ###

### 任务2.2答题区(另起一行时请记得加注释符号#) ###

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 请务必运行下面的代码块先做初步检查.

代码

文本

[88]

X_train, X_test, y_true_train, y_true_test = split_and_normalize(X, y_true)

X_train.shape, X_test.shape, y_true_train.shape, y_true_test.shape, X_train[:5, :]

((5, 9),
 (23, 9),
 (5,),
 (23,),
 array([[-0.34635605,  1.25894362,  1.35941822,  2.        ,  2.        ,
          2.        , -1.06399778,  0.10454167, -0.96283258],
        [-1.53695496, -1.24770305, -0.96537699, -0.5       , -0.5       ,
         -0.5       ,  1.13561301,  0.10454167, -0.06418884],
        [-0.12988352, -0.61823124, -0.62588432, -0.5       , -0.5       ,
         -0.5       , -0.19438421,  1.41131261,  0.23268454],
        [ 0.51953407,  1.10719595,  1.06497965, -0.5       , -0.5       ,
         -0.5       , -1.06399778,  0.10454167, -0.96283258],
        [ 1.49366046, -0.50020528, -0.83313656, -0.5       , -0.5       ,
         -0.5       ,  1.18676675, -1.72493763,  1.75716947]]))

代码

文本

完成初步检查后, 可以运行下面的代码块对该函数进行测试. 测试通过情况将关系到该任务的得分. 请勿修改该代码块中的任何内容.

代码

文本

[89]

n_samples, n_features, test_size = 40, 9, 7

test_X = np.random.randn(n_samples, n_features)

test_y_true = np.random.randn(n_samples)

type(X)

def test_split_and_normalize_shapes(X, y_true):

X_train_, X_test_, y_true_train_, y_true_test_ = split_and_normalize(X, y_true, test_size=test_size)

assert (X_train_.shape, X_test_.shape) == ((n_samples - test_size, 9), (test_size, 9)), "wrong X shape."

assert (*y_true_train_.shape, *y_true_test_.shape) == (n_samples - test_size, test_size), "wrong y_true shape."

def test_split_and_normalize_normalized(X, y_true):

X_train_, X_test_, _, _ = split_and_normalize(X, y_true)

np.testing.assert_allclose(

X_train_.mean(axis=0), np.zeros(n_features), atol=1.0e-5,

err_msg="invalid X_train normalization: non-zero mean."

)

np.testing.assert_allclose(

X_train_.std(axis=0), np.ones(n_features), rtol=1.0e-3,

err_msg="invalid X_train normalization: non-unit std."

)

try:

test_split_and_normalize_shapes(test_X, test_y_true)

print("shapes after split: PASSED")

except Exception as err:

print(err)

try:

test_split_and_normalize_normalized(test_X, test_y_true)

print("normalization: PASSED")

except Exception as err:

print(err)

wrong X shape.
normalization: PASSED

代码

文本

模型训练与评估

我们在9个特征与23个训练样本上训练一个线性回归模型, 并评估其在5个测试样本上的预测表现.

任务3: 搭建模型, 完成训练与评估.

任务3.1: 完成函数train_model()的编写. 输入:
- 训练集X_train, y_true_train.
返回一个训练好的线性回归模型sklearn.linear_model.LinearRegression.
任务3.2: 编写函数evaluate_model(), 实现预测图、 $R^{2}$ 与RMSE值的报告(作图用的函数plot_prediction()已经给出). 输入:
- 模型model与数据集X, y_true;
- 评估模式mode, 可在plot和metrics中二选一, 如果传入plot则作预测图(不返回任何内容), 如果传入metrics则返回RMSE与 $R^{2}$ 值.
任务3.3: 调用你编写的函数evaluate_model(), 完成下述任务:
- 将训练集上的预测图保存为文件pred.png, 以备提交;
- 将测试集上的RMSE值与 $R^{2}$ 值打印出来.

代码

文本

[104]

from sklearn.linear_model import LinearRegression as LR

def train_model(X_train, y_true_train) -> object:

### BEGIN YOUR SOLUTION ###

reg = LR().fit(X_train, y_true_train)

y_pred_test = reg.predict(X_test)

return(reg)

train_model(X_train, y_true_train)

### END YOUR SOLUTION ###

LinearRegression()

代码

文本

完成任务3.1后, 请务必运行下面的代码块做初步检查: 输出的在训练集上的RMSE值应当与0很接近.

代码

文本

[138]

from sklearn.metrics import mean_squared_error

lr = train_model(X_train, y_true_train)

RMSE = mean_squared_error(lr.predict(X_train), y_true_train, squared=False)

r2 = r2_score(y_true, y_pred_test)

print(RMSE)

print(r2)

3.7666698642845686e-11
0.8079282296215671

代码

文本

检查完成后, 你可以继续完成任务3.2: 结合给出的函数plot_prediction()(已经给出, 取自上机实习1), 编写函数evaluate_model().

代码

文本

[110]

from matplotlib import pyplot as plt

from sklearn.metrics import r2_score

def plot_prediction(y_true: np.array, y_pred: np.array):

r2 = r2_score(y_true, y_pred)

# text annotation setup

plt.title(r"True values vs predicted values ($R^2$ = " + f"{r2:.4f}" + ")")

plt.xlabel("True values")

plt.ylabel("Predicted values")

# plot the scatter and line

plt.scatter(y_true, y_pred, c="red", marker="o")

plt.plot(y_true, y_true, "b--")

# show the plot!

plt.show()

代码

文本

[143]

def evaluate_model(model, X, y_true, mode: str):

if mode=="plot":

y_pred = model.predict(X)

plt.scatter(y_true, y_pred)

elif mode=="metrics":

y_pred = model.predict(X)

RMSE = mean_squared_error(y_true, y_pred, squared=False)

r2 = r2_score(y_true, y_pred)

return RMSE, r2

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 你可以继续完成任务3.3: 运行下面的代码块. 我们先在训练集上作出预测图, 再报告测试集上的RMSE与 $R^{2}$ .

代码

文本

[144]

evaluate_model(lr, X_train, y_true_train, mode="plot")

RMSE, r2 = evaluate_model(lr, X_test, y_true_test, mode="metrics")

print(f"RMSE: {RMSE}")

print(f"r2: {r2}")

RMSE: 67650.02114874308
r2: -3.9105580896759387

代码

文本

选做任务: 所选参数是足够合理的吗? 一方面, 实验数据本身可能就带有噪声; 另一方面, 所选取的几何参数数目还是过多了. 事实上, 若想进一步对特征进行降维(或挑选), 还可以有很多策略, 例如:

主元分析法(PCA, 将在后续课程学到).
$L_{1}$ -正则化. 理论课已经讲过, $L_{1}$ -范数约束相比于 $L_{2}$ -范数约束更易产生稀疏解, 我们借助正则化强度的控制, 可以让很多参数变为0, 对应的那些分量就相应被舍弃, 达到特征选择(feature selection)的目的.
根据各个特征可以解释的数据方差大小进行特征选择, 详情可阅读sklearn.feature_selection模块的官方文档.

你可以选择在下面的代码块中自由实践你感兴趣的特征选择方法, 并在注释的答题区中给出你的讨论.

代码

文本

[ ]

### 选做任务答题区(另起一行时请记得加注释符号#) ###

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

结果分析

代码

文本

任务4: 运行下面的代码块, 根据参数正负, 考察各个因素对对映选择性的影响, 在注释的答题区中写下你的分析.

$π$ - $π$ 堆积作用如何影响对映选择性? (从几何结构与能量两方面考虑) 这符合你的直觉吗?
烷基的几何参数如何影响对映选择性? 和芳基几何参数相比, 总体看, 谁的影响更大?

代码

文本

提示

$ΔΔ G^{‡}$ 越大, 表明两种对映异构体的动力学活性相差越大, 于是对映选择性越好.
你可能会得到一些反直觉的结果, 但这是模型过拟合(特征数目过多)导致的问题, 请如实地反映输出结果, 并写出你的疑惑.

代码

文本

[147]

lr.coef_

### 任务4答题区(另起一行时请记得加注释符号#) ###

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

array([ 50635.06035626,   7791.7273241 ,  -3562.31001067,   9919.29876713,
         9919.29876713,   9919.29876713, -26275.45201376,  30486.84862747,
        31987.33731434])

代码

文本

机器学习及其在化学中的应用

点个赞吧