1、基于深度學習的第三代基因測序一致性序列生成摘 要：繼人類基因組計劃開展以來，基因測序已經(jīng)廣泛影響了生命科學的研究方式。基因組組裝是從大量隨機測序獲得的 短片段中重建出基因組長序列的過程，其最終目標是生成完整、準確的一致性序列，為后續(xù)多種研究提供可靠的參考基因組。 第三代基因測序技術(shù)可以產(chǎn)生讀長達幾十kb的片段，其應(yīng)用極大提高了基因組組裝的完整性，但測序的高錯誤率卻限制了最 終一致性序列的準確性。本研究提出基于深度學習的一致性序列生成模型，利用人工神經(jīng)網(wǎng)絡(luò)提取基因多序列比對結(jié)果的 結(jié)構(gòu)特征，生成準確率更高的一致性序列。實驗表明，該模型針對第三代測序數(shù)據(jù)可以生成質(zhì)量較高的一致性序列，并且無 需讀取

2、測序時的質(zhì)量值，也不用一次讀取超長序列，可以更加靈活地處理小數(shù)據(jù)塊。關(guān)鍵詞：第三代基因測序；一致性序列；深度學習；人工神經(jīng)網(wǎng)絡(luò)DLCC： Deep-Learning-Based Consensus Construction from long error-prone readsAbstract】Since the launch of the Human Genome Project, genome sequencing has widely influenced the way life sciences are studied. Genome assembly is the process

3、of reconstructing long genome sequences from a large number of short fragments obtained by random sequencing. Its ultimate goal is to generate complete, accurate consensus, providing reliable reference genome for subsequent analyses. The application of the third-generation sequencing technology, who

4、 can generate the read as long as dozens of Kb, has greatly improved the integrity of genome assembly, but the high error rate of the long read limits the accuracy of the final consensus sequences. This study proposes a consistent sequence generation model based on deep learning, which uses artifici

5、al neural network to extract the structural characteristics of gene multiple sequence alignment results to generate consensus with higher accuracy. Experiments show that the model can generate high-quality consensus for the third-generation sequencing data, while it does not need to read the quality

6、 value during sequencing, nor to read ultra-long sequences at a time, so it can process small data blocks more flexibly .Key words 】the third -generatio n sequencing; consensus; deep learning; artificial neural networko引言基因測序是一種新型基因檢測技術(shù),能夠從血 液或唾液中分析測定基因全序列，來預(yù)測罹患多種 疾病的可能性、個體的行為特征及行為?；驕y序 技術(shù)能鎖定個人病變基因，

7、提前預(yù)防和治療?；?測序相關(guān)產(chǎn)品和技術(shù)已由實驗室研究演變到臨床應(yīng) 用,可以說基因測序技術(shù)，是下一個改變世界的技 術(shù)。2020年初，隨著新型冠狀病毒肺炎疫情的爆 發(fā)，針對病毒基因組分析得到了更多的重視，快速獲 得新冠病毒基因組的參考序列是病毒核酸檢測和生 產(chǎn)疫苗的基礎(chǔ)。對多個病毒基因組和中間宿主 動物基因組的測序和比對,可以了解病毒的來源、感 染物種和在自然界中的變異情況，以此估計病毒的 傳播難易程度，從而為控制疫情的蔓延提供理論指 導。利用第三代測序片段進行基因組組裝、變異檢 測等已經(jīng)成為基因組學領(lǐng)域的基本分析手段。其測 序基因組覆蓋均勻、長讀長的優(yōu)勢極大提高了基因 組的組裝的完整度和連續(xù)性

8、。然而其讀長(1 kb- 50 kb)和高錯誤率(-15%)對組裝過程中的序列比 對產(chǎn)生了極大挑戰(zhàn)，也影響了最終組裝出的參考基 因組序列的準確性目前基于第三代測序序列 得到的組裝、分析結(jié)果仍然需要利用準確率較高的 第二代測序數(shù)據(jù)進行校正。因此,對第三代測序 數(shù)據(jù)進行拋光的重要性不言而喻，這樣能夠快速高 效地對其進行部分錯誤糾正，并生成準確率較高的 一致性序列，對長序列的組裝以及后續(xù)的基因組學 研究都有巨大的意義面。1第三代測序數(shù)據(jù)預(yù)處理引入深度學習方法進行一致性序列生成任務(wù), 首先要對原始數(shù)據(jù)預(yù)處理,使其適合作為人工神經(jīng) 網(wǎng)絡(luò)的輸入。本研究中采用的數(shù)據(jù)集為公共數(shù)據(jù)庫 中的Oxford Nano

9、pore Technology公司的納米孔測序 數(shù)據(jù)(即第三代基因測序ONT數(shù)據(jù))，選取的模式物 種包括大腸桿菌、酵母菌以及果蠅,數(shù)據(jù)集的具體信 息見表1。表1實驗所用ONT數(shù)據(jù)統(tǒng)計Tab. 1 The statistics of the ONT data used in the experimentsONT序列統(tǒng)計E.coli.YeastFly序列數(shù)量34 438146 967663 784平均長度6 71011 03866 956最大長度28 385414 134446 050最小長度2766 0005堿基總數(shù)231 078 8301 622 148 1844 617 586 601測序深

10、度54X134X34Xa數(shù)據(jù)來源：E.coli 和 Fly 數(shù)據(jù)下載自 NCBI( https: / HYPERLINK http:/www.ncbi www.ncbi. /sra)，索引號分別是 ERR1147230 和 SRR6702603 , Yeast 數(shù)據(jù)下載自 http: / / HYPERLINK s.fr/ s.fr/ externe/Download/Projets/ yeast / datasets / raw_data/S288C/3個模式物種的參考基因組長度分別為4,641, 652 個堿基、12,071,326 個堿基、137,547,960 個堿 基。由于這3個模式

11、物種的全基因組長度、復(fù)雜程 度及雜合情況均有較大差異，可以較好地評估本文 提出方法的完整性和全面性。酵母菌參考基因組序 列的片段示例如圖1所示。chr1 tpg|BK006935.2|organism=Saccharomyces cerevisiae S288cstrain=S288c moltype=genomic chromosome=l note=R64-1-l CCACACCACACCCACACACCCACACACCACA CCACACACCACACCACACCCACACACACACA TCCTAACACTACCCTAACACAGCCCTAATCTA ACCCTGGCCAACCTGTC

12、TCTCAACTTACCCTC CATTACCCTGCCTCCACTCGTTACCCTGTCCC 圖1酵母參考基因組染色體1序列片段Fig. 1 The chromosome 1 fragment of Yeast reference genome經(jīng)過不斷試驗，對第三代測序原始數(shù)據(jù)采用如 下處理方式:首先，測序過程不同通量數(shù)據(jù)并非是完 全對齊，而是呈階梯狀排布，相鄰兩條之間都有一定 堿基數(shù)的錯位,因此需要將測序序列比對到初始參 考基因組骨架上，確定測序序列之間的排列順序，此 步驟通過快速比對軟件miniasm2完成;其次，裁剪 比對后的結(jié)果，使大部分序列可以頭尾對齊。以深 度為40X的數(shù)據(jù)塊為

13、例，具體做法是取正向第十五 條序列的尾部向前50個位點作為塊截止位點，反向 第十五條序列頭部向后50個位點作為塊起始位點； 最后，將裁剪后的比對結(jié)果分割為深度為序列乘數(shù)、 寬度為12個位點的小塊，每一個小塊用于預(yù)測塊中 心4個位點的堿基種類，將每個小塊的預(yù)測結(jié)果最 終拼接成完整的一致性序列。2 致性序列生成模型一致性序列生成任務(wù)是通過高通量的序列計算 得到的,并且基因序列中存在一定的結(jié)構(gòu)相關(guān)性,即 序列某一位點前后的堿基對此位點堿基的預(yù)測會起 到影響作用，這些特征從直覺上符合人工神經(jīng)網(wǎng)絡(luò) 適用的場景。最終的實驗結(jié)果論證了在一條序列中 以及多條序列間均存在相關(guān)性，也證明了采用深度 學習方法的可靠

14、性。具體采用的網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。主要分為卷 積神經(jīng)網(wǎng)絡(luò)模塊、通道注意力機制模塊、循環(huán)神經(jīng)網(wǎng) 絡(luò)模塊、多任務(wù)學習模塊四部分。四個結(jié)構(gòu)相結(jié)合, 可以有效地提取出一條序列內(nèi)部的結(jié)構(gòu)相關(guān)性以及 不同通量序列之間的測序時序上的相關(guān)性，還可以 賦予對一致性序列影響較大的部分序列更大的權(quán) 重，將提取到的不同層次相關(guān)性充分利用,以得到更 優(yōu)秀的結(jié)果。圖2 一致性序列生成網(wǎng)絡(luò)結(jié)構(gòu)圖Fig. 2 Consensus construction network網(wǎng)絡(luò)中第一部分是卷積模塊，卷積神經(jīng)網(wǎng)絡(luò) (Convolutional Neural Network , CNN)是深度學習中 一種被廣泛應(yīng)用的網(wǎng)絡(luò)結(jié)構(gòu)，在深度

15、學習技術(shù)快速 發(fā)展的背景之下，首先在計算機視覺領(lǐng)域收獲了卓 越的成績。在一致性序列生成任務(wù)中，由于輸入 到卷積模塊中的數(shù)據(jù)是切割成固定尺寸的第三代測 序數(shù)據(jù)比對塊，其通量和每條序列包含的位點個數(shù) 都是確定的。分析比對塊的特點可知，這是行數(shù)代 表序列通量、列數(shù)代表每條序列位點數(shù)的堿基矩陣。 通過合理推測以及實驗證明，不同序列可以看作不 同時序的信號，用卷積結(jié)構(gòu)對其進行特征提取的效 果并不理想。結(jié)合上述的分析，由于高通量的基因 序列數(shù)據(jù)不同通量之間有時序關(guān)系，經(jīng)過對比發(fā)現(xiàn) 更適合使用循環(huán)網(wǎng)絡(luò)來利用其相關(guān)性。因此，卷積 模塊中采用1x3的一維卷積核，僅僅用作橫向的特 征提取，即提取一條序列中堿基之間

16、的相關(guān)性特征， 此做法不僅更合理地利用了數(shù)據(jù)的特性進行模型設(shè) 計，還可以減少參數(shù)量,使得運算效率大大增加，在 更短時間獲得質(zhì)量更高的結(jié)果，具體卷積的形式如 圖3所示。同樣的，也采用1x2的平均池化，用于 減少參數(shù)量，有效地壓縮數(shù)據(jù)以及參數(shù)的規(guī)模，減小 計算的復(fù)雜性，降低時間代價。此外,考慮到如果采 用平均池化的方式，可能會使得對某個位點堿基種 類影響最大的堿基信息變得模糊，受到其他重要程 度更低的堿基信號影響，使得最終結(jié)果不能達到令 人滿意的準確率，因此決定在模型中采用最大池化 的形式。Fig. 1 One dimensional convolution to extract lateral

17、features在一致性序列預(yù)測任務(wù)中，預(yù)測的堿基在塊中 的位置基本確定，但是需要更多地利用周圍堿基的 信息，以糾正測序過程中產(chǎn)生的錯誤,而非僅僅考慮 對應(yīng)位置的堿基種類。通道注意力機制是通過引入 通道之間的相關(guān)性，而并不會只著重關(guān)注部分堿基 情況。模型中采用 Squeeze-and - Excitation Block 的 結(jié)構(gòu),首先是Squeeze操作，順著空間維度特征壓 縮，將每個二維的特征通道變成一個實數(shù)，這個實數(shù) 某種程度上具有全局的感受野，并且輸出的維度和 輸入的特征通道數(shù)相匹配，表征著在特征通道上響 應(yīng)的全局分布，而且使得靠近輸入的層也可以獲得 全局的感受野；其次是Excitat

18、ion操作，是一個類 似于循環(huán)神經(jīng)網(wǎng)絡(luò)中門的機制，通過參數(shù)來為每個 特征通道生成權(quán)重，其中參數(shù)被學習用來顯式地建 模特征通道間的相關(guān)性。在全局平均池化(Global- Average- Pooling)之后有兩個全連接層(Fully - Connected)，分別應(yīng)用ReLU以及Sigmoid激活函 數(shù)，具有更多的非線性，可以更好地擬合通道間復(fù)雜 的相關(guān)性，極大地減少了參數(shù)量和計算量。通過一 個Sigmoid的門，獲得01之間歸一化的權(quán)重，最后 將歸一化后的權(quán)重加權(quán)到每個通道的特征上。提取出橫向特征并賦予通道注意力后,將會輸 入循環(huán)神經(jīng)網(wǎng)絡(luò)模塊進行不同通量間的特征提取。 循環(huán)神經(jīng)網(wǎng)絡(luò)(Recu

19、rrent Neural Network，RNN)同卷 積神經(jīng)網(wǎng)絡(luò)一樣，是深度學習中一種重要的網(wǎng)絡(luò)結(jié) 構(gòu),其目的是對時序信號有效的處理。循環(huán)神經(jīng)網(wǎng) 絡(luò)的特點在于考慮了原始數(shù)據(jù)中的時序相關(guān)性，在 這個網(wǎng)絡(luò)結(jié)構(gòu)中每一時刻的輸出不僅僅與當前的輸 入有關(guān)，還與之前所有時刻的輸入有著很大的聯(lián)系, 這樣側(cè)重時間相關(guān)性特點的網(wǎng)絡(luò)結(jié)構(gòu)，對于處理有 著明顯時序性特征的信號處理任務(wù)有著巨大的幫 助3。一個經(jīng)典的循環(huán)神經(jīng)網(wǎng)絡(luò)基本結(jié)構(gòu)如圖4所在一致性序列生成模型中，循環(huán)神經(jīng)網(wǎng)絡(luò)部分 采用三層的雙向GRU結(jié)構(gòu),GRU是LSTM的一種變 體,相比于 LSTM 的 3 個門(forget, input, output),

20、GRU只有兩個門(update和reset)。綜合來看兩者 的性能在很多任務(wù)上不分伯仲，但是GRU參數(shù)相對 少，更容易收斂，同時可以一定程度地避免梯度消失 的問題，更適用于本課題的任務(wù)。而雙向GRU是將 前向GRU與后向GRU結(jié)合,用于解決單向的結(jié)構(gòu) 無法編碼從后向前序列信息的問題，這樣可以捕捉 更全面的不同通量序列之間的語義依賴問。此模塊主要分為三部分:通量信息提取層、通量 信息表示層、特定位點預(yù)測層。通量信息提取層：由卷積提取的各個通道 特征進行列維度的拼接后，得到通量信息矩陣，作為 下一步的輸入；通量信息表示層：由于不同通量的序列可 以視作是不同時序的結(jié)果，考慮到測序過程前后兩 個方向可

21、能均包含時序上的信息，因此使用三層雙 向的GRU結(jié)構(gòu)對通量信息建模，從前向后的編碼信 息與從后到前的編碼信息相結(jié)合，此雙向的編碼結(jié) 構(gòu)作為通量信息的表示；特定位點預(yù)測層：由一層感知器結(jié)構(gòu)的全 連接層完成,作為某個位點堿基種類的預(yù)測方式。在循環(huán)網(wǎng)絡(luò)之后，加入多任務(wù)的結(jié)構(gòu)，即對一個 輸入塊最終將輸出4個位點的預(yù)測結(jié)果，多個任務(wù) 同時預(yù)測時，減少了數(shù)據(jù)來源的數(shù)量以及整體模型 參數(shù)的規(guī)模，使預(yù)測更加高效。4個任務(wù)相關(guān)性較 強,經(jīng)過實驗驗證，可以有效提升一致性序列預(yù)測的 準確率。3實驗設(shè)置及結(jié)果分析本研究在酵母菌、大腸桿菌、果蠅3個模式物種 的納米孔單分子測序數(shù)據(jù)上從多個指標來評價所提 出一致性序列生成

22、模型的性能，這3個物種的ONT 序列長度、染色體條數(shù)、復(fù)雜程度以及雜合程度均有 不同,因此可以不同角度考量模型的可靠程度，使得 評價結(jié)果不具有特殊性或偏好性。對最終一致性序列質(zhì)量的評價指標,最直觀也 是最重要的一項就是通過與已發(fā)布的“金標準”參 考基因組進行序列比較，計算與整個參考基因組匹 配的堿基數(shù)目和預(yù)測準確率。但是由于組裝過程中 可能存在組裝不完整、大片段缺失、倒位等問題，引 起一致性序列與參考基因組比對時的困難，因此僅 依靠匹配數(shù)來衡量準確率不能完全展示拋光后的一 致性序列的質(zhì)量。參考本領(lǐng)域內(nèi)其他軟件的評價方 法，本文考慮以如下四點標準衡量該模型：(1 )覆蓋參考基因組的位點數(shù)(cov

23、er reference)，用于評價生成一致性序列的完整度，即 從第一個匹配位點到最后一個匹配位點之間的數(shù) 目，不包含長段的缺失，但是包含插入、小段缺失以 及替換等錯誤。產(chǎn)生的替換數(shù)目，即發(fā)生在某個位點上一 致性序列中的堿基與參考基因組同一位置的堿基不 同，但是產(chǎn)生不同的原因也可能是此位點存在雜合 情況，由于出現(xiàn)的概率不大，因此在統(tǒng)計中均計入替 換數(shù)目。產(chǎn)生的缺失數(shù)目，這也是評判一致性序列 質(zhì)量時主要關(guān)注的問題。小數(shù)目的缺失主要是由第 三代測序原理的特性造成的。當序列中出現(xiàn)連續(xù)的 重復(fù)堿基時，測序儀對熒光信號或者電信號的解讀 精確性不夠，導致轉(zhuǎn)換為堿基信號時容易缺失堿基， 這也是第三代測序技術(shù)

24、準確率提升的最大瓶頸。這 個缺陷在拋光過程中可以被部分地解決，因此在評 價方法有效性時,對不通長度區(qū)間的缺失進行統(tǒng)計， 長度單位為堿基對(bp),主要分為1 bp、2 bp、3-5O bp、50-1 000 bp這4個長度范圍。特別地，由于果 蠅的全基因組包含堿基數(shù)目更為龐大，產(chǎn)生的錯誤 也更嚴重，因此還需統(tǒng)計長度大于1 000 bp的缺失。產(chǎn)生的插入數(shù)目，與缺失數(shù)目的統(tǒng)計方法 類似，分為 1 bp、2 bp、3-50 bp、50-1,000 bp 這 4 個 長度范圍進行計算。同樣地，對于果蠅數(shù)據(jù),還將統(tǒng) 計其長度大于1 000 bp的插入。本文的實驗覆蓋了酵母、大腸桿菌、果蠅這3個 模式物

25、種的全基因組納米孔測序數(shù)據(jù)，可以對模型 在不同結(jié)構(gòu)特征染色體上的性能全面地評估。針對 酵母基因組,增加了一組在不同數(shù)據(jù)量上的實驗，以 研究該模型在不同測序深度上的性能。通過以上實 驗可以清晰地得知本文提出方法在不同物種、不同 規(guī)模數(shù)據(jù)、不同裁剪尺寸上的實際效果。作為比較, 對比實驗將采用目前主流的具有一致性序列生成功 能且結(jié)果優(yōu)秀的軟件進行，包括Wtdbg、Canu、Flye 以及Racon。相比于目前主流的基于第三代測序數(shù)據(jù)的一致 性序列生成或者基因組裝拋光軟件，本文提出的基 于深度學習方法的模型在準確率上有明顯提高，在 參考基因組覆蓋長度、替換數(shù)目、缺失數(shù)目、插入數(shù) 目四項指標上均一定程度

26、地表現(xiàn)出優(yōu)勢，也證明了 算法的有效性及可靠性。表2記錄了采用不同方法對酵母納米孔測序數(shù) 據(jù)拋光的效果對比。與其他軟件生成的一致性序列 相比較，本文提出的方法在覆蓋參考基因組的位點 數(shù)上有明顯的優(yōu)勢，同時替換以及短缺失問題均有 明顯減少。然而，由于在處理原始數(shù)據(jù)時為了盡量 減少大量的空位(gap)以及減少含有重復(fù)片段的序 列，因此導致結(jié)果中大片段的插入數(shù)量稍多。但從 插入的位點數(shù)而言，卻依舊保持了較好的效果，所生 成的一致性序列質(zhì)量有明顯提升。表3記錄了采用不同方法對大腸桿菌納米孔測 序數(shù)據(jù)拋光的效果對比。大腸桿菌基因組結(jié)構(gòu)簡 單，基因組數(shù)據(jù)長度較短，因此測序和組裝過程中得 到的序列質(zhì)量較高，最

27、終生成的一致性序列也都有 很高的準確率。在具體實驗中，本文提出方法在參 考基因組覆蓋率上已達到百分之百,并且在替換、短 缺失兩項依然保持著優(yōu)秀的效果。此外，數(shù)據(jù)預(yù)處 理的方式導致長段插入略微增加，但各項指標依舊 有著不錯的結(jié)果。表4記錄了采用不同方法對果蠅納米孔測序數(shù) 據(jù)拋光的效果對比。果蠅染色體蘊含的信息相比于 前兩個模式物種復(fù)雜程度有極大的提升，不僅長度 的量級有所增加，也有更多雜合位點以及更難處理 的結(jié)構(gòu)信息。不同組裝軟件構(gòu)建的果蠅的基因組結(jié) 果質(zhì)量參差不齊,大片段錯誤較多，為后續(xù)一致性序 列的生成、拋光帶來了非常大的挑戰(zhàn)。本文提出的 模型相比于其他方法依然保持了不錯的效果，在參 考基因

28、組覆蓋率、缺失、替換、插入等各項指標均名 列前茅，但不可避免地由于數(shù)據(jù)本身以及數(shù)據(jù)預(yù)處 理的特點，長段缺失較多，且第四條染色體的堿基發(fā) 生了大量的缺失。總的來說,雖然最終結(jié)果有一些 不足，但是基于深度學習的模型仍然提供了令人滿 意的一致性序列生成結(jié)果。表2評估不同方法對酵母納米孔測序數(shù)據(jù)拋光的效果對比Tab. 2 Assessment of consensus sequences of yeast genomeYeast 130_WtdbgCanuFlyeRaconDLCC參考基因組長度12 071 32612 071 32612 071 32612 071 32612 071 326Cove

29、r Ref11 652 07511 871 96111 803 67011 647 84012 041 917替換11 52111 15021 72431 07711 1641 bp缺失66 97590 13666 30148 58643 8722 bp缺失5293127437 1144 7324 9603-50 bp缺失1 4357 6404 8703 2932 48950-1 000 bp 缺失331351 bp插入10 0751 4258 91910 1785 8692 bp插入4 7316991 4529 7051 0363-50 bp插入5 98462663626 5682 1735

30、0-1 000 bp 插入4787164表3評估不同方法對大腸桿菌納米孔測序數(shù)據(jù)拋光的效果對比Tab. 3 Assessment of consensus sequences of E.coli genomeE.coli 54_WtdbgCanuFlyeRaconDLCC參考基因組長度4 641 6524 641 6524 641 6524 641 6524 641 652Cover Ref4 619 8754 621 4104 633 6814 640 2364 641 652替換1 2102483768474331 bp缺失19 30315 02811 79312 72510 0432 b

31、p缺失4 8734 6422 8823 0462 8263-50 bp缺失92290743845253850-1000 bp 缺失111111 bp插入928531 3521 4831 1892 bp插入48372575553-50 bp插入91194891150-1 000 bp 插入22223表4評估不同方法對果蠅納米孔測序數(shù)據(jù)拋光的效果對比Tab. 4 Assessment of consensus sequences of fly genomeFly 32_Wtdbg mCanu幻Flye19Racon EDLCCRef長度137 547 960137 547 960137 547 960137 547 960137 547 960Cover Ref122 684 667128 546 930130 271 254128 374 616128 445 616替換179 64696 143110 768199 50398 4221 bp缺失1 147 336755 774314 651384 080307 3492 bp缺失164 971122 63224 66743 91732 7683-50 bp缺