2021华为杯D题详细讲解：抗乳腺癌候选药物的优化建模

本题全部代码为python编写~每一个问题都编写了，写了我一晚上，也算是为2022的比赛练习一下吧。祝今年研究生能拿到国奖！

文章目录

D题

一、背景介绍

乳腺癌是目前世界上最常见，致死率较高的癌症之一。乳腺癌的发展与雌激素受体密切相关，有研究发现，雌激素受体α亚型（Estrogen receptors alpha, ERα）在不超过10%的正常乳腺上皮细胞中表达，但大约在50%-80%的乳腺肿瘤细胞中表达；而对ERα基因缺失小鼠的实验结果表明，ERα确实在乳腺发育过程中扮演了十分重要的角色。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。因此，ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。比如，临床治疗乳腺癌的经典药物他莫昔芬和雷诺昔芬就是ERα拮抗剂。
目前，在药物研发中，为了节约时间和成本，通常采用建立化合物活性预测模型的方法来筛选潜在活性化合物。具体做法是：针对与疾病相关的某个靶标（此处为ERα），收集一系列作用于该靶标的化合物及其生物活性数据，然后以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的定量结构-活性关系（Quantitative Structure-Activity Relationship, QSAR）模型，然后使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化。
一个化合物想要成为候选药物，除了需要具备良好的生物活性（此处指抗乳腺癌活性）外，还需要在人体内具备良好的药代动力学性质和安全性，合称为ADMET（Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性）性质。其中，ADME主要指化合物的药代动力学性质，描述了化合物在生物体内的浓度随时间变化的规律，T主要指化合物可能在人体内产生的毒副作用。一个化合物的活性再好，如果其ADMET性质不佳，比如很难被人体吸收，或者体内代谢速度太快，或者具有某种毒性，那么其仍然难以成为药物，因而还需要进行ADMET性质优化。为了方便建模，本试题仅考虑化合物的5种ADMET性质，分别是：1）小肠上皮细胞渗透性（Caco-2），可度量化合物被人体吸收的能力；2）细胞色素P450酶（Cytochrome P450, CYP）3A4亚型（CYP3A4），这是人体内的主要代谢酶，可度量化合物的代谢稳定性；3）化合物心脏安全性评价（human Ether-a-go-go Related Gene, hERG），可度量化合物的心脏毒性；4）人体口服生物利用度（Human Oral Bioavailability, HOB），可度量药物进入人体后被吸收进入人体血液循环的药量比例；5）微核试验（Micronucleus，MN），是检测化合物是否具有遗传毒性的一种方法。

二、数据集介绍及建模目标

本试题针对乳腺癌治疗靶标ERα，首先提供了1974个化合物对ERα的生物活性数据。这些数据包含在文件“ERα_activity.xlsx”的training表（训练集）中。training表包含3列，第一列提供了1974个化合物的结构式，用一维线性表达式SMILES（Simplified Molecular Input Line Entry System）表示；第二列是化合物对ERα的生物活性值（用IC50表示，为实验测定值，单位是nM，值越小代表生物活性越大，对抑制ERα活性越有效）；第三列是将第二列IC50值转化而得的pIC50（即IC50值的负对数，该值通常与生物活性具有正相关性，即pIC50值越大表明生物活性越高；实际QSAR建模中，一般采用pIC50来表示生物活性值）。该文件另有一个test表（测试集），里面提供有50个化合物的SMILES式。
其次，在文件“Molecular_Descriptor.xlsx”的training表（训练集）中，给出了上述1974个化合物的729个分子描述符信息（即自变量）。其中第一列也是化合物的SMILES式（编号顺序与上表一样），其后共有729列，每列代表化合物的一个分子描述符（即一个自变量）。化合物的分子描述符是一系列用于描述化合物的结构和性质特征的参数，包括物理化学性质（如分子量，LogP等），拓扑结构特征（如氢键供体数量，氢键受体数量等），等等。关于每个分子描述符的具体含义，请参见文件“分子描述符含义解释.xlsx”。同样地，该文件也有一个test表，里面给出了上述50个测试集化合物的729个分子描述符。
最后，在关注化合物生物活性的同时，还需要考虑其ADMET性质。因此，在文件“ADMET.xlsx”的training表（训练集）中，提供了上述1974个化合物的5种ADMET性质的数据。其中第一列也是表示化合物结构的SMILES式（编号顺序与前面一样），其后5列分别对应每个化合物的ADMET性质，采用二分类法提供相应的取值。Caco-2：‘1’代表该化合物的小肠上皮细胞渗透性较好，‘0’代表该化合物的小肠上皮细胞渗透性较差；CYP3A4：‘1’代表该化合物能够被CYP3A4代谢，‘0’代表该化合物不能被CYP3A4代谢；hERG：‘1’代表该化合物具有心脏毒性，‘0’代表该化合物不具有心脏毒性；HOB：‘1’代表该化合物的口服生物利用度较好，‘0’代表该化合物的口服生物利用度较差；MN：‘1’代表该化合物具有遗传毒性，‘0’代表该化合物不具有遗传毒性。同样地，该文件也有一个test表，里面提供有上述50个化合物的SMILES式（编号顺序同上）。
建模目标：根据提供的ERα拮抗剂信息（1974个化合物样本，每个样本都有729个分子描述符变量，1个生物活性数据，5个ADMET性质数据），构建化合物生物活性的定量预测模型和ADMET性质的分类预测模型，从而为同时优化ERα拮抗剂的生物活性和ADMET性质提供预测服务。

三、需解决问题

问题1. 根据文件“Molecular_Descriptor.xlsx”和“ERα_activity.xlsx”提供的数据，针对1974个化合物的729个分子描述符进行变量选择，根据变量对生物活性影响的重要性进行排序，并给出前20个对生物活性最具有显著影响的分子描述符（即变量），并请详细说明分子描述符筛选过程及其合理性。
问题2. 请结合问题1，选择不超过20个分子描述符变量，构建化合物对ERα生物活性的定量预测模型，请叙述建模过程。然后使用构建的预测模型，对文件“ERα_activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值预测，并将结果分别填入“ERα_activity.xlsx”的test表中的IC50_nM列及对应的pIC50列。
问题3. 请利用文件“Molecular_Descriptor.xlsx”提供的729个分子描述符，针对文件“ADMET.xlsx”中提供的1974个化合物的ADMET数据，分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型，并简要叙述建模过程。然后使用所构建的5个分类预测模型，对文件“ADMET.xlsx”的test表中的50个化合物进行相应的预测，并将结果填入“ADMET.xlsx”的test表中对应的Caco-2、CYP3A4、hERG、HOB、MN列。
问题4. 寻找并阐述化合物的哪些分子描述符，以及这些分子描述符在什么取值或者处于什么取值范围时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的ADMET性质（给定的五个ADMET性质中，至少三个性质较好）。

附件：
附件一：ERα_activity.xlsx
附件二：Molecular_Descriptor.xlsx
附件三：分子描述符含义解释.xlsx
附件四：ADMET.xlsx

简单描述题目

1.给出了1974个训练样本和50个测试样本，每个样本有729个特征
2. 每个训练样本有7个标签，分别是IC50值、pIC50值、和ADMET性质（包含5个标签）
3. IC50值、pIC50值是两个相关的连续变量。pIC50是IC50的负对数
4.ADMET性质的五个变标签都是布尔值

问题简化

问题一

根据文件“Molecular_Descriptor.xlsx”和“ERα_activity.xlsx”提供的数据，针对1974个化合物的729个分子描述符进行变量选择，根据变量对生物活性影响的重要性进行排序，并给出前20个对生物活性最具有显著影响的分子描述符（即变量），并请详细说明分子描述符筛选过程及其合理性。

简化描述问题：根据特征对IC50值和pIC50值影响的重要性进行排序，并给出前20个对IC50值和pIC50值最具有显著影响的特征

方案：特征选择。

问题二

请结合问题1，选择不超过20个分子描述符变量，构建化合物对ERα生物活性的定量预测模型，请叙述建模过程。然后使用构建的预测模型，对文件“ERα_activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值预测，并将结果分别填入“ERα_activity.xlsx”的test表中的IC50_nM列及对应的pIC50列。

简化描述： 选择不超过20个特征，构建IC50值和pIC50值的定量预测模型，并计算测试样本的IC50值和pIC50值。

方案： 机器学习预测算法构建

问题三

请利用文件“Molecular_Descriptor.xlsx”提供的729个分子描述符，针对文件“ADMET.xlsx”中提供的1974个化合物的ADMET数据，分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型，并简要叙述建模过程。然后使用所构建的5个分类预测模型，对文件“ADMET.xlsx”的test表中的50个化合物进行相应的预测，并将结果填入“ADMET.xlsx”的test表中对应的Caco-2、CYP3A4、hERG、HOB、MN列。

方案：机器学习分类模型中五分类。

问题四

寻找并阐述化合物的哪些分子描述符，以及这些分子描述符在什么取值或者处于什么取值范围时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的ADMET性质（给定的五个ADMET性质中，至少三个性质较好）。

方案： 最优化问题

求解

读取合并数据

读取数据：

import pandas as pd
data=pd.read_excel('ERα_activity.xlsx')

如下：

data2=pd.read_excel('Molecular_Descriptor.xlsx')

如下：

合并它俩：

#两个数据集做列合并
data_ER_X = data.drop(["SMILES"], axis=1) # 删除重复特征
data_concat = pd.concat([data,data2], axis=1) # 列合并
data_concat.head()

如下：

数据处理

提取特征和目标：

#读取特征和目标
X = data_concat.drop(['SMILES','IC50_nM','pIC50'],axis=1) # 特征
y = data_concat.loc[:, ['IC50_nM','pIC50']]

拆分数据集：

#拆分数据集
from sklearn.model_selection import train_test_split
X_train,X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=100)

分别提取两个特征的训练和测试数据：

y_IC50_train = y_train.loc[:,"IC50_nM"]
y_pIC50_train = y_train.loc[:,"pIC50"]
y_IC50_test = y_test.loc[:,"IC50_nM"]
y_pIC50_test = y_test.loc[:,"pIC50"]

构建模型

使用随机森林构建先预测IC50_nM

#使用随机森林的方法构建模型
from sklearn.ensemble import RandomForestRegressor
rf_model = RandomForestRegressor(n_estimators=100)
rf_model.fit(X_train,y_pIC50_train)

评估：

# 评估模型
from sklearn.metrics import r2_score
predict = rf_model.predict(X_train)
print("训练准确度为：:",r2_score(y_pIC50_train,predict))
predict = rf_model.predict(X_test)
print("测试准确度为:",r2_score(y_pIC50_test,predict))

结果为（有一些过拟合实际上）：

训练准确度为：: 0.960482848975931
测试准确度为: 0.7672674598150729

获取重要性排名：

features = X.columns # 读取列命
feature_importances = rf_model.feature_importances_ # 获取每个特征重要性
features_df = pd.DataFrame({'特征':features,'重要性':feature_importances}) # 
features_df.sort_values('重要性',inplace=True,ascending=False)

如下：

获取重要性最高的前20名：

features_df[:20]

如下：

特征重要性可视化：

# 可视化重要性排名
import seaborn as sns
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体设置-黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
sns.set(font='SimHei',font_scale=1.5)  # 解决Seaborn中文显示问题并调整字体大小


sns.set(rc={"figure.figsize":(21,8)}) # 设置大小
sns.barplot(features_df['特征'][:20],features_df['重要性'][:20]) # 柱形图
plt.ylabel('importants') # 纵轴标签
plt.xlabel('features') 
sns.despine(bottom=True)
plt.show()

如下：

问题一到此已经做完了。

选择不超过20个特征

使用前面选择好的特征再次训练：

# 读取最重要的20个特征
X_train = X_train.loc[:,features_df[:20]["特征"]]
X_test = X_test.loc[:,features_df[:20]["特征"]]
X_train.head()

如下：

选择不超过20个分子描述符变量，那么主成成分分析：

from sklearn.decomposition import PCA

#对20个特征进行pca主成分分析
pca = PCA(n_components=20)
pca.fit(X_train)

#查看各个特征的方差
pca.explained_variance_ratio_

如下：

array([3.97726755e-01, 2.49503335e-01, 1.80354983e-01, 6.99190239e-02,
       4.13153733e-02, 2.45000808e-02, 1.50078949e-02, 1.00870492e-02,
       5.65101235e-03, 2.22681992e-03, 1.83584025e-03, 1.27162646e-03,
       2.09948336e-04, 1.81292697e-04, 9.49557335e-05, 5.27544290e-05,
       3.52048327e-05, 2.13253115e-05, 3.80211709e-06, 9.22149350e-07])

可视化：

# 可视化
import numpy as np
plt.xlabel("Demensions")
plt.ylabel("explained_variance_ratio")
plt.plot([i for i in range(X_train.shape[1])],[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(X_train.shape[1])])
plt.legend()

如下：

选取贡献率99%的特征：

#选取贡献率99%方差的特征
pca = PCA(0.99)
pca.fit(X_train)
pca.n_components_ # 特征个数

如下（9个）：

筛选到九个特征后标准化：

X_train_pca = pca.transform(X_train) # 筛选到九个特征后标准化
X_test_pca = pca.transform(X_test)
print(X_train_pca.shape)

构建定量预测模型

标准化数据：

# 标准化数据
from sklearn.preprocessing import StandardScaler
std = StandardScaler()
std.fit(X_train)
X_train_std = std.transform(X_train)
X_test_std = std.transform(X_test)
print(X_train_std.shape)

如下：

建立模型并评估：

#建立SVM模型并对模型的R2值进行评估   rbf内核
from sklearn.svm import SVR
svr = SVR(tol=1e-5, kernel='rbf',C=1e1)
svr.fit(X_train_std,y_pIC50_train)
# 判断是否过拟合
print("支持向量机的训练精度：",svr.score(X_train_std, y_pIC50_train))
print("支持向量机的测试精度：",svr.score(X_test_std, y_pIC50_test))

如下：

使用网格搜索对决策树进行寻优并评估：

# 建立决策树模型并对模型的R2值进行评估
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import GridSearchCV # 网格搜索
param_grid = [
    {
        "criterion": ["squared_error","absolute_error","mae"],
        "max_depth": [10,20,30,40,50,100],
        'min_samples_leaf': [1,2,3,5,7,10],
        'min_impurity_decrease': [0.1,0.2,0.3]
    }
]
dtr = DecisionTreeRegressor()
dtr_grid = GridSearchCV(dtr, param_grid, n_jobs=-1, verbose=1)
dtr_grid.fit(X_train,y_pIC50_train)

print("最优超参数为：",dtr_grid.best_params_)
print("训练精度：",dtr_grid.best_estimator_.score(X_train, y_pIC50_train))
print("测试精度：",dtr_grid.best_estimator_.score(X_test, y_pIC50_test))

如下（说明决策树不适合）：

使用KNN算法试试：

from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import GridSearchCV
param_grid = [
    {
        "weights": ["uniform"],
        "n_neighbors": [i for i in range(1, 11)]
    },
    {
        "weights": ["distance"],
        "n_neighbors": [i for i in range(1, 11)],
        "p": [i for i in range(1,6)]
    }
]
# KNN模型
knn = KNeighborsRegressor()
grid_search = GridSearchCV(knn, param_grid, n_jobs=-1, verbose=1)

grid_search.fit(X_train_std, y_pIC50_train)

print("最优超参数为：",grid_search.best_params_)
print("KNN训练精度：",grid_search.best_estimator_.score(X_train_std, y_pIC50_train))
print("KNN测试精度：",grid_search.best_estimator_.score(X_test_std, y_pIC50_test))

结果如下：

使用随机森林模型:

# 随机森林
param_grid_rf = [
    {
        "max_depth": [30,40,50],
        'min_samples_leaf': [1,2,3],
        'n_estimators': [300,400,500],
    }
]
rf_model = RandomForestRegressor()
rf_grid = GridSearchCV(rf_model, param_grid_rf, n_jobs=-1, verbose=1)
rf_grid.fit(X_train,y_pIC50_train)


print("最优超参数为：",rf_grid.best_params_)
print("RandomForest训练精度",rf_grid.best_estimator_.score(X_train,y_pIC50_train))
print("RandomForest测试精度:",rf_grid.best_estimator_.score(X_test,y_pIC50_test))

如下：

预测并填入表格

读取需要预测的数据：

data3=pd.read_excel('Molecular_Descriptor.xlsx',sheet_name='test')

如下：

再读取另一个活性表格测试：

data4=pd.read_excel("./ERα_activity.xlsx",sheet_name='test')

如下：

读取前二十个特征：

#测试集的特征选择
X_final = data3.drop(['SMILES'],axis=1) # 删除描述
X_final = X_final.loc[:,features_df[:20]["特征"]] # 前二十个特征
X_final.head()

如下：

使用模型预测并保存：

# 使用训练好的随机森林模型进行预测
y_final_pIC50 = rf_grid.predict(X_final) # 预测数据
data4.loc[:,["pIC50"]] = y_final_pIC50 
y_final_IC50 = np.power(10,-y_final_pIC50)/(10**-9)
data4.loc[:,["IC50_nM"]] = y_final_IC50

data4.to_excel("ERα_activity_predict.xlsx",index=False)

查看预测结果：

data4 # 预测结果

如下：

构建五分类模型预测

分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型，并简要叙述建模过程。然后使用所构建的5个分类预测模型

读取数据并分割：

data_ADMET = pd.read_excel("./ADMET.xlsx") # 读取训练数据
X = pd.read_excel("./Molecular_Descriptor.xlsx") # 读取数据
X = X.drop(["SMILES"], axis=1) # 删除描述
y = data_ADMET.drop(["SMILES"], axis=1) #删除描述
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=666)

数据标准化：

# 数据标准化
std = StandardScaler()
std.fit(X_train)
X_train_std = std.transform(X_train)
X_test_std = std.transform(X_test)

提取五个特征的训练和测试：

# 提取五个特征的训练和测试
y_Caco_2_train = y_train.loc[:,"Caco-2"] 
y_CYP3A4_train = y_train.loc[:,"CYP3A4"] 
y_hERG_train = y_train.loc[:,"hERG"] 
y_HOB_train = y_train.loc[:,"HOB"] 
y_MN_train = y_train.loc[:,"MN"]

y_Caco_2_test = y_test.loc[:,"Caco-2"] 
y_CYP3A4_test = y_test.loc[:,"CYP3A4"] 
y_hERG_test = y_test.loc[:,"hERG"] 
y_HOB_test = y_test.loc[:,"HOB"] 
y_MN_test = y_test.loc[:,"MN"]

搭建模型：

#搭建神经网络模型
from keras.models import Sequential
from keras.layers import Dense,Activation
mlp_1 = Sequential()
mlp_1.add(Dense(units=365,activation="relu", input_dim=729))
mlp_1.add(Dense(units=365,activation="softmax"))
mlp_1.add(Dense(units=1,activation="sigmoid"))
mlp_1.summary()

如下：

对Caco-2进行训练

#对Caco-2进行训练
mlp_1.fit(X_train_std,y_Caco_2_train,epochs=400)

如下：

预测Caco_2：

from sklearn.metrics import accuracy_score
y_Caco_2_test_predict = mlp_1.predict(X_test_std) # 使用模型预测
for i in range(y_Caco_2_test_predict.shape[0]):
    if y_Caco_2_test_predict[i,:]>0.5:
        y_Caco_2_test_predict[i,:]=1
    else:
        y_Caco_2_test_predict[i,:]=0
y_Caco_2_test_predict = np.array(y_Caco_2_test_predict,dtype=int)
Caco_2_accuracy = accuracy_score(y_Caco_2_test,y_Caco_2_test_predict)
print("Caco_2精度:",Caco_2_accuracy)

如下：

同理预测CYP3A4：

#对CYP3A4进行训练并预测
mlp_2 = Sequential()
mlp_2.add(Dense(units=365,activation="relu", input_dim=729))
mlp_2.add(Dense(units=365,activation="softmax"))
mlp_2.add(Dense(units=1,activation="sigmoid"))
mlp_2.compile(optimizer='adam',loss='binary_crossentropy')
mlp_2.fit(X_train_std,y_CYP3A4_train,epochs=400)

y_CYP3A4_test_predict = mlp_2.predict(X_test_std)
for i in range(y_CYP3A4_test_predict.shape[0]):
    if y_CYP3A4_test_predict[i,:]>0.5:
        y_CYP3A4_test_predict[i,:]=1
    else:
        y_CYP3A4_test_predict[i,:]=0
y_CYP3A4_test_predict = np.array(y_CYP3A4_test_predict,dtype=int)
CYP3A4_accuracy = accuracy_score(y_CYP3A4_test,y_CYP3A4_test_predict)
print("CYP3A4精度:",CYP3A4_accuracy)

如下：

同理预测hERG：

mlp_3 = Sequential()
mlp_3.add(Dense(units=365,activation="relu", input_dim=729))
mlp_3.add(Dense(units=365,activation="softmax"))
mlp_3.add(Dense(units=1,activation="sigmoid"))
mlp_3.compile(optimizer='adam',loss='binary_crossentropy')
mlp_3.fit(X_train_std,y_hERG_train,epochs=400)

y_hERG_test_predict = mlp_3.predict(X_test_std)
for i in range(y_hERG_test_predict.shape[0]):
    if y_hERG_test_predict[i,:]>0.5:
        y_hERG_test_predict[i,:]=1
    else:
        y_hERG_test_predict[i,:]=0
y_hERG_test_predict = np.array(y_hERG_test_predict,dtype=int)
hERG_accuracy = accuracy_score(y_hERG_test,y_hERG_test_predict)
print("hERG精度:",hERG_accuracy)

如下：

同理HOB：

#对HOB进行训练并预测
mlp_4 = Sequential()
mlp_4.add(Dense(units=365,activation="relu", input_dim=729))
mlp_4.add(Dense(units=365,activation="softmax"))
mlp_4.add(Dense(units=1,activation="sigmoid"))
mlp_4.compile(optimizer='adam',loss='binary_crossentropy')
mlp_4.fit(X_train_std,y_HOB_train,epochs=400)

y_HOB_test_predict = mlp_4.predict(X_test_std)
for i in range(y_HOB_test_predict.shape[0]):
    if y_HOB_test_predict[i,:]>0.5:
        y_HOB_test_predict[i,:]=1
    else:
        y_HOB_test_predict[i,:]=0
y_HOB_test_predict = np.array(y_HOB_test_predict,dtype=int)
HOB_accuracy = accuracy_score(y_HOB_test,y_HOB_test_predict)
print("HOB精度:",HOB_accuracy)

如下：

同理MN：

#对MN进行训练并预测
mlp_5 = Sequential()
mlp_5.add(Dense(units=365,activation="relu", input_dim=729))
mlp_5.add(Dense(units=365,activation="softmax"))
mlp_5.add(Dense(units=1,activation="sigmoid"))
mlp_5.compile(optimizer='adam',loss='binary_crossentropy')
mlp_5.fit(X_train_std,y_MN_train,epochs=400)

y_MN_test_predict = mlp_5.predict(X_test_std)
for i in range(y_MN_test_predict.shape[0]):
    if y_MN_test_predict[i,:]>0.5:
        y_MN_test_predict[i,:]=1
    else:
        y_MN_test_predict[i,:]=0
y_MN_test_predict = np.array(y_MN_test_predict,dtype=int)
MN_accuracy = accuracy_score(y_MN_test,y_MN_test_predict)
print("MN精度:",MN_accuracy)

如下：

读取需要填充的表格：

test = pd.read_excel("./ADMET.xlsx",sheet_name='test')

如下：

使用模型预测并填充到表格：

data_ADMET_test = pd.read_excel("./ADMET.xlsx",sheet_name="test")
origin_molecular = pd.read_excel("./Molecular_Descriptor.xlsx",sheet_name='test')
X_final = origin_molecular.drop(['SMILES'],axis=1)

#对预测值做均值方差归一化
std = StandardScaler()
std.fit(X_final)
X_final_std = std.transform(X_final)

#预测Caco_2指标并填入excel表
Caco_2_predict = mlp_1.predict(X_final_std)
for i in range(Caco_2_predict.shape[0]):
    if Caco_2_predict[i,:]>0.5:
        Caco_2_predict[i,:]=1
    else:
        Caco_2_predict[i,:]=0
Caco_2_predict = np.array(Caco_2_predict,dtype=int)
data_ADMET_test.loc[:,["Caco-2"]] = Caco_2_predict


#预测CYP3A4指标并填入excel表
CYP3A4_predict = mlp_2.predict(X_final_std)
for i in range(CYP3A4_predict.shape[0]):
    if CYP3A4_predict[i,:]>0.5:
        CYP3A4_predict[i,:]=1
    else:
        CYP3A4_predict[i,:]=0
CYP3A4_predict = np.array(CYP3A4_predict,dtype=int)
data_ADMET_test.loc[:,["CYP3A4"]] = CYP3A4_predict

#预测hERG指标并填入excel表
hERG_predict = mlp_3.predict(X_final_std)
for i in range(hERG_predict.shape[0]):
    if hERG_predict[i,:]>0.5:
        hERG_predict[i,:]=1
    else:
        hERG_predict[i,:]=0
hERG_predict = np.array(hERG_predict,dtype=int)
data_ADMET_test.loc[:,["hERG"]] = hERG_predict


#预测HOB指标并填入excel表
HOB_predict = mlp_4.predict(X_final_std)
for i in range(HOB_predict.shape[0]):
    if HOB_predict[i,:]>0.5:
        HOB_predict[i,:]=1
    else:
        HOB_predict[i,:]=0
HOB_predict = np.array(HOB_predict,dtype=int)
data_ADMET_test.loc[:,["HOB"]] = HOB_predict

#预测MN指标并填入excel表
MN_predict = mlp_5.predict(X_final_std)
for i in range(MN_predict.shape[0]):
    if MN_predict[i,:]>0.5:
        MN_predict[i,:]=1
    else:
        MN_predict[i,:]=0
MN_predict = np.array(MN_predict,dtype=int)
data_ADMET_test.loc[:,["MN"]] =

预测并保存：

data_ADMET_test.to_excel("ADMET_predict.xlsx",index=False)

如下：

寻找使化合物对抑制ERα具有更好的生物活性…

用线性回归来求解各个特征的系数，寻找对抑制生物活性具有正性影响的分子描述符即可。

筛选出ADMET中五项指标其中有三项指标值为1的样本

# #筛选出ADMET中五项指标其中有三项指标值为1的样本
ADMET_data = pd.read_excel("ADMET.xlsx").drop(["SMILES"],axis=1)
ER_data = pd.read_excel("ERα_activity.xlsx").drop(["SMILES","IC50_nM"],axis=1)
ADMET_data.loc[:,["hERG","MN"]] = ADMET_data.loc[:,["hERG","MN"]].replace({0:1, 1:0})
comprehensive_value = []
for i in range(ADMET_data.shape[0]):
    if ADMET_data.loc[:,"Caco-2"][i]+ADMET_data.loc[:,"CYP3A4"][i]+ADMET_data.loc[:,"hERG"][i]+ADMET_data.loc[:,"HOB"][i]+ADMET_data.loc[:,"MN"][i]>2:
        comprehensive_value.append(1)
    else:
        comprehensive_value.append(0)

读取：

ADMET_data["comprehensive_value"]=comprehensive_value
ADMET_data = ADMET_data.loc[:,"comprehensive_value"]
Molecular_data = pd.read_excel("Molecular_Descriptor.xlsx").loc[:,features_df[:20]["特征"]]
Data_Comprehensive = pd.concat([ADMET_data,Molecular_data,ER_data],axis=1)
Data_Comprehensive.head()

如下：

读取特征和目标：

from collections import Counter
Counter(Data_Comprehensive.loc[:,"comprehensive_value"])

X = Data_Comprehensive.drop(["pIC50"],axis=1)
y = Data_Comprehensive.loc[:,"pIC50"]
print(X.shape,y.shape)
X.head()

如下：

训练：

#用线性回归来求解各个特征的系数
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X,y)
print(np.sort(-lr.coef_))
print(X.columns[np.argsort(lr.coef_)])
y_columns = X.columns[np.argsort(lr.coef_)][np.sort(lr.coef_)>0]
print(y_columns)

如下：

排序：

X_select = X.loc[:,y_columns]
data = pd.concat([X_select,y],axis=1)
data.head()

如下：

取值范围：

print("MDEC-23_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"MDEC-23"]),np.max(data.loc[:,"MDEC-23"])))
print("SHBint10_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"SHBint10"]),np.max(data.loc[:,"SHBint10"]))) 
# print("LipoaffinityIndex_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"LipoaffinityIndex"]),np.max(data.loc[:,"LipoaffinityIndex"])))
print("minHBint5_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"minHBint5"]),np.max(data.loc[:,"minHBint5"])))
print("minsssN_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"minsssN"]),np.max(data.loc[:,"minsssN"])))
print("nC_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"nC"]),np.max(data.loc[:,"nC"])))
# print("minHsOH_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"minHsOH"]),np.max(data.loc[:,"minHsOH"])))
print("VC-5_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"VC-5"]),np.max(data.loc[:,"VC-5"]))) 
print("MLFER_A_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"MLFER_A"]),np.max(data.loc[:,"MLFER_A"])))
print("maxHsOH_range:{:.3f}~{:.3f}".format(np.min(data.loc[:,"maxHsOH"]),np.max(data.loc[:,"maxHsOH"])))

如下：

做完~

总结

这个题其实并不难，理解题意就很好做。多尝试一些模型，选择一个更佳的模型就好了。

数模问题请私聊联系~