Protein Sequence Classification1

Import Dataset, drop NaN's, select Proteins
這個筆記本顯示瞭如何根據氨基酸序列對蛋白質家族進行分類。 這項工作基於當前自然語言處理（NLP）中深度學習模型的成功，並假設蛋白質序列可以被視為一種語言。 請注意，此任務有一些值得注意的搜索引擎，如BLAST。

在我們深入研究之前，先對數據進行一些預處理：

1. 在structureId上合併
2. 刪除沒有標籤的行
3. 刪除沒有序列的行
4. 選擇蛋白質

import numpy as np
import pandas as pd

# Merge the two Data set together
df = pd.read_csv('../input/pdb_data_no_dups.csv').merge(pd.read_csv('../input/pdb_data_seq.csv'), how='inner', on='structureId')
# Drop rows with missing labels
df = df[[type(c) == type('') for c in df.classification.values]]
df = df[[type(c) == type('') for c in df.sequence.values]]
# select proteins
df = df[df.macromoleculeType_x == 'Protein']
df.reset_index()
df.shape

Preprocessing and visualization of dataset
理想情況下：為了進行比較，我還決定只關注那些實例數大於1000的類（如Akil的這個內核中），它對應於43個最常見的類別。

但是：對於這樣的大數據集4個CPU跑一個小時是不夠的，而只考慮10個最常見的類別。

import matplotlib.pyplot as plt
from collections import Counter

# count numbers of instances per class
cnt = Counter(df.classification)
# select only 10 most common classes!
top_classes = 10
# sort classes
sorted_classes = cnt.most_common()[:top_classes]
classes = [c[0] for c in sorted_classes]
counts = [c[1] for c in sorted_classes]
print("at least " + str(counts[-1]) + " instances per class")

# apply to dataframe
print(str(df.shape[0]) + " instances before")
df = df[[c in classes for c in df.classification]]
print(str(df.shape[0]) + " instances after")

seqs = df.sequence.values
lengths = [len(s) for s in seqs]

# visualize
fig, axarr = plt.subplots(1,2, figsize=(20,5))
axarr[0].bar(range(len(classes)), counts)
plt.sca(axarr[0])
plt.xticks(range(len(classes)), classes, rotation='vertical')
axarr[0].set_ylabel('frequency')

axarr[1].hist(lengths, bins=100, normed=None)
axarr[1].set_xlabel('sequence length')
axarr[1].set_ylabel('# sequences')
plt.show()

ps. most_common([n])
Return a list of the n most common elements and their counts from the most common to the least. If n is omitted or None, most_common() returns all elements in the counter. Elements with equal counts are ordered arbitrarily:

>>> print(sorted_classes)
[('HYDROLASE', 46336), ('TRANSFERASE', 36424), ('OXIDOREDUCTASE', 34321), ('IMMUNE SYSTEM', 15615), ('LYASE', 11682), ('HYDROLASE/HYDROLASE INHIBITOR', 11218), ('TRANSCRIPTION', 8919), ('VIRAL PROTEIN', 8495), ('TRANSPORT PROTEIN', 8371), ('VIRUS', 6972)]

>>> classes = [c[0] for c in sorted_classes]
>>> print(classes)
['HYDROLASE', 'TRANSFERASE', 'OXIDOREDUCTASE', 'IMMUNE SYSTEM', 'LYASE', 'HYDROLASE/HYDROLASE INHIBITOR', 'TRANSCRIPTION', 'VIRAL PROTEIN', 'TRANSPORT PROTEIN', 'VIRUS']

>>> counts = [c[1] for c in sorted_classes]
>>> print(counts)
[46336, 36424, 34321, 15615, 11682, 11218, 8919, 8495, 8371, 6972]

參考
https://docs.python.org/2/library/collections.html

Transform labels
首先，數據集被簡化為十個最常見類別之一的樣本。 序列的長度範圍從極少數氨基酸到幾千個氨基酸（參見下圖）。

其次，使用sklearn.preprocessing中的LabelBinarizer將字符串中的標籤轉換為one hot表示。

from sklearn.preprocessing import LabelBinarizer

# Transform labels to one-hot
lb = LabelBinarizer()
Y = lb.fit_transform(df.classification)

Further preprocessing of sequences with keras
使用keras函式庫進行文本處理，

Tokenizer：將序列的每個字符轉換為數字
pad_sequences：確保每個序列具有相同的長度（max_length）。 我決定使用最大長度512，這對大多數序列來說應該足夠了。 注意：由於計算時間的原因，減少到256。
train_test_split：來自sklearn將數據分成訓練和測試樣本。

from keras.preprocessing import text, sequence
from keras.preprocessing.text import Tokenizer
from sklearn.model_selection import train_test_split

# maximum length of sequence, everything afterwards is discarded!
max_length = 256

#create and fit tokenizer
tokenizer = Tokenizer(char_level=True)
tokenizer.fit_on_texts(seqs)
#represent input data as word rank number sequences
X = tokenizer.texts_to_sequences(seqs)
X = sequence.pad_sequences(X, maxlen=max_length)

print(X.shape)
print(type(X))
print(type(X[0]))
print(type(X[0][0]))

(188353, 256)
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
<class 'numpy.int32'>

Build keras model and fit
以前的方法：

hnike的核心依賴於數據的數值特徵（分子量，殘留計數等）90％，僅考慮三種最常見的類別。 在沒有任何一般性聲明的情況下，使用此方案的一些快速測試顯示出優異的結果

abharg的核心使用來自sklearn的CountVectorizer，使用4-gram，然後使用簡單的分類模型，考慮到前43個類別（超過1000個樣本的類別），產生的準確率約為76％。 請注意，此方法也基於序列，並將與我的結果進行比較。

Further improvements
現在這個內核最重要的部分：

1. 最近在成功的NLP中建議使用word embedding操作在keras層(Embedding)。請注意，在我們的例子中只有20個不同的單詞(每個氨基酸)，可以嘗試將序列的one hot表示取代word embedding來進行2D卷積，但是這種方法側重於將NLP理論應用於蛋白質序列。

2. 考慮在嵌入序列上使用1D卷積，而不是使用每個n-gram。convolutional kernel的大小可以看作n-gram的大小，過濾器的數量可以看作單字的數量。

3. 為了提高性能，還可以使用深度結構（後續的捲積和池化層），這裡有兩層，其中第一層有64個濾波器，convolutional kernel大小為6，第二層有32個濾波器，大小為3。

4. 平化並傳遞激活函數到完全連接層，其中最後一層是softmax激活，並且大小對應於類別的數量。

使用標準adam優化器的最小化分類交叉熵(categorical crossentropy)來優化此結構。

from keras.models import Sequential
from keras.layers import Dense, Conv1D, MaxPooling1D, Flatten
from keras.layers import LSTM
from keras.layers.embeddings import Embedding

embedding_dim = 8

# create the model
model = Sequential()
model.add(Embedding(len(tokenizer.word_index)+1, embedding_dim, input_length=max_length))
model.add(Conv1D(filters=64, kernel_size=6, padding='same', activation='relu'))
model.add(MaxPooling1D(pool_size=2))
model.add(Conv1D(filters=32, kernel_size=3, padding='same', activation='relu'))
model.add(MaxPooling1D(pool_size=2))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dense(top_classes, activation='softmax'))
model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
print(model.summary())

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=.2)
model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=15, batch_size=128)

Evaluate model
評估僅基於準確性(accuracy)和混淆矩陣(confusion matrix)。

from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report
import itertools

train_pred = model.predict(X_train)
test_pred = model.predict(X_test)
print("train-acc = " + str(accuracy_score(np.argmax(y_train, axis=1), np.argmax(train_pred, axis=1))))
print("test-acc = " + str(accuracy_score(np.argmax(y_test, axis=1), np.argmax(test_pred, axis=1))))

# Compute confusion matrix
cm = confusion_matrix(np.argmax(y_test, axis=1), np.argmax(test_pred, axis=1))

# Plot normalized confusion matrix
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
np.set_printoptions(precision=2)
plt.figure(figsize=(10,10))
plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
plt.title('Confusion matrix')
plt.colorbar()
tick_marks = np.arange(len(lb.classes_))
plt.xticks(tick_marks, lb.classes_, rotation=90)
plt.yticks(tick_marks, lb.classes_)
#for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
#    plt.text(j, i, format(cm[i, j], '.2f'), horizontalalignment="center", color="white" if cm[i, j] > cm.max() / 2. else "black")
plt.ylabel('True label')
plt.xlabel('Predicted label')
plt.show()

print(classification_report(np.argmax(y_test, axis=1), np.argmax(test_pred, axis=1), target_names=lb.classes_))

該結果表明NLP理論適用於蛋白質序列分類。 可以進一步改進（更深的模型，更小的序列，使用n-gram作為word embedding的單字來加速學習），並且應該通過交叉驗證進行驗證。 另一個主要問題是關於overfitting小型訓練數據，但這可以通過諸如Dropout或正規化(regularizing)等著名技術通過增加大權重(kernel_regularizer)或大型激活(activation_regularizer)的懲罰來克服。



參考
https://www.kaggle.com/helmehelmuto/cnn-keras-and-innvestigate