基于動態規劃的DNA堿基識別峰匹配方法的設計與實現

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

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

DNA測序技術自1980年發展至今取得了令人矚目的進展，基于熒光誘導熒光毛細管電泳檢測的DNA測序技術因其高精度DAN序列識別準確率而在個體識別、親緣鑒定、災難身份識別、法醫DNA鑒定等領域依然發揮著無可替代的重要作用[1-3]。

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

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

1 峰匹配動態規劃方法理論分析與設計

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

圖1 DNA測序熒光光譜譜圖信號

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

與DNA片段STR分析不同，儀器進行DNA測序時沒有相應的計量手段。STR分析時，待檢DNA樣品與內標同時檢測，根據內標確定DNA片段的長度，并結合Ladder進行等位基因質量匹配[14，15]。因此，針對DNA測序必須提供一個客觀的測量方法，對DNA測序數據進行正確的堿基識別。這里，不妨假設，識別峰與預測峰的特征信息已經確定，從而將DNA序列預測峰作為計量標準，根據預測峰的峰位置及峰周期、識別峰的峰位置及峰相對面積等特征信息，進行識別峰與預測峰的動態匹配。

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

本研究的目的是利用動態規劃方法，給每一個預測峰合理分配一個識別峰，其關鍵是如何設置合理的得分函數，從而得到各種匹配方案的總得分，從中選擇得分最高的作為最優峰匹配結果。

1.1 峰匹配仿射得分函數的設計

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

在動態規劃方法應用設計時，常見的得分函數有簡單線性函數，還有凸函數、凹函數等多種形式，需結合實際以簡單、準確、有效為最佳選取準則[16-18]。針對DNA測序熒光光譜信號，通過全面綜合分析信號堿基峰與各種偽峰的峰特征信息特點，在研究分析相關文獻資料的基礎上[11-12，16-18]，研究設計了改進的峰匹配得分函數如下：

(1)

其中，[shift]是對shift取整，shift為該識別峰idePeak與預測峰perdPeak的峰位偏移值：

(2)

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

該匹配得分函數與識別峰的偏移值、實際面積、峰型有關。當偏移絕對值小、面積大、峰型好時得分高，反之得分低。此得分函數在偏移值為零的左右兩側都為下凹函數，在零點處得分最高。函數在零點的兩側的下降速度取決于懲罰參數的設置。

1.2 峰匹配方法參數選項設計

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

表1 峰匹配參數選項表

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

在峰匹配方法設計時，根據DNA測序熒光光譜信號各種偽峰及真峰的特征信息，通過設置合理的匹配參數閾值可優化峰匹配的實際效果。

1.3 預測峰識別峰動態匹配的設計與實現

峰匹配包括如下三個階段：

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

峰匹配算法的實現，如流程圖2。

圖2 峰匹配算法實現流程圖

2 數值實驗及結果分析

根據上述算法設計和程序流程圖，在Visual Studio 2005下實現了峰匹配，因限于篇幅具體程序省略。

本研究示例對圖1所示的DNA測序數據，設置峰匹配參數值：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所示為仿真示例峰匹配結果，其中4條曲線分別代表堿基‘T’、‘G’、‘C’、‘A’通道的熒光光譜信號，峰匹配對應的堿基位置用‘○’來標示，偽峰對應的峰頂點位置用‘*’來標示。

圖3 峰匹配堿基識別偽峰排除結果顯示圖(星標所示)

圖4 初始信號段峰匹配堿基識別結果顯示圖

圖5 信號質量較好段峰匹配堿基識別結果顯示圖

圖6 信號質量衰減段峰匹配堿基識別結果顯示圖

圖7 低信號質量尾段峰匹配堿基識別結果顯示圖

如圖3所示，采用本研究所設計的方法，摒棄了因試劑非正常干擾所引起的偽峰，有效提取了真實堿基信號峰，優化了文獻[11]所提出的方法，所獲得的峰匹配結果更加準確合理。

如圖4所示，在信號質量較低的初始段，雖然噪聲較大，峰匹配結果準確度很高。除最初的幾個寬峰外均與實際堿基序列結果相吻合。

如圖5、圖6所示，排除了因熒光串擾、噪聲干擾所引起的各種偽峰，并提取了這真實的堿基信號峰。即使在信號質量衰減段，與實際堿基序列結果相吻合的程度極高。

如圖7所示，在信號質量低的尾段，因堿基信號周期處理方式，堿基識別的個數較多，本研究識別峰與預測峰匹配算法在信號尾段的處理結果受預測峰識別算法的影響。但因尾段信號質量太差不可靠，在后續DNA序列分析實際應用時并不用這段的堿基排序結果。

通過多組峰匹配仿真實驗與結果比對，采用本研究所設計的方法可獲得與實際堿基序列吻合程度極高的結果，尤其在中段信號質量較好階段。

3 結論

針對DNA測序熒光光譜信號所設計的峰匹配動態規劃算法，充分考慮了信號堿基峰、及各種偽峰的峰特征信息的不同特點，通過設計改進的匹配得分標準，動態規劃可以最大限度的將識別峰和預測峰進行匹配，盡量做到不錯配、不漏配，能確定高準確度的待檢DNA序列的堿基排序結果，為下一步DNA序列各堿基質量評分中的參數估計提供可靠的依據，從而確保后續DNA圖譜分析的準確有效性。