基因組重測(cè)序的應(yīng)用領(lǐng)域有哪些

背景介紹

1.1

基因測(cè)序在癌癥領(lǐng)域上的應(yīng)用

癌癥是目前人類(lèi)所面臨的最大敵人，其發(fā)病率隨著年齡增長(zhǎng)而升高，在人口老齡化加劇的當(dāng)下，全社會(huì)的癌癥負(fù)擔(dān)也將愈發(fā)嚴(yán)峻。癌癥難以治愈的原因之一是腫瘤具有強(qiáng)異質(zhì)性，即相同癥狀、相同病理改變卻可能由完全不同的基因變化而造成，以至于同類(lèi)型癌癥患者對(duì)于相同藥物的藥效反應(yīng)有很大的差別 。這些基因上的變化包含了單堿基突變（SNV）、小片段插入缺失（InDel）、結(jié)構(gòu)與拷貝數(shù)變異（SV & CNV）和甲基化等。隨著近幾年基因測(cè)序成本如圖 1所示不斷下降，在萬(wàn)元內(nèi)即可完成人類(lèi)的全基因組測(cè)序，GPU的技術(shù)發(fā)展也帶來(lái)分析成本與時(shí)間的下降，于是用于檢測(cè)基因組變化的重測(cè)序技術(shù)在癌癥治療中起到了越來(lái)越重要的作用?；蚪M重測(cè)序的應(yīng)用主要有三個(gè)方向：

01

基因測(cè)序技術(shù)提供腫瘤內(nèi)基因分子水平的變化，為研究提供預(yù)見(jiàn)性，支持腫瘤新靶向藥物研發(fā)。

02

提高腫瘤早期診斷準(zhǔn)確率，降低高危風(fēng)險(xiǎn)，從而提高癌癥存活率。以在兒童中高發(fā)的兒童神經(jīng)母細(xì)胞瘤（NB）為例，目前國(guó)內(nèi)低、中危組NB患兒，長(zhǎng)期無(wú)瘤生存率可達(dá)90%，但高危組NB患兒的5年總生存率不足30%。而如果能夠做到早期診斷的話就可以通過(guò)手術(shù)切除聯(lián)合化療的方式進(jìn)行治療，從而提高治愈率。

03

基于易感基因分子水平的基因變化檢測(cè)，可以根據(jù)變化情況為患者提供個(gè)體化和預(yù)見(jiàn)性的治療。

圖 1 每兆數(shù)據(jù)DNA測(cè)序開(kāi)銷(xiāo)

圖中數(shù)據(jù)來(lái)源于Wetterstrand KA. DNA Sequencing Costs: Data from the NHGRI Genome Sequencing Program (GSP) 和 Next-Generation-Sequencing.v1.10.25

1.2

基因測(cè)序手段介紹

當(dāng)下基因組重測(cè)序常見(jiàn)的測(cè)序手段分為以Illumina雙端150bp為代表的短讀長(zhǎng)測(cè)序與以PacBio HiFi或Oxford Nanopore Technology（ONT）為代表的長(zhǎng)讀長(zhǎng)測(cè)序兩種。如圖1所示，Illumina雙端測(cè)序?qū)NA打斷兩端加上不同接頭與流動(dòng)池中分布的接頭結(jié)合不斷合成、退火與擴(kuò)增，擴(kuò)增后利用帶熒光的dNTP與流動(dòng)池中的DNA鏈結(jié)合，由于不同的dNTP發(fā)出不同熒光，于是可以通過(guò)掃描流動(dòng)池中的熒光信號(hào)，得到序列信息；PacBio 將單分子與DNA聚合酶固定在零模波導(dǎo)孔上，激光從孔的底部照入，激光不能穿過(guò)小孔且只能照亮孔底的一小部分區(qū)域。當(dāng)DNA聚合酶使用游離的dNTP進(jìn)行聚合反應(yīng)時(shí)，dNTP上的熒光基團(tuán)被激光激發(fā)，通過(guò)傳感器將光信號(hào)轉(zhuǎn)化為電信號(hào)。ONT則是通過(guò)讓 DNA 單分子通過(guò)帶電的納米孔，不同的堿基使得電流發(fā)生改變，傳感器記錄電流變化，測(cè)序儀將電流變化轉(zhuǎn)化為堿基值。

圖 2 測(cè)序原理示意圖 (a) Illumina 雙端測(cè)序，

圖片來(lái)源于DMLapato, Illumina dye sequencing, https://en.wikipedia.org/w/index.php?title=Illumina_dye_sequencing&oldid=1159305769#cite_note-Illumina,_Inc_2013-7

(b) PacBio 零模波導(dǎo)孔單分子實(shí)時(shí)測(cè)序，

圖片來(lái)源于百邁客生物, pacbio三代測(cè)序儀, http://www.biomarker.com.cn/technology-services/pacbio

(c) ONT 納米孔測(cè)序，

圖片來(lái)源于 Kraft, Long-read sequencing in human genetics, DOI: 10.1007/s11825-019-0249-z

短讀長(zhǎng)測(cè)序的優(yōu)勢(shì)在于價(jià)格低廉與測(cè)序堿基準(zhǔn)確度高，長(zhǎng)讀長(zhǎng)測(cè)序擁有以下幾個(gè)優(yōu)點(diǎn)：

01

長(zhǎng)讀長(zhǎng)測(cè)序使用單分子測(cè)序，不需要擴(kuò)增反應(yīng)。一方面節(jié)省了時(shí)間，另一方面也避免了擴(kuò)增中引入的錯(cuò)誤，比如擴(kuò)增效率不同帶來(lái)的偏好性問(wèn)題。

02

長(zhǎng)讀長(zhǎng)測(cè)序相比于短讀長(zhǎng)測(cè)序更容易跨過(guò)基因組的重復(fù)區(qū)域。在短讀長(zhǎng)測(cè)序結(jié)果中，至少748個(gè)與人類(lèi)相關(guān)的基因存在測(cè)序暗黑區(qū)域，這個(gè)區(qū)域序列的Read會(huì)出現(xiàn)低深度、低覆蓋度及低比對(duì)質(zhì)量。而造成這個(gè)現(xiàn)象的主要原因是在參考基因組中該區(qū)域的序列多次重復(fù)。而測(cè)序序列長(zhǎng)到能夠跨過(guò)重復(fù)區(qū)域的長(zhǎng)讀長(zhǎng)測(cè)序正好解決了這個(gè)問(wèn)題。

03

長(zhǎng)讀長(zhǎng)測(cè)序相比于短序列可以更容易得到核酸的甲基化水平。短讀長(zhǎng)測(cè)序基于熒光信號(hào)，通常需要將樣本分為兩份，一份利用亞硫酸氫鹽處理將未甲基化的C轉(zhuǎn)化為U，另一份直接測(cè)序進(jìn)行比較獲得基因組甲基化位點(diǎn)。而基于電信號(hào)的三代測(cè)序可以直接區(qū)分出甲基化的C與未甲基化的C，避免了亞硫酸氫鹽等藥品處理和多次測(cè)序引入的錯(cuò)誤。

而近幾年來(lái)測(cè)序的技術(shù)有兩個(gè)明顯的發(fā)展趨勢(shì)：

01

短讀長(zhǎng)測(cè)序和長(zhǎng)讀長(zhǎng)測(cè)序成本均逐步降低。短讀長(zhǎng)測(cè)序成本降低幅度變緩，長(zhǎng)讀長(zhǎng)測(cè)序的價(jià)格隨著儀器的升級(jí)已經(jīng)降到了三年前短讀長(zhǎng)測(cè)序成本。在一些應(yīng)用場(chǎng)景中，比如結(jié)構(gòu)變異的檢測(cè)，由于長(zhǎng)讀長(zhǎng)測(cè)序只需要三分之一短讀長(zhǎng)測(cè)序的數(shù)據(jù)量，長(zhǎng)讀長(zhǎng)測(cè)序的成本已經(jīng)低于短讀長(zhǎng)測(cè)序。

02

隨著試劑與算法的研發(fā)，長(zhǎng)讀長(zhǎng)測(cè)序的準(zhǔn)確度也逐漸接近短讀長(zhǎng)測(cè)序。在單堿基質(zhì)量和測(cè)序價(jià)格的保證下，長(zhǎng)序列能保留更原始的基因組特征。

于是目前研究所選擇的測(cè)序手段逐漸從短讀長(zhǎng)測(cè)序向長(zhǎng)讀長(zhǎng)測(cè)序轉(zhuǎn)移。

1.3

基因組重測(cè)序的流程介紹

如圖3所示，重測(cè)序的流程主要包含以下大步驟：

圖 3 基因組重測(cè)序流程與軟件示意圖。

圖中紫色矩形部分為主要分析步驟，黃色矩形部分為分析類(lèi)別，綠色矩形部分為軟件，藍(lán)色上標(biāo)的軟件為利用到了機(jī)器學(xué)習(xí)的軟件。

01

實(shí)驗(yàn)設(shè)計(jì)、分子實(shí)驗(yàn)與上機(jī)測(cè)序：這一步包括了需求分析和選擇測(cè)序手段。分子實(shí)驗(yàn)通常包括核酸提取、核酸分子打斷和樣品的制備。上機(jī)測(cè)序?qū)悠芳尤霚y(cè)序儀中，測(cè)序儀產(chǎn)生序列信息，完成測(cè)序。

02

數(shù)據(jù)預(yù)處理：通常測(cè)序儀會(huì)利用傳感器捕獲電信號(hào)，測(cè)序儀的Basecalling步驟和甲基化檢測(cè)步驟將電信號(hào)轉(zhuǎn)化為文本形式的序列?；蚪M測(cè)序的結(jié)果保存檢測(cè)得到每一個(gè)核酸分子的堿基和該堿基的質(zhì)量值，ONT測(cè)序和PacBio測(cè)序也會(huì)保存原始電信號(hào)。甲基化檢測(cè)的結(jié)果會(huì)同時(shí)輸出每一個(gè)堿基的甲基化概率。長(zhǎng)序列檢測(cè)會(huì)發(fā)生隨機(jī)的測(cè)序錯(cuò)誤，被認(rèn)為準(zhǔn)確度不高，需要進(jìn)行堿基修正，修正測(cè)序錯(cuò)誤來(lái)提升堿基質(zhì)量值。

03

比對(duì)：這一步驟將測(cè)序得到的序列與參考基因組比對(duì)，確定該序列來(lái)自于參考基因組的哪個(gè)部分。RNA測(cè)序得到的序列相比于DNA缺失了內(nèi)含子部分，長(zhǎng)讀長(zhǎng)測(cè)序與短讀長(zhǎng)測(cè)序在錯(cuò)誤率、復(fù)雜度上有區(qū)別。所以軟件上會(huì)根據(jù)樣本的類(lèi)型和測(cè)序序列長(zhǎng)度做區(qū)分。

04

變異檢測(cè)：變異檢測(cè)將獲得樣本與參考基因組的差異，這些差異包括了單堿基的變異（SNV）、小的插入缺失（InDel）以及大的結(jié)構(gòu)變異（SV）。短讀長(zhǎng)測(cè)序由于單堿基的準(zhǔn)確度高，在SNV和InDel上有優(yōu)勢(shì)，而對(duì)于長(zhǎng)度大于短讀長(zhǎng)測(cè)序序列長(zhǎng)度的結(jié)構(gòu)變異而言，長(zhǎng)讀長(zhǎng)測(cè)序更有優(yōu)勢(shì)。所以軟件上會(huì)根據(jù)需要檢測(cè)的變異長(zhǎng)度和測(cè)序序列長(zhǎng)度做區(qū)分。

05

后續(xù)分析：后續(xù)分析則包括了變異的過(guò)濾、注釋、比較以及關(guān)聯(lián)分析等過(guò)程，這一步將以萬(wàn)計(jì)的變異中篩選出與研究課題息息相關(guān)的變異，會(huì)跟隨著課題不同而不同。

而這五個(gè)大步驟之中，第一步的時(shí)間改進(jìn)則依賴(lài)于實(shí)驗(yàn)設(shè)計(jì)的改良以及試劑的研發(fā)。而后續(xù)步驟則需要依賴(lài)于算法的改進(jìn)和計(jì)算設(shè)備的提升，在本文中將會(huì)介紹這些分析的常用軟件和算法，并指出這些軟件均可以通過(guò)對(duì)GPU和并行算法的優(yōu)化進(jìn)行提速。比如Carpi等人利用合理實(shí)驗(yàn)設(shè)計(jì)和GPU，使得惡性瘧原蟲(chóng)在計(jì)算成本減少4.6倍的同時(shí)整體分析速度提升27倍。

二

測(cè)序數(shù)據(jù)預(yù)處理

2.1

利用GPU獲得堿基序列

將測(cè)序過(guò)程中產(chǎn)生的電信號(hào)轉(zhuǎn)化為堿基的步驟被稱(chēng)為 Basecalling。以 ONT測(cè)序?yàn)槔?，如圖 4所示，當(dāng)核酸分子通過(guò)芯片上納米孔時(shí)，生成不斷變化的電信號(hào)，測(cè)序儀記錄下電信號(hào)，并轉(zhuǎn)化為堿基。

圖 4 ONT Basecalling 過(guò)程示意圖

圖片來(lái)源于What is Oxford Nanopore Technology (ONT)sequencing?,https://www.yourgenome.org/facts/what-is-oxford-nanopore-technology-ont-sequencing/

由于堿基通過(guò)納米孔的速度并不勻速，電流信號(hào)與預(yù)測(cè)得到的堿基單調(diào)但非對(duì)齊，為了解決這個(gè)問(wèn)題，ONT判斷堿基的序列則是通過(guò)語(yǔ)音識(shí)別中常用的 CTC 作為損失函數(shù)來(lái)實(shí)現(xiàn)。ONT開(kāi)源Basecalling 軟件Dorado 的r10.4.1系列模型中，模型分為編碼器和解碼器兩個(gè)部分，解碼器的模型如圖 3所示，其接續(xù)線性條件隨機(jī)場(chǎng)（CRF）；解碼器則是用Viterbi算法進(jìn)行的CRF解碼。Dorado對(duì)每一個(gè)測(cè)序儀器構(gòu)建Fast、Hac 和 Sup 三種模型，模型依次增大且得到的堿基結(jié)果逐漸準(zhǔn)確。r10.4.1模型中，在Fast 模式下，解碼器進(jìn)入到線性CRF層的特征向量為96，狀態(tài)長(zhǎng)度為3；Hac模式為特征向量128，狀態(tài)長(zhǎng)度為4，而Sup模式特征向量達(dá)到了256，狀態(tài)長(zhǎng)度為5。在 Fast模型下GPU加速效果可以達(dá)到60倍以上。

圖 5 Dorado r10.4.1 Fast模式的模型示意圖

2.2

DeepConsensus深度學(xué)習(xí)修正測(cè)序錯(cuò)誤

由于PacBio會(huì)有完全隨機(jī)的錯(cuò)誤發(fā)生（主要為短片段的插入與缺失），PacBio使用同一個(gè)零模波導(dǎo)孔內(nèi)SMRT 圈多次測(cè)序得到多條Subreads序列的方法去除測(cè)序錯(cuò)誤，矯正得到一個(gè)低錯(cuò)誤率高質(zhì)量值的序列，平均SMRT 圈的測(cè)序圈數(shù)決定了測(cè)序的時(shí)間。PacBio Revio測(cè)序儀使用CCS 和DeepConsensus 模型來(lái)矯正測(cè)序錯(cuò)誤獲得高質(zhì)量的序列。CCS基于隱馬爾可夫模型，其GPU版本可加速至12倍。DeepConsensus是一個(gè)編碼器-only的Transformer模型，與使用CCS的結(jié)果相比，DeepConsensus能夠得到更準(zhǔn)確的多聚核苷酸位點(diǎn)。

DeepConsensus將測(cè)序得到的Reads切割，每份 100bp，每一份向量包含了序列、鏈方向、脈沖寬度、脈沖間隔以及信噪比等信息。通過(guò)一個(gè)6個(gè)自注意力層的編碼器后，使用以Softmax作為激活函數(shù)的前饋神經(jīng)網(wǎng)絡(luò)解碼得到真實(shí)的序列。由于GPU在神經(jīng)網(wǎng)絡(luò)訓(xùn)練和推理上的優(yōu)勢(shì)，相比于CPU版本的DeepConsensus，GPU加速效果能夠達(dá)到 3.3倍以上。

圖 6 DeepConsensus 流程示意圖

圖片來(lái)源于Baid et al. DeepConsensus improves the accuracy of sequences with a gap-aware sequence transformer. DOI: 10.1038/s41587-022-01435-7

三

測(cè)序序列比對(duì)

兩條序列比對(duì)問(wèn)題是序列比對(duì)的基礎(chǔ)，以圖7所示，序列1 GACTTTCC 與 序列2 TAGACTACC 進(jìn)行比對(duì)得到比對(duì)表。兩條DNA序列的比對(duì)通常有四種狀態(tài)：插入、缺失、匹配和錯(cuò)配，而匹配和錯(cuò)配往往會(huì)在比對(duì)問(wèn)題中合并為同一個(gè)狀態(tài)。

圖 7 seq1與seq2的序列比對(duì)示意圖

Smith-Waterman算法和 Needleman-Wunsch算法是序列比對(duì)中常用的算法，其中Smith-Waterman算法用于局部比對(duì)而Needleman-Wunsch算法應(yīng)用于全局序列比對(duì)。而實(shí)際案例中，由于參考基因組的序列長(zhǎng)，比如人的參考序列總長(zhǎng)度在3Gb，一次測(cè)序深度10x測(cè)序所得Reads數(shù)目在2百萬(wàn)左右，直接使用局部比對(duì)的Smith-Waterman算法在時(shí)間成本上是不可接受的，于是有兩種加速方案：第一種為各測(cè)序所得序列間并行運(yùn)算，另一種則是對(duì)現(xiàn)有算法進(jìn)行改進(jìn)，Minimap2等基Needleman-Wunsch的比對(duì)軟件便應(yīng)運(yùn)而生。在本節(jié)中首先介紹Smith-Waterman與Needleman-Wunsch的動(dòng)態(tài)規(guī)劃算法與其GPU加速情況，然后介紹Minimap2特殊的GPU加速。

3.1

Smith-Waterman與Needleman-Wunsch 算法與GPU加速

兩個(gè)算法擁有相同的步驟：首先初始化得分矩陣，根據(jù)相似性矩陣和罰分矩陣填滿表格，后根據(jù)分?jǐn)?shù)回溯得到比對(duì)結(jié)果。如圖 8所示，Gotoh提出的比對(duì)方式具體如下：

01

確定罰分矩陣，罰分矩陣定義了堿基匹配的獎(jiǎng)勵(lì)分s、開(kāi)啟間隔區(qū)域（插入或者缺失）的懲罰分q和持續(xù)間隔e的懲罰分。這一步往往根據(jù)測(cè)序的特性進(jìn)行確定。

02

根據(jù)罰分矩陣來(lái)初始化打分矩陣，通常會(huì)將分?jǐn)?shù)存儲(chǔ)在三個(gè)打分矩陣 H、E和F中，分別代表著匹配分?jǐn)?shù)、插入分?jǐn)?shù)和缺失分?jǐn)?shù)。Smith-Waterman算法中，打分矩陣H、E 和F的首列與首行初始化為0；而Needleman-Wunsch算法中，打分矩陣F首列初始化公式為F0,j=F0,j-1-e，打分矩陣E首行初始化公式為Ei,0=Ei-1,0-e，打分矩陣H的首列和首行初始化為Hi,j=max{Ei,j,Fi,j}

03

計(jì)算打分矩陣H算法中的每一個(gè)位置的分?jǐn)?shù)Hi,j 。Smith-Waterman 算法中Hi,j=max{Hi-1,j-1+s,Ei,j,Fi,j,0} ，而Needleman-Wunsch算法中，Hi,j=max?{Hi-1,j-1+s,Ei,j,Fi,j}，其中 Ei,j=max?{Hi-1,j-q,Ei-1,j}-e,Fi,j=max?{Hi,j-1-q,Fi,j-1}-e。

當(dāng)打分矩陣H、E和F全部填寫(xiě)完成后，需要回溯來(lái)獲取最優(yōu)的比對(duì)序列。Smith-Waterman算法從打分矩陣中的最大值開(kāi)始回溯直到Hi,j=0為止；而Needleman-Wunsch算法從打分矩陣H右下角開(kāi)始回溯，直到左上角為止?；厮萋窂郊礊楸葘?duì)路徑。此時(shí)對(duì)于長(zhǎng)度M和長(zhǎng)度N的序列對(duì)而言，時(shí)間復(fù)雜度為O(MN)。

Myers提出了Bit-Vector加速該算法，打分矩陣Hi,j中的依賴(lài)于Hi-1,j-1、Hi-1,j和Hi,j-1，如果以列向量的角度看的話，如圖 9所示，那么Hj僅依賴(lài)于Hj-1。該算法將計(jì)算打分矩陣Hi,j變?yōu)橛?jì)算行差值Δhi,j與列差值Δvi,j。將時(shí)間復(fù)雜度降為O(N)。

圖 9 Bit-Vector 加速示意圖

Siriwardena等人提出了另一個(gè)加速方案，如圖 10所示，次對(duì)角線上的各個(gè)元素沒(méi)有相互依賴(lài)，故拓展到次對(duì)角線的Block沒(méi)有依賴(lài)，因此通過(guò)各個(gè)Block之間的并行計(jì)算進(jìn)行GPU加速。該方案運(yùn)行時(shí)間與比對(duì)序列長(zhǎng)度成非線性正相關(guān)，在如人染色體的數(shù)量級(jí)上能夠達(dá)到40倍的加速效果。

圖 10 次對(duì)角線分塊加速示意圖

與CPU 和 GPU執(zhí)行時(shí)間比較圖

除了上述提到的方案以外，使用動(dòng)態(tài)規(guī)劃算法的時(shí)候需要頻繁用到max、min等計(jì)算，GPU同時(shí)采用特殊的硬件指令來(lái)加速動(dòng)態(tài)規(guī)劃算法，預(yù)計(jì)可以加速7倍。

3.2

比對(duì)軟件的新星Minimap2的流程與GPU加速

為了解決Smith-Waterman局部比對(duì)在大規(guī)模數(shù)據(jù)中耗時(shí)長(zhǎng)的問(wèn)題，Li開(kāi)發(fā)了Minimap2 軟件極大減少了比對(duì)的時(shí)間。Minimap2 是在比對(duì)領(lǐng)域較新且廣泛使用的軟件，相比于其它軟件如BWA-mem、STAR等有特定使用場(chǎng)景的限制，Minimap2的使用范圍比較廣泛，其適用于二代的DNA短讀長(zhǎng)測(cè)序數(shù)據(jù)以及三代DNA長(zhǎng)讀長(zhǎng)測(cè)序數(shù)據(jù)，也可以被應(yīng)用RNA-seq這種有大間隔的數(shù)據(jù)。如圖11所示，Minimap2將測(cè)序得到的 Reads與參考基因組的比對(duì)的過(guò)程主要分為 Seeding、Chaining 和 Alignment 三個(gè)步驟。

圖 11 Minimap2 流程示意圖

01

Seeding這一步快速定位參考基因組與Read之間完全匹配，無(wú)插入與缺失的短序列，得到若干有匹配位置信息的錨點(diǎn)（Anchors）。錨點(diǎn)使用A(x,y,w)表示，意味著參考基因組的坐標(biāo)[x-w+1, x]與Read的坐標(biāo)[y-w+1, y]完全匹配，每一條Read會(huì)得到若干個(gè)錨點(diǎn)。

02

Chaining這一步則是將Seeding 這步得到的錨點(diǎn)連起來(lái)，確定這條Read在參考基因組上的真實(shí)位置。在 Minimap2中，Chaining是用一個(gè)一維的動(dòng)態(tài)規(guī)劃算法計(jì)算f 解決的，將所有錨點(diǎn)根據(jù)參考基因組的位置坐標(biāo)x進(jìn)行排序，第 i個(gè)anchor的分?jǐn)?shù)f(i)可以由公式 f(i)=max?{max?{f(j)+α(j,i)-β(j,i),wi}(i>j≥1)}, 計(jì)算得到，其中 α(j,i)=min?{min??{yi-yj,xi-xj},wi} 作為獎(jiǎng)勵(lì)分。為兩個(gè)錨點(diǎn)之間相同的長(zhǎng)度；β(j,i)=γc ((yi-yj )-(xi-xj ))，作為兩個(gè)錨點(diǎn)之間發(fā)生插入或者缺失間隔的罰分，γc 是一個(gè)與間隔長(zhǎng)度相關(guān)的函數(shù)，可以調(diào)整 γc 以適應(yīng)不同測(cè)序類(lèi)型。Minimap2也會(huì)設(shè)置一個(gè)最長(zhǎng)間隔G，當(dāng)max?{yi-yj,xi-xj}>G時(shí)，罰分β(j,i)=∞。若有N個(gè)錨點(diǎn)，那么直接計(jì)算最大f時(shí)間復(fù)雜度為O(N2)，基于大部分f(i) 的最大值在i附近，故只會(huì)計(jì)算h次，于是Minimap2將時(shí)間控制在O(hN)內(nèi)。對(duì)于i每次得到f后，會(huì)回溯到 j 直到f(j)=wj ，此時(shí)將回溯過(guò)程中的錨點(diǎn)記為“已用”，得到一條鏈，每一個(gè)鏈根據(jù)匹配長(zhǎng)度和間隔長(zhǎng)度計(jì)算比對(duì)分?jǐn)?shù)S。當(dāng)我們將所有的錨點(diǎn)標(biāo)記為“已用”后，可以得到所有鏈和其分?jǐn)?shù)S，分?jǐn)?shù)最高的鏈便是得到的最優(yōu)鏈。

圖 12 Minimap2不同步驟所花時(shí)間比例圖。圖片來(lái)源于Sadasivan et al. Accelerating Minimap2 for Accurate Long Read Alignment on GPUs. DOI: 10.26502/jbb.2642-

03

Alignment 這一步，則是將測(cè)序得到的 Read 與 Chaining 這一步得到的鏈序列進(jìn)行 base to base 的全局比對(duì)。這一步則可以用上一小節(jié)所提到的Needleman-Wunsch算法完成。

Minimap2的耗時(shí)時(shí)間如圖 12所示，Chaining和 Alignment 是耗時(shí)最長(zhǎng)的時(shí)間。因此提升 Chaining 和 Alignment步驟的速度是提升 Miniamp2 運(yùn)行速度的關(guān)鍵。

Sadasivan 等人提出了Minimap2 Chaining步驟的并行方案。如圖 13所示，他們首先將 Minimap2 的Chaining步驟的逆向搜索變成順向搜索，并且對(duì)于每一個(gè)i錨點(diǎn)并行計(jì)算其順向的 >i 錨點(diǎn)是否為最大分?jǐn)?shù)。GPU通過(guò)這種方式并行即可以將Minimap2提速3.8倍左右。

圖 13 mm2-ax 加速 Minimap2 Chaining 的方案

圖片來(lái)源于Sadasivan et al. Accelerating Minimap2 for Accurate Long Read Alignment on GPUs. DOI: 10.26502/jbb.2642-91280067

四

變異檢測(cè)

4.1

SNV與InDel的檢測(cè)

在生物學(xué)上，通常需要關(guān)注的是測(cè)序得到的變異。單堿基變異（SNV）是常見(jiàn)的一種變異方式，SNV 是單堿基的變異，堿基的突變有時(shí)候會(huì)改變氨基酸，或者使得轉(zhuǎn)錄提前終止。SNV的檢測(cè)主要分為兩個(gè)步驟：1. 確定并過(guò)濾變異位置； 2. 確定基因型。在本節(jié)中將介紹兩種軟件，一種是以 HaplotypeCaller和MuTect2為代表的基于傳統(tǒng)概率學(xué)模型建模的軟件，另一種則是以DeepVariant為代表的利用深度學(xué)習(xí)進(jìn)行檢測(cè)的軟件。

4.1.1

概率模型 HaplotypeCaller / MuTect2檢測(cè)流程

HaplotypeCaller 與 MuTect2 是 GATK 中重要的 SNV 與 InDel Calling 的軟件。兩款軟件在一定程度上擁有相同的原理，在概率模型上的區(qū)別，導(dǎo)致HaplotypeCaller 主要用于性細(xì)胞變異的檢測(cè)，而MuTect2 更適用于體細(xì)胞變異的檢測(cè)。如圖 14兩款軟件都在確定變異區(qū)間后，利用局部組裝得到的單倍型并將測(cè)序得到的Reads用Smith-Waterman算法比對(duì)到單倍型中來(lái)減少測(cè)序、比對(duì)等錯(cuò)誤信息，通過(guò) PairHMM 來(lái)得到每一個(gè)單倍型的似然值，然后根據(jù)各自的概率模型獲得最終的基因型。

圖 14 HaplotypeCaller與MuTect2 的通用流程示意圖

每一個(gè)單倍型的似然值時(shí)HaplotypeCaller和MuTect2的關(guān)鍵步驟，即求單倍型的集合H對(duì)于所有測(cè)序得到的Reads的集合R的比對(duì)結(jié)果集合A的概率，即 P(R│H)=∑AP(R,A|H)，那么只需要求解每一個(gè)比對(duì)的概率即可。其中我們已經(jīng)看到了重新比對(duì)中Smith-Waterman的GPU加速方案，HaplotypeCaller中的PairHMM這一步亦可以用GPU并行加速。

PairHMM 的過(guò)程與前述的比對(duì)算法相似，定義了三個(gè)狀態(tài): 匹配M, 插入I和缺失D。如圖15所示，三個(gè)狀態(tài)會(huì)發(fā)生轉(zhuǎn)移, M→I或者 M→D，即從匹配狀態(tài)轉(zhuǎn)化為插入狀態(tài)或者缺失狀態(tài)概率為δ、 I→M 從插入狀態(tài)結(jié)束間隔變回匹配狀態(tài)概率為ε、D→M從缺失狀態(tài)結(jié)束間隔變回匹配狀態(tài)概率為ε，也會(huì)發(fā)生持續(xù) M→M 即延續(xù)匹配狀態(tài)概率為1-2δ、I→I 延續(xù)插入狀態(tài)概率為1-ε和D→D 延續(xù)缺失狀態(tài)，概率為1-ε。基于生物的序列特征，往往1-ε<δ?1-2δ，即維持間隔的概率小于發(fā)生間隔的概率，遠(yuǎn)小于匹配的概率。

圖 15 PairHMM 狀態(tài)轉(zhuǎn)化示意圖

計(jì)算時(shí)定義了三個(gè)狀態(tài)分別定義了一個(gè)矩陣 M、I和D。Mi,j表示測(cè)序序列i與參考基因組j位置匹配的概率，同理 Ii,j和Di,j分別代表著發(fā)生插入的概率和缺失的概率。矩陣M可以通過(guò)公式 Mi,j=S(Mi-1,j-1 (1-2δ)+Ii-1,j ε+Di,j-1 ε) 計(jì)算而得，其中S與堿基的質(zhì)量值錯(cuò)配概率等有關(guān)；矩陣I可以通過(guò)公式Ii,j=Mi-1,j δ+Ii-1,j (1-ε) 得到；矩陣D可以通過(guò)公式Di,j=Mi,j-1 δ+Di,j-1 (1-ε) 得到；當(dāng)計(jì)算完矩陣M、I后，將M與I的最后一列相加即可得到似然值。

PairHMM是檢測(cè)流程中最密集的部分，因此提升 PairHMM的速度就是提升 HaplotypeCaller/MuTect2 的關(guān)鍵，可以發(fā)現(xiàn)PairHMM的分?jǐn)?shù)矩陣計(jì)算過(guò)程與Smith-Waterman的打分矩陣過(guò)程相似，于是Li 等人參考Smith-Waterman的GPU加速方案提出的基于 GPU的 PairHMM 比CPU版本加速 783 倍。

4.1.2

DeepVariant深度學(xué)習(xí)檢測(cè)SNV

HaplotypeCaller/Mutect2 傳統(tǒng)的概率模型算法存在兩個(gè)問(wèn)題：

1

由于不同的測(cè)序技術(shù)或者物種不同，HaplotypeCaller所需要的往往先驗(yàn)參數(shù)不同。比如PacBio測(cè)序中，插入發(fā)生的概率大于缺失發(fā)生的概率，遠(yuǎn)大于傳統(tǒng)的 Illumina 二代的概率，于是 HaplotypeCaller 需要根據(jù)測(cè)序技術(shù)的更迭不斷嘗試獲得最優(yōu)的先驗(yàn)參數(shù)。

2

傳統(tǒng)概率模型的建模中假設(shè)了SNV在基因組發(fā)生的概率相同。而由于自然選擇的存在，SNV的發(fā)生概率在基因組中并非均勻。在生物中，往往位于內(nèi)含子區(qū)域的變異發(fā)生概率會(huì)大于外顯子區(qū)域概率，而管家基因內(nèi)變異發(fā)生概率會(huì)遠(yuǎn)小于其它基因變異發(fā)生概率。

為了解決這個(gè)問(wèn)題，DeepVariant與 HaplotypeCaller 不同，利用深度學(xué)習(xí)的方法代替?zhèn)鹘y(tǒng)的概率模型。如圖 16所示，在確定變異的位置后，DeepVariant將覆蓋該位置的 Reads 提取出來(lái)，組成一個(gè) PileUp 的圖像，圖像包含了該位點(diǎn)上下游100bp的序列、質(zhì)量值、Reads 鏈、比對(duì)質(zhì)量、變異支持及覆蓋情況等信息，這些信息將起到與 HaplotypeCaller 組裝單倍型的相似作用。

圖 16 DeepVariant MakeExample 所以生成的典型圖片，分別為純合SNV變異、雜合SNV變異和純合未變異。圖片來(lái)源于Looking Through DeepVariant’s Eyes. https://google.github.io/deepvariant/posts/2020-02-20-looking-through-deepvariants-eyes/

得到圖片后，DeepVariant 使用 CNN 將 PileUp 得到的圖片推理直接獲得該位點(diǎn)基因型的似然值，從而判斷是否發(fā)生變異以及變異的基因型。

圖 17 DeepVariant流程示意圖

圖片來(lái)源于Poplin et al. A universal SNP and small-indel variant caller using deep neural networks. DOI: 10.1038/nbt.4235

DeepVariant 在測(cè)序技術(shù)和物種都具有泛化性，用戶可以通過(guò)兩個(gè)個(gè)體的測(cè)序數(shù)據(jù)構(gòu)建純合和雜合的“銀標(biāo)準(zhǔn)”對(duì)DeepVariant的模型進(jìn)行微調(diào)，使DeepVariant適合不同的測(cè)序技術(shù)（如 Illumina、PacBio HiFi、 ONT R10.4.1等等）。另一方面，DeepVariant 也通過(guò)對(duì)輸入或者輸出改造，已經(jīng)衍生出應(yīng)用于其他場(chǎng)景的模型，比如適用于家系的 DeepTrio，或者利用腫瘤樣本-正常細(xì)胞樣本對(duì)檢測(cè)體細(xì)胞變異的 VarNet。如圖 18所示，DeepVariant 在 GPU 的幫助下可以顯著減少 Variants Calling 這個(gè)步驟。

圖 18 DeepVariant不同硬件不同環(huán)節(jié)時(shí)間統(tǒng)計(jì)圖。

圖片來(lái)源于 Huang et al. DeepVariant-on-Spark: Small-Scale Genome Analysis Using a Cloud-Based Computing Framework. DOI: 10.1155/2020/7231205

4.2

結(jié)構(gòu)變異檢測(cè)

結(jié)構(gòu)變異通常包括了插入、刪除、串聯(lián)重復(fù)、染色體倒位、染色體易位、拷貝數(shù)變異和復(fù)雜的復(fù)合結(jié)構(gòu)變異等變異。結(jié)構(gòu)變異涉及到>50bp的序列，與SNV或小片段的InDel相比，一旦變異發(fā)生通常更不可逆。結(jié)構(gòu)變異的檢測(cè)有著重要的意義：

1

由于不可逆的特性，結(jié)構(gòu)變異往往是人種區(qū)別、物種分化的關(guān)鍵。

2

稀有的結(jié)構(gòu)變異往往與疾病或者癌癥發(fā)生關(guān)聯(lián)。

3

在腫瘤細(xì)胞中變異會(huì)出現(xiàn)大量的高密度的異常的結(jié)構(gòu)變異，甚至?xí)霈F(xiàn)染色體碎裂的現(xiàn)象。

結(jié)構(gòu)變異的檢測(cè)的方法如圖 19所示，包括了基于雙端測(cè)序方向的Read Pair、基于測(cè)序深度變化的 Read Depth、基于比對(duì)時(shí)異常比對(duì)的Split Reads 和將測(cè)序得到的 Reads組裝與參考基因組比對(duì)的 Assembly。目前短序列有Manta、Delly，而長(zhǎng)序列有 Sniffles2、CuteSV等軟件，基于比對(duì)時(shí)異常比對(duì)的Split Reads是主流的方法。

圖 19 結(jié)構(gòu)變異檢測(cè)的常用方法

圖片來(lái)源于Alkan et al. Genome structural variation discovery and genotyping. DOI: 10.1038/nrg2958

4.2.1

SVision深度學(xué)習(xí)模型

如 Sniffles2、Manta等結(jié)構(gòu)變異軟件，過(guò)濾得到比對(duì)異常的Reads，將比對(duì)異常點(diǎn)根據(jù)其在參考基因組的位置、異常的長(zhǎng)度等方式進(jìn)行聚類(lèi)，每一個(gè)類(lèi)即獲得一個(gè)結(jié)構(gòu)變異位點(diǎn)。但是這缺少了發(fā)現(xiàn)復(fù)合結(jié)構(gòu)變異的能力。通常復(fù)合結(jié)構(gòu)變異檢測(cè)都依賴(lài)于現(xiàn)有的模式，如易位后缺失等，缺乏一個(gè)統(tǒng)一的結(jié)構(gòu)。

SVision 便是為了解決這個(gè)問(wèn)題而出現(xiàn)的模型31。如圖 20所示，SVision包含了編碼器、tMOR和解釋器。編碼器將測(cè)序的序列轉(zhuǎn)化為圖片，其識(shí)別支持變異序列和參考基因組序列的堿基通過(guò)VAR-REF和REF-REF圖像，從VAR-REF中去除REF-REF中已有的重復(fù)序列，以提高后續(xù)分析的準(zhǔn)確度。在tMOR步驟中，SVison將包含多個(gè)復(fù)合結(jié)構(gòu)變異斷點(diǎn)的圖像通過(guò)5層卷積層、1層展平層和2層全連接層的卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測(cè)發(fā)生變異的概率。最后，解釋器描述了不同的復(fù)合結(jié)構(gòu)變異。由于這是一個(gè)卷積神經(jīng)網(wǎng)絡(luò)的模型，GPU的矩陣乘法上的優(yōu)勢(shì)能夠比CPU版本的SVision快1.8倍。

圖 20 SVision 流程示意圖

圖片來(lái)源于 Lin et al. SVision: a deep learning approach to resolve complex structural variants. DOI: 10.1038/s41592-022-01609-w

四

總結(jié)與展望

變異檢測(cè)是定位功能基因或疾病基因的基礎(chǔ)。如表格 1所示，從測(cè)序數(shù)據(jù)的產(chǎn)生、預(yù)處理、基因組的比對(duì)、變異檢測(cè)都可以使用GPU加速來(lái)提升檢測(cè)的速度。高速分析能夠節(jié)省測(cè)序的時(shí)間成本，從而降低測(cè)序價(jià)格和分析成本。這可以讓研究人員更愿意進(jìn)行測(cè)序來(lái)建立與完善疾病腫瘤相關(guān)的數(shù)據(jù)庫(kù)，為將來(lái)攻克疾病與腫瘤打下基礎(chǔ)。

表格 1 文中所提的軟件GPU與CPU運(yùn)行時(shí)間比

可以看到，在目前基因測(cè)序領(lǐng)域，傳統(tǒng)軟件算法因?yàn)樵O(shè)備局限性，在面對(duì)復(fù)雜的生物系統(tǒng)時(shí)，選擇使用簡(jiǎn)化的概率模型建模來(lái)犧牲準(zhǔn)確度而獲得更快的速度。而隨著GPU算力的進(jìn)步，越來(lái)越多研究人員選擇深度學(xué)習(xí)或者大模型來(lái)得到更精確的結(jié)果。除了文中所提到的數(shù)據(jù)預(yù)處理、基因組比對(duì)和變異檢測(cè)步驟以外，在后續(xù)分析中，機(jī)器學(xué)習(xí)也起到了更重要的作用，

比如用線性模型建模關(guān)聯(lián)樣本表型與基因型的GWAS可以使用 DeepNull 得到更好的關(guān)聯(lián)結(jié)果。

我們預(yù)計(jì)，隨著機(jī)器學(xué)習(xí)算法和GPU設(shè)備的發(fā)展，快速而準(zhǔn)確的基因測(cè)序可以幫助研發(fā)人員設(shè)計(jì)靶向藥物，幫助普通民眾更早發(fā)現(xiàn)潛在風(fēng)險(xiǎn)，幫助醫(yī)生更準(zhǔn)確地對(duì)病情進(jìn)行診斷治療，精準(zhǔn)醫(yī)療未來(lái)可期。

名詞列表

責(zé)任編輯：彭菁