Protein Sequence Classification2

Introduction

Background
該數據集由許多不同類型的大分子組成，這些大分子具有生物意義。大多數數據記錄都是蛋白質。 DNA是RNA的前體，當轉錄時，蛋白質是在生物途徑和循環中直接相互作用的生物分子。蛋白質通常以一個或幾個工作為中心，這個工作由其家族類型定義。例如，我們可以使用來自水解酶組的蛋白質，其專注於催化水解（通過添加水來破壞鍵）以幫助促進蛋白質鍊或其他分子的鏈的破壞。另一個例子是作為載體蛋白的蛋白質，其允許其他分子如蔗糖，果糖或甚至水進入細胞內外。

Goals
由於這些蛋白質具有不同的家族類型，如果可以基於序列確定蛋白質的家族類型，則會出現問題。有一些著名的搜索引擎如BLAST具有這種能力，但如果機器學習方法可以根據蛋白質序列很好地對蛋白質家族進行分類是很有趣的。

1). Import Dataset
首先，加載函式庫和數據集。由於序列是一種字符串，這似乎是嘗試使用CountVectorizer的好機會。CountVectorizer是與NLP機器學習模型一起使用的特徵提取器(feature extractor)。

# 1). ----- Import Libraries and Datasets ------

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report


# Import Datasets
df_seq = pd.read_csv('../input/pdb_data_seq.csv')
df_char = pd.read_csv('../input/pdb_data_no_dups.csv')

print('Datasets have been loaded...')

2) Filter and Process Data
將數據加載到兩個單獨的pandas dataframes中，必須執行過濾器，項目和連接以將數據組合在一起。 這裡有一種好的函式庫dfply用來檢查，類似於R的dplyr函式庫。儘管如此，pandas提供了一種執行這些類似SQL的命令的方法。 首先過濾分類等於'蛋白質'的數據集，然後除去所有其他變量，除了pdb_data_seq.csv的structureId和classification以及pdb_data_no_dups.csv數據集中的structureId和sequence。

# 2). ----- Filter and Process Dataset ------

# Filter for only proteins
protein_char = df_char[df_char.macromoleculeType == 'Protein']
protein_seq = df_seq[df_seq.macromoleculeType == 'Protein']

# Select only necessary variables to join
protein_char = protein_char[['structureId','classification']]
protein_seq = protein_seq[['structureId','sequence']]
protein_seq.head()

protein_char.head()

protein_seq正確包含蛋白質，因為不同的序列以Methoinine開頭，以及包含比'ACTG'更多的字母的序列。 我們現在可以使用structureId作為索引來執行連接。 我們將使用pandas 'join' 方法。 為此，我們必須將每個dataframe的索引設置為“structureId”。

# Join two datasets on structureId
model_f = protein_char.set_index('structureId').join(protein_seq.set_index('structureId'))
model_f.head()

print('%d is the number of rows in the joined dataset' %model_f.shape[0])

這兩個dataframes正式加入了一個擁有346,325個蛋白質。 數據處理尚未完成，因為查看與columns關聯的mising是很重要的。

# Check NA counts
model_f.isnull().sum()

使用346,325種蛋白質，似乎只需刪除缺失值即可。

# Drop rows with missing values
model_f = model_f.dropna()
print('%d is the number of proteins that have a classification and sequence' %model_f.shape[0])

最後，重要的是要看一下分類的家庭群體的類型。

# Look at classification type counts
counts = model_f.classification.value_counts()
print(counts)

#plot counts
plt.figure()
sns.distplot(counts, hist = False, color = 'purple')
plt.title('Count Distribution for Family Types')
plt.ylabel('% of records')
plt.show()

家庭類型的數量似乎有很多分佈。 過濾具有特定家庭類型的一定數量的記錄可能是個好主意。 1,000似乎是一個可靠的數字，它允許機器學習模型學習特定類的模式。

# Get classification types where counts are over 1000
types = np.asarray(counts[(counts > 1000)].index)

# Filter dataset's records for classification types > 1000
data = model_f[model_f.classification.isin(types)]

print(types)
print('%d is the number of records in the final filtered dataset' %data.shape[0])

3). Train Test Split
在最終過濾數據集之後，必須對數據進行拆分以創建訓練和測試集。分割數據後，利用CountVectorizer創建由訓練數據集組成的字典非常重要。這將提取單個字符或字符子集以獲得特徵。關於使用CountVectorizer的一個重要注意事項是ngram_range的規範。

在蛋白質中，蛋白質中不是單個氨基酸可以識別它的目的。它是經由序列中氨基酸鍵形成的二級和三級結構。此外，鏈以及氨基酸的不同指出使用大於一個單元的特徵是重要的。因此，使用（4,4）的ngram_range似乎是特徵提取的合法選擇。這將提取長度為4的不同子集，允許氨基酸使用其鄰居來幫助我們進行分類。

# 3). ----- Train Test Split -----

# Split Data
X_train, X_test,y_train,y_test = train_test_split(data['sequence'], data['classification'], test_size = 0.2, random_state = 1)

# Create a Count Vectorizer to gather the unique elements in sequence
vect = CountVectorizer(analyzer = 'char_wb', ngram_range = (4,4))

# Fit and Transform CountVectorizer
vect.fit(X_train)
X_train_df = vect.transform(X_train)
X_test_df = vect.transform(X_test)

#Print a few of the features
print(vect.get_feature_names()[-20:])

4). Machine Learning Models
通過提取的功能，是時候使用機器學習模型了。 傳統上，樸素貝葉斯(Naive Bayes)方法適用於這些類型的計數向量化特徵(count vectorized features)。 Adaboost也將用來比較。

# 4). ------ Machine Learning Models ------

# Make a prediction dictionary to store accuracys
prediction = dict()

# Naive Bayes Model
from sklearn.naive_bayes import MultinomialNB
model = MultinomialNB()
model.fit(X_train_df, y_train)
NB_pred = model.predict(X_test_df)
prediction["MultinomialNB"] = accuracy_score(NB_pred, y_test)
print( prediction['MultinomialNB'])

# Adaboost
from sklearn.ensemble import AdaBoostClassifier
model = AdaBoostClassifier()
model.fit(X_train_df,y_train)
ADA_pred = model.predict(X_test_df)
prediction["Adaboost"] = accuracy_score(ADA_pred , y_test)
print(prediction["Adaboost"])

5). Visualize Metrics
似乎Naive Bayes在分類方面比Adaboost更好。 混淆矩陣的可視化和Navie Bayes預測的分類可以幫忙告訴我們模型表現是否不佳。

# 5). ----- Plot Confusion Matrix for NB -----

# Plot confusion matrix
conf_mat = confusion_matrix(y_test, NB_pred, labels = types)

#Normalize confusion_matrix
conf_mat = conf_mat.astype('float')/ conf_mat.sum(axis=1)[:, np.newaxis]

# Plot Heat Map
fig , ax = plt.subplots()
fig.set_size_inches(13, 8)
sns.heatmap(conf_mat)

混淆矩陣顯示標籤索引3被錯誤地分類為索引38。 根據下面列出的名稱，將這兩者混淆是有道理的。

print(types[3])
print(types[38])

最後，分類報告的矩陣顯示每個類的案例指標應該是富有洞察力的

#Print F1 score metrics
print(classification_report(y_test, NB_pred, target_names = types))

Reasons for Model Error
蛋白質通常可以是一種酶，或信號蛋白，結構和各種其他選擇。一小部分蛋白質傾向於具有非常相似的特徵，因為一些蛋白質與其他蛋白質結合在相似的區域中。例如，水解酶和水解酶抑製劑蛋白將具有相似的結構，因為它們將靶向非常相似的區域。這反映在混淆矩陣和熱圖(heat map)中。基因調節蛋白(Gene regulator proteins)與RNA結合蛋白，DNA結合蛋白以及轉錄蛋白具有相似性。最重要的是要注意，因為該模型最多只使用4種氨基酸的特徵。在理論上使用更高程度氨基酸的可能性應該能夠產生更高的準確度。

Future Work
該模型肯定有改進的餘地。利用諸如pH，分子量和其他組成的因素可能能夠產生關於家族群的更多信息。此外，如果可能的話，增加ngram_range的長度以包含超過4個字符，以允許氨基酸之間更高的相互作用，如現實中所反映的



參考
https://www.kaggle.com/abharg16/predicting-protein-classification