一种基于BLAST的多种遗传标记序列比对方法＊

赵钊，王彤，孙洁，周勤虎，刘亚楠，胡蓉

（公安部物证鉴定中心，北京100038；法医遗传学公安部重点实验室，北京100038）

一种基于BLAST的多种遗传标记序列比对方法＊

赵钊，王彤，孙洁，周勤虎，刘亚楠，胡蓉

（公安部物证鉴定中心，北京100038；法医遗传学公安部重点实验室，北京100038）

设计并实现了一种基于BLAST（Basic Local Alignment Search Tool）的生物序列比对方法，它可以应用于法庭科学多种遗传标记的检索比对。引入BLAST算法并设计合理的打分函数，使其在兼容序列多态性比如线粒体DNA比对的同时，可以有效地进行长度多态性比如常染色体STR的比对。实验表明，采用这种方法可以实现不同序列的有效比对，输出准确的有序结果集，具有较高的比对效率。

法庭科学DNA数据库；检索比对；局部序列比对算法；序列多态性

近年来，法庭科学DNA检验在方法学层面上已经逐渐趋于成熟和稳定，在应用层面上已为社会普遍接受。法庭科学DNA数据库的建设应用也取得了令人瞩目的成效，数据量增长趋势明显。法庭科学DNA数据库作为集成现代DNA检验技术、信息技术和科学管理等要素，实现跨时空多元化应用、精确打击犯罪的武器，破案整体效益日益突出[1-4]。随着DNA检验技术的发展，越来越多的序列多态性数据逐渐被应用于实践[5-6]，尤其是针对重大疑难案件的深度研判，需要设计一种适用于多种遗传标记的通用比对算法。这样既可以作为现有DNA数据库检索的有益补充，也符合DNA检验与应用的整体发展趋势。

1 相关技术

1.1 序列比对

序列比对主要应用于DNA或蛋白质的相似性比对上，是生物信息学中的重要内容，也是生物数据库的基础。进行序列比对的算法主要有全局比对算法和局部比对算法以及在二者基础上的改进算法[7-8]。1990年由Altschul等人发布的基于局部序列比对的BLAST（Basic Local Alignment Search Tool）序列检索工具，采用短片段匹配算法和一种有效的统计模型，可以找出目的序列与数据库之间的最佳局部比对效果[9-10]。该算法具有较快的比对速度和较高的比对精度，在速度上比只使用动态规划大约快50倍左右，逐渐成为现今最为流行的数据库搜索算法之一，被广泛应用于解决蛋白质DNA序列的分析问题[11]。例如最著名的NCBI（美国国立生物技术信息中心）重点开发的GenBank DNA序列数据库便利用BLAST进行相似搜索。BLAST能够在15 s的时间内对整个DNA数据库执行序列搜索，同时，它也是EBI（欧洲生物信息研究中心）、Sanger（英国桑格研究院）等多个主要研究机构常用的生物信息数据库的搜索工具[12]。

1.2 法庭科学DNA数据库

从20世纪90年代初开始，美国、英国、加拿大、法国和日本等许多发达国家和地区已经开始建设DNA数据库。随着时间的推移，DNA数据库在刑事案件的侦破方面发挥了巨大的作用，而且发展迅猛[13-15]。我国公安机关从1999年开始研究和应用DNA数据库，近几年来取得了规模化的发展。我国的法庭科学DNA数据库已经容纳了五千万以上的短串联重复序列（short tandem repeat，STR）数据，是目前世界上最大的DNA数据库。近年来，国内外陆续建立了Y-STR数据库和线粒体DNA（mtDNA）数据库，但是，与常染色体STR数据库比对相互独立，并没有多种遗传标记的通用比对模式[16-17]。目前，中国法庭科学DNA数据库检索主要针对常染色体STR进行逐个数字的容差式精确匹配，通常多于1～2个等位基因值的容差即跳出，比对结果为有比中或无比中，比对模式单一，虽有较高的比对效率，但结果信息不够全面。本文实现的方法支持不同形式的遗传标记数据的检索、比对，支持相似性模糊匹配，在一些重大案件的线索排查中可以充分发挥DNA信息的作用。

2 基于BLAST的多种遗传标记检索

2.1 数据准备

不同遗传标记呈现的形式不同，比如常染色体STR、Y-STR是数字形式，而mtDNA是彼此不一的字母形式。本文设计的方法首先支持常染色体STR、Y-STR、mtDNA等不同遗传标记的统一字符串格式化存储，以此作为检索比对的基础。实现的比对系统基于Oracle 10g完成了样本数据库的搭建，并实现了基于高效的C/C++文件流进行序列文件的导入，数据举例如表1所示。其中，Y-STR、STR各基因座之间以空格符相隔，而同一基因座杂合子等位基因值之间以“/”相隔。

表1 数据格式化存储示例

2.2 BLAST算法步骤

BLAST的核心算法步骤是将需查询的序列数据读入内存，然后从数据库的第一条记录开始到最后一条记录逐条与待查询的序列进行比较。具体步骤如下：①设置一个最小片断长度T，对查询序列建立长度为T的片段索引表。②从数据库中查找与表中匹配的无空位片段（如果遇到空格或“/”忽略不计入）。③设定片断匹配最小分值S，对这些片段利用动态规划方法进行无空位延伸，匹配得正分，不匹配得负分，直到得到一个局部最高分，也就是高分序列片段HSP。如果得到的HSP值小于S，则舍弃此序列。同时，为了避免无限制延伸，算法规定如果片段在进行一段延伸后发现此时的分值与当前的HSP值相差超过一个设定的值D，则停止延伸。比较此时的HSP与S，决定此片断的取舍。查询序列的相似度由这些局部最高分片断决定。

其中，利用动态规划进行无空位延伸为BLAST的核心，相关参数设置为，最小片段长度T值为11，HSP相差阈值D为4，正确匹配得分为1，错误匹配得分为-2.

以mtDNA为例，长度为11的初始序列片段对如下：

（初始序列片段延伸位置22，延伸分数11）

此时，延伸位置为22，延伸分数为11，HSP值为11，位置为22.

基于初始序列片段对向右进行延伸，从位置22向右延伸，发现从第23位到第27位的延伸均匹配，延伸分数随着向右的扩展在逐步增长。与此同时，HSP值和位置也发生了相应的变化，具体情况如下：

（延伸位置27，延伸分数16）

此时，HSP值相应变为16，HSP位置变为27.当延伸至第28位时，T与C不匹配，延伸分数变为14.当延伸到第28位时，HSP值保持不变为16（位置为27），具体如下：

（延伸位置28，延伸分数14）

此时的HSP值和位置仍然保持不变，分别为16和27；延伸分数与HSP值的差值为16－14＝2，没有大于HSP相差阈值D＝4，因此，继续向右延伸。

从第29位到第31位的延伸均匹配，当延伸到第29位时，延伸分数为15.此时的延伸分数小于HSP值16，并未改变HSP。同样，当延伸到第30位时也未改变HSP。当延伸到第31位时，此时的延伸分数为17，延伸分数大于HSP值16，HSP的值和位置变为17和31.第29位到第31位的匹配和衍生情况如下：

（延伸位置31，延伸分数17）

此时，HSP值相应变为17，HSP位置变为31.继续向右延伸，32位到34位均不同，当延伸到第34位时，具体延伸情况如下：

（延伸位置34，延伸分数11）

此时的延伸分数为11，HSP值与延伸分数的差值为17－11＝6，差值大于HSP相差阈值D，因此，算法停止向右延伸。这时，即使后续的3个位置（第35到第37位）均匹配，也不再进行向右的延伸。在以上延伸过程中，延伸分数、延伸位置和HSP的变化情况如表2所示。

至此，算法完成了向右延伸。最终的延伸位置为31，HSP值为17.向左的延伸同向右的步骤。

表2 BLAST延伸过程中的数据变化

2.3 打分机制

根据以上算法步骤逐一比对后，需要设立打分机制实现最终的检索结果排序。在匹配算法的过程中，用打分函数控制整个匹配过程。本文通过分析mtDNA、STR和YSTR的特征，设计了相关的打分函数。

针对mtDNA，对于序列Si和Sj，打分函数Score（Si，Sj）定义如下：

式（1）中：α为序列中匹配字符的数量；R为匹配字符奖励的分数（即R≥0）；β为序列中不匹配字符的数量；P为不匹配字符惩罚的分数（P≤0）。

对于小节2.2中BLAST的动态规划无空位延伸中的例子，当延伸位置为31时，R＝1，P＝－2，α＝19，β＝1，此时的分数为α×R＋β×P＝19×1＋1×（－2）＝17.

针对STR和YSTR，根据其数字型的类型，项目定制了如下打分函数。对于长度为n的序列，Si和Sj分别为：Si＝｛Si1，Si2…Sik…Sin｝，Sj＝｛Sj1，Sj2…Sjk…Sjn｝.

2.4 实验结果

实现基于BLAST的多种遗传标记序列比对系统，搭建实验数据库，以mtDNA、YSTR、STR真实数据为基础模拟批量数据进行算法验证。系统比对结果按照匹配度从高到低返回TOP序列匹配表，同时，也可以输出序列结果集的匹配详情，以*标记不同，-标记相同。以mtDNA、STR（YSTR与STR类似）为例，结果形式见图1、图2、图3和图4.

图1 mtDNA序列匹配表

图2 mtDNA序列匹配详情

图3 STR序列匹配表

图4 STR序列匹配详情

为了提高比对效率，系统使用加载数据到内存后进行比对，以不同数据量进行性能测试，结果见表3、表4、表5.

表3 mtDNA序列比对性能数据

表4 STR序列比对性能数据

表5 YSTR序列比对性能数据

3 结束语

应用BLAST算法与不同打分函数实现对法庭科学领域常见遗传标记的序列检索和结果筛选，从综合应用角度看，数据处理、检索比对的性能较为理想。本文的检索方式相比于目前中国法庭科学DNA数据库容差比对方式，更能体现生物序列的匹配程度，能够提供更为灵活和全面的有序结果集，而不是简单的有无比中。中国法庭科学DNA数据国家库常染色体STR数据已臻数千万，随机匹配率比较高[18]，如果将本文方法应用于STR数据的全库比对，从效率上来说并不适合。其应用价值体现在以下2个方面：①可以就个案在小范围（比如某地区或某年龄段）内进行相似度匹配检索，以求得隐藏生物线索；②为实现包含多种遗传标记的综合性法庭科学DNA数据库提供算法基础和检索对照。

［1］陈连康，周怀谷，顾丽华，等.上海市公安局“法庭科学DNA数据库”介绍［J］.刑事技术，2003（5）：3-6.

［2］侯一平，王保捷，丛斌，等.中国法医学会物证专业委员会法医DNA分析的若干建议［J］.中国法医学杂志，2006，21（5）：257-259.

［3］张怀才.试述我国公安DNA数据库在侦查中的应用与展望［D］.上海：华东政法大学，2014.

［4］姜先华.中国DNA数据库建设应用技术现状及发展趋势［J］.中国法医学杂志，2011，26（5）：383-386.

［5］王乐，叶健，白雪，等.二代测序技术及其在法医遗传学中的应用［J］.刑事技术，2015，40（5）：353-358.

［6］张素华，边英男，赵琪，等.二代测序技术在法医学中的应用进展［J］.法医学杂志，2016，32（4）：282-289.

［7］Xiong，J.Essential Bioinformatics［M］.New York：Cambridge University Press，2006.

［8］Niranjan P.Hidden Markov Models Incorporaitiong Fuzzy Measures and Integrals for Protein Sequence Identification and Alignment.Genomics Proteomics&bioinformatics，2008，6（2）：98-110.

［9］PAlexander，FJW Iii.Having a BLAST with bioinformatics（and avoiding BLAST phemy）.Genome Biology，2001，2（10）：1-10.

［10］李红燕.基于BLAST算法的序列分析软件开发［D］.长沙：中南大学，2009.

［11］SF Altschul，TL Madden，AA Schäffer，et al.Gapped BLAST and PSI BLAST：a new generation of protein database search programs.Nucl.Acids Res，1997（25）：3389-3402.

［12］吕军，张颖，冯立芹，等.生物信息学工具BLAST的使用简介［J］.内蒙古大学学报，2003，34（2）：179-186.

［13］焦文慧，宋辉.英美国家犯罪DNA数据库建设及应用［J］.上海公安高等专科学校学报，2013，23（2）：86-91.

［14］葛百川，王海鸥，陈连康，等.赴美考察DNA数据库及DNA实验室的情况介绍［J］.刑事技术，2010（3）：3-6.

［15］赵兴春，李路平，高洵.英国DNA技术应用与国家DNA数据库［J］.公安大学学报：自然科学版，1999（2）：17-21.

［16］刘冰.现阶段我国DNA数据库发展的几个关键问题［J］.刑事技术，2015，40（4）：318-323.

［17］S Montague.Global Governance of Forensic DNA Profiling and Databasing.Yale Journal of Biology& Medicine，2011，84（3）：326.

［18］葛建业，严江伟，Bruce Budowle，等.关于法庭科学DNA数据库若干问题的探讨［J］.中国法医学杂志，2011，26（3）：252-255.

D919.2

A

10.15913/j.cnki.kjycx.2017.18.001

2095－6835（2017）18－0001－04

赵钊（1984—），女，河北宣化人，助理研究员，硕士，主要研究领域是法庭科学DNA数据库应用。

〔编辑：白洁〕

公安部科技强警基础工作专项项目“法医DNA二代测序数据批量比对关键技术研究“（编号：2016GABJC18）