網頁

2019年7月8日 星期一

Working with DNA sequence data for machine learning part 2 - clasification of gene function

接下來將應用scikit learn的ordinal vectors, one-hot encoded vectors or as counts of k-mers來建立一個分類模型，該模型可以根據編碼序列的DNA序列單獨預測基因的功能。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

開啟數據

human = pd.read_table('../input/humandata/human_data.txt')
human.head()


除了有一些人類DNA序列編碼區和類別標籤的數據以外，還有黑猩猩的數據和更加不同的物種，狗。

chimp = pd.read_table('../input/chimp-data/chimp_data.txt')
dog = pd.read_table('../input/dog-data/dog_data.txt')
chimp.head()
dog.head()


以下是7個類別中每個類的定義以及人類訓練數據中的數量。 它們是基因序列功能組。

from IPython.display import Image
Image("../input/image2/Capture1.PNG")


讓我們定義一個函數來從任何序列字符串中收集指定長度的所有可能的重疊k-mers。

# function to convert sequence strings into k-mer words, default size = 6 (hexamer words)
def getKmers(sequence, size=6):
    return [sequence[x:x+size].lower() for x in range(len(sequence) - size + 1)]

現在我們可以將訓練數據序列轉換為長度為6的短重疊k-mers。讓我們使用getKmers函數對每種數據進行處理。

human['words'] = human.apply(lambda x: getKmers(x['sequence']), axis=1)
human = human.drop('sequence', axis=1)
chimp['words'] = chimp.apply(lambda x: getKmers(x['sequence']), axis=1)
chimp = chimp.drop('sequence', axis=1)
dog['words'] = dog.apply(lambda x: getKmers(x['sequence']), axis=1)
dog = dog.drop('sequence', axis=1)

現在，我們的編碼序列數據被改為小寫，分成長度為6的所有可能的k-mer字，並為下一步做好準備。

human.head()


由於我們將使用scikit-learn自然語言處理工具進行k-mer計數，我們現在需要將每個基因的k-mers list轉換為計數向量化器可以使用的字符串句子。 我們還可以創建一個y變量來保存類別標籤。

human_texts = list(human['words'])
for item in range(len(human_texts)):
    human_texts[item] = ' '.join(human_texts[item])
y_h = human.iloc[:, 0].values                         #y_h for human

human_texts[0]


原先human_texts = list(human_data['words'])的type跟shape為list巢狀結構
print(type(human_texts))
print(type(human_texts[0]))
print(type(human_texts[0][0]))
print(np.array(human_texts).shape)
print(np.array(human_texts[0]).shape)
print(np.array(human_texts[0][0]).shape)

<class 'list'>
<class 'list'>
<class 'str'>
(4380,)
(202,)
()

執行完
for item in range(len(human_texts)):
    human_texts[item] = ' '.join(human_texts[item])

human_texts第二層list的202個list會變成一個字串

print(type(human_texts))
print(type(human_texts[0]))
print(type(human_texts[0][0]))
print(np.array(human_texts).shape)
print(len(human_texts[0]))
print(len(human_texts[0][0]))

<class 'list'>
<class 'str'>
<class 'str'>
(4380,)
1413
1

y_h


現在讓我們為黑猩猩和狗做同樣的事情。

chimp_texts = list(chimp['words'])
for item in range(len(chimp_texts)):
    chimp_texts[item] = ' '.join(chimp_texts[item])
y_c = chimp.iloc[:, 0].values                       # y_c for chimp

dog_texts = list(dog['words'])
for item in range(len(dog_texts)):
    dog_texts[item] = ' '.join(dog_texts[item])
y_d = dog.iloc[:, 0].values                         # y_d for dog

現在讓我們回顧一下如何使用sklearn的“自然語言”處理工具將我們的k-mer單詞轉換為統一長度的數字向量，這些向量表示詞彙表中每個k-mer的計數。

# Creating the Bag of Words model using CountVectorizer()
# This is equivalent to k-mer counting
# The n-gram size of 4 was previously determined by testing
from sklearn.feature_extraction.text import CountVectorizer
cv = CountVectorizer(ngram_range=(4,4))
X = cv.fit_transform(human_texts)
X_chimp = cv.transform(chimp_texts)
X_dog = cv.transform(dog_texts)

ps. 要提取的不同n-gram的n-values範圍的下邊界和上邊界。 n的所有值都被使用在min_n <= n <= max_n。以這個範例來說，human_texts[0]有1413個長度為6的短重疊k-mers構成，當CountVectorizer的ngram_range設定為(4,4)表示1413個k-mers會取4個k-mers形成一個4-gram序列然後去計算有重複的個數。例如list human_texts[0]本來有[atgtgt tgtgtg gtgtgg tgtggc gtggca tggcat ggcatt gcattt....]1413個k-mers，執行CountVectorizer(ngram_range=(4,4))後，會產生4-gram序列['atgtgt tgtgtg gtgtgg tgtggc', 'atgtgt tgtgtg gtgtgg gtggca',.....]。

參考
https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.text.CountVectorizer.html

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

(4380, 232414)
[[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]]
<class 'scipy.sparse.csr.csr_matrix'>
<class 'scipy.sparse.csr.csr_matrix'>
<class 'scipy.sparse.csr.csr_matrix'>

讓我們看看我們擁有什麼......對於人類，我們將4380個基因轉換成4-gram k-mer（長度6）計數的均勻長度特徵向量。 對於黑猩猩和狗，我們分別具有1682和820個基因的預期相同數量的特徵。

print(X.shape)
print(X_chimp.shape)
print(X_dog.shape)


如果我們看看類別平衡，我們可以看到我們有相對平衡的數據集

human['class'].value_counts().sort_index().plot.bar()


chimp['class'].value_counts().sort_index().plot.bar()


dog['class'].value_counts().sort_index().plot.bar()


現在，我們知道如何將我們的DNA序列轉換為k-mer計數和ngrams形式的均勻長度數值向量，我們現在可以繼續構建一個分類模型，該模型僅基於序列本身來預測DNA序列功能。

這裡將使用人類數據來訓練模型，拿出20％的人體數據來測試模型。 然後我們可以通過嘗試預測其他物種（黑猩猩和狗）的序列功能來挑戰模型的通則性(generalizability)。

所以下面我們將 1：訓練/測試分化。 2：構建簡單的多項式樸素貝葉斯分類器(multinomial naive Bayes classifier )和3：測試模型性能。

# Splitting the human dataset into the training set and test set
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, 
                                                    y_h, 
                                                    test_size = 0.20, 
                                                    random_state=42)

print(X_train.shape)
print(X_test.shape)


將創建一個多項式樸素貝葉斯分類器。 之前做了一些參數調整發現ngram大小為4(反映在Countvectorizer()實例中)，模型alpha為0.1做得最好。

### Multinomial Naive Bayes Classifier ###
# The alpha parameter was determined by grid search previously
from sklearn.naive_bayes import MultinomialNB
classifier = MultinomialNB(alpha=0.1)
classifier.fit(X_train, y_train)


現在讓我們對人類測試集進行預測，看看它如何在未看見數據上進行預測。

y_pred = classifier.predict(X_test)

讓我們來看看一些模型性能指標，如混淆矩陣(confusion matrix)，準確度(accuracy)，精確度(precision)，召回率(at)和f1 score。 我們在未看見數據上獲得了非常好的結果，所以看起來我們的模型沒有overfit訓練數據。 在一個真實的項目中，由於我們的數據集相對較小，我會回去再取樣更多的訓練測試分化。

from sklearn.metrics import accuracy_score, f1_score, precision_score, recall_score
print("Confusion matrix\n")
print(pd.crosstab(pd.Series(y_test, name='Actual'), pd.Series(y_pred, name='Predicted')))
def get_metrics(y_test, y_predicted):
    accuracy = accuracy_score(y_test, y_predicted)
    precision = precision_score(y_test, y_predicted, average='weighted')
    recall = recall_score(y_test, y_predicted, average='weighted')
    f1 = f1_score(y_test, y_predicted, average='weighted')
    return accuracy, precision, recall, f1
accuracy, precision, recall, f1 = get_metrics(y_test, y_pred)
print("accuracy = %.3f \nprecision = %.3f \nrecall = %.3f \nf1 = %.3f" % (accuracy, precision, recall, f1))


現在進行真正的測試。 讓我們看看我們的模型對來自其他物種的DNA序列的性能如何。 首先，我們將嘗試黑猩猩，我們希望它與人類非常相似。 然後將嘗試狗DNA序列。

預測黑猩猩和狗的序列

# Predicting the chimp, dog and worm sequences
y_pred_chimp = classifier.predict(X_chimp)
y_pred_dog = classifier.predict(X_dog)

現在，讓我們來檢查每個物種的典型性能指標。

# performance on chimp genes
print("Confusion matrix\n")
print(pd.crosstab(pd.Series(y_c, name='Actual'), pd.Series(y_pred_chimp, name='Predicted')))
accuracy, precision, recall, f1 = get_metrics(y_c, y_pred_chimp)
print("accuracy = %.3f \nprecision = %.3f \nrecall = %.3f \nf1 = %.3f" % (accuracy, precision, recall, f1))


# performance on dog genes
print("Confusion matrix\n")
print(pd.crosstab(pd.Series(y_d, name='Actual'), pd.Series(y_pred_dog, name='Predicted')))
accuracy, precision, recall, f1 = get_metrics(y_d, y_pred_dog)
print("accuracy = %.3f \nprecision = %.3f \nrecall = %.3f \nf1 = %.3f" % (accuracy, precision, recall, f1))


該模型似乎在人類數據上表現良好。 它也適用於黑猩猩。 這可能不是一個驚喜，因為黑猩猩和人類在遺傳上非常相似。 狗的表現不太好。 我們會期待這一點，因為狗比黑猩猩更偏離人類。



參考
https://www.kaggle.com/thomasnelson/working-with-dna-sequence-data-for-ml-part-2
張貼者： 於
標籤：

沒有留言:

張貼留言

訂閱： 張貼留言 (Atom)