抗乳腺癌候选药物的优化建模

抗乳腺癌候选药物的优化建模

背景介绍

  乳腺癌是目前世界上最常见，致死率较高的癌症之一。乳腺癌的发展与雌激素受体密切相关，有研究发现，雌激素受体α亚型（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），是检测化合物是否具有遗传毒性的一种方法。

需解决问题

  问题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性质中，至少三个性质较好）。

解题

问题一

  针对1974个化合物的729个分子描述符信息特征数据，通过python对该数据采用主成分分析方法，得出的特征值及特征向量结果，并计算肯德尔相关系数：

  计算主成分贡献率及累计贡献率,得到前20个主成分的贡献率已达88%。
  给出IC50和pIC50之间的关系方程式。

问题二

  通过E和M两个表中smiles字段进行关联，用trainning数据集训练bagging回归模型，然后应用到测试集上，利用问题一求解的二十个变量，使用M数据文件中的训练集trainning数据与题中给好的预测结果E数据trainning文件训练模型，应用到M和E的测试集test上，问题一因为ic50和plc50成负对数关系，所以只要预测一个即可，另一个就出来了。

问题三

  使用trainning中的数据训练，将训练好准确率高的模型应用到“ADMET.xlsx”的test表中的50个化合物中，得出50个化合物的预测结果，该文采用机器学习中无监督学习下的分类算法，众多分类算法中我们首先采用KNN算法，后续继续挑选几个分类算法进行调参后找到最优的分类算法。

问题四

AHP法
  层次分析法(Analytic Hierarehy Process，AHP)是20世纪70年代由美国学者A.L.saaty创立的一种定性与定量分析相结合的多目标评价决策方法。它基于多目标的层次结构，根据主观判断计算一系列备选方案的相对重要程度，由顶层依次向下计算，通过决策者为每一水平和子水平提供主观两两相对重要性的判断，为每一单元创立向下的成对比较矩阵。通过计算比较矩阵的特征向量得到同层次各元素对上一层次同一单元的相对重要性，然后再按照从底层依次向上的顺序，计算综合重要度，最后得到各备选方案(决策单元)的排序值。因而AHP法被广泛应用于质量控制系统、优先级评价、企业发展规划的选择方面，它适用于长期合作伙伴的评价、选择。
  AHP首先把问题层次化，然后根据问题的性质和需要达到的总目标，将问题分解为不同的组成因素，并按照因素间的相互关联影响以及隶属关系将因素按不同层次聚集组合，形成一个层次分析模型，并最终把系统分析归结为最低层相对最高层的相对重要性权值的确定或相对优劣次序的排序问题。其模型的建立如下:

程序代码

getSortdata.py

import openpyxl
import numpy as np
import scipy.stats as ss

# build a excel file obj and active it
excel = openpyxl.load_workbook('肯德尔相关系数.xlsx').active
# set list to store all variables name
var_name = []

# get variables names
for line in excel.iter_rows(min_row=1, max_row=1, min_col=2, max_col=730):
    temp = []
    for cell in line:
        temp.append(cell.value)
    var_name.append(temp)

# get names
var_name = var_name[0]

# set list to store kendall matrix
kendall_data = []

# get kendall attr data
for line in excel.iter_rows(min_row=2, max_row=2, min_col=2, max_col=730):
    temp = []
    for cell in line:
        temp.append(cell.value)
    kendall_data.append(temp)

# store data after processing
data_deal = []

# data process
for item in kendall_data[0]:
    if item == None:
        data_deal.append(0)
    else:
        data_deal.append(abs(item))

# data normalize
data = np.loadtxt('dataAttr.txt')
norm_data = ss.zscore(data)

# data deal
norm_data_new = []
for item in norm_data:
    norm_data_new.append([round(d, 2) for d in item])
norm_data_new = np.array(norm_data_new)

# replace nan data
nan = np.isnan(norm_data_new)
norm_data_new[nan] = 0

# correlation_coefficient value
correlation_coefficient = np.corrcoef(norm_data_new.T)
nan = np.isnan(correlation_coefficient)
correlation_coefficient[nan] = 0

# eigenvalue and eigenvector
eigenvalue, eigenvector = np.linalg.eig(correlation_coefficient)

# contribution_rate
contribution_rate = eigenvalue / eigenvalue.sum()


def looper(limit):
    cols = ['x1', 'x2', 'x3', 'x5', 'x6', 'x7', 'x8', 'x9']
    for i in range(len(cols)):
        data1 = data[cols]
        x = sm.add_constant(data1) #生成自变量
        y = data['y'] #生成因变量
        model = sm.OLS(y, x) #生成模型
        result = model.fit() #模型拟合
        pvalues = result.pvalues #得到结果中所有P值
        pvalues.drop('const',inplace=True) #把const取得
        pmax = max(pvalues) #选出最大的P值
        if pmax>limit:
            ind = pvalues.idxmax() #找出最大P值的index
            cols.remove(ind) #把这个index从cols中删除
        else:
            return result
#多元线性的效果偏差极大对比实际数据之后
def linearRegression(data_X,data_Y,learningRate,loopNum):

W = np.zeros(shape=[1, data_X.shape[1]])

# W的shape取决于特征个数，而x的行是样本个数，x的列是特征值个数

# 所需要的W的形式为 行=特征个数，列=1 这样的矩阵。但也可以用1行，再进行转置：W.T

# X.shape[0]取X的行数，X.shape[1]取X的列数

b = 0

#梯度下降

for i in range(loopNum):

W_derivative = np.zeros(shape=[1, data_X.shape[1]])

b_derivative, cost = 0, 0

WXPlusb = np.dot(data_X, W.T) + b # W.T：W的转置

W_derivative += np.dot((WXPlusb - data_Y).T, data_X) # np.dot:矩阵乘法

b_derivative += np.dot(np.ones(shape=[1, data_X.shape[0]]), WXPlusb - data_Y)

cost += (WXPlusb - data_Y)*(WXPlusb - data_Y)

W_derivative = W_derivative / data_X.shape[0] # data_X.shape[0]:data_X矩阵的行数，即样本个数

b_derivative = b_derivative / data_X.shape[0]

W = W - learningRate*W_derivative

b = b - learningRate*b_derivative

cost = cost/(2*data_X.shape[0])

if i % 100 == 0:

print(cost)

print(W)

print(b)

result = looper(0.05)
result.summary()

# build dict, as variable name: value
var_imp = dict(list(zip(var_name, data_deal)))

# dict sort for get 20 variables
var_sort = sorted(var_imp.items(), key = lambda item:(item[1]), reverse=False)
var_sort = [item for item in var_sort if item[1] != 0]
print('all var sort：\n', var_sort)
# write 20 variables to file
with open('var_sort_20.txt', 'a', encoding='utf-8') as f:
    for step,data in enumerate(var_sort):
        if step < 20:
            f.write(data[0] + '  ' + str(data[1]) + '\n')

getAttr.py

import openpyxl

# build a excel file obj and active it
excel = openpyxl.load_workbook('Molecular_Descriptor.xlsx').active

# set list to store 729 variables attr
data = []

# read 729 variables attr data
for line in excel.iter_rows(min_row=2, max_row=1975, min_col=2, max_col=730):
    temp = []
    for cell in line:
        temp.append(cell.value)
    data.append(temp)

# save attr data to txt file for next deal
with open('dataAttr.txt', 'a', encoding='utf-8') as f:
    for step,item in enumerate(data):
        for d in item:
            f.write(str(d)+' ')
        f.write('\n')