簡介

在本文中，我們將 DNA 測序糾錯(cuò)表述為多級(jí)分類問題，并提出兩種深度學(xué)習(xí)解決方案。第一種方法是在單次讀取中糾錯(cuò)，而第二種方法（如圖 1 所示）則通過多次讀取來達(dá)成共識(shí)，以預(yù)測正確的 DNA 序列。我們的Colab 筆記教程使用Nucleus和TensorFlow庫實(shí)現(xiàn)第二種方法。本文旨在向您展示如何同時(shí)使用 Nucleus 與 TensorFlow 解決基因組學(xué)領(lǐng)域的機(jī)器學(xué)習(xí)問題。

問題概覽

盡管 DNA 測序日漸快捷和便宜，其過程仍容易出錯(cuò)。使用Illumina等公司開發(fā)的新一代測序 (NGS) 技術(shù)處理原始數(shù)據(jù)時(shí)，錯(cuò)誤率約為 1%。第三代技術(shù)，例如Pacific BioSciences(PacBio) 公司開發(fā)的技術(shù)，正日益普及，其錯(cuò)誤率約為 15%。測序錯(cuò)誤可分為替換、插入和缺失，后兩者通常稱為 indel。所有這些錯(cuò)誤均不利于下游的分析步驟，例如變異檢測和基因組組裝。

如要獲取較高質(zhì)量的數(shù)據(jù)集，一個(gè)簡單的方法是舍棄可能包含錯(cuò)誤的數(shù)據(jù)，丟棄全部讀取內(nèi)容或去除低質(zhì)量區(qū)域皆可。該方法并非理想之選，因?yàn)檫@會(huì)導(dǎo)致最終的數(shù)據(jù)集會(huì)變小。此外，某些序列上下文本來就有較高的錯(cuò)誤率，進(jìn)而導(dǎo)致采樣出現(xiàn)偏差。因此，大量研究都側(cè)重于開發(fā)更成熟的糾錯(cuò)方法。大多數(shù)已開發(fā)的方法均可歸類為以下兩組之一：

對(duì)單次讀取進(jìn)行操作的方法，旨在確定正確的讀取序列

對(duì)多次讀取進(jìn)行操作的方法，以共識(shí)為基礎(chǔ)，旨在確定正確的基礎(chǔ) DNA 序列

深度學(xué)習(xí)概覽

本文中闡述的兩種方法均使用深度神經(jīng)網(wǎng)絡(luò)，學(xué)習(xí)將輸入映射至輸出的函數(shù)。神經(jīng)網(wǎng)絡(luò)由若干層線性與非線性運(yùn)算構(gòu)成，而這些運(yùn)算會(huì)依次應(yīng)用至輸入。神經(jīng)網(wǎng)絡(luò)已成功應(yīng)用于包括圖像分類和自然語言翻譯在內(nèi)的多個(gè)問題領(lǐng)域。最近，神經(jīng)網(wǎng)絡(luò)也被用于解決基因組學(xué)問題，例如蛋白質(zhì)結(jié)構(gòu)預(yù)測和變異檢測。

方法

Nucleus

我們的實(shí)現(xiàn)需要Nucleus，這是由 Google Brain 的 Genomics 團(tuán)隊(duì)開發(fā)的內(nèi)容庫，用于處理基因組學(xué)數(shù)據(jù)。Nucleus 使用專門的 reader 對(duì)象與 writer 對(duì)象，可輕松讀取、寫入和分析常見基因組文件格式（如 BAM、FASTA 和 VCF）中的數(shù)據(jù)。Nucleus 讓我們能夠：

針對(duì)指定基因組區(qū)域中的所有變異查詢 VCF 文件

針對(duì)映射至指定基因組范圍的所有讀取內(nèi)容查詢 BAM 文件

針對(duì)從指定位置開始的參考序列查詢 FASTA 文件

我們還能使用 Nucleus 將數(shù)據(jù)寫入TFRecords，這種二進(jìn)制文件格式由協(xié)議緩沖區(qū)消息構(gòu)成，可由 TensorFlow 輕松讀取。讀取 TFRecords 文件后，我們會(huì)使用Estimator API訓(xùn)練和評(píng)估卷積神經(jīng)網(wǎng)絡(luò)。

數(shù)據(jù)

以下是我們?cè)趯?shí)現(xiàn)過程中所使用的文件列表。所有數(shù)據(jù)均公開提供，且此 教程包含下載鏈接與說明。

NA12878_sliced.bam — 從 20 號(hào)染色體（位置 10,000,000–10,100,000）獲得的 Illumina HiSeq 讀取內(nèi)容，降采樣至 30x 的覆蓋度

NA12878_sliced.bam.bai — NA12878_sliced.bam 的索引

NA12878_calls.vcf.gz — 瓶中基因組 NA12878 變異的真值集合

NA12878_calls.vcf.gz.tbi — NA12878_calls.vcf.gz 的索引

hs37d5.fa.gz — hs37d5 的參考基因組

hs37d5.fa.gz.fai 和 hs37d5.fa.gz.gzi — hs37d5.fa.gz 的索引文件

注：教程 鏈接

https://colab.research.google.com/github/google/nucleus/blob/master/nucleus/examples/dna_sequencing_error_correction.ipynb

網(wǎng)絡(luò)架構(gòu)

卷積神經(jīng)網(wǎng)絡(luò)通常用于處理計(jì)算機(jī)視覺任務(wù)，但也非常適用于基因組學(xué)。每個(gè)卷積層都會(huì)反復(fù)將學(xué)習(xí)后的過濾器應(yīng)用于輸入數(shù)據(jù)。在網(wǎng)絡(luò)中早期出現(xiàn)的卷積過濾器會(huì)學(xué)習(xí)識(shí)別輸入數(shù)據(jù)的低級(jí)特征（如圖像中的邊緣及色彩梯度），而后期出現(xiàn)的過濾器則會(huì)學(xué)習(xí)識(shí)別更復(fù)雜的低級(jí)特征組合。對(duì)于 DNA 序列輸入，低級(jí)卷積過濾器會(huì)充當(dāng) motif 檢測器，這類似于序列標(biāo)識(shí)圖的位置權(quán)重矩陣。

在實(shí)現(xiàn)過程中，我們使用的標(biāo)準(zhǔn)卷積架構(gòu)依次由兩個(gè)卷積層及三個(gè)全連接層組成。我們使用非線性 ReLU 層提升模型的表現(xiàn)能力。當(dāng)卷積層減少輸入量后，我們會(huì)進(jìn)行最大池化，并會(huì)在全連接層充當(dāng)正則化矩陣后退出此過程。請(qǐng)注意，在得到最終的全連接層后，我們不會(huì)加入 softmax 層，因?yàn)槲覀兪褂玫膿p失函數(shù)是在內(nèi)部應(yīng)用 softmax。如需了解每層的詳情，請(qǐng)參閱此教程。

注：教程 鏈接

https://colab.research.google.com/github/google/nucleus/blob/master/nucleus/examples/dna_sequencing_error_correction.ipynb

方法 1：單次讀取的糾錯(cuò)

為了糾正序列讀取中的錯(cuò)誤，我們使用深度學(xué)習(xí)來訓(xùn)練神經(jīng)網(wǎng)絡(luò)，以解決一個(gè)較為普遍的問題：填充 DNA 序列中缺失的堿基。此方法旨在開發(fā)一種可理解 DNA 序列語法的模型。若僅靠真實(shí)序列的語法，我們可能無法獲取充足的信息來開發(fā)可用于生產(chǎn)環(huán)境的解決方案。盡管如此，這依然是一個(gè)簡單明了的示例應(yīng)用。

出于指導(dǎo)目的，我們通過以下方法簡化此問題：

僅考慮存在替換錯(cuò)誤的區(qū)域，并忽略 indel 錯(cuò)誤

僅考慮未存在已知變異的區(qū)域

我們可以在參考基因組的區(qū)域中訓(xùn)練該神經(jīng)網(wǎng)絡(luò)。此網(wǎng)絡(luò)的輸入是定長的 DNA 序列，其核心是我們希望預(yù)測的堿基。此網(wǎng)絡(luò)的輸出是可能出現(xiàn)的堿基分布，且最終預(yù)測結(jié)果為可能性最高的堿基。我們使用在參考基因組中觀測到的堿基產(chǎn)生標(biāo)簽集。由于我們僅使用映射至未存在已知真值變異之區(qū)域的讀取內(nèi)容，因此可以將參考基因組中存在的堿基明確標(biāo)記為標(biāo)簽。

我們將參考基因組分割成非重疊的定長片段，以產(chǎn)生輸入序列。在訓(xùn)練、評(píng)估和測試時(shí)，我們將參考序列中的一個(gè)堿基置零，以模擬缺失的堿基，如圖 3 所示（位置 5）。除了使用參考基因組來模擬缺失的數(shù)據(jù)之外，我們還可將此類模型應(yīng)用于序列讀取的數(shù)據(jù)，特別是質(zhì)量評(píng)分低于閾值的堿基。

方法 2：基于共識(shí)的糾錯(cuò)

糾錯(cuò)的最終目的是確定基礎(chǔ) DNA 序列，而非為了糾正單次讀取的錯(cuò)誤。在本部分，我們通過匯總序列堆疊來使用多次讀取達(dá)成的共識(shí)。如此一來，無需糾正單次讀取的中間步驟即可直接確定 DNA 序列。有關(guān)序列堆疊的示例如下方圖 4 所示。請(qǐng)注意，下圖僅展示了此窗口中存在的讀取部分。

出于指導(dǎo)目的，我們?cè)俅瓮ㄟ^以下方法簡化此問題：

僅考慮存在替換錯(cuò)誤的區(qū)域，并忽略 indel 錯(cuò)誤

僅考慮未存在已知變異的區(qū)域

與第一種方法不同，我們并未在參考基因組中訓(xùn)練此模型。相反，我們的訓(xùn)練數(shù)據(jù)來自所映射的 Illumina HiSeq 讀取內(nèi)容。此網(wǎng)絡(luò)的輸入是在所映射的讀取內(nèi)容中觀測到的標(biāo)準(zhǔn)化堿基數(shù)矩陣，其核心是我們希望預(yù)測的正確堿基的位置。Clairvoyante（一種用于變異檢測的神經(jīng)網(wǎng)絡(luò)）的作者以及Jason Chin 的示例方法中皆使用了類似的特征化方式。此網(wǎng)絡(luò)的輸出是可能出現(xiàn)的堿基分布，且最終預(yù)測結(jié)果為可能性最高的堿基。與第一種方法類似，我們使用在參考基因組中觀測到的堿基來產(chǎn)生標(biāo)簽集。我們將包含錯(cuò)誤（在堆疊中至少有一次讀取與中心位置中的參考序列不符）的示例及未包含錯(cuò)誤（堆疊中的所有讀取均與中心位置中的參考序列相符）的示例結(jié)合使用。

結(jié)論

此隨附教程演示了本文所述的第二種方法。盡管我們分析的示例較為簡單，不適合在生產(chǎn)環(huán)境中部署，但我們希望它們能幫助開發(fā)者學(xué)會(huì)高效利用 Nucleus 和深度學(xué)習(xí)解決基因組學(xué)領(lǐng)域的問題。