劉夢雅, 徐 云

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

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

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

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

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

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

1 相關工作

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

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

亞硫酸氫鹽測序技術通過對基因片段進行預處理,使得甲基化的胞嘧啶(C)保持不變, 未發(fā)生甲基化的C先轉(zhuǎn)換成尿嘧啶(U), 再轉(zhuǎn)換為胸腺嘧啶(T), 如圖1所示[10].因此, 在DNA甲基化序列比對的過程中, 序列中的T有可能比對到參考基因組上的T或C, 但反之不行, 導致比對的難度增加[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不對稱比對的特點, 適應序列長度短(40 bp-400 bp)、數(shù)量多、規(guī)模大的特性, 實現(xiàn)將甲基化序列快速比對到參考基因組上, 使得全基因組甲基化分析成為可能.但亞硫酸氫鹽預處理將未發(fā)生甲基化的C轉(zhuǎn)化成T, 在大部分序列比對中, 字符集從4字符集(A、T、C、G)變成了3字符集(A、T、G),降低了序列的復雜性, 增加了BS-reads唯一比對位置確定的難度, 同時使參考基因組中重復區(qū)域的甲基化狀態(tài)分析更為艱難.

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

導向定位測序是一種新的全基因組DNA甲基化檢測的方法.每條DNA鏈是由磷酸和脫氧核糖構(gòu)成,3’端和5’端表示DNA鏈的兩端, 其中連接磷酸基團的一端為5’端, 另一端是3’端.DNA的復制方向是從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比對到參考基因組的范圍也相應確定.其中Read2對Read1位置的確定起到定位作用, 為后續(xù)全基因組甲基化的分析奠定了基礎.

圖2 導向定位測序過程

現(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的不對稱性比對問題.為受亞硫酸氫鹽預處理影響較大的基因片段和部分物種, 提供了甲基化分析的新方法, 使得這部分序列甲基化信息的精確檢測成為可能.但使用Smith-Waterman算法對多個候選比對位置進行動態(tài)規(guī)劃驗證, 需要大量的時間.且未考慮僅允許Read1中的C比對到參考基因組中的C, 有可能導致Read1的錯誤比對, 從而影響到后續(xù)全基因組甲基化的分析.且現(xiàn)有GPS數(shù)據(jù)的唯一比對比例為79.8%-82.3%, 仍有提升的空間.

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

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

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

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

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

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

2.2 定位候選位置

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

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

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

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

2.3 過濾偏離區(qū)域

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

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

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

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

2.4 確定最佳位置

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

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

其中,R1的第一個堿基比對到參考基因組對應位置的概率為s1, 依次類推得第K個堿基比對到參考基因組對應位置的概率為sK.如圖6所示,s1-sK的計算使用工具BAM-ABS[20], 該工具使用貝葉斯模型, 以Multireads和參考基因組之間的錯配信息和對應甲基化區(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 評價指標

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

(1)唯一比對比率

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

(2)準確率

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

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)有方法相比, 準確率相差不大, 最多相差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é)論與展望

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