基于动态规划的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图谱分析的准确有效性。