「脫氧核糖核酸（DNA）是一種分子，其中包含每個物種獨特的生物學指令。DNA及其包含的說明在繁殖過程中從成年生物傳給其後代。「
簡介
基因組是生物體中DNA的完整集合。所有生物物種都有一個基因組，但是它們的差異很大。例如，人類基因組被排列成23條染色體，這有點像百科全書被編輯成23卷。如果算上所有字符（單個DNA「鹼基對」），每個人類基因組中將有超過60億個字符。所以這是一個巨大的工程。
人類基因組大約有60億個字符。如果您認為基因組（完整的DNA序列）就像一本書，那就是一本由大約60億個「 A」，「 C」，「 G」和「 T」字母組成的書。每個人都有獨特的基因組。儘管如此，科學家發現人類基因組的大部分彼此相似。
作為數據驅動的科學，基因組學廣泛地利用機器學習來捕獲數據中的關係並推斷出新的生物學假設。但是，要想擁有從不斷增長的基因組學數據中提取新發現的能力，就需要更強大的機器學習模型。通過有效利用大型數據集，深度學習已給了計算機視覺和自然語言處理等領域帶來了很大進步。它已成為許多基因組建模任務的首選方法，包括預測遺傳變異對基因調控機制（如DNA接受性和DNA剪接）的影響。
在本文中，我們將了解如何解釋DNA結構以及如何使用機器學習算法來建立DNA序列數據的預測模型。
DNA序列如何表示？
該圖顯示了DNA雙螺旋結構的一小部分。
DNA的雙螺旋結構
雙螺旋是DNA的化學表示。但是DNA很特別。它是由四種類型的鹼基組成的核苷酸：腺嘌呤（A），胸腺嘧啶（T），鳥嘌呤（G）和胞嘧啶（C）。我們總是稱它們為A，C，G和T。
這四種化學物質通過氫鍵以任何可能的順序連接在一起，形成一條鏈，這形成了DNA雙螺旋的一條線。雙螺旋的第二條線平衡了第一根。因此，如果第一條線上有A，則第二條線上對應位置必須為T。此外，C和G始終保持平衡。因此，一旦您確定了螺旋的一個螺紋，就可以隨時拼寫另一個螺紋。
單個DNA線（螺旋）的序列
這些鹼基的順序決定了DNA鏈中包含哪些生物學指令。例如，序列ATCGTT可能指示藍眼睛，而ATCGCT可能指示棕色眼睛。
使用Python處理DNA序列數據
熟悉諸如Biopython和squiggle之類的Python包將在處理Python中的生物序列數據時為您提供幫助。
Biopython是python模塊的集合，這些模塊提供處理DNA，RNA和蛋白質序列操作的功能，例如DNA字符串的反向互補，尋找蛋白質序列中的基序列等。它提供了很多解析器，可以讀取所有主要的遺傳資料庫，例如 GenBank，SwissPort，FASTA等，
安裝Biopython
pip install biopython
**Squiggle：**這是一個工具，可以自動生成交互式網頁下的原始DNA序列的二維圖形表示。考慮到易用性，Squiggle實現了幾種序列可視化算法，並引入了「為方便人類使用而做」的新穎可視化方法。
安裝Squiggle
pip install Squiggle
DNA序列數據通常以「 fasta」格式的文件格式儲存。Fasta格式通過包含注釋的大於號和包含序列的另一行作為組成自己的單行數據:
「AAGGTGAGTGAAATCTCAACACGAGTATGGTTCTGAGAGTAGCTCTGTAACTCTGAGG」*
文件可以包含一個或多個DNA序列。還有許多其他格式，但是fasta是最常見的格式。
這是使用Biopython處理Fasta格式的DNA序列的簡要示例。序列對象將包含諸如序列ID和sequence等屬性以及可以直接使用的序列長度。
我們將使用Biopython的Bio.SeqIO來解析DNA序列數據（fasta）。它提供了一個簡單的統一界面來輸入和輸出各種文件格式。
from Bio import SeqIOfor for sequence in SeqIO.parse('./drive/My Drive/example.fa', "fasta"): print(sequence.id) print(sequence.seq)print(len(sequence))
這樣就產生了序列ID，序列本身和序列長度。
ENST00000435737.5ATGTTTCGCATCACCAACATTGAGTTTCTTCCCGAATACCGACAAAAGGAGTCCAGGGAATTTCTTTCAGTGTCACGGACTGTGCAGCAAGTGATAAACCTGGTTTATACAACATCTGCCTTCTCCAAATTTTATGAGCAGTCTGTTGTTGCAGATGTCAGCAACAACAAAGGCGGCCTCCTTGTCCACTTTTGGATTGTTTTTGTCATGCCACGTGCCAAAGGCCACATCTTCTGTGAAGACTGTGTTGCCGCCATCTTGAAGGACTCCATCCAGACAAGCATCATAAACCGGACCTCTGTGGGGAGCTTGCAGGGACTGGCTGTGGACATGGACTCTGTGGTACTAAATGAAGTCCTGGGGCTGACTCTCATTGTCTGGATTGACTGA390
我們可以可視化這些DNA序列嗎？
是的，我們可以使用Squiggle python庫將這些DNA序列可視化，進行操作並運行：
Squiggle example.fa --method=gates
DNA序列被轉換為2D圖像，其中T，A，C和G分別在上，下，左和右方位。這給每個序列一個「形狀」。
現在，我們來可視化另一個包含6個DNA序列的fasta數據。
Squiggle example.fasta
在此，首先使用2位編碼方案將DNA序列轉換為二進位序列，該方案將T映射為00，C映射為01，A映射為10，G映射為11。
現在我們可以輕鬆加載和操作生物序列數據，那麼怎麼將數據用於機器學習或深度學習？
由於機器學習或深度學習模型要求輸入必須是特徵矩陣或數字值，但目前我們仍然以字符或字符串格式存儲數據。因此，下一步是將這些字符編碼為矩陣。
編碼序列數據有3種通用方法：
順序編碼DNA序列獨熱（one-hot）編碼DNA序列DNA序列作為獨立「語言」，稱為k-mer計數讓我們分別進行實現，看看哪一個為我們提供了完美的輸入。
順序編碼DNA序列
在這種方法中，我們需要將每個鹼基編碼為序數值。例如，「 ATGC」變為[0.25、0.5、0.75、1.0]。任何其他字符（例如「 N」）都可以為0。
因此，讓我們創建一些函數，例如從序列字符串創建NumPy數組對象，以及帶有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數組，其中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序列。這在深度學習方法中得到了廣泛使用，非常適合卷積神經網絡之類的算法。在此示例中，「 ATGC」將變為[0,0,0,1]，[0,0,1,0]，[0,1,0,0]，[1,0,0,0]。這些編碼的矢量可以連接起來，也可以變成二維數組。
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計數
仍然存在的問題是，以上所有方法均不能產生長度一致的向量，這是將數據導入分類或回歸算法的必要條件。因此，使用上述方法，您必須輔助諸如截斷序列或用「 n」/「 0」填充的方法，以獲取長度一致的向量。
DNA和蛋白質序列可以看作是生命的語言。該語言對所有生命形式中存在的分子的指令和功能進行編碼。基因組與序列語言和書是相似的，子序列（基因和基因家族）是句子和章節，k-mers和肽是單詞，核苷酸鹼基和胺基酸是字母。自然語言處理（NLP）也應採用和DNA及蛋白質序列相似的處理方式是有理由的。
我們在這裡使用的方法是易於管理的。我們首先採用較長的生物學序列，並將其分解為k-mer長度重疊的「單詞」。例如，如果我們使用長度為6（六進位）的「單詞」，則「 ATGCATGCA」將變為：「 ATGCAT」，「 TGCATG」，「 GCATGC」，「 CATGCA」。因此，我們的示例序列分為4個六聚體字（hexamer words）。
在基因組學中，我們將這種類型的操作稱為「 k-mer計數」，或者對每種可能出現的k-mer序列進行計數，而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'
您可以調整單詞長度和重疊量。這使您可以確定DNA序列信息和詞彙量在您的應用程式中的重要程度。例如，如果您使用長度為6的單詞，並且有4個字母，則詞彙量為4096個可能的單詞。然後，您可以像在NLP中一樣繼續創建單詞庫（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
單詞袋
下面開始進行機器學習
現在，我們已經學習了如何從DNA序列中提取特徵矩陣，讓我們將我們新獲得的知識應用於機器學習用例。
用例：建立一個在人類DNA序列上受訓的分類模型，並可以根據編碼序列的DNA序列預測基因家族。為了測試該模型，我們將使用人，狗和黑猩猩的DNA序列進行訓練，並測試其準確性。
基因家族是一組具有共同祖先的相關基因。基因家族的成員可以是旁系同源物或直系同源物。基因旁系同源物是來自相同物種的具有相似序列的基因，而基因直系同源物是在不同物種中具有相似序列的基因。
數據集包含人類的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個類別中每個類別的定義，以及人工訓練數據中的類別數：
人類DNA數據集中存在帶有類別標籤的基因家族
現在我們已經加載了所有數據，下一步是將字符序列轉換為k-mer詞，默認大小為6（六進位）。函數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字
現在，我們需要將每個基因的k-mers列錶轉換為可用於創建單詞袋模型的字符串句子。我們將創建一個目標變量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
因此，目標變量包含一個保存類別的數組。
array([4, 4, 3, …, 6, 6, 6])
將我們的k-mer單詞轉換為均等長度的數字矢量，這些矢量代表詞彙中每個k-mer的計數：
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)
您可能要檢查每個訓練集數據的形狀。
print(X.shape)print(X_chimp.shape)print(X_dog.shape)#shapes(4380, 232414) (1682, 232414) (820, 232414)
因此，對於人類，我們已經將4380個基因轉換為4-gram的k-mer（長度6）計數的均勻長度特徵向量。對於黑猩猩和狗，我們分別具有1682和820個基因的相同形狀的特徵。
既然我們知道如何將我們的DNA序列轉換為k-mer計數和n-gram形式的均勻長度的數字矢量，那麼我們現在就可以繼續構建一個分類模型，該模型可以僅基於序列本身來預測DNA序列功能 。
在這裡，我將使用人類數據來訓練模型，並拿出20％的人類數據來測試模型。然後，通過嘗試預測其他物種（黑猩猩和狗）的序列功能，我們可以測試模型的可推廣性。
接下來，拆分用來訓練/測試的人類數據集並構建簡單的多項樸素貝葉斯分類器。
您可能需要進行一些參數調整，並構建具有不同n-gram大小的模型，在這裡，我將繼續使用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)
現在，讓我們對測試集進行預測，看看它的性能如何。
y_pred = classifier.predict(X_test)
下面讓我們來看一些模型性能指標，例如混淆矩陣，準確性，召回率和f1得分。我們在數據上獲得了非常好的結果，因此看來我們的模型並未過擬合訓練數據。
from sklearn.metrics import accuracy_score, f1_score, precision_score, recall_scoreprint("Confusion matrix for predictions on human test DNA sequence\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, f1accuracy, 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 sequencesy_pred_chimp = classifier.predict(X_chimp)
檢查精度矩陣：
# performance on chimpanzee genesprint("Confusion matrix for predictions on Chimpanzee test DNA sequence\n")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 \nprecision = %.3f \nrecall = %.3f \nf1 = %.3f" % (accuracy, precision, recall, f1))
現在讓我們對狗的測試DNA序列進行預測：
y_pred_dog = classifier.predict(X_dog)
檢查精度矩陣：
# performance on dog genesprint("Confusion matrix for predictions on Dog test DNA sequence\n")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 \nprecision = %.3f \nrecall = %.3f \nf1 = %.3f" % (accuracy, precision, recall, f1))
該模型似乎在人類數據上產生了良好的結果。黑猩猩也是如此，這是因為黑猩猩和人類具有相同的遺傳層次。狗的表現不太好，這是因為狗比黑猩猩與人類的差異更大。
結論
在本文中，我們學習了如何分析DNA序列數據，如何對其進行可視化，以及如何使用不同的編碼技術將這些序列表示為矩陣。最後，我們創建了一個Naive Byes模型，可以在人，狗和黑猩猩的測試數據中檢測基因家族。