基于机器学习的化合物活性预测模型

利用化合物的结构与活性数据，基于RDKit和Python3的机器学习活性预测模型小示例。

代码示例：

---

#导入必须的包
#!/usr/bin/env python3
from rdkit.Chem import Descriptors
from rdkit.Chem import AllChem as ch
from rdkit.Chem import Draw as d
from rdkit import DataStructs
import pandas as pd
from rdkit.Chem import rdMolDescriptors as rdescriptors
from matplotlib.mlab import PCA
import matplotlib.pyplot as plt
import csv
from rdkit.SimDivFilters.rdSimDivPickers import MaxMinPicker
import sklearn
from rdkit.Chem import PandasTools, Descriptors, MolFromSmiles
from pandas import DataFrame
from sklearn.model_selection import train_test_split
from sklearn import svm
import numpy as np
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, GradientBoostingRegressor

#载入数据
with open ('receptor-bioactivity.txt', 'r') as f:
    content_raw = list((csv.reader(f, delimiter = '\t')))
len(content_raw)

#移除重复数据
content=[]
for i in range(0,len(content_raw)):
    if i == 0:
        chembl_id=content_raw[i][0]
        content.append(content_raw[i])
    elif content_raw[i][0]!=chembl_id:
        chembl_id=content_raw[i][0]
        content.append(content_raw[i])

#提取有用信息
names = []
smiles = []
activity = []
mols = []
for i in range(1,len(content)):
    names.append(content[i][0])
    smiles.append(content[i][5])
    activity.append(content[i][7])
    mols.append(ch.MolFromSmiles(content[i][5]))

#为深入分析创建数据集
def HBA(mol):
    return rdescriptors.CalcNumLipinskiHBA(mol)
def HBD(mol):
    return rdescriptors.CalcNumLipinskiHBD(mol)
def logP(mol):
    return Descriptors.MolLogP(mol)
def TPSA(mol):
    return Descriptors.TPSA(mol)
def num_rotatable_bonds(mol):
    return Descriptors.NumRotatableBonds(mol)
def num_heavy_atoms(mol):
    return mol.GetNumHeavyAtoms()
def MW(mol):
    return Descriptors.MolWt(mol)
data=[]
for i,mol in enumerate(mols):
    data.append([names[i], float(activity[i]), HBA(mol), HBD(mol), float(MW(mol)), logP(mol),float(TPSA(mol)),num_rotatable_bonds(mol),num_heavy_atoms(mol),smiles[i]])
dataframe=pd.DataFrame(data,columns=["CHEMBL_ID", "activity", "HBA", "HBD", "MW", "logP", "TPSA","rb",'smiles'])
dataframe.set_index("CHEMBL_ID",inplace=True)

PCA分析，数据降维也称主成分分析

#PCA分析
pca1=PCA(dataframe.drop(['smiles'],axis=1))
plt.rcParams["figure.figsize"] = [15, 15]
plt.style.use('ggplot')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('Ghrelin Receptor Homo sapiens PCA')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
X = [x[0] for x in pca1.Y]
Y = [y[1] for y in pca1.Y]
plt.scatter(X,Y)
plt.show()

#运用随机森林模型，并为其选择有用数据
model=dataframe.loc[:,["smiles", "activity"]]
desc_list = Descriptors.descList
model["pic50"] = model.activity.apply(lambda x : -1.0 * np.log10(x / 1.0e9))
for desc_name, function in desc_list:
    values = []
    for smiles in model["smiles"]:
        mol = MolFromSmiles(smiles)
        values.append(function(mol))
    model[desc_name] = values 
columns = [x[0] for x in desc_list[:30]]

#划分数据集，训练模型
train_data, test_data = train_test_split(model, test_size=0.3)
model2 = RandomForestRegressor(n_estimators=15)
model2.fit(train_data[columns], train_data["pic50"])

#测试模型，绘图
plt.rcParams["figure.figsize"] = [15, 15]
span = (1,12)
axes = plt.gca()
axes.set_xlim(span)
axes.set_ylim(span)
plt.plot((span[0],span[1]), (span[0],span[1]), linestyle='--')
plt.scatter(
    test_data["pic50"]
    , model2.predict(test_data[columns])
    , c='blue'
    , s=20
)
plt.show()

---