劉夢雅, 徐 云

1(中國科學(xué)技術(shù)大學(xué) 計算機科學(xué)與技術(shù)學(xué)院, 合肥 230027)

2(安徽省高性能計算重點實驗室, 合肥 230026)

DNA甲基化是指在DNA甲基化轉(zhuǎn)移酶的作用下, 將甲基化基團選擇性的添加到胞嘧啶(C)的過程.因此, 在人類基因組中一部分C被甲基化, 另一部分C未被甲基化, 在未改變基因序列的前提下, 控制基因表達[1,2].大量研究表明, 基因組中甲基化C的比例和所在區(qū)域, 能夠為疾病的預(yù)測提供幫助, 同時也被證實在包括癌癥在內(nèi)的諸多疾病的治療中發(fā)揮著重要的作用[3-6].

圍繞全基因組甲基化的分析是近年研究的熱點內(nèi)容, 其中最重要的一步是將測序所產(chǎn)生的序列, 比對到參考基因組上, 獲取整個基因組的甲基化狀態(tài)[7].目前常用的測序技術(shù)是亞硫酸氫鹽的全基因組甲基化測序,但由于此類測序技術(shù)需要用亞硫酸氫鹽對原始DNA片段進行預(yù)處理, 降低了序列的復(fù)雜性, 增加了后續(xù)比對的難度.2019年出現(xiàn)的導(dǎo)向定位測序數(shù)據(jù)很好地解決這一問題, 其利用雙端測序的優(yōu)勢: 一端是亞硫酸氫鹽處理后的序列Read1, 另一端是原始序列Read2, 沒有經(jīng)過處理的原始序列更容易確定在參考基因組上的位置, 通過雙端測序序列的位置關(guān)系, 實現(xiàn)對甲基化的精確檢測[8].

然而, 現(xiàn)有導(dǎo)向定位測序數(shù)據(jù)(GPS)的比對方法先確定原始序列R2的20個候選比對位置, 時間消耗大; 之后再用動態(tài)規(guī)劃算法確定甲基化序列Read1的比對位置, 算法本身的時間成本高, 且需對多個候選位置進行動態(tài)規(guī)劃驗證.同時, 根據(jù)Read2確定Read1的比對位置過于絕對, 可能會產(chǎn)生誤判.現(xiàn)有亞硫酸氫鹽測序(BS)中的比對方法能將70%-90%的序列確定到唯一的位置, 比對的準(zhǔn)確率高達99%, GPS數(shù)據(jù)的現(xiàn)有比對方法, 相比之下仍有較大改進空間[9].

因此, 本文提出一種新的導(dǎo)向定位測序數(shù)據(jù)的比對算法.由于亞硫酸氫鹽序列比對精度高達99%, 對于能確定唯一位置的甲基化序列不再用常規(guī)序列進行定位, 保證高精度的同時節(jié)約了時間.首先確定導(dǎo)向定位測序數(shù)據(jù)中的甲基化序列的候選比對位置; 然后根據(jù)甲基化序列和常規(guī)序列在參考基因組上對應(yīng)的位置關(guān)系過濾偏離區(qū)域; 最后使用唯一比對序列的信息確定最佳比對位置.充分利用輔助信息, 實現(xiàn)以時間高效的方式將更多的甲基化序列比對到參考基因組上.

1 相關(guān)工作

目前, 對DNA甲基化進行檢測的金標(biāo)準(zhǔn)是亞硫酸氫鹽測序的全基因組甲基化測序, 隨著導(dǎo)向定位測序數(shù)據(jù)的出現(xiàn), 在實現(xiàn)對全基因組甲基化位點高度覆蓋的同時, 帶來了新的研究問題.接下來根據(jù)全基因組DNA甲基化測序數(shù)據(jù)的類型, 分別介紹數(shù)據(jù)的特點和相應(yīng)比對方法, 分析其優(yōu)缺點.

1.1 亞硫酸氫鹽測序(BS)及其比對方法

亞硫酸氫鹽測序技術(shù)通過對基因片段進行預(yù)處理,使得甲基化的胞嘧啶(C)保持不變, 未發(fā)生甲基化的C先轉(zhuǎn)換成尿嘧啶(U), 再轉(zhuǎn)換為胸腺嘧啶(T), 如圖1所示[10].因此, 在DNA甲基化序列比對的過程中, 序列中的T有可能比對到參考基因組上的T或C, 但反之不行, 導(dǎo)致比對的難度增加[11].這是甲基化序列比對, 同常規(guī)DNA序列比對的不同之處.測序得到的基因序列, 稱為BS-reads.甲基化分析中很重要的一步就是將BS-reads比對到參考基因組上, 確定其位置.

圖1 亞硫酸氫鹽測序過程

現(xiàn)有亞硫酸氫鹽測序的比對方法分成兩種, 分別是基于三字符集和基于通配符的比對方法.基于三字符集方法的特性, 是把BS-reads和參考基因組中的C都轉(zhuǎn)化成T, 將問題轉(zhuǎn)化成常規(guī)的DNA序列比對,在候選位置確定后, 再根據(jù)轉(zhuǎn)化前的相似性對候選位置進行過濾, 代表方法有Bismark[12], GEMBS[13], BSSeeker3[14], BatMeth2[15].基于通配符方法的特性是BS-reads中的C轉(zhuǎn)化成一個通配符, 同時允許通配符比對到參考基因組上的C和T, 代表方法有BSMAP[16]和RMAP[17].

隨著越來越多的甲基化數(shù)據(jù)被測出, 這兩類比對方法針對BS-reads不對稱比對的特點, 適應(yīng)序列長度短(40 bp-400 bp)、數(shù)量多、規(guī)模大的特性, 實現(xiàn)將甲基化序列快速比對到參考基因組上, 使得全基因組甲基化分析成為可能.但亞硫酸氫鹽預(yù)處理將未發(fā)生甲基化的C轉(zhuǎn)化成T, 在大部分序列比對中, 字符集從4字符集(A、T、C、G)變成了3字符集(A、T、G),降低了序列的復(fù)雜性, 增加了BS-reads唯一比對位置確定的難度, 同時使參考基因組中重復(fù)區(qū)域的甲基化狀態(tài)分析更為艱難.

1.2 導(dǎo)向定位測序(GPS)及其比對方法

導(dǎo)向定位測序是一種新的全基因組DNA甲基化檢測的方法.每條DNA鏈?zhǔn)怯闪姿岷兔撗鹾颂菢?gòu)成,3’端和5’端表示DNA鏈的兩端, 其中連接磷酸基團的一端為5’端, 另一端是3’端.DNA的復(fù)制方向是從5’端到3’端.測序數(shù)據(jù)中3’端的序列保持不變, 5’端的未甲基化的C轉(zhuǎn)化成T, 甲基化的C保持不變[8].獲得的兩條DNA序列(Read1和Read2), 其中Read1中未甲基化的C轉(zhuǎn)化成T, 和亞硫酸氫鹽測序方法處理后的序列特性一致; Read2是原始DNA序列, 更容易比對到參考基因組上, 如圖2所示.在Read2比對到參考基因組之后, Read1比對到參考基因組的范圍也相應(yīng)確定.其中Read2對Read1位置的確定起到定位作用, 為后續(xù)全基因組甲基化的分析奠定了基礎(chǔ).

圖2 導(dǎo)向定位測序過程

現(xiàn)有方法調(diào)用bowtie2[18]將Read2比對到參考基因組上, 獲取Read2在參考基因組上的20個候選比對位置.由于Illumina測序原理可得, Read1位于Read2下游的相反鏈上, 且由于GPS測序庫中的碎片大小是400 bp-500 bp, 可以確定Read2的比對范圍.通過使用Smith-Waterman算法[19], 將Read1比對到Read2所在參考基因組下游1 kb的反鏈上, 獲取Read1和參考基因組上局部相似性最高的位置.

新型測序數(shù)據(jù)的出現(xiàn), 實現(xiàn)以較低的覆蓋率(5X)獲得甲基化序列, 降低了甲基化測序的成本, 檢測甲基化沒有序列偏好.同時, 比對過程中使用Smith-Waterman算法, 允許Read1中的T比對到參考基因組上的T或C, 以解決甲基化序列比對中C/T的不對稱性比對問題.為受亞硫酸氫鹽預(yù)處理影響較大的基因片段和部分物種, 提供了甲基化分析的新方法, 使得這部分序列甲基化信息的精確檢測成為可能.但使用Smith-Waterman算法對多個候選比對位置進行動態(tài)規(guī)劃驗證, 需要大量的時間.且未考慮僅允許Read1中的C比對到參考基因組中的C, 有可能導(dǎo)致Read1的錯誤比對, 從而影響到后續(xù)全基因組甲基化的分析.且現(xiàn)有GPS數(shù)據(jù)的唯一比對比例為79.8%-82.3%, 仍有提升的空間.

2 比對算法設(shè)計和優(yōu)化

本文首先將導(dǎo)向定位數(shù)據(jù)中的甲基化序列比對到參考基因組上, 隨后利用和常規(guī)序列間的位置關(guān)系對候選位置進行過濾, 最后對仍不能確定位置的甲基化序列, 利用唯一比對位置的信息進行定位, 該方法主要包括4個步驟: (1)數(shù)據(jù)預(yù)處理; (2)定位候選位置;(3)過濾偏離區(qū)域; (4)確定最佳位置.

2.1 數(shù)據(jù)預(yù)處理

由于GPS庫的建立, 需要用到T4 DNA聚合酶處理基因片段, 從而保證Read2中的序列和原始DNA片段一致, 最后獲取雙端測序序列(Read1和Read2).但T4 DNA聚合酶可能產(chǎn)生處理不足或過度處理的現(xiàn)象,直接影響獲取數(shù)據(jù)的準(zhǔn)確性, 影響比對的效率.所以,需要找到Read1和Read2處理的邊界, 進而對數(shù)據(jù)進行預(yù)處理[8].

參考基因組中CH的甲基化水平較低, 若序列中出現(xiàn)CH, 則說明酶處理充分.Read2位于參考基因組的反鏈上, 根據(jù)堿基互補配對原理, 可知CH在Read2上的表現(xiàn)形式是[A/G/T]G.通過尋找[A/G/T]G確定酶處理邊界, 對Read2進行預(yù)處理.如圖3所示, 最靠近右端, 且滿足要求的處理邊界是TG.確定處理邊界后,保留邊界右邊的序列作為處理后的Read2序列.

圖3 數(shù)據(jù)處理示意圖

2.2 定位候選位置

本文直接將甲基化序列比對到參考基因組上.一方面由于亞硫酸氫鹽測序序列比對準(zhǔn)確率較高, 另一方面易比對到多個位置的比例約為20%-30%, 直接比對甲基化序列在保證準(zhǔn)確率的同時, 減少了后續(xù)的計算成本.本文使用基于三字符集方法和種子擴展策略的亞硫酸氫鹽比對工具BitmapperBS[9]進行修改, 其包含高效的數(shù)據(jù)結(jié)構(gòu)FM-tree, 針對數(shù)據(jù)三字符集特性對傳統(tǒng)FM-index索引進行優(yōu)化, 能夠獲得高達99.36%的準(zhǔn)確率.

首先Read1比對到參考基因組后, 分成兩部分.如圖4所示, 將能夠確定唯一位置的序列稱為Unique Reads; 比對到多個位置的序列稱為Multireads, 這部分序列比對到參考基因組的多個相似度較高的位置, 或者比對到了參考基因組的重復(fù)區(qū)域.

圖4 唯一比對序列和多比對序列舉例

后續(xù)處理主要針對Multireads, 找到其至多20個候選比對位置.將Multireads符號化表示為集合M, 設(shè)R1為集合M中的一條序列, 候選比對位置的個數(shù)為n,其候選比對位置集合P(R1)表示為:

2.3 過濾偏離區(qū)域

針對GPS數(shù)據(jù)中的測序原理可得, Read1位于Read2下游的相反鏈上, 且距離相差不大于1000 bp.通過利用Read2的位置信息作為輔助信息, 對Read1的候選比對位置進行限制, 過濾位于偏離區(qū)域的候選位置.

設(shè)與R1相對應(yīng)的另一端序列是R2, 首先使用bowtie2將R2比對到參考基因組上, 其候選比對的個數(shù)為m, 得到候選比對集合:

對Read1和Read2的候選位置進行兩兩比較, 過濾掉Read1候選比對集合中不能與Read2成對的位置.如圖5所示,R1的候選比對位置集合P(R1)中只有p1和p3存在與之相對應(yīng)的pos1和pos3, 所以對其余位置進行過濾, 此時P(R1)={p1,p3}.若此時R1的候選比對位置個數(shù)為1, 則轉(zhuǎn)化為Unique Reads, 否則其仍在Multireads的集合M中.

圖5 過濾偏離區(qū)域

2.4 確定最佳位置

過濾偏離區(qū)域后, 使用與Multireads重疊的Unique Reads信息, 以及對應(yīng)參考基因組之間的信息, 逐個堿基計算相應(yīng)位置的可能性, 最后對候選集合中每個位置得到一個總的得分, 確定最有可能的比對位置.

設(shè)甲基化序列R1的長度為K, 比對到候選比對位置的概率S為:

其中,R1的第一個堿基比對到參考基因組對應(yīng)位置的概率為s1, 依次類推得第K個堿基比對到參考基因組對應(yīng)位置的概率為sK.如圖6所示,s1-sK的計算使用工具BAM-ABS[20], 該工具使用貝葉斯模型, 以Multireads和參考基因組之間的錯配信息和對應(yīng)甲基化區(qū)域信息;以及重疊Unique Reads中獲得的SNP和甲基化區(qū)域信息作為先驗概率, 計算比對到每個位置的可能性.最后選取候選比對集合中得分最高的位置為最佳比對位置.

圖6 找最佳位置的過程

3 實驗分析

3.1 數(shù)據(jù)集和實驗環(huán)境

本文分別在真實和模擬數(shù)據(jù)集中對兩種方法進行比較, 真實數(shù)據(jù)集來自GSE92328, 在文獻[8]中提出并被證實有利于甲基化信息的分析.本文使用其中的GPS數(shù)據(jù)SRR6443657和SRR6443658進行后續(xù)實驗.模擬數(shù)據(jù)集使用模擬工具ART和Astair獲得, 先用ART生成常規(guī)DNA數(shù)據(jù), 再通過Astair對其中一條序列進行甲基化模擬.在未特殊聲明時, 本文使用工具的默認參數(shù)進行比較.

本文的實驗平臺包括2個14核Intel Xeon Gold 5120處理器和512 GB內(nèi)存, 操作系統(tǒng)為64位的Ubuntu 18.04.

3.2 評價指標(biāo)

分別使用時間、唯一比對比率和準(zhǔn)確率與現(xiàn)有方法進行比較.其中時間包括數(shù)據(jù)預(yù)處理和得到最終結(jié)果的時間, 建索引的時間不包括在內(nèi), 因為索引只需建造一次, 在后續(xù)實驗中通用.

(1)唯一比對比率

該評價指標(biāo)表示比對到唯一位置的甲基化序列占全部甲基化序列的比例.如式(4)所示,U表示唯一比對序列集合,n(U)表示唯一比對序列集合中序列的條數(shù),N表示全部甲基化序列的條數(shù).

(2)準(zhǔn)確率

準(zhǔn)確率這里表示唯一比對序列中, 比對到正確的位置所占的比例.如式(5)所示,n(R)是唯一比對集合U中比對到正確位置的序列個數(shù).模擬數(shù)據(jù)集中序列在參考基因組上的位置是已知的, 當(dāng)真實位置和比對結(jié)果相差200 bp以內(nèi), 則認為比對正確.真實數(shù)據(jù)集中序列在參考基因組上的位置是未知的, 故不進行準(zhǔn)確率的驗證.

3.3 實驗結(jié)果

分別使用模擬數(shù)據(jù)集和真實數(shù)據(jù)集探究本文方法和現(xiàn)有方法[8]的性能優(yōu)劣.模擬數(shù)據(jù)集和真實數(shù)據(jù)的序列長度均為100 bp.數(shù)據(jù)規(guī)模分別為甲基化序列1w條、5w條、10w條, 常規(guī)DNA序列1w條、5w條、10w條.

(1)模擬數(shù)據(jù)集實驗結(jié)果

如表1所示, 本文改進方法和現(xiàn)有方法相比, 準(zhǔn)確率相差不大, 最多相差0.7%.而本文方法獲得3-30倍時間性能的提升, 隨著數(shù)據(jù)規(guī)模的增大, 對時間性能的提升越明顯.同時本文方法獲得6%-10%唯一比對比率的提升, 將更多的序列比對到唯一位置, 有利于后續(xù)甲基化信息的分析.因模擬數(shù)據(jù)集不能完全模擬真實數(shù)據(jù)中插入、刪除, 以及發(fā)生測序錯誤、結(jié)構(gòu)變異的情況, 更容易比對到參考基因組上, 唯一比對比率相比真實數(shù)據(jù)更高.

表1 模擬數(shù)據(jù)集實驗結(jié)果

(2)真實數(shù)據(jù)集實驗結(jié)果

通過實驗探究了不同方法在運行時間、唯一比對比率方面的優(yōu)劣.如表2所示, 在SRR6443657數(shù)據(jù)集中, GPS方法的運行時間從4 min到38 min, 受數(shù)據(jù)規(guī)模影響較大; 本文的改進方法在這3種數(shù)據(jù)規(guī)模下運行時間相差不大, 為56-67 s, 小數(shù)據(jù)集時比現(xiàn)有方法快約3倍, 大數(shù)據(jù)集時比現(xiàn)有方法快30倍, 對大規(guī)模數(shù)據(jù)集的提升效果更為明顯.同時, GPS比對方法的唯一比對比率為79.32%-80.09%, 數(shù)據(jù)集規(guī)模對唯一比對比率的影響不大; 本文的改進方法唯一比對比率從85.37%到89.32%, 比之前方法提升了5%-10%, 且本文方法隨著數(shù)據(jù)集規(guī)模越大, 唯一比對比率越來越大,因獲取比對到唯一位置的序列信息越多, 更容易比對到唯一位置.第2個數(shù)據(jù)集整體結(jié)果和第1個數(shù)據(jù)集相似, 但唯一比對比率提升約為2%-6%, 較上一個數(shù)據(jù)集提升不明顯.實驗中發(fā)現(xiàn)部分甲基化序列未能找到與之配對的常規(guī)DNA序列, 使得該數(shù)據(jù)集比對難度增加.

表2 真實數(shù)據(jù)集實驗結(jié)果

4 結(jié)論與展望

本文提出了一種高效的導(dǎo)向定位測序數(shù)據(jù)的比對算法, 首先對數(shù)據(jù)進行預(yù)處理, 將甲基化序列定位到參考基因組上; 再利用雙端測序中兩端序列的位置關(guān)系,對甲基化序列的候選比對位置集合進行過濾; 最后通過比對到唯一位置的序列包含的信息, 找到最佳比對位置.實驗結(jié)果表明, 本文方法能夠加速比對過程, 將更多的甲基化序列比對到唯一位置, 且對大規(guī)模數(shù)據(jù)集的性能提升效果更為明顯.下一步的研究工作是提出啟發(fā)式的算法, 探究影響準(zhǔn)確率的因素, 在比對精度上取得更好的效果, 并探究比對性能的提升對后續(xù)甲基化信息的影響.