戴道成

（西安歐亞學(xué)院 金融學(xué)院，陜西 西安 710065）

0 引 言

眾所周知，桑格（Sanger）測(cè)序法是人類歷史上第一代DNA測(cè)序技術(shù)，因其成本高、通量低的缺點(diǎn)，而被以Roche公司的454技術(shù)、illumina公司的Solexa技術(shù)和ABI公司的Solid技術(shù)為代表的二代測(cè)序技術(shù)（Next Generation Sequencing, NGS）所取代。相較于Sanger測(cè)序，NGS測(cè)序的成本大大降低，同時(shí)實(shí)現(xiàn)了較高的測(cè)序精度，但read長(zhǎng)度的極大縮短限制了NGS的廣泛應(yīng)用。以Pacific Biosciences公司SMRT技術(shù)和Oxford Nanopore Technologies公司納米孔單分子技術(shù)為代表的三代測(cè)序技術(shù)（Third Generation Sequencing, TGS），不僅繼承了NGS的優(yōu)點(diǎn)，而且能產(chǎn)生長(zhǎng)度大于10 kbp的長(zhǎng)read，因而在序列組裝、基因突變鑒定以及疾病診斷等諸多領(lǐng)域得到廣泛的應(yīng)用。但是錯(cuò)誤隨機(jī)分布的特性和15%的測(cè)序錯(cuò)誤率是限制TGS大范圍應(yīng)用的主要瓶頸，由此產(chǎn)生兩大類針對(duì)三代測(cè)序數(shù)據(jù)進(jìn)行錯(cuò)誤校正的方法。第一類是三代測(cè)序數(shù)據(jù)（long reads, LRs）的自校正方法，以HGAP和LoRMA為例，這些方法需要較高水平的覆蓋度來(lái)確保質(zhì)量，這無(wú)疑是增加了基因組計(jì)劃的成本。第二類是混合校正方法，此方法需要同時(shí)比較LRs和二代測(cè)序數(shù)據(jù)（short reads, SRs），如proovread、LoRDEC和Jabba等，其經(jīng)常需要計(jì)算和處理數(shù)百萬(wàn)個(gè)SRs到LRs的比對(duì)結(jié)果，這也是一個(gè)十分艱巨的挑戰(zhàn)，需要消耗較多的資源和內(nèi)存。為了避免以上問(wèn)題，充分發(fā)揮已有糾錯(cuò)方法的優(yōu)點(diǎn)，基于機(jī)器學(xué)習(xí)中的眾包標(biāo)注思想，通過(guò)計(jì)算高質(zhì)量SRs集合的一致度、能力和準(zhǔn)確度來(lái)對(duì)低質(zhì)量LRs的位點(diǎn)進(jìn)行標(biāo)注，從而實(shí)現(xiàn)三代測(cè)序數(shù)據(jù)的精準(zhǔn)糾錯(cuò)。

1 方法設(shè)計(jì)

1.1 概述

在基于眾包標(biāo)注的三代測(cè)序數(shù)據(jù)的糾錯(cuò)方法（Error Correction method of Third generation sequencing Data based on Crowdsourcing annotation, CTDC）中，將通過(guò)基因組仿真得到的三代測(cè)序數(shù)據(jù)視為待標(biāo)注任務(wù)（Long Reads Assignment,LRA），將通過(guò)相同基因組仿真得到的二代測(cè)序數(shù)據(jù)視為標(biāo)注工作者（Short Reads Worker, SRW），利用SRW對(duì)LRA進(jìn)行標(biāo)注，從而完成三代測(cè)序數(shù)據(jù)的糾錯(cuò)過(guò)程。

在CTDC中，對(duì)LRA中的每個(gè)位點(diǎn)進(jìn)行眾包標(biāo)注的準(zhǔn)確性主要取決于參與標(biāo)注任務(wù)的SRW的能力，SRW能力越強(qiáng)，說(shuō)明其可信度越大，相應(yīng)位點(diǎn)的標(biāo)簽準(zhǔn)確率越高。而SRW的能力受其對(duì)應(yīng)LRA的位點(diǎn)任務(wù)難度以及比對(duì)質(zhì)量等的影響，所以為了從眾多標(biāo)簽中選出最準(zhǔn)確的一個(gè)，需要將所有可能影響SRW能力的因素都考慮進(jìn)去，包括LRA中每個(gè)位點(diǎn)的任務(wù)難度水平、SRW的準(zhǔn)確度以及覆蓋在此位點(diǎn)下的SRW之間的一致性。整個(gè)算法的設(shè)計(jì)過(guò)程如圖1所示。

圖1 眾包標(biāo)注的算法設(shè)計(jì)過(guò)程

1.2 SRs的數(shù)據(jù)預(yù)處理

通常來(lái)講，對(duì)于一個(gè)待標(biāo)注的LRs任務(wù)，初始采集到的來(lái)自同一樣本的SRs包括三類：完全正確的SRs、具有測(cè)序誤差噪聲的SRs、比對(duì)錯(cuò)誤而導(dǎo)致錯(cuò)位的SRs。因此，有必要對(duì)初始采集到的SRs進(jìn)行數(shù)據(jù)預(yù)處理，過(guò)濾掉相關(guān)噪聲數(shù)據(jù)，從而盡可能提高SRs的精度。

首先，由于二代測(cè)序數(shù)據(jù)的錯(cuò)誤通常位于read的末端，在使用BLASR進(jìn)行SRs到LRs的比對(duì)時(shí)將read尾部的5bp截掉，因此BLASR的參數(shù)為：blasr Sk Li --header --out result -m 5。

其次，保證SR與LR之間的相似度Similarity≥0.7，其中Similarity=numMatch/(numMatch + numSub + numIns +numDel)，numMatch、numSub、numIns、numDel分別代表正確比對(duì)的位點(diǎn)數(shù)目、錯(cuò)誤比對(duì)的位點(diǎn)數(shù)目、插入的位點(diǎn)數(shù)目和刪除的位點(diǎn)數(shù)目。

同時(shí)，由于BLASR是一種高容錯(cuò)性的比對(duì)工具，每條read在M5格式中會(huì)產(chǎn)生多條比對(duì)結(jié)果，并且這些結(jié)果已經(jīng)按照比對(duì)質(zhì)量score值進(jìn)行降序排序，因此選擇最上面的比對(duì)結(jié)果。

1.3 SRW對(duì)LRA的眾包標(biāo)注

1.3.1 LRA難度的計(jì)算

根據(jù)信息熵值衡量標(biāo)簽結(jié)果分布的平衡度，熵值越大，說(shuō)明標(biāo)簽越分散，平衡度越高，任務(wù)難度越大；反之，標(biāo)簽越集中，平衡度越低，任務(wù)難度越小。

1.3.2 SRW能力的計(jì)算

在眾包標(biāo)注中，需要考慮SRW的三個(gè)要素，自身能力、標(biāo)注一致度和標(biāo)注準(zhǔn)確度。下文分別進(jìn)行分析和計(jì)算。

1.3.2.1 SRW一致度

1.3.2.2 SRW準(zhǔn)確度

SRW準(zhǔn)確度通過(guò)其所標(biāo)注任務(wù)的位點(diǎn)中正確標(biāo)注的位點(diǎn)所占的比例來(lái)衡量。在SRW中，對(duì)于任意S，令A(yù)代表其準(zhǔn)確度，計(jì)算公式為：

（1）根據(jù)LRs的測(cè)序錯(cuò)誤率，為SRW中的每個(gè)S設(shè)置一個(gè)準(zhǔn)確度初始值0.85；

（2）對(duì)于第1位堿基，以S的準(zhǔn)確度初始值0.85計(jì)算相應(yīng)的支持率，完成眾包標(biāo)注；

（3）對(duì)于第k位堿基，以前k-1位更新過(guò)的準(zhǔn)確度值作為第k位的準(zhǔn)確度初值，計(jì)算位點(diǎn)對(duì)象的支持率，從而實(shí)現(xiàn)位點(diǎn)的校正。同時(shí)根據(jù)校正結(jié)果重新計(jì)算第k位的準(zhǔn)確度值，作為第k+1位的準(zhǔn)確度初值。

（4）重復(fù)以上步驟，直至所有堿基都得到處理。

1.3.2.3 SRW能力

基于眾包標(biāo)注的已有研究，對(duì)于任意S，令B代表其在某一堿基 上的能力，計(jì)算公式為：

其中，λ為一個(gè)調(diào)和參數(shù)，用作權(quán)衡SRW準(zhǔn)確度與一致度的權(quán)重比例，默認(rèn)取值為0.5。

1.4 算法集成

綜上，當(dāng)算法對(duì)最終標(biāo)簽進(jìn)行整合時(shí)，通常采用SRW的加權(quán)能力值進(jìn)行以上過(guò)程，并且具有更高能力值的工作者對(duì)標(biāo)簽選擇的貢獻(xiàn)權(quán)重更大。因此，對(duì)式（1）的支持率M進(jìn)行如下改進(jìn)。

最大期望算法的計(jì)算過(guò)程為：

（2）根據(jù)式（6）求得所有待標(biāo)注對(duì)象的各個(gè)標(biāo)簽的支持率，然后根據(jù)最大支持率更新相應(yīng)標(biāo)簽；

（3）根據(jù)式（3）更新S的一致度G，同時(shí)根據(jù)式（4）更新S的準(zhǔn)確度A，然后根據(jù)式（5）更新S的能力值B；

（5）迭代次數(shù)加1，如果迭代次數(shù)超過(guò)預(yù)設(shè)閾值，則算法轉(zhuǎn)到（6），否則，算法轉(zhuǎn)到（5）；

（6）重復(fù)以上步驟，直至LRA中的所有堿基都得到標(biāo)注。

2 實(shí)驗(yàn)分析

2.1 實(shí)驗(yàn)設(shè)計(jì)

2.1.1 數(shù)據(jù)說(shuō)明

考慮到真實(shí)數(shù)據(jù)的巨大和難處理性，CTDC算法的性能驗(yàn)證實(shí)驗(yàn)中均使用仿真數(shù)據(jù)。首先，將人類基因組的chr19染色體作為參考序列，然后通過(guò)二代仿真工具TNsim生成長(zhǎng)度為100 bp的SRs，并隨機(jī)植入0.1%的變異位點(diǎn)，通過(guò)三代仿真工具PBSIM生成長(zhǎng)度為15 000 bp、測(cè)序錯(cuò)誤率為15%的LRs。

2.1.2 實(shí)驗(yàn)說(shuō)明

本次實(shí)驗(yàn)主要驗(yàn)證CTDC算法對(duì)于三代測(cè)序數(shù)據(jù)的糾錯(cuò)能力。考慮到CTDC算法的糾錯(cuò)過(guò)程以及已有混合校正方法在測(cè)序覆蓋度、資源消耗和運(yùn)行內(nèi)存上的不足，設(shè)計(jì)以下三組實(shí)驗(yàn)進(jìn)行CTDC算法的性能驗(yàn)證。第一組實(shí)驗(yàn)：CTDC在SRs不同測(cè)序覆蓋度下的精確度分析；第二組實(shí)驗(yàn)：CTDC與LoRDEC在SRs不同測(cè)序覆蓋度下的糾錯(cuò)結(jié)果比較；第三組實(shí)驗(yàn)：CTDC、LoRDEC與proovread在SRs不同測(cè)序覆蓋度下的運(yùn)行時(shí)間和內(nèi)存比較。

2.2 評(píng)估指標(biāo)

TP：原位點(diǎn)有測(cè)序錯(cuò)誤，且通過(guò)CTDC算法進(jìn)行糾錯(cuò)；

TN：原位點(diǎn)無(wú)測(cè)序錯(cuò)誤，且通過(guò)CTDC算法未進(jìn)行糾錯(cuò)；

FP：原位點(diǎn)無(wú)測(cè)序錯(cuò)誤，且通過(guò)CTDC算法進(jìn)行糾錯(cuò)；

FN：原位點(diǎn)有測(cè)序錯(cuò)誤，且通過(guò)CTDC算法未進(jìn)行糾錯(cuò)；

Accuracy=(TP+TN)/(TP+TN+FN+FP)，即通過(guò)CTDC算法成功糾錯(cuò)的位點(diǎn)數(shù)占位點(diǎn)總數(shù)的比例，通常來(lái)說(shuō)，精確度越高，糾錯(cuò)效果越好。

2.3 結(jié)果分析

2.3.1 CTDC在SRs不同測(cè)序覆蓋度下的精確度分析

對(duì)于SRs，設(shè)置生成的測(cè)序覆蓋度分別為5×、10×、20×、30×、50×和100×，考慮到運(yùn)行時(shí)間，將LRs的測(cè)序覆蓋度設(shè)置為5×。經(jīng)過(guò)CTDC糾錯(cuò)，對(duì)LRs的每一個(gè)位點(diǎn)進(jìn)行統(tǒng)計(jì)分析，得到如圖2所示的結(jié)果。

圖2 SRs的不同測(cè)序覆蓋度下CTDC對(duì)于三代測(cè)序數(shù)據(jù)的校正精確度

在圖2中，隨著SRs測(cè)序覆蓋度的增加，越來(lái)越多的SRs參與眾包標(biāo)注，因此CTDC對(duì)于三代測(cè)序數(shù)據(jù)的校正性能越來(lái)越好。同時(shí)，當(dāng)SRs的測(cè)序覆蓋度達(dá)到50×以上時(shí)，CTDC的精確度達(dá)到最佳97.6%。

2.3.2 CTDC與LoRDEC在SRs不同測(cè)序覆蓋度下的糾錯(cuò)結(jié)果比較

基于以上實(shí)驗(yàn)分析，設(shè)置如上的SRs與LRs，同時(shí)選擇LoRDEC進(jìn)行比對(duì)實(shí)驗(yàn)。經(jīng)過(guò)糾錯(cuò)后的位點(diǎn)統(tǒng)計(jì)，得到如圖3所示的結(jié)果。

圖3 SRs的不同測(cè)序覆蓋度下CTDC與LoRDEC對(duì)于三代測(cè)序數(shù)據(jù)的校正結(jié)果

由圖3可知，隨著SRs測(cè)序覆蓋度的增加，CTDC和LoRDEC對(duì)于三代測(cè)序數(shù)據(jù)的校正性能逐漸增強(qiáng)，但是在測(cè)序覆蓋度相同時(shí)，CTDC的性能優(yōu)于LoRDEC。為了便于比較CTDC和LoRDEC，繪制圖4進(jìn)行分析。

圖4 SRs的不同測(cè)序覆蓋度下CTDC與LoRDEC對(duì)于三代測(cè)序數(shù)據(jù)的校正精確度

總而言之，對(duì)于三代測(cè)序數(shù)據(jù)的校正，CTDC的性能優(yōu)于LoRDEC。

2.3.3 CTDC、LoRDEC與proovread在SRs不同測(cè)序覆蓋度下的運(yùn)行時(shí)間和內(nèi)存比較

為了驗(yàn)證CTDC在資源消耗和運(yùn)行內(nèi)存上的性能，同樣設(shè)置如上的SRs與LRs，同時(shí)選擇LoRDEC和proovread進(jìn)行比對(duì)實(shí)驗(yàn)。經(jīng)過(guò)對(duì)三代測(cè)序數(shù)據(jù)的糾錯(cuò)，所統(tǒng)計(jì)的實(shí)驗(yàn)結(jié)果如表1所示。

表1 SRs的不同測(cè)序覆蓋度下CTDC、LoRDEC和Proovread的內(nèi)存消耗比較

由表1可知，在低于20×的測(cè)序覆蓋度下，對(duì)于三代測(cè)序數(shù)據(jù)的校正，CTDC的運(yùn)行時(shí)間顯著低于LoRDEC和proovread。隨著SRs測(cè)序覆蓋度的增加，LoRDEC與proovread的運(yùn)行時(shí)間也明顯增加，但均低于CTDC。同時(shí)，在運(yùn)行內(nèi)存上，CTDC由于計(jì)算上的優(yōu)勢(shì)，資源消耗明顯低于LoRDEC和proovread。

3 結(jié) 論

針對(duì)測(cè)序數(shù)據(jù)進(jìn)行有效的校正是獲得高精度的基因序列的關(guān)鍵技術(shù)。針對(duì)現(xiàn)有三代測(cè)序數(shù)據(jù)糾錯(cuò)方法的不足，本文提出一種基于眾包標(biāo)注的三代測(cè)序數(shù)據(jù)糾錯(cuò)方法CTDC。CTDC使用經(jīng)處理的高質(zhì)量SRs對(duì)LRs的任務(wù)難度進(jìn)行計(jì)算，結(jié)合其一致度、能力和準(zhǔn)確度對(duì)當(dāng)前低質(zhì)量的各個(gè)位點(diǎn)對(duì)象分別進(jìn)行眾包標(biāo)注，并使用最大期望算法完成糾錯(cuò)。根據(jù)不同測(cè)序覆蓋度下的實(shí)驗(yàn)結(jié)果，當(dāng)SRs的測(cè)序覆蓋度為50×?xí)r，CTDC對(duì)于三代測(cè)序數(shù)據(jù)的校正精確度達(dá)到97.6%。同時(shí)，在運(yùn)行時(shí)間和占用內(nèi)存上，CTDC的性能明顯優(yōu)于已知的混合校正算法LORDEC和proovread。后續(xù)需要在眾包標(biāo)注算法中考慮LRs的重疊片段，以進(jìn)一步提高三代測(cè)序數(shù)據(jù)的糾錯(cuò)性能，從而獲取更高精度的測(cè)序數(shù)據(jù)，實(shí)現(xiàn)后續(xù)的基因分析乃至疾病的精準(zhǔn)治療。