目录
1 知识入门
2 LinearRegression实现
2.1 语法
LinearRegression（fit_intercept,normalize,copy_X,n_jobs）2.2 参数讲解
fit_intercept :
布尔类型，初始为True。决定在这个模型中是否有intercept，即偏移量，即类似于线性函数y = w1x1 + w0 中的w0。 如果False则无。
normalize:
布尔类型，初始为False。如果fit_intercept设置为False，那这个参数就会被忽略。反之，设为True，则模型在回归之前，会对特征集X减去平均数并除以L2范式(没懂)，理解为一种标准化的规则。如果设为了False，而你又想标准化特征集合，则需要使用 sklearn.preprocessing.StandardScaler类来进行预处理。
copy_X：
布尔类型，初始化为True。True则，特征集合不变，反之会被复写。
n_jobs:
The number of jobs to use for the computation
初始为None，表示用1个处理器计算；-1代表所有处理器，只用于多个目标集问题的提速和大型问题。
2.3 属性
coef_：权重矩阵，理解为线性函数y = w1x1 + w0 中的W1
intercept_ :偏移量，理解为线性函数y = w1x1 + w0 中的W0
rank_: 特征矩阵的秩
singular_:特征矩阵的奇异值
3 算例及Python代码实现
3.1 算例
这个算例是2021年“华为杯”比赛的题目：
乳腺癌是目前世界上最常见，致死率较高的癌症之一。乳腺癌的发展与雌激素受体密切相关，有研究发现，雌激素受体α亚型（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性质提供预测服务。
3.2 问题
本次分析第四题
问题4. 寻找并阐述化合物的哪些分子描述符，以及这些分子描述符在什么取值或者处于什么取值范围时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的ADMET性质（给定的五个ADMET性质中，至少三个性质较好）。
3.3 Python代码实现
'''导入相关库'''
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split #拆分数据集
from sklearn.ensemble import RandomForestRegressor #随机森林
from sklearn.metrics import r2_score #评估模型
import seaborn as sns
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文字体设置-黑体
plt.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
sns.set(font='SimHei',font_scale=1.5) # 解决Seaborn中文显示问题并调整字体大小
'''1 获取重要性排名：'''
'''(1)读取数据'''
data=pd.read_excel('ERα_activity.xlsx')
data2=pd.read_excel('Molecular_Descriptor.xlsx')
#====两个数据集做列合并==========
data_ER_X = data.drop(["SMILES"], axis=1) # 删除重复特征
data3 = pd.concat([data,data2], axis=1) # 列合并
#res=data3
#res.to_excel('合并数据.xlsx')
#data3=pd.read_excel('合并数据.xlsx')
'''(2)数据处理'''
#===1)提取特征和目标：============
X = data3.drop(['SMILES','IC50_nM','pIC50'],axis=1) # 特征
y = data3.loc[:, ['IC50_nM','pIC50']] # 目标
#===2)拆分数据集==============
X_train,X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=100)
#====3）分别提取两个特征的训练和测试数据：===========
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"]
'''(3)构建模型'''
#使用随机森林构建先预测IC50_nM
#使用随机森林的方法构建模型
rf_model = RandomForestRegressor(n_estimators=100)
rf_model.fit(X_train,y_pIC50_train) # 预测IC50_nM的模型训练
# 评估模型
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))
'''(4)获取重要性排名：'''
features = X.columns # 读取列命
feature_importances = rf_model.feature_importances_ # 获取每个特征重要性
te_zheng_df = pd.DataFrame({'特征':features,'重要性':feature_importances}) #
te_zheng_df.sort_values('重要性',inplace=True,ascending=False) # 排序
print(te_zheng_df)
#获取重要性最高的前20名：
print(te_zheng_df[:20])
'''2 筛选出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}) #互换0和1值
man_zu_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:
man_zu_value.append(1)
else:
man_zu_value.append(0)
ADMET_data["man_zu_value"]=man_zu_value
ADMET_data = ADMET_data.loc[:,"man_zu_value"]
Molecular_data = pd.read_excel("Molecular_Descriptor.xlsx").loc[:,te_zheng_df[:20]["特征"]]
Data_Comprehensive = pd.concat([ADMET_data,Molecular_data,ER_data],axis=1)
#print(Data_Comprehensive.head())
'''3 读取特征和目标：'''
from collections import Counter #python自带计数器
Counter(Data_Comprehensive.loc[:,"man_zu_value"])
X = Data_Comprehensive.drop(["pIC50"],axis=1)
y = Data_Comprehensive.loc[:,"pIC50"]
print(X.shape,y.shape)
'''4 用线性回归来求解各个特征的系数'''
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X,y)
print(np.sort(-lr.coef_))
print(X.columns[np.argsort(lr.coef_)])
XX_columns = X.columns[np.argsort(lr.coef_)][np.sort(lr.coef_)>0]
print(XX_columns)
'''5 排序：'''
X_select = X.loc[:,XX_columns]
data = pd.concat([X_select,y],axis=1)
print(data.head())
data.head().to_csv('排序.csv')
'''6 取值范围'''
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"]))) #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"])))3.4 结果
D:\编程软件\Lib\project\Scripts\python.exe "C:/Users/20796/Downloads/2 每日数学建模学习/练习题程序/华为杯/2021年D题/第五题.py"
训练准确度为：: 0.9605400522521138
测试准确度为: 0.7662024651056416
特征 重要性
659 MDEC-23 0.099459
406 minsssN 0.076583
587 LipoaffinityIndex 0.071537
357 minHsOH 0.035872
39 BCUTc-1l 0.032639
.. ... ...
522 maxsssNHp 0.000000
314 SssssssS 0.000000
315 SSm 0.000000
317 SsGeH3 0.000000
364 minHssNH2p 0.000000
[729 rows x 2 columns]
特征 重要性
659 MDEC-23 0.099459
406 minsssN 0.076583
587 LipoaffinityIndex 0.071537
357 minHsOH 0.035872
39 BCUTc-1l 0.032639
476 maxHsOH 0.032489
639 nHBAcc 0.031353
652 MLogP 0.027206
11 nC 0.025845
56 C1SP2 0.025653
673 MLFER_A 0.023056
351 minHBint5 0.010904
291 SsOH 0.010305
238 SHsOH 0.010104
233 SHBint6 0.009305
716 TopoPSA 0.009223
79 VC-5 0.009090
237 SHBint10 0.007814
661 MDEC-33 0.007308
531 maxssO 0.007132
(1974, 21) (1974,)
[-2.27738763e+00 -1.75667314e+00 -1.14019313e+00 -8.31994737e-01
-2.77316882e-01 -1.20447439e-01 -8.55746980e-02 -8.03030264e-02
-5.57762963e-02 -2.54642345e-02 -1.11323525e-02 -5.19576386e-03
-7.77450168e-04 7.84512111e-03 4.56472495e-02 5.83745717e-02
1.80026275e-01 2.20208593e-01 2.65991765e-01 2.95358155e-01
3.37575048e-01]
Index(['MLogP', 'nHBAcc', 'SHsOH', 'minHsOH', 'C1SP2', 'SsOH', 'maxssO',
'SHBint6', 'MDEC-33', 'SHBint10', 'TopoPSA', 'MDEC-23', 'nC',
'LipoaffinityIndex', 'minHBint5', 'man_zu_value', 'minsssN', 'VC-5',
'MLFER_A', 'maxHsOH', 'BCUTc-1l'],
dtype='object')
Index(['MDEC-33', 'SHBint10', 'TopoPSA', 'MDEC-23', 'nC', 'LipoaffinityIndex',
'minHBint5', 'man_zu_value', 'minsssN', 'VC-5', 'MLFER_A', 'maxHsOH',
'BCUTc-1l'],
dtype='object')
MDEC-33 SHBint10 TopoPSA ... maxHsOH BCUTc-1l pIC50
0 9.238227 0.000000 67.23 ... 0.469126 -0.360525 8.602060
1 9.238227 0.000000 67.23 ... 0.449126 -0.360530 8.124939
2 10.328977 9.842059 87.46 ... 0.516534 -0.361379 8.508638
3 8.529910 0.000000 67.23 ... 0.456486 -0.360530 8.408935
4 8.529910 0.000000 67.23 ... 0.473631 -0.360530 8.130768
[5 rows x 14 columns]
MDEC-23_range:0.000~54.020
minHBint5_range:-1.198~11.555
minsssN_range:0.000~2.735
nC_range:7.000~95.000
VC-5_range:0.000~1.481
MLFER_A_range:-0.856~7.754
maxHsOH_range:0.000~0.853
Process finished with exit code 0
699
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
