一种短串联重复序列重复数的检测和分型方法

摘要

本发明涉及测序数据的生信分析领域，具体提供一种基于二代测序的短串联重复序列重复数的检测和分型方法。所述检测方法能够跨过中间重复区域，对两个侧翼序列实现同时比对，并且同时考虑侧翼序列的错配、插入和缺失；所述分型方法是基于检测有效峰值，结合重复数差异进行动态阈值设定的STR分型方法，可应用于不同二代测序平台的重复数的检测。

权利要求

1.一种基于二代测序的STR侧翼序列搜索比对方法，其特征在于，包括如下步骤：步骤1）原始数据处理步骤：将fastq文件和STR参考序列库进行比对；

步骤2）比对算法构建步骤：

a、已知左侧侧翼序列flank_left、右侧侧翼序列flank_right和待分析序列S1，将flank_left、flank_right按顺序拼接成序列S2；

b、初始化得分矩阵：以S1为列、S2为行构建得分矩阵，设置第1列惩罚得分都为0；设置第1行S2中左右侧翼序列拼接处的两个位置惩罚得分分别为‑10，其余位置惩罚得分为‑1；

c、填充得分矩阵：设置匹配得分match_score = 1，错配得分mismatch_score = ‑1，gap得分gap_score= ‑2；对于矩阵中每个单元格（i,j），其中的i代表矩阵行的索引，j代表矩阵列的索引，得分计算公式为：

若j ≠x 且 j ≠ y : F_ij=max(F_(i‑1,j‑1)+S(A_i,B_j),F_(i,j‑1)+d,F_(i‑1,j)+d) ；

若 j = x : F_ij=max(F_(i‑1,j‑1)+2S(A_i,B_j),F_(i,j‑1)+0,F_(i‑1,j)+2×d)；

若 j = y : F_ij=max(F_(i‑1,j‑1)+2S(A_i,B_j),F_(i,j‑1)+2×d,F_(i‑1,j)+0)；

其中x是flank_left的长度，y为x+1，d = gap_score，S(A_i,B_j) = match_score 或者 mismatch_score；

步骤3）STR位点重复区域锚定：读取步骤1）比对到的每条read，根据是否能够比对上STR参考序列分为两种情况：若比对上STR位点，则使用该STR位点的侧翼序列和上述比对算法搜索STR重复区域；若未比对上任何STR位点，则遍历所有STR位点的侧翼序列，并且使用上述比对算法搜索STR重复区域；若侧翼序列比对的结果总共不超过2个错配、插入和缺失，则输出比对上的STR位点以及重复区域序列，否则舍弃这条read；最终获得每条read锚定的重复区域序列长度。

2.权利要求1所述的基于二代测序的STR侧翼序列搜索比对方法，其特征在于，所述方法进一步包括：

步骤4）：根据每条read锚定的重复区域序列长度，计算每条read比对上的STR位点的重复数。

3.权利要求1‑2任一所述的基于二代测序的STR侧翼序列搜索比对方法，其特征在于，所述步骤2）进一步包括：

d、获取重复区域序列：在矩阵中最后一列由下往上选取该列第一个最大分值s作为侧翼序列的比对得分，假设序列S2的长度为n，若n ‑ s ≤ 4，则在允许2个错配、插入或缺失的情况下，S1和S2可以比对上；以s的坐标为源头往左上回溯，得出j=x和j=y时行的坐标i1、i2,取S1上i1、i2位置中间的序列为STR位点重复区域序列。

4.一种基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型方法，其特征在于：

1）基于二代测序的单条read比对上的STR位点的重复数，统计样本每个STR位点上每种重复数检测到的reads数量（RC）；

2）选取第一分型：根据峰值位置来选取分型，对于每个STR位点，RC最大的重复数A为峰,将该重复数判断为第一分型，反之则判断该STR位点未检测到任何分型；

3）遍历选取第二分型：每次找到一个峰之后将该峰RC值赋0，首先将第一分型的RC赋0，即令RC（A）= 0；假设重复数A、B的RC分别为RC(A)、RC（B），需同时满足以下3个条件则判定重复数B为峰：

a. RC(B)/RC（A）≥ 待检测分型B的梯度阈值T，所述阈值T=初始梯度 + 0.25 × （左侧梯度或者右侧梯度）×（A ‑ B）；

b.RC（B）≥ 3；

c.RC（B）高于相邻重复数的RC;

4）输出分型结果，若同时找到满足要求的第一分型和第二分型，则分型结果为（A，B）；

若找到满足要求的第一分型但没有找到满足要求的第二分型，则分型结果为（A）；若没有找到满足要求的第一分型和第二分型，则未检测到分型；

所述步骤1）基于二代测序的单条read比对上的STR位点的重复数采用权利要求1‑3任一所述方法。

5.权利要求4所述的基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型方法，其特征在于，所述步骤2）中，所述最大的重复数A的RC ≥ 20。

6.权利要求5所述的基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型方法，其特征在于，所述步骤3）中，所述待检测分型B的梯度阈值T的计算如下：T = 0.28 + 0.25 × step ×（A ‑ B），其中，若A ‑ B > 0，step= 0.19；

若A ‑ B < 0 ，step = 0.33。

7.权利要求6所述的基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型方法，其特征在于，所述T=0.02。

8.一种基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型系统，其特征在于，所述系统包括如下模块，所述模块执行权利要求4‑7任一所述方法的步骤。

9.一种计算机可读介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现权利要求1‑7任一所述方法。

10.一种电子设备，包括处理器以及存储器，所述存储器上存储一条或多条可读指令，所述一条或多条可读指令被所述处理器执行时，实现权利要求1‑7任一所述方法。

说明书

一种短串联重复序列重复数的检测和分型方法

技术领域

[0001] 本发明属于生信分析领域，具体涉及一种短串联重复序列重复数的检测和分型方法。

背景技术

[0002] 短串联重复序列(short  tandem  repeat,STR)，也称做微卫星DNA(micrositellite DNA)，是由2～6个碱基作为核心串联合成的DNA序列，STR具有变异率高，
多态性，易检测等特点，因此广泛地应用与法医相关的检测。从1985年开始，STR检测就已经
被应用到法医领域，通常是通过毛细管电泳配合荧光标记来检测。针对不同的STR位点进行
特异性引物设计，通过扩增获得不同长度和不同的荧光标记的扩增产物，通过毛细管电泳
区分不同的STR位点。

[0003] 然而，这种基于毛细管电泳检测的方法通量有限，且无法有效区分长度超过1,000bp的片段。此外，使用的荧光通道的颜色一般不超过6种，否则也会导致不同波长荧光的
渗漏问题。因此，基于毛细管电泳检测的STR的数量一般不超过60个。随着高通量测序的发
展，越来越多的研究开始转向使用高通量测序来检测和识别STR等分子标记物。二代测序具
有通量大、检测位点多等优点，对于使用STR检测来进行身份鉴定和亲缘推断非常合适，也
越来越被广泛地应用于法医领域。例如近年美国破获的金州杀人案就是应用基于二代测序
的亲缘推断找到犯罪嫌疑人的亲属，从而确定犯罪嫌疑人。

[0004] STR是法医检测的重要对象，一方面STR的多态性使得其包含的信息更加丰富，平均一个STR位点的信息量大约是6个SNP位点的信息量；另一方面，司法部门在过去的几十年
中积累了数千万个个体的STR信息，而且其中很多是犯罪嫌疑人。使用二代测序来对STR位
点进行检测，相比于毛细管电泳有更多的优势。二代测序除了能够获得STR重复区域长度的
信息，也能够得到重复区域DNA序列的信息，这两种信息加在一起更加强了STR个体识别的
功能。已经有研究表明在人群中存在长度相同而序列不同的分型，而基于毛细管电泳的STR
检测不能够区分这种序列差异导致的不同分型；二代测序需要的扩增子片段要小于一代测
序，更加适合于一些降解DNA的STR检测；二代测序检测分型不依赖于bin值，因此能够检测
到新的和较为稀有的分型，而且有研究表明有的样本的STR重复区域发生较少bp的插入或
者缺失时可能会被毛细管电泳划分到同一个bin中，带来STR分型错误；通过给不同的样本
加上不同的序列标签，二代测序可以同时检测几百个，甚至上千个样本，检测通量大；二代
测序技术可以把STR、SNP、线粒体DNA等靶分子标记物融合到一个panel中，能够在一个
panel中得到大量有用的信息。目前美国已经开始在部分州和城市的法医系统中使用二代
测序代替毛细管电泳检测STR的测试，预计在近几年内开始在部分州和城市中实施。

[0005] 尽管中国也在相关领域进行积极研究和调查，并且部分公司也推出了自己的panel产品，目前已有的基于二代测序检测STR的方法仍然存在不足之处。目前主流的测序
仪是美国illumina公司的测序仪，其测序长度一般不超过250bp，这对于长重复序列STR的
准确分型提出了挑战。另一方面，二代测序检测STR的方法依赖于其两侧的侧翼序列的匹
配，如果侧翼序列发生错配、插入或者缺失，则无法准确定位到STR区域，影响对其重复数的
检测。此外，目前主流STR分型的生信分析方法是通过查看第二个分型的reads数是否为第
一个分型的40％及以上，如果达到则认为该位点是杂合型，有两个不同的分型(即该STR位
点有两种重复数)。由于二代测序的数据量随片段的长度会发生偏差，一般片段越长数据量
越少，所以在检测杂合型STR位点时，如果两个分型重复数相差较大，两个allele上read的
数量可能会差很多(有时两者的比率会超过1：10)，按照40％的阈值就会造成重复数多的
STR分型漏检。

[0006] 由此可见，目前二代测序对于STR分型检测仍然存在多处分型不准的问题，有鉴于此，提出本发明。

发明内容

[0007] 本发明的目的是针对前面现有技术的问题，为克服上述问题，本发明至少做了如下技术创新：

[0008] 1.针对侧翼序列有可能发生的错配、插入和缺失，本发明提出了一种改进的Needleman‑Wunsch算法，该算法能够在O(MN)时间内对两个侧翼序列实现同时比对，并且同
时考虑错配、插入和缺失，而且跨过了中间的重复区域。

[0009] 2.针对二代测序检测STR分型不准的问题，本发明提出一种基于检测有效峰值(时间复杂度为O(T^2))，结合重复数差异进行动态阈值设定的STR分型的新方法。

[0010] 具体的，本发明提出如下技术方案：

[0011] 本发明首先提供一种基于二代测序的STR侧翼序列搜索比对方法，包括如下步骤：

[0012] 步骤1)原始数据处理步骤：将fastq文件和STR参考序列库进行比对；

[0013] 步骤2)比对算法构建步骤：

[0014] a、已知左侧侧翼序列flank_left、右侧侧翼序列flank_right和待分析序列S1，将flank_left、flank_right按顺序拼接成序列S2；

[0015] b、初始化得分矩阵：以S1为列、S2为行构建得分矩阵，设置第1列惩罚得分都为0；设置第1行S2中左右侧翼序列拼接处的两个位置惩罚得分分别为‑10，其余位置惩罚得分
为‑1；

[0016] c、填充得分矩阵：设置匹配得分match_score＝1，错配得分mismatch_score＝‑1，gap得分gap_score＝‑2；对于矩阵中每个单元格(i,j)，其中的i代表矩阵行的索引，j代表
矩阵列的索引，得分计算公式为：

[0017] 若j≠x且j≠y:F_ij＝max(F_(i‑1,j‑1)+S(A_i,B_j),F_(i,j‑1)+d,F_(i‑1,j)+d)；

[0018] 若j＝x:F_ij＝max(F_(i‑1,j‑1)+2S(A_i,B_j),F_(i,j‑1)+0,F_(i‑1,j)+2×d)；

[0019] 若j＝y:F_ij＝max(F_(i‑1,j‑1)+2S(A_i,B_j),F_(i,j‑1)+2×d,F_(i‑1,j)+0)；

[0020] 其中x是flank_left的长度，y为x+1，d＝gap_score，S(A_i,B_j)＝match_score或者mismatch_score；

[0021] 步骤3)STR位点重复区域锚定：读取步骤1)比对到的每条read，根据是否能够比对上STR参考序列分为两种情况：若比对上STR位点，则使用该STR位点的侧翼序列和上述比对
算法搜索STR重复区域；若未比对上任何STR位点，则遍历所有STR位点的侧翼序列，并且使
用上述比对算法搜索STR重复区域；若侧翼序列比对的结果总共不超过2个错配、插入和缺
失，则输出比对上的STR位点以及重复区域序列，否则舍弃这条read；最终获得每条read锚
定的重复区域序列长度。

[0022] 进一步的，所述方法还包括：

[0023] 步骤4)：根据每条read锚定的重复区域序列长度，计算每条read比对上的STR位点的重复数。

[0024] 进一步的，所述步骤2)包括：

[0025] d、获取重复区域序列：在矩阵中最后一列由下往上选取该列第一个最大分值s作为侧翼序列的比对得分，假设序列S2的长度为n，若n‑s≤4，则在允许2个错配、插入或缺失
的情况下，S1和S2可以比对上；以s的坐标为源头往左上回溯，得出j＝x和j＝y时行的坐标
i1、i2,取S1上i1、i2位置中间的序列为STR位点重复区域序列。

[0026] 本发明还提供一种基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型方法，其特征在于：

[0027] 1)基于二代测序的单条read比对上的STR位点的重复数，统计样本每个STR位点上每种重复数检测到的reads数量(RC)；

[0028] 2)选取第一分型：根据峰值位置来选取分型，对于每个STR位点，RC最大的重复数A为峰,将该重复数判断为第一分型，反之则判断该STR位点未检测到任何分型；

[0029] 3)遍历选取第二分型：每次找到一个峰之后将该峰RC值赋0，首先将第一分型的RC赋0，即令RC(A)＝0；假设重复数A、B的RC分别为RC(A)、RC(B)，需同时满足以下3个条件则判
定重复数B为峰：

[0030] a.RC(B)/RC(A)≥待检测分型B的梯度阈值T，所述阈值T＝初始梯度+0.25×(左侧梯度或者右侧梯度)×(A‑B)；

[0031] b.RC(B)≥3；

[0032] c.RC(B)高于相邻重复数的RC；

[0033] 4)输出分型结果，若同时找到满足要求的第一分型和第二分型，则分型结果为(A，B)；若找到满足要求的第一分型但没有找到满足要求的第二分型，则分型结果为(A)；若没
有找到满足要求的第一分型和第二分型，则未检测到分型。

[0034] 进一步的，所述步骤2)中，所述最大的重复数A的RC≥20。

[0035] 进一步的，所述待检测分型B的梯度阈值T的计算如下：

[0036] T＝0.28+0.25×step×(A‑B)，其中，

[0037] 若A‑B>0，step＝0.19；

[0038] 若A‑B<0，step＝0.33；

[0039] 优选的，所述T＝0.02。

[0040] 进一步的，所述步骤1)基于二代测序的单条read比对上的STR位点的重复数具体可采用上述任一方法。

[0041] 本发明还提供一种基于检测有效峰值且结合重复数差异进行动态阈值设定的STR分型系统，所述系统包括如下模块，所述模块执行上述任一方法的步骤。

[0042] 本发明还提供一种计算机可读介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现上述任一所述方法。

[0043] 本发明还提供一种电子设备，包括处理器以及存储器，所述存储器上存储一条或多条可读指令，所述一条或多条可读指令被所述处理器执行时，实现上述任一所述方法。

[0044] 与现有技术相比，本发明至少具有如下优势：

[0045] 1)本发明具有准确率高、速度快等优点。本发明的搜索侧翼序列的时间复杂度为O2
(MN)，而搜索分型的时间复杂度为O(T)，其中的M是两个侧翼序列的长度和，N为read的长
度，T为对于某一个STR位点的检测到的总的分型的个数。当侧翼序列发生1‑2个错配、插入
或者缺失时，本发明改进的Needleman‑Wunsch算法依然能根据侧翼序列检测出STR重复区
域，有效地避免了常规检测方法漏掉含错配、插入或者缺失的侧翼序列的read的情况，数据
利用率更高，检测结果更准确。

[0046] 2)本发明提出检测峰的算法，对于STR位点的第二分型设置不同的reads数量阈值，阈值随着第二分型重复数的增大而减小，随着第二分型重复数的减小而增大。通过设置
动态阈值的方式，能够有效地校正测序长度对分型检出的影响，提高杂合型STR两个位点同
时检出的概率。

[0047] 3)本发明方法显著改善了二代测序技术用于STR的分型检测，检出结果客观性强、准确度高，具有推广使用价值。

附图说明

[0048] 为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的
附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前
提下，还可以根据这些附图获得其他的附图。

[0049] 图1本发明二代测序检测STR的流程；

[0050] 图2改进的NW算法识别侧翼序列的方法示意图；

[0051] 图3改进的NW算法详细以及打分矩阵实例；

[0052] 图4传统测序长度对较长的分型的read数量的影响；

[0053] 图5本发明二代测序通过峰值检测STR分型的算法示意图；

[0054] 图6本发明read数量阈值参数的调整；

[0055] 图7本发明峰值和梯度阈值的参数调整。

具体实施方式

[0056] 下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术
人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0057] 以下术语或定义仅仅是为了帮助理解本发明而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。

[0058] 除非在下文中另有定义，本发明具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员
很好理解，但仍然阐述以下定义以更好地解释本发明。

[0059] 如本发明中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为
是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方
案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

[0060] 在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。

[0061] 本发明中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。

[0062] 此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在
适当的环境下可互换，并且本发明描述的实施方案能以不同于本发明描述或举例说明的其
它顺序实施。

[0063] 本发明所述的方法，总体思路包括如下两方面：

[0064] 一、本发明首先提出了一种STR侧翼序列搜索比对方法。这种方法能够考虑错配、插入和缺失，并且同时搜索STR位点两个侧翼序列的算法，具体包括如下步骤：

[0065] 1、原始数据处理

[0066] 1)数据收集，二代测序仪输出的BCL格式的原始文件存放在目标路径，使用bcl2fastq将BCL文件拆分为fastq格式的文件并存放在指定路径。

[0067] 2)序列比对，整理所有STR基因座位置往上下游各延伸300bp作为参考序列库，使用BWA将fastq文件和参考序列库进行比对获取SAM格式的输出文件。

[0068] 2、在read中搜索侧翼序列：

[0069] 3)比对算法构建(分以下5步)：

[0070] a)已知左侧侧翼序列flank_left、右侧侧翼序列flank_right和待分析read序列S1，令x＝左侧侧翼序列length(即x为左侧侧翼序列的长度)，令y＝x+1,将flank_left、
flank_right按顺序拼接成序列S2。这里得到的x、y代表的就是侧翼序列连接处的位置。

[0071] b)初始化得分矩阵：以S1为列、S2为行构建得分矩阵，设置第1列惩罚得分都为0(即矩阵第一列的值都为0)；设置第1行侧翼序列连接处(即x和y的位置)惩罚得分为‑10，其
余位置惩罚得分为‑1(即矩阵第一行的值除了侧翼序列连接处为‑10，其余都为‑1)。

[0072] c)填充得分矩阵：假设匹配得分match_score＝1，错配得分mismatch_score＝‑1，gap得分gap_score＝‑2；对于矩阵中每个单元格(i,j)，其中的i代表矩阵行的索引，j代表
矩阵列的索引，得分计算公式为：

[0073] 若j≠x且j≠y:Fij＝max(Fi‑1,j‑1+S(Ai,Bj),Fi,j‑1+d,Fi‑1,j+d)；

[0074] 若j＝x:Fij＝max(Fi‑1，j‑1+2S(Ai,Bj),Fi,j‑1+0,Fi‑1,j+2×d)；

[0075] 若j＝y:Fij＝max(Fi‑1,j‑1+2S(Ai,Bj),Fi,j‑1+2×d,Fi‑1,j+0)；

[0076] 其中d＝gap_score，S(Ai,Bj)＝match_score或者mismatch_score。

[0077] d)获取重复区域序列：在最后一列由下往上选取该列第一个最大分值s作为比对得分，假设序列S2的长度为n，若n‑s≤4，则在允许2个错配、插入或缺失的情况下，S1和S2可
以比对上；以s的坐标为源头往左上回溯，可以得出j＝x和j＝y时行的坐标i1、i2,最后取S1
上i1、i2位置中间的序列即为STR位点重复区域序列。

[0078] 4)STR位点重复区域锚定：读取SAM文件的每条read，根据它是否能够比对上STR参考序列分为两种情况：1、若比对上STR位点，则使用该STR位点的侧翼序列和上述算法找到
STR重复区域；2、若未比对上任何STR位点，则遍历所有STR位点的侧翼序列，并且使用上述
算法找到STR重复区域；若比对的结果总共不超过2个错配、插入或者缺失，则输出比对上的
STR位点以及重复区域序列，否则舍弃这条read；最终获得每条read锚定的重复区域序列长
度。

[0079] 5)根据每条read锚定的重复区域序列长度，计算每条read比对上的STR位点的重复数。

[0080] 二、进一步的针对一个STR位点的所有的read及其重复数，本发明还涉及一种基于检测有效峰值，并且结合重复数差异进行动态阈值设定的STR分型的新方法，具体的：

[0081] 在侧翼序列的搜索比对结束后，以下是基于峰值检测和动态阈值的分型新方法：

[0082] 6)基于二代测序的单条read比对上的STR位点的重复数，统计样本每个STR位点上每种重复数检测到的reads数量(RC)。

[0083] 7)选取第一分型：根据峰值位置来选取分型，对于每个STR位点，RC最大的重复数A为峰，若RC(A)的值过低，则认为该位点没有检测出任何分型。

[0084] 8)遍历选取第二分型：为避免已经找到的峰对后续找峰的干扰，每次找到一个峰之后就将该峰的RC值赋0，即令RC(A)＝0；若同时满足以下两个条件则判定重复数B为峰：

[0085] a.计算待检测分型B的梯度阈值T，T＝初始梯度+0.25×(左侧梯度或者右侧梯度)×(A‑B)。其中的A代表检测到的第一分型，B代表正在检测的分型，ACR值代表B分型和A分型
的RC之间的比值。重复数B的RC要高于该分型的梯度阈值T。

[0086] b.重复数B的RC需要高于相邻重复数的RC。

[0087] 9)输出分型结果，若同时找到满足要求的第一分型和第二分型A、B，则分型结果为(A，B)；若找到满足要求的第一分型A但没有找到满足要求的第二分型，则分型结果为(A)；
若没有找到满足要求的第一分型和第二分型，则未检测到分型。

[0088] 下面为具体的实施方法。

[0089] 实施例1本发明方法体系构建

[0090] 本发明二代测序STR整体分析的流程如图1，流程上，首先对下机数据BCL文件拆分成fastq，接着进行read比对，得到每个read所属的STR位点。针对比对到某一个位点上的所
有的read，使用侧翼序列匹配，如果匹配则使用重复数计算公式计算重复区域包含的重复
数。随后统计每一个位点中的各个分型的read的数目，最后使用分型算法得到正确的分型。
其中关键的步骤是重复区域的检测和分型的确定，具体如下。

[0091] 1、STR侧翼序列的搜索/比对方法构建

[0092] 在使用侧翼序列匹配read的时候，侧翼序列可能包含错配、插入和缺失，例如在D13S317附近就有2个SNP。为了考虑这些因素，本发明改进了Needleman‑Wunsch算法(NW算
法)，使之能够同时搜索左侧和右侧侧翼序列，并且能够跨过中间的未知重复数的重复区
域。识别侧翼序列的方式的图2所示，图中的最上面是理想情况下read的结构和侧翼序列所
在的位置，中间部分为改进点，最下面是两个侧翼序列和read比对的一个例子。

[0093] 具体的，改进的Needleman‑Wunsch比对算法如下：

[0094] 已知左侧侧翼序列flank_left、右侧侧翼序列flank_right和待分析read序列S1，令x＝左侧侧翼序列length(即x为左侧侧翼序列的长度)，令y＝x+1,将flank_left、flank_
right按顺序拼接成序列S2。这里得到的x、y代表的就是侧翼序列连接处的位置。

[0095] 以一个例子演示修改后的Needleman‑Wunsch比对算法，假设flank_left＝TCA，flank_right＝TC,S1＝TCATTTC,则S2＝TCATC，x＝3,y＝4(具体参见图3，其表明了改进的
NW算法详细以及打分矩阵，图3左侧是事例打分矩阵，右侧是具体的对NW算法的改进)。

[0096] 初始化得分矩阵：第1列的惩罚得分为0；第1行侧翼序列连接处惩罚得分为‑10，其余位置惩罚得分为‑1；

[0097] 填充得分矩阵：假设匹配得分match_score＝1，错配得分mismatch_score＝‑1，gap得分gap_score＝‑2；对于每个单元格(i,j)，得分计算公式为：

[0098] 若j！＝x且j！＝y:Fij＝max(Fi‑1,j‑1+S(Ai,Bj),Fi,j‑1+d,Fi‑1,j+d)

[0099] 若j＝x:Fij＝max(Fi‑1,j‑1+2S(Ai,Bj),Fi,j‑1+0,Fi‑1,j+2×d)

[0100] 若j＝y:Fij＝max(Fi‑1,j‑1+2S(Ai,Bj),Fi,j‑1+2×d,Fi‑1,j+0)

[0101] 其中d＝gap_score，S(Ai,Bj)＝match_score or mismatch_score。具体填充出的得分矩阵如图3所示。

[0102] 获取重复区域序列：在最后一列由下往上选取该列第一个最大分值s作为比对得分，假设序列S2的长度n＝length(S2)，若n‑s≤4，则在允许2个错配、插入或缺失的情况下，
S1和S2可以比对上；以s的坐标为源头往左上回溯，可以得出j＝x和j＝y时i的坐标i1、i2,
最后取S1上i1、i2中间的序列(TT)即为STR位点重复区域序列。

[0103] STR位点重复区域锚定：读取SAM文件的每条read，根据它是否能够比对上STR参考序列分为两种情况：1、若比对上STR位点，则使用该STR位点的侧翼序列和上述算法找到STR
重复区域；2、若未比对上任何STR位点，则遍历所有STR位点的侧翼序列，并且使用上述算法
找到STR重复区域；若比对的结果总共不超过2个错配、插入或缺失，则输出比对上的STR位
点以及重复区域序列，否则舍弃这条read；最终获得每条read锚定的重复区域序列长度。

[0104] 根据每条read锚定的重复区域序列长度，计算每条read比对上的STR位点的重复数。

[0105] 2、STR分型方法构建

[0106] 在检测完该STR位点中每一个read的重复区域的长度后，下一步是找到正确的分型。传统的方法是查看read数量第二高分型占第一高的分型的read数量的百分比(ACR值)，
如果ACR值高于40％或者50％，一般就认为第二高的分型是一个正确的分型。然而，在毛细
管电泳和二代测序中，重复区域越长的分型其read数目或者荧光强度通常越低。而且，这一
趋势随着read测序长度的变化而变化，如图4所示。图4为测序读长对较长的分型的read数
量的影响，其中，测序读长和两个分型之间read数量的比值(ACR值)已经标记到图中。对于
同一STR位点的两个allele，测序读长越短，则read数量相差越大。

[0107] 为此，本发明开发了一种基于峰值的检测分型的算法，算法如图5所示，图5为二代测序通过峰值检测STR分型的算法示意图，其中，第二行图重点强调了当两个分型相邻的时
候，去掉之前找到的峰的作用。具体的，read最高的分型必为峰，下一步将最高峰的read数
赋值为0，这是为了避免第二高的正确分型和最高分型相邻的情况。下一步再重新寻找峰值
的位置，依次循环。

[0108] 但是由于测序和PCR的错误，通常在不正确的分型位置上也会有read，因此本发明加入read的数量作为阈值，对于某一个分型除了要为峰，而且read数量必须超过阈值。为了
考虑分型之间长度越大read数量相差越大的情况，本发明使用了梯度阈值，即阈值随着待
检测的分型和最高read数量的分型的重复数差异的变化而变化，如果待检测的分型的重复
数比最高分型的重复数越大，则阈值越低，反之则阈值越高。

[0109] 在二代测序中，由于PCR和测序的不均衡，经常导致某个样本或者某个位点的read数量不足，因此要考虑read数量小于某个阈值的时候排除掉这个位点。本发明使用实际样
本，找到比较理想的阈值为20。如图6所示，当read数量的阈值小于20的时候，不一致率变化
很大，而当read数量的阈值大于20的时候，不一致率的变化很小。其中检测到的分型数量是
随着read数量阈值的增大而逐渐的减小。

[0110] 除了read数量的阈值以外，另外的关键参数是左侧阈值梯度、右侧阈值梯度和初始梯度。假设read数量最高的分型为A，则某个分型B的阈值为：如果B>A，则阈值＝初始阈值
+(A‑B)乘以右侧阈值梯度；如果B使用几十个真实样本和三重循环遍历的方式确定了三个参数的最优值，并且使用地形图来
查看得到的三个阈值是否确实为最优值，结果如图7所示。图7中左图和右图中箭头所指的
位置即为左侧梯度阈值和右侧梯度阈值的最优值。图中的左侧梯度阈值和右侧梯度阈值都
乘以4来让三个值在大小上更相近。可以看出，其中得到的初始梯度的最优值为0.28，左侧
梯度的最优值为0.19，右侧梯度的最优值为0.33，图中箭头的位置即为所述的值的位置。

[0111] 在得到参数后，将具体的改进的分型方法描述如下，基于R的伪代码见表1：

[0112] 选取第一分型：根据峰值位置来选取分型，对于每个STR位点，RC最大的重复数A为峰，若该重复数的RC≥20,则将该重复数判断为第一分型，反之则判断该STR位点未检测到
任何分型。

[0113] 遍历选取第二分型：为避免已经找到的峰对后续找峰的干扰，每次找到一个峰之后就将该峰的RC值赋0，首先将第一分型的RC赋0，即令RC(A)＝0；假设重复数A、B的RC分别
为RC(A)、RC(B)，若同时满足以下两个条件则判定重复数B为峰：

[0114] a.计算待检测分型B的梯度阈值T，T＝初始梯度+0.25×(左侧梯度或者右侧梯度)×(A‑B)。通过优化参数得到的具体的公式为T＝0.28+0.25×step×(A‑B)，其中若A‑B>0，
step＝0.19(左侧梯度)；若A‑B<0，step＝0.33(右侧梯度)。已知ACR＝RC(B)/RC(A)，需同时
满足ACR≥T、ACR≥0.02和RC(B)≥3；其中的A代表检测到的第一分型，B代表正在检测的分
型，ACR值代表B分型和A分型的RC之间的比值。

[0115] b.重复数B的RC需要高于相邻重复数的RC。假设重复数B左右相邻的整数重复数分别为B_left、B_right，相邻重复数对应的RC分别为RC(left)、RC(right)，其中0＜B‑B_left
≤1,0RC(left)和
RC(B)>RC(right)。

[0116] 如果找到第二分型，则输出分型为(A，B)；如果只找到一个分型，则最终的分型为A；若未找到任何分型，则没有找到任何分型。

[0117] 表1：二代测序检测STR分型的具体算法的伪代码，基于R语言。

[0118]

[0119]

[0120] 实施例2样本检测测试

[0121] 在构建完算法以及得到最优的参数体系后，本发明分别使用标准样本和真实样本来对本发明的方法进行测试。具体的，使用标准阳性样本9948和9947细胞系，以及70多组的
真实样本进行测试。真实样本测试的阳性对照为毛细管电泳(CE)的结果。具体的步骤包括
首先是STR位点的二代建库以及高通量测序；下机后的数据拆分、序列比对；最后应用本发
明提出的方法进行重复数的检测和STR的分型，以旧方法(通过ACR值直接进行分型)作对
照。

[0122] 对于标准样本的检测结果：

[0123] 使用新方法针对9948细胞系的66个位点检测STR分型，NGS和CE的一致率为100％，而使用旧方法检测的一致率为98％。针对9947细胞系的37个位点，NGS和CE的一致率为
99.9％，使用旧的方法为99.1％。

[0124] 对于真实样本的检测结果：

[0125] 使用新的侧翼序列搜索方法和新的检测STR分型方法针对73个真实样本66个位点的一致率为99.6％，其中少量不一致的部分有NGS的原因，也有少部分是CE的原因。而使用
旧方法检测的一致率为93.8％。具体结果如表2和表3所示。

[0126] 表2：真实样本使用不同方法在部分STR位点上得到的结果的区别。NGS和CE一致的结果用灰色框表示。

[0127]

[0128] 表3：真实样本中使用不同的方法得到的NGS和CE的分型一致率。(a/b)的格式为：

[0129] a是一致的样本‑位点对，b为总的样本‑位点对。

[0130]

[0131] 由此可见，本发明提出的方法相比于传统方法具有很强优势。无论是单独采用本发明单独的新重复数检测方法或新STR的分型，还是联用本发明两个方法，分析一致率都显
著提高。尤其联用本发明两个方法以后，一致率从93.8％增加到了99.6％。而且在一些STR
位点上，例如D13S317，同样看到了明显的效果的改善。

[0132] 前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变
和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应
用，从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及
各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。