机器学习从蛋白序列预测蛋白分类（二）

半路出家学机器学习，先在自己熟悉的领域尝试，每天进步一点点，记录自己成长过程（正式开始前的小叨叨）
接（[机器学习从蛋白序列预测蛋白分类（一）]）继续分析
三，特征提取
其实在（机器学习从蛋白序列预测蛋白分类（一）中已经做过一次特征筛选，即从两个原始文件中提取出sequence和classification信息并整合至一个文件，这里特征提取主要针对蛋白序列。为何要对蛋白序列进行二次特征提取呢？是因为蛋白序列中包含了大量信息，比如20种必需氨基酸的组成顺序、每种氨基酸其实还携带着氨基酸残基大小、电性等信息。蛋白序列可以看成是以文本存储，所以处理方式类似于文本处理。
1，对数据data分训练集和测试集
采用常用的8：2分法， random_state = 1表示随机分配

X=data['sequence']
Y=data['classification']
X_train,X_test,Y_train,Y_test=train_test_split(X,Y,test_size = 0.2, random_state = 1)

2，设置特征提取方法，这里借鉴NLP中词袋模型（用CountVectorizer()函数）,这里用char_wb取词方式（只在句子内部分割，一个序列看作一个句子），ngram_range=(4,4)表示以4个氨基酸长度作为扫描窗口，产生词袋组合

vect=CountVectorizer(analyzer='char_wb',ngram_range=(4,4))

3，用上述特征提取方式对训练集进行提取：

vect.fit(X_train)
X_train_df=vect.transform(X_train)
X_test_df=vect.transform(X_test)
print(X_test_df)

输出，括号内第一个数字表示首个序列（从0开始编号），第二个数字表示该词在词袋模型中的索引，括号外的1表示索引为2969的词在首个序列中出现1次：

 (0, 2969) 1
  (0, 5775) 1
  ...
  (55773, 173528) 1
  (55773, 173553) 1

看一下CountVectorizer从蛋白序列提取出哪些基本单元

print(vect.get_feature_names()[-20:])

输出（倒数20个）：

['zhhh', 'ziar', 'zigi', 'ziwz', 'zkal', 'zkky', 'zknt', 'zkyh', 'zlik', 'zlzk', 'zpvm', 'zrgd', 'zrvi', 'ztvl', 'ztzk', 'zvbd', 'zvib', 'zvka', 'zwdl', 'zzvb']

四，模型构建
用上述提取好的特征训练模型，并计算准确度打分
选择朴素贝叶斯模型Naive Bayes

from sklearn.naive_bayes import MultinomialNB
model=MultinomialNB()
model.fit(X_train_df,Y_train)
NB_pred=model.predict(X_test_df)
accuracy_score=accuracy_score(Y_test, NB_pred)
print(accuracy_score)

输出：

0.7638863986803887

五，模型评估
1，混淆矩阵
由于这个项目的label有多个，属于多分类的混淆矩阵，分类标准以types中存储的为准

#confusion matrix
conf_mat=confusion_matrix(Y_test,NB_pred,label=types)
print(conf_mat)

输出：

[[7215  262   55 ...    2   13   11]
 [ 230 5563   69 ...    5    9   20]
 [ 240  166 5464 ...    1    2   13]
 ...
 [   0    1    0 ...  203    0    0]
 [   8   10    0 ...    0  165    1]
 [   6    5    6 ...    0    1  174]]

print(conf_mat.sum(axis=1))  #按行计算总和，用于归一化处理

[9274 7262 6895 3204 2324 2228 1652 1703 1641 1351 1284 1312  963  991
  953  884  859  710  690  622  686  601  614  607  589  583  526  541
  509  396  326  293  309  308  250  251  269  237  224  224  209  209
  211]

np.newaxis 表示加一个维度，加在行上还是列上，取决于np.newaxis在行的位置还是在列的位置，比如下面的例子就是加在列上具体用法（有实例，讲述清晰）

print(conf_mat_sum[:, np.newaxis])

输出：

[[9274]
 [7262]
 [6895]
 [3204]
 [2324]
 ...
  [ 224]
 [ 224]
 [ 209]
 [ 209]
 [ 211]]

归一化处理，以概率形式方便作图，并不影响混淆矩阵结果

conf_mat_Nor = conf_mat.astype('float') / conf_mat.sum(axis=1)[:, np.newaxis]
print(conf_mat_Nor)

输出,这里除法运算规则是[7215 262 55 … 2 13 11]中每一个元素除以9274，以下以此类推：

[[7.77981454e-01 2.82510244e-02 5.93055855e-03 ... 2.15656675e-04
  1.40176838e-03 1.18611171e-03]
 [3.16717158e-02 7.66042413e-01 9.50151473e-03 ... 6.88515560e-04
  1.23932801e-03 2.75406224e-03]
 [3.48078318e-02 2.40754170e-02 7.92458303e-01 ... 1.45032632e-04
  2.90065265e-04 1.88542422e-03]
 ...
 [0.00000000e+00 4.78468900e-03 0.00000000e+00 ... 9.71291866e-01
  0.00000000e+00 0.00000000e+00]
 [3.82775120e-02 4.78468900e-02 0.00000000e+00 ... 0.00000000e+00
  7.89473684e-01 4.78468900e-03]
 [2.84360190e-02 2.36966825e-02 2.84360190e-02 ... 0.00000000e+00
  4.73933649e-03 8.24644550e-01]]

来个直观点的，绘制热图看一下

plt.figure(figsize=(13,8))
sns.heatmap(conf_mat)
plt.show()

出图：

上图可见，索引为3和38的两个出现了错判

2，F1打分（F1-score是精确度（precision）和召回率（recall）的平衡）
在二分类问题中，recall也成为sensitivity，precision也成为specificity。
sklearn中F1_score位于metrics.classification_report模块中

from sklearn.metrics import classification_report
print(classification_report(Y_test,NB_pred,target_names=types))

输出：

                                                      precision    recall    f1-score   support
                              HYDROLASE       0.45      0.74      0.56       250
                          TRANSFERASE       0.60      0.82      0.70       211
                   OXIDOREDUCTASE       0.75      0.77      0.76       589
                         ...
            OSYNTHETIC PROTEIN       0.87      0.86      0.86      1351
                                    avg / total       0.79      0.76      0.77     55774                                                                                                                                               

最后唠叨一下最近的学习心得，也算是给自己的提醒：刚学的东西一定要自己敲代码实现，这是最直观的感触；另一点就是，代码实现之后短期内要整体再复盘一下，理解更上一个层面。