🌞欢迎来到AI+生物医疗的世界
🌈博客主页：卿云阁💌欢迎关注🎉点赞👍收藏⭐️留言📝
🌟本文由卿云阁原创！
📆首发时间：🌹2024年3月19日🌹
✉️希望可以和大家一起完成进阶之路！
🙏作者水平很有限，如果发现错误，请留言轰炸哦！万分感谢！
理论介绍
研究背景
微管是由 α- 和 β- 微管蛋白异二聚体以头尾相连的方式形成，是真核细胞中细胞骨架的关键元素。微管动态平衡的破坏可以诱导细胞周期停滞在 G2/M 期，最终导致细胞凋亡。因此，微管蛋白靶向剂已被广泛应用于癌症的治疗。秋水仙素 (Colchicine)是一种生物碱，属于 秋水仙碱位点抑制剂(Colchicine Binding Site Iinhibitors，CBSIs)的一种， 但对正常细胞系也具有较强的细胞毒性。因此，发现新型的 CBSIs ， 并用于癌症的治疗仍是抗癌药物研究领域的热点之一。CBSIs 活性预测模型的构建
本研究以 1076 个骨架新颖的秋水仙碱位点抑制剂为数据集，以 4:1 的比例将数据集随机分为训练集和测试集，基于分子指纹并采用机器学习方法构建了468 个理论预测模型。微管蛋白秋水仙碱位点抑制剂
研究成果模型构建和评估基于 CBSIs 活性预测模型先导化合物的发现
先导化合物 23g 抗癌作用机制的研究
代码实现
数据集准备
数据原始的样子
- 删除空值
- 去金属离子
- 去除不确定值，只保留=号的值；
- 所有单位要统一
- IC50/EC50/Ki/Kd求平均值代表最终活性
- 添加数据标签，活性数值>10uM的代表无活性，活性数据<或=10uM代表有活性。
读取数据import pandas as pd import rdkit from rdkit import Chem from standardiser import standardise import logging from rdkit.Chem import Descriptors METAL_ELEMENTS = ['Li', 'Be', 'Na', 'Mg', 'Al', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Cs', 'Ba', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'Fr', 'Ra', 'Lr', 'Ho'] # 加载数据 df = pd.read_csv('tubulin.csv')删除Smiles为空值的行
# 删除Smiles为空值的行 df.dropna(axis=0, subset=["Smiles"], inplace=True)去金属离子
#去金属离子 for metal in METAL_ELEMENTS: df = df[~df['Smiles'].str.contains(metal)] # df['Smiles'].str.contains(metal)取出包含这个离子的 print("no {}".format(metal),df.shape)分子标准化
for i in df.index: try: smi = df.loc[i, 'Smiles'] # print(smi) mol = Chem.MolFromSmiles(smi) mol = Chem.AddHs(mol) parent = standardise.run(mol) mol_ok_smi = Chem.MolToSmiles(parent) df.loc[i, 'Smiles'] = mol_ok_smi # print(i, 'done') except standardise.StandardiseException as e: logging.warning(e.message)删除重复值
df.drop_duplicates(keep='first', inplace=True) print('删除重复值df:', df.shape)删除Standard type中非'IC50 EC50 Ki Kd'的数据
df1=df[df['Standard Type'].isin(["IC50"])] # df2=df[df['Standard Type'].isin(["EC50"])] df3=df[df['Standard Type'].isin(["Ki"])] df4=df[df['Standard Type'].isin(["Kd"])] # result = df1.append(df2).append(df3).append(df4) # frames = [df1, df2, df3, df4] frames = [df1, df3, df4] df = pd.concat(frames)#删除Standard Relation中非'='的数据 df=df[df['Standard Relation'].isin(["'='"])] print('删除非=的df:', df.shape)#删除assay_type中非"B"的数据 df=df[df['Assay Type'].isin(["B"])] print('删除非B的df:', df.shape)#删除行中有空值的行 df=df.dropna(how='any',subset=(['Smiles','Standard Units','Molecule ChEMBL ID'])) print('删除存在空值的df:',df.shape)# 计算分子MW molweight = [] for smi in list(df['Smiles']): molweight.append(Descriptors.MolWt(Chem.MolFromSmiles(smi))) df['molecular_weight'] = molweight # 删除分子量大于1000的 df = df[ df['molecular_weight']<=1000 ] # 计算分子logP logP = [] for smi in list(df['Smiles']): logP.append(Descriptors.MolLogP(Chem.MolFromSmiles(smi))) df['logP'] = logP# 单位换算 print('去掉Da大于1000前', df.shape[0]) df = df[df['Molecular Weight'] <= 1000] print('去掉Da大于1000后', df.shape[0]) df1 = df[df['standard_units'].isin(["nM"])] print('单位为nM的df:', df1.shape) df2 = df[df['standard_units'].isin(["ug.mL-1"])] print('单位为ug/ml的df:', df2.shape) # df3 = df[df['standard_units'].isin(["uM"])] # print('单位为uM的df:',df3.shape) df2['standard_value'] = df2['standard_value'].astype(float) / df2['Molecular Weight'].astype(float) * 1000000 # g*e-6/L*e-3 / (g/mol) = e-9*mol/L Molecular Weight是M分子质量 df2['standard_units'] = "nM" # df3['translate'] = df3['standard_value']*1000 # df3['standard_value'] = df3['translate'] df = df1.append(df2) # df = df3.append(df4) print('单位换算后的df:', df.shape)# 根据相同SMILES的取standard_value的平均值 if df1.shape[0]!=0: df1_mean = df1.groupby('SMILES')['standard_value'].mean() print('df', df1_mean.shape) # 将df_mean转换成表 # df_mean = pd.DataFrame({'standard_value':df_mean}) df1_mean_dict = df1_mean.to_dict() print(df1_mean) df1 = df1.copy() df1['standard_value_mean'] = df1['SMILES'].apply(lambda x: df1_mean_dict[x]) # 添加 每个化合物的平均值 df1.drop('standard_value', axis=1, inplace=True) # 未平均的值已经求完，可舍弃 df1.drop_duplicates(subset=['SMILES'], inplace=True) # 对SMILES去重，因为现在有很多求过平均没删掉的数 df1 = df1.copy() df1.loc[df1['standard_value_mean'] < 10000, 'standard_value_mean'] = 1 df1.loc[df1['standard_value_mean'] >= 10000, 'standard_value_mean'] = 0 df1=df1[['SMILES','Molecular Weight','AlogP','standard_value_mean']]# 标签转换 df.loc[df['standard_value_mean']<=10000,'standard_value_mean']=1 df.loc[df['standard_value_mean']>10000,'standard_value_mean']=0 print('转换平均值后的df：',df.shape) df = pd.DataFrame(df,columns = ['Smiles','standard_value_mean']) # df = pd.DataFrame(df,columns = ['Smiles','standard_value_mean'，'molecular_weight'，'logp']) # 修改列名 df.columns = ['SMILES','LABEL'] # 保存清洗后的数据 df.to_csv('data.csv', index=None)建模及评估
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem import AllChem from rdkit.Chem import MACCSkeys from sklearn.neighbors import KNeighborsClassifier from sklearn.svm import SVC from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.metrics import roc_auc_score, confusion_matrix import math # 加载数据集 # 绝对路径 df = pd.read_csv('/home/hdpaii01/Day3/1data_extract/data.csv') # 相对路径 # df = pd.read_csv('../1data_extract/data.csv') # 生成分子对应的指纹 X = np.array( [AllChem.GetMorganFingerprintAsBitVect(Chem.MolFromSmiles(smi), 2, nBits=1024) for smi in list(df.iloc[:,0])]) y = df['LABEL'].values X_train,X_test,y_train,y_test = train_test_split(X,y,test_size = 0.2,random_state=0) knn = KNeighborsClassifier() # 定义超参数优化词典 param = {'n_neighbors':[3,5,7,9]} # param = {'n_neighbors':[3,5,7,9],'weights':['uniform', 'distance']} # 超参数优化 gc_knn = GridSearchCV(knn, param_grid=param, cv=5, scoring='roc_auc') gc_knn.fit(X_train,y_train)XGBoost
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem import AllChem from rdkit.Chem import MACCSkeys from sklearn.neighbors import KNeighborsClassifier from sklearn.svm import SVC from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.metrics import roc_auc_score, confusion_matrix import math import xgboost # from rdkit.Chem import Descriptors # # # def get_rdkit_des(mol): # features = [] # for desc_name, function in Descriptors.descList: # if desc_name == 'Ipc': # feature = function(mol, avg=True) # else: # feature = function(mol) # features.append(feature) # return np.asarray(features) # 加载数据集 df = pd.read_csv('C:/Users/Xsc/Desktop/backup/data.csv') # 生成分子对应的指纹 # Morgan # X = np.array( # [AllChem.GetMorganFingerprintAsBitVect(Chem.MolFromSmiles(smi), 2, nBits=1024) for smi in list(df.iloc[:, 1])]) # MACCS # from rdkit.Chem import MACCSkeys X = np.array([MACCSkeys.GenMACCSKeys(Chem.MolFromSmiles(smi)) for smi in list(df.iloc[:, 0])]) # Rdkit描述符 # X = np.array([get_rdkit_des(Chem.MolFromSmiles(smi)) for smi in list(df.iloc[:, 1])]) y = df['Label'].values X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0) # 构建模型 xgb = xgboost.XGBClassifier(use_label_encoder=False, random_state=123) # 定义超参数优化词典 param = {"max_depth": [3, 5, 10]} # 超参数优化 gc_knn = GridSearchCV(xgb, param_grid=param, cv=2, scoring='roc_auc') gc_knn.fit(X_train, y_train) # 得到最优参数 print(gc_knn.best_params_) # 在测试集上评估 y_pred = gc_knn.predict_proba(X_test)[:, 1] # 计算评估指标 auc_roc_score = roc_auc_score(y_test, y_pred) y_pred_print = [round(y, 0) for y in y_pred] tn, fp, fn, tp = confusion_matrix(y_test, y_pred_print).ravel() se = tp / (tp + fn) sp = tn / (tn + fp) # 也是R q = (tp + tn) / (tp + fn + tn + fp) mcc = (tp * tn - fn * fp) / math.sqrt((tp + fn) * (tp + fp) * (tn + fn) * (tn + fp)) P = tp / (tp + fp) F1 = (P * se * 2) / (P + se) BA = (se + sp) / 2 print(tp, tn, fn, fp, se, sp, mcc, q, auc_roc_score, F1, BA)原始分子，分子骨架，分子图结构 骨架分析
Y随机验证
打乱标签，重新训练
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem import AllChem from rdkit.Chem import MACCSkeys from sklearn.neighbors import KNeighborsClassifier from sklearn.svm import SVC from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.metrics import roc_auc_score, confusion_matrix import math import xgboost import pickle def score(y_true,y_pre): auc_roc_score = roc_auc_score(y_test, y_pre) y_pred_print = [round(y, 0) for y in y_pre] tn, fp, fn, tp = confusion_matrix(y_test, y_pred_print).ravel() se = tp / (tp + fn) sp = tn / (tn + fp) # 也是R q = (tp + tn) / (tp + fn + tn + fp) mcc = (tp * tn - fn * fp) / math.sqrt((tp + fn) * (tp + fp) * (tn + fn) * (tn + fp)) P = tp / (tp + fp) F1 = (P * se * 2) / (P + se) BA = (se + sp) / 2 return tp, tn, fn, fp, se, sp, mcc, q, auc_roc_score, F1, BA # 加载数据集 df = pd.read_csv('../1数据收集和分析/data.csv') # 生成分子对应的指纹 # Morgan X = np.array([AllChem.GetMorganFingerprintAsBitVect(Chem.MolFromSmiles(smi), 2, nBits=1024) for smi in list(df.iloc[:, 0])]) # 得到对应标签 y = df['LABEL'].values # 数据集拆分 X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0) # 构建模型 xgb = xgboost.XGBClassifier(use_label_encoder=False, random_state=123,max_depth = 5) # 模型训练 xgb.fit(X_train, y_train) # 模型保存 with open('xgb.pickle','wb') as f: pickle.dump(xgb,f)子结构
Draw.DrawMorganBit(mol, 1873, bi, useSVG=True)
使用shap统计每一个特征的贡献值
# 调用解释器的方法 计算shap shap_values = explainer.shap_values(X_train)比如预测这个分子有0.8的可能是有活性的，计算每个特征对这0.8的贡献。
某个样本的某个特征shap值为正数，会让最终预测值变大，即对于最终预测为活性（0.5为阈值，大于0.5为有活性，小于0.5为非活性）的贡献是正的，反之亦然，某个样本的某个特征shap值为负数，会让最终预测值变小，即对于最终预测为活性（0.5为阈值，大于0.5为有活性，小于0.5为非活性）的贡献是负的，上面的图只反映了对于一个，每个特征对于最终预测的贡献值，我们还想知道特征的大小和贡献之间的关系
这个图展示了每个样本（图上的每个点代表一个样本）的特征大小，以及对当前样本的预测值的贡献，横坐标代表特征的shap值（对预测值的贡献）每一行代表某个特征，图上展示的是对模型预测值影响最大的前20个特征，颜色代表特征值的大小，蓝色为0，红色为1（因为特征是指纹，只存在0和1，因此只有两种颜色）因此，对于上面这个图的解读，以Morgan746这个特征为例，我们只关注最上面一行，可以看到以shap值的0轴线为分界，蓝色点集中在左边，红色点集中在右边，说明：当分子中不存在这个子结构片段（蓝色代表0，指纹上的0代表分子不存在这个子结构）时，对预测值贡献为负值（使预测值减小），即分子不存在这个片段时更有可能对该靶标无活性；当分子中存在这个子结构片段（红色代表1，指纹上的1代表分子存在这个子结构）时，对预测值贡献为正值（使预测值变大），即分子存在这个片段时更有可能对该靶标有活性。
取绝对值
上图是shap_values，行是样本，列是特征，我们关注的是特征对最终预测值的影响大小，正负只代表该特征对预测值的影响是正向还是负向的，因此，需要取绝对值来避免正负抵消（对于每一个样本，该特征可能存在或者不存在，还是以Morgan746为例子，该特征存在的时候贡献为正，不存在的时候贡献为负，要计算一个特征的贡献大小，可以采用两种方式
- 单独取出存在该特征的样本（包含这个子结构的那些分子，该位上为1）的贡献求平均，不存在该特征（不包含这个子结构的那些分子，该位上为0）的样本求平均，然后和其他特征比较
- 直接取绝对值求平均
否则，因为这些样本包含存在和不存在这个结构，因此正负之间会抵消（直接求平均存在该特征子结构的样本和不存在该特征子结构的样本的贡献会抵消），最终得到的就不是特征对于预测值的影响
除了可以使用XGB之外还可以使用KNN，RF等算法。
扫一扫
10-29
1657
1657
10-31
149
149
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
