连不上机器判断机器状态_实战机器学习:使用朴素贝叶斯判断分子药物致突变性...

"埃姆斯试验"(引用自维基百科)

埃姆斯试验（英语：Ames test，又称为细菌回复突变试验或安姆氏突变试验）为美国人布鲁斯·埃姆斯博士于1983年所提出的突变物测试方法。此法检测出175已知之致癌物，并发现其中有超过90%者皆显示为阳性(具突变性)，因此，埃姆斯试验也被做为先期测试方法之一。而现在正以埃姆斯试验为模式系统，进行各式化学物可能的毒性及致突变性之安全评估，但仍需其它的测试辅助测试而进一步确认。

我们也可以使用机器学习的手法在已有的数据上建立QSAR(Quantitative Structure-Activity(or Affinity) Relationship )模型，通过药物的分子式来预测一个未知的药物是否具有突变性。这里显然就变成了一个分类问题。最近在实装实战机器学习里面的算法，今天就拿朴素贝叶斯来试试手。

AMES的数据集可以在这篇文章的附加内容里面找到:

Benchmark Data Set for in Silico Prediction of Ames Mutagenicity​pubs.acs.org

Benchmark Data Set for in Silico Prediction of Ames Mutagenicity:

J. Chem. Inf. Model. 2009, 49, 2077. (DOI: 10.1021/ci900161g

或者以下地址下载(提取码q6h9):

https://pan.baidu.com/s/1gBG_84z9ArilB-i3-5dKaw

导入所需的library:

from rdkit import *
from rdkit import Chem
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
from math import *

读取数据集

df = pd.read_csv('/smiles_cas_N6512.smi',header=None, sep='t')
df.columns = ['smiles', 'CAS_NO', 'activity']
df.head()

可以得到以下结果:

第一列是化合物的smiles，第二个是CAS代号，第三个是关于该化合物的突变性信息，有突变性=1,没有突变性性=0。使用pandas工具可以可视化化合物的分子式:

这里，我们发现有一些分子是无法直接从smiles导出的(None)。其中的原因可能是smiles出现了一些语法错误。图中的smiles一部分构造混入了符号，这代表这些样本属于混合物。

我们可以先去除掉这部分化合物。

df['ROMol'] = df['ROMol'].map(lambda x: False if x == None else True)
del_index = df[df['ROMol'] == False].index
df2 = df.drop(del_index)
#判断是否可以导出为mol文件，如果出错的smiles会被判定为false然后drop掉(去除掉)

之后，将剩下的分子数据从smiles导出为mol文件之后我们使用maccskey来提取分子指纹。

mols = [Chem.MolFromSmiles(smile) for smile in df2['smiles']]
activity=df2['activity']
maccskeys = []
for m in mols:
    maccskey = [x for x in AllChem.GetMACCSKeysFingerprint(m)]
    maccskeys.append(maccskey)
maccskeys = np.array(maccskeys)
maccskeys.shape
#(6506, 167)有6506个分子，maccskey把分子分成了167个部分结构。

之后，划分数据集，以7:3的比例将全部的数据分为训练和测试数据集:

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(maccskeys, activity, test_size=0.3, random_state=0)

朴素贝叶斯与条件概率

已知w是由多个独立事件组成的。这里的独立事件是w1,w2,w3，由于我们使用分子指纹作为记述子，因此w1,w2,w3...这些独立事件表示的是某个分子的部分结构存在与否。

里面的i有两个取值1和0，也就是有或者没有突变性。因此这里的贝叶斯公式里的

的意思就是已知某个分子的部分结构存在的前提下，分子具有突变性的概率是多少。这个概率可以通过贝叶斯公式计算:

同样的，对于

就是分子具有药突变性的概率，反之

就是分子不具有突变性的概率。p(w)是某个独立事件发生的几率，也就是数据集中所有分子具有该部分结构的概率。

是反过来，知道分子有突变性的时候，某个部分结构存在的概率。

通过条件概率，如果对于一个分子来说，计算出来的拥有突变性的概率大于没有突变性的概率，那么我们判断它属于有突变性的分子，反之就没有突变性, 感觉很是神奇。

训练模型

那么我们的目的就转换成训练一个模型来计算分子具有突变性的概率:

def trainNaiveBayes(descriptor, label):
    #数据中的分子数目
    num_data=len(descriptor)
    #分子指纹的长度
    num_fingerprint=len(descriptor[0])
    #分子中具有突变性数据所占比例p(c=1)
    pAbusive=sum(label)/float(num_data)
    #创建一个长度为分子指纹等长的列表
    p0Num=np.ones(num_fingerprint);p1Num=np.ones(num_fingerprint)
    p0Denom=2.0;p1Denom=2.0
    #遍历数据中的所有分子
    for i in range(num_data):
        #如果某个分子具有突变性，也就是标签为1
        if label[i]==1:
            p1Num+=descriptor[i]
            #有突变性的分子中，具有某个部分结构的分子数量
            p1Denom+=sum(descriptor[i])
            #具有突变性的分子总数 
        else:
            #没有突变性的分子中，具有某个部分结构的分子数量
            p0Num+=descriptor[i]
            #不具有突变性的分子总数
            p0Denom+=sum(descriptor[i])
    #计算p(w|c1)
    p1Vect=p1Num/p1Denom #之后会改为np.log(p1Num/p1Denom)，便于计算 
    p0Vect=p0Num/p0Denom #之后会改为np.log(p0Num/p0Denom)
    return p0Vect,p1Vect,pActivity#pActivity也就是p(c)

之后计算所有分子具有突变性的概率:

p1=[]
for i in range(len(X_train)):
    p1.append(sum(X_train[i]*p1V)+pAV)
len(p1)

分子分子具有突变性的概率分布如下:

import plotly.offline as offline
offline.init_notebook_mode()
data = go.Histogram(x=p1, 
                    xbins=dict(start=0, end=1, size=0.01))
fig = dict(data=[data])
offline.iplot(fig)

首先，我们使用主成分分析PCA计算出我们数据集的第一主成分和第二主成分，把他们分别作为x和y轴，将我们的化学空间可视化。同时，使用计算出的分子具有突变性的概率概率作为代表每一个化合物的点的颜色:

import plotly
import plotly.graph_objs as go
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
pca.fit(X_train)
pca_result = pca.transform(X_train)
PC1=pca_result[:,0]
PC2=pca_result[:,1]

train = go.Scatter( 
    x = PC1,
    y = PC2,
    mode = 'markers',
    name = 'Test',
    marker = dict( 
        color = p1,
        size = '3',
        colorscale = 'Earth',
    )
)
data=[train]

layout = go.Layout(
    autosize = False,     
    width = 600,
    height = 600)
fig = go.Figure(data=data, layout=layout)
plotly.offline.init_notebook_mode()
plotly.offline.iplot(fig,filename='basic-scatter')

可以看到我们的化学空间如下:

可以看到越往左下，分布的分子颜色越深，象征着计算出的具有突变性的概率相对较低。

之后，计算出p0之后，得到我们的分类器如下:

def classifyNaiveBayes(Classify,p0Vec,p1Vec,pActivity):
    #使用朴素贝叶斯计算分子具有突变性的概率
    p1=sum(Classify*p1Vec)+pActivity#之后改为np.log(pActivity)，便于计算 
    p0=sum(Classify*p0Vec)+(1.0-pActivity)#之后改为np.log(1.0-pActivity)，便于计算 
    if p1>p0:
        return 1
    else:
        return 0

使用测试数据测试我们的朴素贝叶斯分类器:

y_pred=[]
for i in range(len(X_test)):
    y_pred.append(classifyNB(X_test[i].tolist(),p0V,p1V,pAb))
y_pred[0]==y_test.tolist()[0]
True    

运气很好，第一个分子分类正确~

我们可以调用一些sklearn的包来评价我们的分类结果:

from sklearn.metrics import accuracy_score
accuracy_score(y_test.tolist(), y_pred)
0.6956967213114754

达到及格分数。。。

混淆行列如下:

from sklearn.metrics import confusion_matrix
import seaborn as sns
cm=confusion_matrix(y_test.tolist(), y_pred, labels=[1, 0])
sns.heatmap(cm, annot=True, cmap='Blues')

思考

朴素贝叶斯的前提条件是"假设特征之间强（朴素）独立", 但是我们使用的分子指纹是将单一分子切分成许多个部分构造，各个部分构造之间是有很强的相互关系的，所以客观条件使得这次选用的记述子不太适用于朴素贝叶斯，导致预测精度没有达到非常满意的结果。