基于動(dòng)態(tài)規(guī)劃的DNA堿基識(shí)別峰匹配方法的設(shè)計(jì)與實(shí)現(xiàn)

賈二惠 李 曉2 張 濤 李 彬 趙麗華2 常海龍 金 川

(1.公安部第一研究所，北京 102200；2.太原理工大學(xué)數(shù)學(xué)院，太原 030024)

DNA測(cè)序技術(shù)自1980年發(fā)展至今取得了令人矚目的進(jìn)展，基于熒光誘導(dǎo)熒光毛細(xì)管電泳檢測(cè)的DNA測(cè)序技術(shù)因其高精度DAN序列識(shí)別準(zhǔn)確率而在個(gè)體識(shí)別、親緣鑒定、災(zāi)難身份識(shí)別、法醫(yī)DNA鑒定等領(lǐng)域依然發(fā)揮著無可替代的重要作用[1-3]。

鑒于當(dāng)前DNA測(cè)序儀器、DNA測(cè)序反應(yīng)與堿基合成等技術(shù)本身的限制, 儀器所采集的熒光光譜信號(hào)并不十分理想[4-8]。需對(duì)原始檢測(cè)信號(hào)進(jìn)行數(shù)據(jù)預(yù)處理、遷移率校正、DNA序列峰預(yù)測(cè)、Spacing校正、峰匹配等[9-13]，從而確定待檢DNA序列的堿基排序結(jié)果。但對(duì)一些具體環(huán)節(jié)所需的方法還缺少專門的深入分析，其中對(duì)峰匹配是將各堿基通道信號(hào)的識(shí)別峰與DNA序列的預(yù)測(cè)峰進(jìn)行匹配，它是堿基識(shí)別過程中的最復(fù)雜最關(guān)鍵部分，同時(shí)，峰匹配在下一步DNA序列準(zhǔn)確度質(zhì)量評(píng)估及后續(xù)DNA圖譜分析中起著至關(guān)重要的作用。

截至目前，對(duì)專門的峰匹配數(shù)據(jù)處理環(huán)節(jié)的研究文獻(xiàn)甚少，學(xué)術(shù)價(jià)值較高的當(dāng)屬文獻(xiàn)[11]。本研究通過綜合考慮DNA測(cè)序熒光光譜信號(hào)的實(shí)際特性，分析比較信號(hào)堿基峰、及各種偽峰的峰特征信息特點(diǎn)，根據(jù)預(yù)測(cè)峰的峰位置及峰周期、識(shí)別峰的峰位置及峰相對(duì)面積、峰型等特征信息，將動(dòng)態(tài)規(guī)劃的思想應(yīng)用到實(shí)際DNA測(cè)序數(shù)據(jù)處理中，在文獻(xiàn)[11]的基礎(chǔ)上設(shè)計(jì)了一種改進(jìn)的的峰匹配得分標(biāo)準(zhǔn)。當(dāng)峰匹配參數(shù)即峰限制閾值設(shè)置適中時(shí)，動(dòng)態(tài)規(guī)劃可以最大限度的將識(shí)別峰和預(yù)測(cè)峰進(jìn)行匹配，盡量做到不錯(cuò)配、不漏配。采用該方法所獲得的峰匹配結(jié)果更加準(zhǔn)確合理，可直接為下一步DNA序列各堿基質(zhì)量評(píng)分中的參數(shù)估計(jì)提供可靠的依據(jù)，從而確保后續(xù)DNA圖譜分析的準(zhǔn)確有效性。

1 峰匹配動(dòng)態(tài)規(guī)劃方法理論分析與設(shè)計(jì)

通常情況下，DNA測(cè)序熒光光譜信號(hào)并不完全理想[5-8]。圖1所示為一組DNA測(cè)序熒光光譜譜圖。

圖1 DNA測(cè)序熒光光譜譜圖信號(hào)

研究中觀察到，各通道熒光光譜中的DNA堿基信號(hào)呈現(xiàn)程度不一的重疊混雜現(xiàn)象，除包含待檢DNA序列應(yīng)有的堿基信號(hào)峰外，還混雜了多種偽峰如因運(yùn)行環(huán)境及試劑引起的既寬又高的無用雜峰、因不同熒光串?dāng)_帶來的的非堿基峰、以及儀器正常運(yùn)行存在不可避免的噪聲信號(hào)峰等。因此，需對(duì)各堿基通道信號(hào)的識(shí)別峰進(jìn)行進(jìn)一步的判斷。

與DNA片段STR分析不同，儀器進(jìn)行DNA測(cè)序時(shí)沒有相應(yīng)的計(jì)量手段。STR分析時(shí)，待檢DNA樣品與內(nèi)標(biāo)同時(shí)檢測(cè)，根據(jù)內(nèi)標(biāo)確定DNA片段的長(zhǎng)度，并結(jié)合Ladder進(jìn)行等位基因質(zhì)量匹配[14，15]。因此，針對(duì)DNA測(cè)序必須提供一個(gè)客觀的測(cè)量方法，對(duì)DNA測(cè)序數(shù)據(jù)進(jìn)行正確的堿基識(shí)別。這里，不妨假設(shè)，識(shí)別峰與預(yù)測(cè)峰的特征信息已經(jīng)確定，從而將DNA序列預(yù)測(cè)峰作為計(jì)量標(biāo)準(zhǔn)，根據(jù)預(yù)測(cè)峰的峰位置及峰周期、識(shí)別峰的峰位置及峰相對(duì)面積等特征信息，進(jìn)行識(shí)別峰與預(yù)測(cè)峰的動(dòng)態(tài)匹配。

預(yù)測(cè)峰特征信息是對(duì)四堿基通道的加和信號(hào)通過微分法峰識(shí)別及傅里葉方法估計(jì)信號(hào)峰周期而確定，包括峰位置及峰周期，它基本對(duì)應(yīng)了待檢DNA序列各堿基的理想峰位置，在堿基識(shí)別過程中起著相對(duì)計(jì)量的作用[11]。識(shí)別峰特征信息是對(duì)各堿基通道信號(hào)通過微分法或其它方法而確定，包括峰位置、峰高、峰相對(duì)面積、峰型(是否重疊峰、峰對(duì)稱性等)及所代表的堿基，并根據(jù)其表現(xiàn)值的優(yōu)劣性及時(shí)摒棄了噪聲信號(hào)峰及熒光串?dāng)_峰等偽峰，在堿基識(shí)別過程中去偽存真進(jìn)而提供真實(shí)的峰堿基信息。

本研究的目的是利用動(dòng)態(tài)規(guī)劃方法，給每一個(gè)預(yù)測(cè)峰合理分配一個(gè)識(shí)別峰，其關(guān)鍵是如何設(shè)置合理的得分函數(shù)，從而得到各種匹配方案的總得分，從中選擇得分最高的作為最優(yōu)峰匹配結(jié)果。

1.1 峰匹配仿射得分函數(shù)的設(shè)計(jì)

從直觀上講，分配給每一預(yù)測(cè)峰的識(shí)別峰應(yīng)該是真峰即正確的堿基峰，需在充分把握DNA測(cè)序熒光光譜信號(hào)各種偽峰及真峰實(shí)際特點(diǎn)的基礎(chǔ)上，以預(yù)測(cè)峰的峰位置及峰周期作為相對(duì)計(jì)量標(biāo)尺，根據(jù)識(shí)別峰的峰位置、峰高、峰相對(duì)面積、峰型(是否重疊峰、峰對(duì)稱性)等特征信息表現(xiàn)值進(jìn)行判斷，排除該預(yù)測(cè)峰附近識(shí)別峰中的各種偽峰，提取對(duì)應(yīng)的真實(shí)堿基峰。因此，峰匹配得分函數(shù)值的大小應(yīng)與識(shí)別峰的峰位置、峰高、峰相對(duì)面積、峰型等特征信息表現(xiàn)值的優(yōu)劣性相吻合。

在動(dòng)態(tài)規(guī)劃方法應(yīng)用設(shè)計(jì)時(shí)，常見的得分函數(shù)有簡(jiǎn)單線性函數(shù)，還有凸函數(shù)、凹函數(shù)等多種形式，需結(jié)合實(shí)際以簡(jiǎn)單、準(zhǔn)確、有效為最佳選取準(zhǔn)則[16-18]。針對(duì)DNA測(cè)序熒光光譜信號(hào)，通過全面綜合分析信號(hào)堿基峰與各種偽峰的峰特征信息特點(diǎn)，在研究分析相關(guān)文獻(xiàn)資料的基礎(chǔ)上[11-12，16-18]，研究設(shè)計(jì)了改進(jìn)的峰匹配得分函數(shù)如下：

(1)

其中，[shift]是對(duì)shift取整，shift為該識(shí)別峰idePeak與預(yù)測(cè)峰perdPeak的峰位偏移值：

(2)

idearea為識(shí)別峰的實(shí)際面積；α、β分別為識(shí)別峰在右邊、在左邊的懲罰因子；perdLoc為預(yù)測(cè)峰的位置，ideLoc為識(shí)別峰的位置，perdPeriod為預(yù)測(cè)峰的周期。

該匹配得分函數(shù)與識(shí)別峰的偏移值、實(shí)際面積、峰型有關(guān)。當(dāng)偏移絕對(duì)值小、面積大、峰型好時(shí)得分高，反之得分低。此得分函數(shù)在偏移值為零的左右兩側(cè)都為下凹函數(shù)，在零點(diǎn)處得分最高。函數(shù)在零點(diǎn)的兩側(cè)的下降速度取決于懲罰參數(shù)的設(shè)置。

1.2 峰匹配方法參數(shù)選項(xiàng)設(shè)計(jì)

該方法涉及峰偏移、峰面積、峰型共三種峰匹配閾值參數(shù)。其中，峰偏移閾值參數(shù)包括懲罰因子α,β及峰位偏移MinShift、MaxShift與峰個(gè)數(shù)偏移LeftIndex、RightIndex；峰面積參數(shù)為最小相對(duì)面積MinArea；峰型參數(shù)為最小峰分割面積MinSplitArea與峰對(duì)稱性系數(shù)SymCoeff，詳見公式(1)、(2)和參數(shù)選項(xiàng)表1。

表1 峰匹配參數(shù)選項(xiàng)表

參數(shù)選項(xiàng)描述參數(shù)值選擇準(zhǔn)則峰偏移參數(shù)MinShiftMaxShiftLeftIndexRightIndexα,βMinShift為最小峰位偏移閾值,為負(fù)值;MaxShift為最大峰位偏移閾值,為正值;峰位偏移在允許范圍內(nèi),預(yù)測(cè)峰向左或向右移動(dòng)并尋找匹配方案。LeftIndex為預(yù)測(cè)峰向左偏移最大峰個(gè)數(shù); RightIndex 為預(yù)測(cè)峰向右偏移最大峰個(gè)數(shù);峰偏移個(gè)數(shù)在允許范圍內(nèi),預(yù)測(cè)峰向左或向右移動(dòng)并尋找匹配方案。α,β分別為識(shí)別峰在右邊、在左邊的懲罰因子,直接影響峰匹配得分值的大小。根據(jù)信號(hào)實(shí)際特性而定, Min-Shift、MaxShift、LeftIndex、Right-Index應(yīng)大小適中,以保證其對(duì)應(yīng)的信號(hào)數(shù)據(jù)點(diǎn)的偏移個(gè)數(shù)大約在半個(gè)峰周期到兩個(gè)峰周期范圍內(nèi)。根據(jù)最佳匹配實(shí)際需求及經(jīng)驗(yàn)確定懲罰因子,當(dāng)識(shí)別峰與當(dāng)前預(yù)測(cè)峰峰位一致時(shí)為理想匹配。峰面積參數(shù)MinAreaMinArea為識(shí)別峰最小相對(duì)面積閾值,識(shí)別峰相對(duì)面積大于閾值時(shí)參與匹配。根據(jù)信號(hào)實(shí)際特性而定,以排除偽峰信號(hào),提取可能的真實(shí)信號(hào)。峰型參數(shù)SplitAreaSymCoeffSplitArea為識(shí)別峰可分割相對(duì)面積閾值,識(shí)別峰相對(duì)面積大于閾值時(shí),進(jìn)行峰分割。SymCoeff為峰對(duì)稱性(左、右半峰占全峰面積的最小比例)閾值,僅當(dāng)識(shí)別峰相對(duì)面積較大進(jìn)行峰分割時(shí),對(duì)分割峰附近的其它識(shí)別峰進(jìn)行峰對(duì)稱性判斷,大于閾值時(shí)替代分割峰進(jìn)行最佳匹配。一般設(shè)定SplitArea>1.5,因此時(shí)可能為重疊峰情形。設(shè)定 Sym-Coeff>0.38,以保證排除既寬又高的無用雜峰,選取峰對(duì)稱性較好的真實(shí)堿基峰。

在峰匹配方法設(shè)計(jì)時(shí)，根據(jù)DNA測(cè)序熒光光譜信號(hào)各種偽峰及真峰的特征信息，通過設(shè)置合理的匹配參數(shù)閾值可優(yōu)化峰匹配的實(shí)際效果。

1.3 預(yù)測(cè)峰識(shí)別峰動(dòng)態(tài)匹配的設(shè)計(jì)與實(shí)現(xiàn)

峰匹配包括如下三個(gè)階段：

(1)確定容易匹配的：利用設(shè)定的峰匹配參數(shù)閾值，根據(jù)預(yù)測(cè)峰、識(shí)別峰的峰位偏移值和識(shí)別峰相對(duì)面積，當(dāng)兩者都在允許范圍內(nèi)時(shí)將該識(shí)別峰匹配給當(dāng)前預(yù)測(cè)峰；(2)用動(dòng)態(tài)規(guī)劃算法對(duì)第一階段未匹配的進(jìn)行匹配：對(duì)于每一對(duì)未能成功匹配的預(yù)測(cè)峰和識(shí)別峰，首先按照設(shè)計(jì)的峰匹配得分函數(shù)公式(1)計(jì)算匹配得分值，再對(duì)每一預(yù)測(cè)峰與識(shí)別峰進(jìn)行匹配，從所有匹配方案中找出得分最高的分配方案即為最優(yōu)的，將得分值最高的識(shí)別峰匹配給這一預(yù)測(cè)峰；(3)對(duì)第一、二階段都未匹配但確實(shí)認(rèn)為是堿基峰的進(jìn)行匹配：對(duì)于沒有匹配到預(yù)測(cè)峰的識(shí)別峰、沒有匹配到識(shí)別峰的預(yù)測(cè)峰這兩種情況，分別檢查附近的預(yù)測(cè)峰、附近的識(shí)別峰是否都已匹配，如果是根據(jù)相對(duì)面積大小匹配，如果不是將識(shí)別峰匹配給預(yù)測(cè)峰。最后，對(duì)于以上三階段中沒有任何識(shí)別峰可匹配的預(yù)測(cè)峰，對(duì)應(yīng)的堿基定義為N，這種情況的出現(xiàn)是非常少見的。

峰匹配算法的實(shí)現(xiàn)，如流程圖2。

圖2 峰匹配算法實(shí)現(xiàn)流程圖

2 數(shù)值實(shí)驗(yàn)及結(jié)果分析

根據(jù)上述算法設(shè)計(jì)和程序流程圖，在Visual Studio 2005下實(shí)現(xiàn)了峰匹配，因限于篇幅具體程序省略。

本研究示例對(duì)圖1所示的DNA測(cè)序數(shù)據(jù)，設(shè)置峰匹配參數(shù)值：MinShift=-0.5、MaxShift=2.1、LeftIndex=2、RightIndex=4、α=0.1、β=0.3、MinArea=0.2、SplitArea=1.6、SymCoeff=0.39。

如圖3至圖7所示為仿真示例峰匹配結(jié)果，其中4條曲線分別代表堿基‘T’、‘G’、‘C’、‘A’通道的熒光光譜信號(hào)，峰匹配對(duì)應(yīng)的堿基位置用‘○’來標(biāo)示，偽峰對(duì)應(yīng)的峰頂點(diǎn)位置用‘*’來標(biāo)示。

圖3 峰匹配堿基識(shí)別偽峰排除結(jié)果顯示圖(星標(biāo)所示)

圖4 初始信號(hào)段峰匹配堿基識(shí)別結(jié)果顯示圖

圖5 信號(hào)質(zhì)量較好段峰匹配堿基識(shí)別結(jié)果顯示圖

圖6 信號(hào)質(zhì)量衰減段峰匹配堿基識(shí)別結(jié)果顯示圖

圖7 低信號(hào)質(zhì)量尾段峰匹配堿基識(shí)別結(jié)果顯示圖

如圖3所示，采用本研究所設(shè)計(jì)的方法，摒棄了因試劑非正常干擾所引起的偽峰，有效提取了真實(shí)堿基信號(hào)峰，優(yōu)化了文獻(xiàn)[11]所提出的方法，所獲得的峰匹配結(jié)果更加準(zhǔn)確合理。

如圖4所示，在信號(hào)質(zhì)量較低的初始段，雖然噪聲較大，峰匹配結(jié)果準(zhǔn)確度很高。除最初的幾個(gè)寬峰外均與實(shí)際堿基序列結(jié)果相吻合。

如圖5、圖6所示，排除了因熒光串?dāng)_、噪聲干擾所引起的各種偽峰，并提取了這真實(shí)的堿基信號(hào)峰。即使在信號(hào)質(zhì)量衰減段，與實(shí)際堿基序列結(jié)果相吻合的程度極高。

如圖7所示，在信號(hào)質(zhì)量低的尾段，因堿基信號(hào)周期處理方式，堿基識(shí)別的個(gè)數(shù)較多，本研究識(shí)別峰與預(yù)測(cè)峰匹配算法在信號(hào)尾段的處理結(jié)果受預(yù)測(cè)峰識(shí)別算法的影響。但因尾段信號(hào)質(zhì)量太差不可靠，在后續(xù)DNA序列分析實(shí)際應(yīng)用時(shí)并不用這段的堿基排序結(jié)果。

通過多組峰匹配仿真實(shí)驗(yàn)與結(jié)果比對(duì)，采用本研究所設(shè)計(jì)的方法可獲得與實(shí)際堿基序列吻合程度極高的結(jié)果，尤其在中段信號(hào)質(zhì)量較好階段。

3 結(jié)論

針對(duì)DNA測(cè)序熒光光譜信號(hào)所設(shè)計(jì)的峰匹配動(dòng)態(tài)規(guī)劃算法，充分考慮了信號(hào)堿基峰、及各種偽峰的峰特征信息的不同特點(diǎn)，通過設(shè)計(jì)改進(jìn)的匹配得分標(biāo)準(zhǔn)，動(dòng)態(tài)規(guī)劃可以最大限度的將識(shí)別峰和預(yù)測(cè)峰進(jìn)行匹配，盡量做到不錯(cuò)配、不漏配，能確定高準(zhǔn)確度的待檢DNA序列的堿基排序結(jié)果，為下一步DNA序列各堿基質(zhì)量評(píng)分中的參數(shù)估計(jì)提供可靠的依據(jù)，從而確保后續(xù)DNA圖譜分析的準(zhǔn)確有效性。