“脫氧核糖核酸(DNA)是一種分子，其中包含每個物種獨特的生物學指令。DNA及其包含的說明在繁殖過程中從成年生物傳給其后代。“

簡介

基因組是生物體中DNA的完整集合。所有生物物種都有一個基因組，但是它們的差異很大。例如，人類基因組被排列成23條染色體，這有點像百科全書被編輯成23卷。如果算上所有字符(單個DNA“堿基對”)，每個人類基因組中將有超過60億個字符。所以這是一個巨大的工程。

人類基因組大約有60億個字符。如果您認為基因組(完整的DNA序列)就像一本書，那就是一本由大約60億個“ A”，“ C”，“ G”和“ T”字母組成的書。每個人都有獨特的基因組。盡管如此，科學家發(fā)現(xiàn)人類基因組的大部分彼此相似。

作為數(shù)據(jù)驅(qū)動的科學，基因組學廣泛地利用機器學習來捕獲數(shù)據(jù)中的關(guān)系并推斷出新的生物學假設。但是，要想擁有從不斷增長的基因組學數(shù)據(jù)中提取新發(fā)現(xiàn)的能力，就需要更強大的機器學習模型。通過有效利用大型數(shù)據(jù)集，深度學習已給了計算機視覺和自然語言處理等領(lǐng)域帶來了很大進步。它已成為許多基因組建模任務的首選方法，包括預測遺傳變異對基因調(diào)控機制(如DNA接受性和DNA剪接)的影響。

在本文中，我們將了解如何解釋DNA結(jié)構(gòu)以及如何使用機器學習算法來建立DNA序列數(shù)據(jù)的預測模型。

DNA序列如何表示？

該圖顯示了DNA雙螺旋結(jié)構(gòu)的一小部分。

DNA的雙螺旋結(jié)構(gòu)

雙螺旋是DNA的化學表示。但是DNA很特別。它是由四種類型的堿基組成的核苷酸：腺嘌呤(A)，胸腺嘧啶(T)，鳥嘌呤(G)和胞嘧啶(C)。我們總是稱它們?yōu)锳，C，G和T。

這四種化學物質(zhì)通過氫鍵以任何可能的順序連接在一起，形成一條鏈，這形成了DNA雙螺旋的一條線。雙螺旋的第二條線平衡了第一根。因此，如果第一條線上有A，則第二條線上對應位置必須為T。此外，C和G始終保持平衡。因此，一旦您確定了螺旋的一個螺紋，就可以隨時拼寫另一個螺紋。

單個DNA線(螺旋)的序列

這些堿基的順序決定了DNA鏈中包含哪些生物學指令。例如，序列ATCGTT可能指示藍眼睛，而ATCGCT可能指示棕色眼睛。

使用Python處理DNA序列數(shù)據(jù)

熟悉諸如Biopython和squiggle之類的Python包將在處理Python中的生物序列數(shù)據(jù)時為您提供幫助。

Biopython是python模塊的集合，這些模塊提供處理DNA，RNA和蛋白質(zhì)序列操作的功能，例如DNA字符串的反向互補，尋找蛋白質(zhì)序列中的基序列等。它提供了很多解析器，可以讀取所有主要的遺傳數(shù)據(jù)庫，例如 GenBank，SwissPort，FASTA等，

安裝Biopython

pip install biopython

**Squiggle：**這是一個工具，可以自動生成交互式網(wǎng)頁下的原始DNA序列的二維圖形表示。考慮到易用性，Squiggle實現(xiàn)了幾種序列可視化算法，并引入了“為方便人類使用而做”的新穎可視化方法。

安裝Squiggle

pip install Squiggle

DNA序列數(shù)據(jù)通常以“ fasta”格式的文件格式儲存。Fasta格式通過包含注釋的大于號和包含序列的另一行作為組成自己的單行數(shù)據(jù):

“AAGGTGAGTGAAATCTCAACACGAGTATGGTTCTGAGAGTAGCTCTGTAACTCTGAGG”*

文件可以包含一個或多個DNA序列。還有許多其他格式，但是fasta是最常見的格式。

這是使用Biopython處理Fasta格式的DNA序列的簡要示例。序列對象將包含諸如序列ID和sequence等屬性以及可以直接使用的序列長度。

我們將使用Biopython的Bio.SeqIO來解析DNA序列數(shù)據(jù)(fasta)。它提供了一個簡單的統(tǒng)一界面來輸入和輸出各種文件格式。

from Bio import SeqIOfor for sequence in SeqIO.parse('./drive/My Drive/example.fa', "fasta"): print(sequence.id) print(sequence.seq)print(len(sequence))

這樣就產(chǎn)生了序列ID，序列本身和序列長度。

ENST00000435737.5ATGTTTCGCATCACCAACATTGAGTTTCTTCCCGAATACCGACAAAAGGAGTCCAGGGAATTTCTTTCAGTGTCACGGACTGTGCAGCAAGTGATAAACCTGGTTTATACAACATCTGCCTTCTCCAAATTTTATGAGCAGTCTGTTGTTGCAGATGTCAGCAACAACAAAGGCGGCCTCCTTGTCCACTTTTGGATTGTTTTTGTCATGCCACGTGCCAAAGGCCACATCTTCTGTGAAGACTGTGTTGCCGCCATCTTGAAGGACTCCATCCAGACAAGCATCATAAACCGGACCTCTGTGGGGAGCTTGCAGGGACTGGCTGTGGACATGGACTCTGTGGTACTAAATGAAGTCCTGGGGCTGACTCTCATTGTCTGGATTGACTGA390

我們可以可視化這些DNA序列嗎？

是的，我們可以使用Squiggle python庫將這些DNA序列可視化，進行操作并運行：

Squiggle example.fa --method=gates

DNA序列被轉(zhuǎn)換為2D圖像，其中T，A，C和G分別在上，下，左和右方位。這給每個序列一個“形狀”。

現(xiàn)在，我們來可視化另一個包含6個DNA序列的fasta數(shù)據(jù)。

Squiggle example.fasta

在此，首先使用2位編碼方案將DNA序列轉(zhuǎn)換為二進制序列，該方案將T映射為00，C映射為01，A映射為10，G映射為11。

現(xiàn)在我們可以輕松加載和操作生物序列數(shù)據(jù)，那么怎么將數(shù)據(jù)用于機器學習或深度學習？

由于機器學習或深度學習模型要求輸入必須是特征矩陣或數(shù)字值，但目前我們?nèi)匀灰宰址蜃址袷酱鎯?shù)據(jù)。因此，下一步是將這些字符編碼為矩陣。

編碼序列數(shù)據(jù)有3種通用方法：

順序編碼DNA序列

獨熱(one-hot)編碼DNA序列

DNA序列作為獨立“語言”，稱為k-mer計數(shù)

讓我們分別進行實現(xiàn)，看看哪一個為我們提供了完美的輸入。

順序編碼DNA序列

在這種方法中，我們需要將每個堿基編碼為序數(shù)值。例如，“ ATGC”變?yōu)閇0.25、0.5、0.75、1.0]。任何其他字符(例如“ N”)都可以為0。

因此，讓我們創(chuàng)建一些函數(shù)，例如從序列字符串創(chuàng)建NumPy數(shù)組對象，以及帶有DNA序列字母“ a”，“ c”，“ g”和“ t”的標簽編碼器，以及其他任何字符比如“n”的編碼器。

import numpy as npimport redef string_to_array(seq_string): seq_string = seq_string.lower seq_string = re.sub('[^acgt]', 'n', seq_string) seq_string = np.array(list(seq_string)) return seq_string# create a label encoder with 'acgtn' alphabetfrom sklearn.preprocessing import LabelEncoderlabel_encoder = LabelEncoderlabel_encoder.fit(np.array(['a','c','g','t','z']))

這是一種將DNA序列字符串編碼為有序載體的功能。它返回一個NumPy數(shù)組，其中A = 0.25，C = 0.50，G = 0.75，T = 1.00，n = 0.00。

def ordinal_encoder(my_array): integer_encoded = label_encoder.transform(my_array) float_encoded = integer_encoded.astype(float)float_encoded[float_encoded == 0] = 0.25 # A float_encoded[float_encoded == 1] = 0.50 # C float_encoded[float_encoded == 2] = 0.75 # G float_encoded[float_encoded == 3] = 1.00 # T float_encoded[float_encoded == 4] = 0.00 # anything else, lets say n return float_encoded

讓我們嘗試一個簡單的序列：

seq_test = 'TTCAGCCAGTG'ordinal_encoder(string_to_array(seq_test))

獨熱編碼DNA序列

另一種方法是使用獨熱編碼來表示DNA序列。這在深度學習方法中得到了廣泛使用，非常適合卷積神經(jīng)網(wǎng)絡之類的算法。在此示例中，“ ATGC”將變?yōu)閇0,0,0,1]，[0,0,1,0]，[0,1,0,0]，[1,0,0,0]。這些編碼的矢量可以連接起來，也可以變成二維數(shù)組。

from sklearn.preprocessing import OneHotEncoderdef one_hot_encoder(seq_string):int_encoded = label_encoder.transform(seq_string) onehot_encoder = OneHotEncoder(sparse=False, dtype=int) int_encoded = int_encoded.reshape(len(int_encoded), 1) onehot_encoded = onehot_encoder.fit_transform(int_encoded) onehot_encoded = np.delete(onehot_encoded, -1, 1) return onehot_encoded

讓我們以一個簡單的序列來嘗試一下：

seq_test = 'GAATTCTCGAA'one_hot_encoder(string_to_array(seq_test))

DNA序列作為獨立“語言”，稱為k-mer計數(shù)

仍然存在的問題是，以上所有方法均不能產(chǎn)生長度一致的向量，這是將數(shù)據(jù)導入分類或回歸算法的必要條件。因此，使用上述方法，您必須輔助諸如截斷序列或用“ n”/“ 0”填充的方法，以獲取長度一致的向量。

DNA和蛋白質(zhì)序列可以看作是生命的語言。該語言對所有生命形式中存在的分子的指令和功能進行編碼。基因組與序列語言和書是相似的，子序列(基因和基因家族)是句子和章節(jié)，k-mers和肽是單詞，核苷酸堿基和氨基酸是字母。自然語言處理(NLP)也應采用和DNA及蛋白質(zhì)序列相似的處理方式是有理由的。

我們在這里使用的方法是易于管理的。我們首先采用較長的生物學序列，并將其分解為k-mer長度重疊的“單詞”。例如，如果我們使用長度為6(六進制)的“單詞”，則“ ATGCATGCA”將變?yōu)?#xff1a;“ ATGCAT”，“ TGCATG”，“ GCATGC”，“ CATGCA”。因此，我們的示例序列分為4個六聚體字(hexamer words)。

在基因組學中，我們將這種類型的操作稱為“ k-mer計數(shù)”，或者對每種可能出現(xiàn)的k-mer序列進行計數(shù)，而Python的自然語言處理工具使其變得非常容易。

def Kmers_funct(seq, size): return [seq[x:x+size].lower for x in range(len(seq) - size + 1)]

因此，讓我們以一個簡單的序列嘗試一下：

mySeq = 'GTGCCCAGGTTCAGTGAGTGACACAGGCAG'Kmers_funct(mySeq, size=7)

它返回k-mer“單詞”的列表。然后，您可以將“單詞”加入到“句子”中，然后像往常一樣在“句子”上應用自己喜歡的自然語言處理方法。

words = Kmers_funct(mySeq, size=6)joined_sentence = ' '.join(words)joined_sentence

連接的句子：

'gtgccc tgccca gcccag cccagg ccaggt caggtt aggttc ggttca gttcag ttcagt tcagtg cagtga agtgag gtgagt tgagtg gagtga agtgac gtgaca tgacac gacaca acacag cacagg acaggc caggca aggcag'

您可以調(diào)整單詞長度和重疊量。這使您可以確定DNA序列信息和詞匯量在您的應用程序中的重要程度。例如，如果您使用長度為6的單詞，并且有4個字母，則詞匯量為4096個可能的單詞。然后，您可以像在NLP中一樣繼續(xù)創(chuàng)建單詞庫(bag-of-words)模型。

讓我們增加一些更有趣的“句子”。

mySeq1 = 'TCTCACACATGTGCCAATCACTGTCACCC'mySeq2 = 'GTGCCCAGGTTCAGTGAGTGACACAGGCAG'sentence1 = ' '.join(Kmers_funct(mySeq1, size=6))sentence2 = ' '.join(Kmers_funct(mySeq2, size=6))

建立單詞袋模型：

from sklearn.feature_extraction.text import CountVectorizercv = CountVectorizerX = cv.fit_transform([joined_sentence, sentence1, sentence2]).toarrayX

單詞袋

下面開始進行機器學習

現(xiàn)在，我們已經(jīng)學習了如何從DNA序列中提取特征矩陣，讓我們將我們新獲得的知識應用于機器學習用例。

用例：建立一個在人類DNA序列上受訓的分類模型，并可以根據(jù)編碼序列的DNA序列預測基因家族。為了測試該模型，我們將使用人，狗和黑猩猩的DNA序列進行訓練，并測試其準確性。

基因家族是一組具有共同祖先的相關(guān)基因。基因家族的成員可以是旁系同源物或直系同源物。基因旁系同源物是來自相同物種的具有相似序列的基因，而基因直系同源物是在不同物種中具有相似序列的基因。

數(shù)據(jù)集包含人類的DNA序列，狗的DNA序列和黑猩猩的DNA序列。

加載人類DNA序列。

import numpy as npimport pandas as pdimport matplotlib.pyplot as plt%matplotlib inlinehuman_dna = pd.read_table('./drive/My Drive/human_data.txt')human_dna.head

人類DNA序列和類別標簽

加載黑猩猩和狗的DNA序列。

chimp_dna = pd.read_table('./drive/My Drive/chimp_data.txt')dog_dna = pd.read_table('./drive/My Drive/dog_data.txt')

狗DNA序列和類別標簽

這是7個類別中每個類別的定義，以及人工訓練數(shù)據(jù)中的類別數(shù)：

人類DNA數(shù)據(jù)集中存在帶有類別標簽的基因家族

現(xiàn)在我們已經(jīng)加載了所有數(shù)據(jù)，下一步是將字符序列轉(zhuǎn)換為k-mer詞，默認大小為6(六進制)。函數(shù)Kmers_funct將從序列字符串中收集指定長度的所有可能重疊的k-mers。

def Kmers_funct(seq, size=6): return [seq[x:x+size].lower for x in range(len(seq) - size + 1)]#convert our training data sequences into short overlapping k-mers of length 6. Lets do that for each species of data we have using our Kmers_funct function.human_dna['words'] = human_dna.apply(lambda x: Kmers_funct(x['sequence']), axis=1)human_dna = human_dna.drop('sequence', axis=1)chimp_dna['words'] = chimp_dna.apply(lambda x: Kmers_funct(x['sequence']), axis=1)chimp_dna = chimp_dna.drop('sequence', axis=1)dog_dna['words'] = dog_dna.apply(lambda x: Kmers_funct(x['sequence']), axis=1)dog_dna = dog_dna.drop('sequence', axis=1)

將DNA序列更改為小寫，分為所有可能的長度為6的k-mer字，并準備下一步。

human_dna.head

人類DNA序列中長度為6的k-mer字

現(xiàn)在，我們需要將每個基因的k-mers列表轉(zhuǎn)換為可用于創(chuàng)建單詞袋模型的字符串句子。我們將創(chuàng)建一個目標變量y來保存類標簽。

對黑猩猩和狗也進行一樣的操作。

human_texts = list(human_dna['words'])for item in range(len(human_texts)):human_texts[item] = ' '.join(human_texts[item])#separate labelsy_human = human_dna.iloc[:, 0].values # y_human for human_dna#Now let's do the same for chimp and dog.chimp_texts = list(chimp_dna['words'])for item in range(len(chimp_texts)): chimp_texts[item] = ' '.join(chimp_texts[item])#separate labelsy_chim = chimp_dna.iloc[:, 0].values # y_chim for chimp_dnadog_texts = list(dog_dna['words'])for item in range(len(dog_texts)): dog_texts[item] = ' '.join(dog_texts[item])#separate labelsy_dog = dog_dna.iloc[:, 0].values # y_dog for dog_dna

因此，目標變量包含一個保存類別的數(shù)組。

array([4, 4, 3, …, 6, 6, 6])

將我們的k-mer單詞轉(zhuǎn)換為均等長度的數(shù)字矢量，這些矢量代表詞匯中每個k-mer的計數(shù)：

from sklearn.feature_extraction.text import CountVectorizercv = CountVectorizer(ngram_range=(4,4)) #The n-gram size of 4 is previously determined by testingX = cv.fit_transform(human_texts)X_chimp = cv.transform(chimp_texts)X_dog = cv.transform(dog_texts)

您可能要檢查每個訓練集數(shù)據(jù)的形狀。

print(X.shape)print(X_chimp.shape)print(X_dog.shape)#shapes(4380, 232414) (1682, 232414) (820, 232414)

因此，對于人類，我們已經(jīng)將4380個基因轉(zhuǎn)換為4-gram的k-mer(長度6)計數(shù)的均勻長度特征向量。對于黑猩猩和狗，我們分別具有1682和820個基因的相同形狀的特征。

既然我們知道如何將我們的DNA序列轉(zhuǎn)換為k-mer計數(shù)和n-gram形式的均勻長度的數(shù)字矢量，那么我們現(xiàn)在就可以繼續(xù)構(gòu)建一個分類模型，該模型可以僅基于序列本身來預測DNA序列功能 。

在這里，我將使用人類數(shù)據(jù)來訓練模型，并拿出20％的人類數(shù)據(jù)來測試模型。然后，通過嘗試預測其他物種(黑猩猩和狗)的序列功能，我們可以測試模型的可推廣性。

接下來，拆分用來訓練/測試的人類數(shù)據(jù)集并構(gòu)建簡單的多項樸素貝葉斯分類器。

您可能需要進行一些參數(shù)調(diào)整，并構(gòu)建具有不同n-gram大小的模型，在這里，我將繼續(xù)使用n-gram大小為4和alpha為0.1的模型。

# Splitting the human dataset into the training set and test setfrom sklearn.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(X, y_human, test_size = 0.20, random_state=42)## Multinomial Naive Bayes Classifier ##from sklearn.naive_bayes import MultinomialNBclassifier = MultinomialNB(alpha=0.1)classifier.fit(X_train, y_train)

現(xiàn)在，讓我們對測試集進行預測，看看它的性能如何。

y_pred = classifier.predict(X_test)

下面讓我們來看一些模型性能指標，例如混淆矩陣，準確性，召回率和f1得分。我們在數(shù)據(jù)上獲得了非常好的結(jié)果，因此看來我們的模型并未過擬合訓練數(shù)據(jù)。

from sklearn.metrics import accuracy_score, f1_score, precision_score, recall_scoreprint("Confusion matrix for predictions on human test DNA sequence")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, f1accuracy, precision, recall, f1 = get_metrics(y_test, y_pred)print("accuracy = %.3f precision = %.3f recall = %.3f f1 = %.3f" % (accuracy, precision, recall, f1))

讓我們看看我們的模型如何處理其他物種的DNA序列。首先，我們將嘗試黑猩猩，我們希望它與人類非常相似。然后是狗DNA序列。

# Predicting the chimp, dog and worm sequencesy_pred_chimp = classifier.predict(X_chimp)

檢查精度矩陣：

# performance on chimpanzee genesprint("Confusion matrix for predictions on Chimpanzee test DNA sequence")print(pd.crosstab(pd.Series(y_chim, name='Actual'), pd.Series(y_pred_chimp, name='Predicted')))accuracy, precision, recall, f1 = get_metrics(y_chim, y_pred_chimp)print("accuracy = %.3f precision = %.3f recall = %.3f f1 = %.3f" % (accuracy, precision, recall, f1))

現(xiàn)在讓我們對狗的測試DNA序列進行預測：

y_pred_dog = classifier.predict(X_dog)

檢查精度矩陣：

# performance on dog genesprint("Confusion matrix for predictions on Dog test DNA sequence")print(pd.crosstab(pd.Series(y_dog, name='Actual'), pd.Series(y_pred_dog, name='Predicted')))accuracy, precision, recall, f1 = get_metrics(y_dog, y_pred_dog)print("accuracy = %.3f precision = %.3f recall = %.3f f1 = %.3f" % (accuracy, precision, recall, f1))

該模型似乎在人類數(shù)據(jù)上產(chǎn)生了良好的結(jié)果。黑猩猩也是如此，這是因為黑猩猩和人類具有相同的遺傳層次。狗的表現(xiàn)不太好，這是因為狗比黑猩猩與人類的差異更大。

結(jié)論

在本文中，我們學習了如何分析DNA序列數(shù)據(jù)，如何對其進行可視化，以及如何使用不同的編碼技術(shù)將這些序列表示為矩陣。最后，我們創(chuàng)建了一個Naive Byes模型，可以在人，狗和黑猩猩的測試數(shù)據(jù)中檢測基因家族。

總結(jié)

以上是生活随笔為你收集整理的numpy序列预处理dna序列_使用机器学习和Python揭开DNA测序神秘面纱的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。