我演示了如何将深度学习用于古代DNA，单细胞生物学，OMICs数据集成，临床诊断和显微镜成像。在之前的文章《尼安德特人基因的深度学习》中，我强调了深度学习和自然语言处理对于古代基因组学的巨大潜力， 并演示了如何实际开始使用它来推断现代人类基因组中尼安德特人的基因渗入区域。现在，我们将充分利用长记忆循环神经网络（RNN）的强大功能 为了更好地预测现代人类基因组中尼安德特人的祖先的片段。

针对古代基因组学的HMM与LSTM

DNA是一个具有长记忆的序列，通过沿着该序列的长距离相关性表现出来，称为连锁不平衡。但是，在古代基因组学中，许多分析都是使用隐马尔可夫模型（HMM）进行的，该模型是一种无记忆模型，只能考虑一些先前的步骤。在长短期记忆（LSTM）人工神经网络已被证明，以超越的HMM利用其独特的记忆能力在许多应用，如语音识别，文本生成等，很多以前的步骤序列中 提供足够的数据。

来自Panzner和Cimiano的《机器学习，优化和大数据》，2016年

基因组学和古基因组学代表了一个真正的大数据资源，它提供了数以百万计的序列，可用于训练LSTMs的示例/统计观察。因此，古老的基因组数据是深度学习的礼物！

尽管基因组学/古基因组学的训练资源非常丰富，但大多数分析仍然是在模拟基因组数据上进行的，合并理论可能是模拟人口统计学和种群遗传学选择的最流行的框架。

LSTM + Word嵌入尼安德特人的DNA

在这里，我将在之前的文章《关于尼安德特人基因的深度学习》中，进一步阐述利用NLP进行古代基因组学研究的思路，并论证LSTMs在测序数据分析方面的优越性。您可能已经注意到，上一篇文章中的“深度学习”不是特别“深度”而是“广泛”，此外，随机森林似乎在“语言包”模型中展示了更好的性能。在这里，我将实现一个LSTM模型，该模型可以检测到从尼安德特人那里继承来的DNA片段，其准确性达到99％。在上一篇文章中，我展示了如何提取尼安德特人渐渗和枯竭的序列，所以在这里，我将从读取序列开始，将它们分成200个核苷酸长的子序列，每个子序列代表一个句子，因此我们可以将每个句子进一步拆分为k-mers /句子中的单词。

from Bio import SeqIOfrom Bio.Seq import Seqintr_seqs = []; depl_seqs = []; e = 0intr_f = 'hg19_intr_clean.fa'; depl_f = 'hg19_depl_clean.fa'for intr, depl in zip(SeqIO.parse(intr_f, 'fasta'), SeqIO.parse(depl_f, 'fasta')):    step = 200; jump = 1; a = 0; b = step; n_jumps = 5    for j in range(n_jumps):        s_intr = str(intr.seq)[a:b]; s_depl = str(depl.seq)[a:b]        intr_seqs.append(s_intr); depl_seqs.append(s_depl)        a = a + jump; b = a + step        e = e + 1    if e%20000 == 0:        print('Finished ' + str(e) + ' entries')        def getKmers(sequence, size):    return [sequence[x:x+size].upper() for x in range(len(sequence) - size + 1)]kmer = 10intr_texts = [' '.join(getKmers(i, kmer)) for i in intr_seqs]depl_texts = [' '.join(getKmers(i, kmer)) for i in depl_seqs]

我们最终得到737 340个句子，它们属于两个类：尼安德特人的渗入和枯竭。下一步是通过使用Python中的Tokenizer类将k-mers /word转换为整数，对句子进行独热编码。请注意，在我们的情况下，实际上并不需要填充，因为所有句子的长度都相同，长度为191 k-mers，但是为了通用性，我在这里使用了填充。

import numpy as npfrom keras.preprocessing.text import Tokenizerfrom sklearn.model_selection import train_test_splitfrom keras.preprocessing.sequence import pad_sequencesmerge_texts = intr_texts + depl_textslabels = list(np.ones(len(intr_texts))) + list(np.zeros(len(depl_texts)))tokenizer = Tokenizer()tokenizer.fit_on_texts(merge_texts)encoded_docs = tokenizer.texts_to_sequences(merge_texts)max_length = max([len(s.split()) for s in merge_texts])X = pad_sequences(encoded_docs, maxlen = max_length, padding = 'post')X_train,X_test,y_train,y_test = train_test_split(X,labels,                                                 test_size=0.20,random_state=42)vocab_size = len(tokenizer.word_index) + 1

在我们的情况下，词汇量为964114，小于4 ^ 10 = 1 048 576，这意味着并非所有由4个字符构建的10-mers都存在于序列中。最后，我们定义了一个顺序模型，首先是Keras嵌入层，它在对序列进行分类时学习wod嵌入，然后是双向LSTM和密集层。

from keras.models import Sequentialfrom keras.callbacks import ModelCheckpointfrom keras.optimizers import SGD, Adam, Adadelta, RMSpropfrom keras.layers import Embedding, LSTM, Dense, Bidirectionalmodel = Sequential()model.add(Embedding(vocab_size, 10))model.add(Bidirectional(LSTM(10)))model.add(Dense(10, activation = 'relu'))model.add(Dense(1, activation = 'sigmoid'))epochs = 5model.compile(loss='binary_crossentropy',optimizer='rmsprop',metrics=['accuracy'])checkpoint=ModelCheckpoint("LSTM.weights.best.hdf5", monitor='val_acc',verbose = 1,                           save_best_only = True, mode = 'max')model.summary()

在第一层使用Word嵌入的好处是可以将维度从964 114减少到10个维度。这减少了过度拟合，并通过迫使相似的单词共享拟合参数/权重来提高模型的泛化能力。现在我们可以开始训练我们的模型了。

import matplotlib.pyplot as plthistory = model.fit(X_train, y_train,                     epochs = epochs, verbose = 1, validation_split = 0.2,                     batch_size = 32, shuffle = True,                     callbacks = [checkpoint])plt.figure(figsize=(20,15))plt.plot(history.history['loss']); plt.plot(history.history['val_loss'])plt.title('Model Loss', fontsize = 20)plt.ylabel('Loss', fontsize = 20); plt.xlabel('Epoch', fontsize = 20)plt.legend(['Train', 'Validation'], fontsize = 20)plt.figure(figsize=(20,15))plt.plot(history.history['acc']); plt.plot(history.history['val_acc'])plt.title('Model Accuracy', fontsize = 20)plt.ylabel('Accuracy', fontsize = 20); plt.xlabel('Epoch', fontsize = 20)plt.legend(['Train', 'Validation'], fontsize = 20)

仅仅进行5个epoch的训练就足够了，因为模型对验证数据集的准确率很快达到了惊人的99%。。

对测试数据集上的模型性能进行评估，我们再次获得了99%的尼安德特人渐渗序列与枯竭序列分类的准确性 。

import itertoolsimport matplotlib.pyplot as pltfrom sklearn.metrics import confusion_matrixplt.figure(figsize = (15,10))predicted_labels = model.predict(X_test)cm = confusion_matrix(y_test, [np.round(i[0]) for i in predicted_labels])print('Confusion matrix:\n',cm)cm = cm.astype('float') / cm.sum(axis = 1)[:, np.newaxis]plt.imshow(cm, cmap = plt.cm.Blues)plt.title('Normalized confusion matrix', fontsize = 20)plt.colorbar()plt.xlabel('True label',fontsize=20); plt.ylabel('Predicted label',fontsize=20)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', verticalalignment = 'center',              fontsize = 20, color='white' if cm[i, j] > 0.5 else 'black')scores = model.evaluate(X_test, y_test, verbose = 0)print("Accuracy: %.2f%%" % (scores[1]*100))

现在，我们有了训练好的LSTM模型，我们将在稍后使用该模型来预测从尼安德特人继承的基因序列。现在是时候解释模型了，在是解释模型的时候了，该模型包括词汇的可视化和检测驱动分类的大多数预测性k-mers。

可视化单词嵌入

由于我们在网络的前端应用了嵌入层，因此目前每个k-mer /词都由10维向量表示。为了直观地想象埋层了解到，我们首先需要保存的权重层和单词的词汇。

import ioe = model.layers[0]weights = e.get_weights()[0]words = [i.upper() for i in list(tokenizer.index_word.values())]out_v = io.open('vecs.tsv', 'w', encoding='utf-8')out_m = io.open('meta.tsv', 'w', encoding='utf-8')for num, word in enumerate(words):  vec = weights[num + 1] # skip 0, it's padding.  out_m.write(word + "\n")  out_v.write('\t'.join([str(x) for x in vec]) + "\n")out_v.close()out_m.close()

为了可视化，我们可以使用来自http://projector.tensorflow.org/上的Tensorflow嵌入投影仪，并使用UMAP非线性降维技术研究k-mers的聚类。

通过嵌入层获得的k-mers与前一篇文章中的Word2Vec嵌入有很大的不同，没有发现AT-rich和GC-rich k-mers有明显的聚类。。

确定最具预测性的K-mers

建立神经网络特征重要性的一种方法是对输入数据施加一些干扰，并在网络输出处监视预测准确性的变化。在这里，我们想找到最有用的k-mers，但是，在我们的情况下，输入矩阵X的维度（737340，191），第一个维度代表的数量为神经网络训练的例子，第二个维度对应单词的数量在每个句子/序列。这意味着每个word/ k-mer 的索引（或句子中的位置）是特征，因此，如果我们对737340个样本中的191个特征逐个进行随机混洗，并检查与未扰动的输入矩阵相比，准确性的下降，我们就可以根据特征对最终预测的重要性对它们进行排序。在我们的例子中，排序特征相当于确定所有句子/序列中word/ k-mers的最重要位置。

import matplotlib.pyplot as pltperm_scores = []for i in range(X_test.shape[1]):    X_test_perm = X_test    X_test_perm[:,i] = random.sample(list(X_test[:,i]), len(list(X_test[:,i])))    perm_scores.append(abs(model.evaluate(X_test_perm, y_test, verbose = 0)[1]*100 -                            model.evaluate(X_test, y_test, verbose = 0)[1]*100))plt.figure(figsize=[16,5])barlist = plt.bar(np.arange(len(perm_scores)), perm_scores)plt.xlabel('Bases'); plt.ylabel('Magnitude of saliency values')

如果我们选择精度高于0.01的单词的位置/指数，我们会发现0、4、7、18、34、52、75、81、104、130、146和182是最重要的。也许，我们在句子的开头看到了一些丰富的重要词汇。现在我们简单地检查一下在以上位置的所有句子中最常见的word/ k-mers。

import pandas as pdfrom collections import Counterscores_df = pd.DataFrame({'Base': range(len(perm_scores)),'Score': perm_scores})informative_indices = list(scores_df[scores_df['Score'] > 0.01].index)reverse_word_map = dict(map(reversed, tokenizer.word_index.items()))def sequence_to_text(list_of_indices):    return [reverse_word_map.get(letter) for letter in list_of_indices]my_texts = list(map(sequence_to_text, X_test[informative_indices, :]))fig = plt.figure(figsize = (20, 15))D = dict(Counter([item.upper() for sublist in my_texts                   for item in sublist]).most_common(20))plt.bar(range(len(D)), list(D.values()), align = 'center')plt.title('Most Predictive K-mers', fontsize = 20)plt.ylabel("Counts", fontsize = 20); plt.xticks(rotation = 90)plt.xticks(range(len(D)), list(D.keys()), fontsize = 20)

令人放心的是，我们观察到，AT-rich k-mers在尼安德特人渐渗序列和枯竭序列之间的区别最大，正如先前的文章中所示。

预测尼安德特人基因

现在，我们将考虑人类基因序列和训练好的LSTM模型，以便预测每个人类基因从尼安德特人那里继承的可能性。在上一篇文章中，我解释了如何将基因序列提取到单独的fasta文件中，所以现在我们将读取这些序列，将它们分成k-mers，使用Tokenizer将k-mers转换为整数，并使用训练好的LSTM模型。

import pandas as pdfrom Bio import SeqIOfrom keras.preprocessing.text import Tokenizerfrom keras.preprocessing.sequence import pad_sequencese = 0; gene_seqs = []; gene_ids = []for gene in SeqIO.parse('hg19_gene_clean.fa', 'fasta'):    cutoff = 200    if len(str(gene.seq)) < cutoff:        continue    gene_ids.append(str(gene.id)); s_gene = str(gene.seq)[0:cutoff]    gene_seqs.append(s_gene)    e = e + 1    if e%10000 == 0:        print('Finished ' + str(e) + ' genes')def getKmers(sequence, size):    return [sequence[x:x+size].upper() for x in range(len(sequence) - size + 1)]kmer = 10gene_texts = [' '.join(getKmers(i, kmer)) for i in gene_seqs]tokenizer = Tokenizer(); tokenizer.fit_on_texts(gene_texts)encoded_docs = tokenizer.texts_to_sequences(gene_texts)max_length = max([len(s.split()) for s in gene_texts])X_gene = pad_sequences(encoded_docs, maxlen = max_length, padding = 'post')gene_predictions = model.predict_classes(X_gene)gene_predictions_prob = model.predict_proba(X_gene)gene_pred_df = pd.DataFrame({'Gene': gene_ids, 'Gene_Symbol': gene_symbol,                              'Predict': list(gene_predictions.flatten()),                              'Prob': list(gene_predictions_prob.flatten())})gene_pred_df = gene_pred_df.sort_values(['Prob'], ascending = False)gene_pred_df[(gene_pred_df['Predict'] == 1) & (gene_pred_df['Prob'] > 0.8)]

通过对尼安德特人遗传基因数量的预测，我们证实了之前文章的结论:绝大多数基因，即3.1万个被分析的人类基因中的2.2万个，都与尼安德特人的祖先无关。因此，我们再次得出结论，由于某种原因，进化将尼安德特人的祖先从基因中赶了出来，而基因是人类基因组中最具特征的元素。。仔细研究“幸存的”尼安德特人基因列表，并将其与先前文章中“语言包”模型中随机森林的预测进行比较，我们观察到许多与各种人类特征和疾病相关的基因。例如，LSTM和随机森林均以>99%的概率预测ANK3基因为尼安德特人祖先。已知该基因与躁郁症有关，主要在人脑中表达。

ANK3基因可能起源于尼安德特人与躁郁症有关

ANK3是一种基于GTEX的人类大脑特异性基因

此外，一些研究表明，与精神障碍有关的人脑部分可能是由尼安德特人的基因遗传塑造的。

人脑包藏着尼安德特人基因的“residual echo”

这就是进化科学可以为生物医学提供信息的地方，关于导致现代人类疾病的特征的历史发展和起源。

摘要

在这篇文章中，我们了解到，将LSTM等长记忆模型应用于DNA测序数据可能是有利的，因为已知这些数据具有沿序列的长期相关性。我们证明，LSTM取得了令人印象深刻的准确性，并优于所有其他模型在现代人类基因组中对尼安德特人基因渗入区域的检测。对该模型的解释表明，AT-rich k-mers是尼安德特人和自然人现代DNA序列之间的关键区别。我们预测大多数人类基因都没有尼安德特人的祖先，这意味着在进化过程中，尼安德特人的祖先会遭到强烈的负面选择。据预测，尼安德特人遗传的一些基因与人类有趣的特征有关。如躁郁症。