CN115035957A - 基于粒子群算法的改进最小残差法分析混合str图谱 - Google Patents

Abstract

一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法：S100：输入来自两个供体的混合STR图谱；S200：在(0，0.5]区间内随机生成3个粒子的初始位置，并在[‑0.01，0.01]区间内随机生成所述3个粒子的3个初始速度；S300：分别根据所述3个粒子计算得出3个residual_sum；S400：分别计算3个粒子各自的个体最优位置，获得residual_sum最小时Mx的值；S500：计算整个粒子群体的全局最优位置，获得residual_sum最小时Mx的值；S600：更新所述3个粒子的速度和位置；S700：如果迭代次数达到阈值，则获得residual_sum最小值、混合比例预测值Mx’和拆分结果；否则返回到步骤S300继续执行。本方法改进了最小残差法，重新定义混合比例实现等位基因模型(allele model)的化简；粒子群算法的引入，能够实现混合STR图谱的快速分析。

Description

基于粒子群算法的改进最小残差法分析混合STR图谱

技术领域

本公开属于法医遗传学和法医物证学技术领域，特别涉及一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法。

背景技术

在法医遗传学领域对混合STR图谱分析的研究一直是难点和热点。混合DNA常见于刑事案件，案件现场会采集到两个或者更多个体的混合血迹或脱落上皮细胞的混合物等。对于混合样本的DNA分型，目前法医DNA鉴定常规使用短串联重复(short tandem repeats，STR)检测技术，应用荧光STR图谱辨识出各个供者的DNA分型，对混合检材结果统计的理论方面已经有较成熟的研究。STR等位基因的峰高和面积信息能够用于分析混合检材中的基因分型。

目前，国内大部分实验室对于混合STR图谱的拆分都是根据每个位点下峰个数、峰高等参数进行人工拆分。然而，人工拆分只是单个位点逐个的去拆分，没有考虑整个图谱；并且拆分时没有量化的标准，主观因素较大，导致不确定性较大。而国际上已有专门的拆分软件，主要是利用统计学的方法去解释混合图谱的数据。软件拆分的优势在于：有可量化的标准、可重复性好、减少主观性，结果更加客观公正。

当前国际拆分软件按照分析模型划分主要有三类：一类是二进制法(Binary)，主要思想是设定阈值来处理峰的随机性以及去除一些不合适的数据，但是这种方法在处理低拷贝检材、降解检材以及检材混合比例差异较大时结果不理想；第二类是半连续法(Semi-continuous)，是概率法的一种，该法没有考虑峰值的变异性、混合比例及Stutter峰所占百分比等因素，在拆分的合理性、准确性上较差；第三类是连续法(Continuous)，也称概率法，采用马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo，MCMC)计算方法，利用计算机模拟分析所有图谱数据，以概率的形式给出可能性的组合，但这种方法不能解决维度灾难。

因此，现有技术中的软件存在计算量大或结果不尽如人意的问题。

发明内容

为了解决上述技术问题，本公开揭示了一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法，其包括如下步骤：

所述方法用于法医DNA鉴定的STR图谱；

所述方法包括如下步骤：

S100：输入来自两个供体的混合STR图谱，其包括基因座、等位基因和峰面积；

S200：在(0，0.5)区间内随机生成3个粒子的初始位置，并在[-0.01，0.01]区间内随机生成所述3个粒子的3个初始速度，以此执行粒子群算法，其中，所述随机生成的3个粒子的初始位置还作为混合比例Mx的三个初始取值；

S300：根据粒子群算法计算出3个粒子各自对应的、所有基因座的残差的最小值之和residual_sum；

S400：根据3个粒子所对应的各个residual_sum的最小值，计算3个粒子各自的个体最优位置，并以此作为混合比例Mx的3个个体极值；

S500：根据3个粒子所对应的各个residual_sum的最小值，进一步获得其中最小的一个residual_sum，以此作为整个粒子群体的全局最优位置，并将该全局最优位置作为混合比例Mx的全局极值；

S600：3个粒子中的每一个粒子，通过跟踪个体极值和全局极值来更新自己在解空间中的位置与速度，最终找到全局最优；其中，若粒子群算法的迭代次数已达到阈值，则获得residual_sum最小值、混合比例预测值Mx和建立在二者之上的基因座的拆分结果；

S700：若粒子群算法的迭代次数未达到阈值，则更新所述3个粒子的速度和位置，迭代执行步骤S300至S600。

通过上述技术方案，本方法通过改进后的最小残差法，以及重新定义混合比例实现等位基因模型(allele model)的化简；粒子群算法的引入，能够实现混合STR图谱的快速分析。

附图说明

图1是本公开一个实施例中所提供的一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法流程图；

图2是本公开一个实施例中所提供的粒子群算法residual_sum变化趋势图；

图3是本公开一个实施例中所提供的Mx在[0.01，0.5]内residual_sum的计算结果图；

图4是是本公开一个实施例中所提供的10次迭代的目标函数值变化趋势图。

具体实施方式

为了使本领域技术人员理解本公开所披露的技术方案，下面将结合实施例及有关附图1至图4，对各个实施例的技术方案进行描述，所描述的实施例是本公开的一部分实施例，而不是全部的实施例。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本公开的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其他实施例互斥的独立的或备选的实施例。本领域技术人员可以理解的是，本文所描述的实施例可以与其他实施例相结合。

参见图1，在一个实施例中，本公开揭示了一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法，其包括如下步骤：

所述方法用于法医DNA鉴定的STR图谱；

所述方法包括如下步骤：

S100：输入来自两个供体的混合STR图谱，其包括基因座、等位基因和峰面积；

S200：在(0，0.5)区间内随机生成3个粒子的初始位置，并在[-0.01，0.01]区间内随机生成所述3个粒子的3个初始速度，以此执行粒子群算法，其中，所述随机生成的3个粒子的初始位置还作为混合比例Mx的三个初始取值；

S300：根据粒子群算法计算出3个粒子各自对应的、所有基因座的残差的最小值之和residual_sum；

S400：根据3个粒子所对应的各个residual_sum的最小值，计算3个粒子各自的个体最优位置，并以此作为混合比例Mx的3个个体极值；

S500：根据3个粒子所对应的各个residual_sum的最小值，进一步获得其中最小的一个residual_sum，以此作为整个粒子群体的全局最优位置，并将该全局最优位置作为混合比例Mx的全局极值；

S600：3个粒子中的每一个粒子，通过跟踪个体极值和全局极值来更新自己在解空间中的位置与速度，最终找到全局最优；其中，若粒子群算法的迭代次数已达到阈值，则获得residual_sum最小值、混合比例预测值Mx和建立在二者之上的基因座的拆分结果；

S700：若粒子群算法的迭代次数未达到阈值，则更新所述3个粒子的速度和位置，迭代执行步骤S300至S600。

就该实施例而言，本方法改进了最小残差法，重新定义混合比例实现等位基因模型(allele model)的化简；粒子群算法的引入，能够实现混合STR图谱的快速分析，有较高的准确率且运行时间短。

混合STR图谱是混合DNA样本经过一些专业仪器检测生成的数据，通过混合STR图谱可以获取混合DNA样本的基因信息(检测位点基因座上有哪些等位基因)。混合DNA样本含有的等位基因可以通过STR图谱得到，等位基因的峰面积是DNA模板量的映射，法医实验已证实，DNA的模板量与STR图谱中的等位基因峰面积存在线性关系，即，样本中的DNA模板量越大，该样本生成的STR图谱上的等位基因的峰面积越大。

一般来说，每个STR基因座有两个等位基因，基因型为纯合子时在图谱上表现为一个峰，基因型为杂合子时在图谱上表现为两个峰。但是对于两人混合STR图谱来说，在一个基因座上可以出现1个等位基因，2个等位基因，3个等位基因以及4个等位基因。对于杂合子来说，两个等位基因的峰面积应大致相等，比例接近于1。基于此原理，分析者可以通过观察到的峰面积来推断可能的基因型组合和各组分的比例。对确定是二供体组成的混合样本来说，分型结果可以表现为一带型、二带型、三带型或四带型。除一带型外，其余各种表现均可提供基因型组合的信息。其中，四带型由于直观地表现出2个供体的基因型，所以最易于通过其推断2个供体的基因型组合及其比例。

对于任意两个人的混合物，基因座上有1、2、3和4个等位基因时可能的组合数分别为1、7、12、6，详见表1。

表1

其中，a，b，c，d分别代表基因座上的等位基因。

本发明中，分析混合STR图谱的第一步是估计混合比例(Mx，Mixtureproportion)，混合比例是指在混合样本中含量最小的供体在混合样本中的比例，此时Mx的取值范围为(0，50％)。当基因座上只有一个等位基因时，可能的基因组合只有一个，即两个体的基因型均为纯合子aa，因此一带型基因座可以直接得到拆分结果，后面的拆分方法对一带型基因座不需要做讨论；当基因座上有四等位基因时，估计出混合比例Mx比较容易做到的。然而，如果基因座上包含三个或更少的等位基因，这个过程就不再简单。法医学实验证明混合比例Mx在混合STR图谱的所有基因座上是一致的。所以一旦知道了混合比例，就可以建立等位基因模型来估计表1中基因组合的期望峰面积，因为Mx为含量最小的供体在混合样本中的比例，因此可以表1中的基因组合部分可以排除不考虑。

设在混合样本中含量最小的供体为个体1，在二组分混合样本中，另一个供体为个体2；

然后，对等位基因进行如下编号；

将每个基因座中的等位基因按照其峰面积大小降序排列，对于四带型基因座，等位基因分别标记为a，b，c，d，即

对于三带型基因座，等位基因分别标记为a，b，c，即

对于二带型基因座，等位基因分别标记为a，b，即

等位基因模型如表2至表4所示，表2为四等位基因模型，表3为三等位基因模型，表4为二等位基因模型。

表2

表3

表4

由表2至表4计算出的等位基因峰面积比率，可以得到等位基因峰面积的期望值，如果基因组合推断正确，等位基因的峰面积期望值与观察值应该是近似的。使用期望值与观察值的残差residual来量化这种近似程度(详见下式)，residual越小，当前基因组合的为真实结果的概率越大：

其中，

表示STR图谱中第i个基因座在基因型组合为com时的残差值，n是基因座i的等位基因数量，

分别表示等位基因j的峰面积观察值与期望值。此时基因座i的推测基因型结果combination_i为：

基因座i的residual的最小值residual_i为：

粒子数量也是粒子群的初始化参数，根据经验，3是在多次实验后选择的一个合适的值(数值太大算法运行时间会太长，数值太小不能保证准确率，3是在保证准确率的情况下选择的一个最小的数值)。在本方法中，引入了粒子群算法，那么混合比例Mx(0＜Mx＜50％)的数值可以视为粒子群算法中的粒子的位置，3个粒子的位置，即Mx的三个取值，随机生成方式是：在0到0.5之间随机生成3个随机数赋值给这三个粒子作为它们的初始位置。

在另一个实施例中，步骤S200中速度的初始化具体为：

惯性权值为0.5，学习因子为2，迭代次数为10。

在另一个实施例中，步骤S700中如果达到迭代次数阈值，则获得residual_sum最小值、Mx′和拆分结果具体为：

在Mx的取值范围中，即(0，0.5)，寻找使得整个图谱的所有基因座的最小残差值之和最小时的数值，把这个数值赋值给Mx′，再根据Mx′分别计算每个基因座的残差的最小值residual_i，取残差值为residual_i时的基因组合作为拆分结果。

基因座i的残差的最小值residual_i为：

整个图谱所有基因座的最小残差值之和residual_sum为：

residual_sum取值最小时，混合比例Mx的值记为Mx’，有：

其中，n为STR图谱中基因座的数量，residual_sum为整个STR图谱中所有基因座的残差最小值之和，同样，residual_sum越小，当前Mx′越接近真实的混合比例，混合比例为Mx′时的拆分结果正确的概率越大，

此时基因座i的基因结果combination_i为

本方法在Mx的取值范围中，寻找使得整个图谱的residual最小的Mx，把这个Mx的值赋值给Mx′，并把这个值代入等位基因简化模型中(表2至表3)，分别计算每个基因座的在基因组合com的

取使得

最小时的基因组合com作为拆分结果。

在另一个实施例中，步骤S300中的计算具体为：

其中，

表示STR图谱中第i个基因座在基因型组合为com时的残差值，n是基因座i的等位基因数量，

分别表示等位基因j的峰面积观察值与期望值；

此时基因座i的基因结果combination_i为

基因座i的残差的最小值residual_i为

就该实施例而言，对于公式

中的Mx′，最简单的计算方法是Mx在0到50％的范围依次取1％，2％，...49％分别代入公式

计算residual_sum，找到使得residual_i最小的Mx。如果使用粒子群算法则可以快速计算出结果。

因为0＜Mx＜50％，所以可以将公式

做为粒子群算法的目标函数f((x)。在进行边界条件处理后，此时粒子群算法的目标函数为：

通过粒子群算法可以快速找到使得这个目标函数residual_sum最小时的Mx的值，即Mx′的值。

在另一个实施例中，步骤S600进一步包括：

速度更新公式为：

粒子位置更新公式为：

X_i(t+1)＝X_i(t)+v_i(t+1)

其中，粒子i的取值范围是1到n，n为粒子数量，t为当前迭代次数，t＝0表示初始化阶段，即第1次迭代前，ω为惯性权值，r₁，r₂为0到1之间的随机数，X_i(t)为粒子i在第t次迭代时的位置，X_i(0)表示粒子i初始化位置，

为粒子i当前的个体最优位置，X^gb为粒子群全局最优位置，v_i(t)为粒子i在第t次迭代时的速度，v_i(0)为粒子i初始化速度，c₁，c₂是学习因子。

就该实施例而言，X_i(0)与v_i(0)是随机数，在算法初始化时需分别给出随机数的取值范围。

粒子群优化算法(PSO：Particle swarm optimization)是一种进化计算技术(evolutionary computation)。源于对鸟群捕食的行为研究。粒子群优化算法的核心思想是在有限空间中创建n个粒子，每个粒子单独搜寻最优解并将最优解由整个粒子群共享，从而达到优化的目的。

粒子群算法可以快速找到粒子的最优位置，即粒子在这个位置上对应的目标函数值最优。

设粒子群算法要优化的目标函数为f(x)，算法的目标是找到f(x)最优值f(x)_best以及目标函数取最优值时x的取值x_best，那么

f(x)_best＝f(x_best)

最优值可以选择取最大值或者最小值，选取什么作为最优值需结合具体的问题来决定，其中argbest就是使相应函数达到最优值时的自变量X的取值。

粒子i每次迭代更新的位置均记录在集合P_i中，经过t次迭代后P_i为

P_i＝{X_i(o)，X_i(1)，X_i(2)，…，X_i(t)}

那么粒子i的个体最优位置为

粒子群中的每个粒子每次迭代更新的位置均记录在集合G中，经过t次迭代后G为

其中，n为粒子的数量。

那么粒子群全局最优位置为

其中，arg是变元(即自变量argument)的英文缩写。arg min就是使后面这个式子达到最小值时的变量的取值，arg max就是使后面的函数达到最大值时的变量的取值，其中arg best就是使后面函数达到最优值时的自变量X的取值。

综上，在混合STR图谱分析中使用粒子群算法，可以快速找到使得residual_sum取值最小时的混合比例Mx的值，即Mx’的值，这混合比例的取值接近真实混合比例的概率最大，在这个取值下获得的拆分结果正确的概率也最大。

另一个实施例中，实验数据是来自两个个体的混合样本，使用GeneMarker HIDV3.0.0对21对STR位点进行基因型分析，样本的基因座、等位基因及峰面积如表5所示，由于数据的私密性，只展示了部分的位点信息及分析结果。表5为来自两个体混合的部分位点信息。

表6为改进最小残差法引入粒子群算法前后对比。

表5

使用粒子群算法后，经过10次迭代得到了residual_sum的最小值为0.025792，如图2所示。由表6知，粒子群算法可以在不影响分析结果的情况下，计算时间减少了近三成。

表6

R语言有个用于法医遗传学DNA混合物分类的工具包mixsep，使用mixsep对表8数据进行了分析，用来和本方法进行比较，如表7所示。表7中也列出了表5数据的真实分型结果。可以发现，改进的最小残差法和mixsep，除了基因座TH01个体1拆分有误，其余分型结果均与真实分型一致。由表5可知，基因座D3S1358和TH01均为二带型基因座，它们的真实分析结果均为两个纯合子，它们的杂合均衡比分别为0.4215、0.2892，具有很大差异，因此可以推断样本采集过程和图谱生成过程等外界因素对基因座TH01上的峰面积产生了较大的影响，等位基因模型生成的峰面积期望值难以接近观察值。

表7

利用STR图谱的定量信息峰面积可对混合DNA样本进行分析，实现了通过计算机对混合STR图谱分析。本方法对二组分混合DNA样本生成的STR图谱进行分析，可推断混合物中个组分的比例，以及各组分的分型结果，具有较高的准确率和较快的计算速度，可以实现对大量混合STR图谱的批量分析，协助法医工作者进行图谱分析。

另一个实施例中，表8是一个两人混合的STR图谱数据

表8

对比例一、不使用粒子群算法加速进行分析

若不使用粒子群算法，则可以将Mx依次赋值0.01，0.02，...，0.50，可以计算出50个residual_sum的值。根据公式

找出使得residual_sum最小的Mx的值，即Mx′。

首先，将表8中每个基因座的等位基因峰面积降序排序，并进行归一化(归一化的目的是便于计算)。

如基因座vWA，将它的等位基因按照峰面积降序排列结果为1318(15)、1200(19)、793(18)、621(16)。括号里为对应的等位基因名称。

归一化过程如下：

该基因座的峰面积之和为1318+1200+793+621＝3921。

那么等位基因15的归一化结果为

等位基因19的归一化结果为

等位基因18的归一化结果为

等位基因16的归一化结果为

同理可以得到其他基因座的归一化结果，如表9所示。归一化结果可以视为公式

的

表9

接下来将Mx依次赋值0.01，0.02，...，0.50，可以计算residual_sum的值。以Mx＝0.20为例(即混合样本中含量最少的供体在混合样本中的比例为20％)。

对于四等位基因基因座vWA，将等位基因15、19、18、16分别编号为a，b，c，d，参考表5可以得到vWA的期望值

如表10所示。

表10

将表9中vWA的等位基因观察值和表10中期望值代入公式

中可得：

表：10中只有(cd，ab)一种基因型组合，由公式

可知：

combination_vWA＝(cd，ab)

即，在Mx＝0.2时，基因座vWA的推断分析结果为(cd，ab)，即(18/16，15/19)。

同理，可以得到其他四等位基因的residual值以及推断分析结果，如表11所示。

表11

对于三等位基因基因座D3S1358，将等位基因15、18、16分别编号为a，b，c，参考表3可以得到D3S1358的期望值

进一步，参见表12：

表12

将表9中D3S1358的等位基因观察值和表12中的每一行期望值代入公式

中计算residual_i，如表13所示。

13

由表13可知，residual_D3S1358的最小值为0.00638，由公式

可得：

combination_D3S1358＝(cc，ab)

(cc，ab)即(16/16，15/18)为基因座D3S1358在Mx＝0.20时的推断分析结果。

上表中，以表15第一行基因型组合(bc，aa)即(18/16，15/15)为例，

同理，可以得到其他三等位基因在Mx＝0.20时的最小residual值及推断分析结果，如表14所示。

表14

对于二等位基因基因座D5S818，将等位基因12、13分别编号为a、b，参考表4可以得到D5S818的期望值

如表15所示。

表15

将表9中D5S818的等位基因观察值和表15中的每一行期望值代入公式

中计算residual_i，如表16所示。

其中，以表18第一行基因型组合(ab，aa)即(12/13，12/12)为例，

表16

由表16可知，residual_D5S818的最小值为0.004398，由公式

可得

combination_D5S818＝(bb，aa)

(bb，aa)即(13/13，12/12)作为基因座D5S818在Mx＝0.20时的推断分析结果。

同理，可以得到其他二等位基因在Mx＝0.20时的最小residual值及推断分析结果，如表17所示。

表17

分析完所有的等位基因可以得到Mx＝0.2时的推断分析结果，如表18所示。

表18

表18中最后一行为所有基因座的最小residual之和，即当Mx＝0.2时，residual_sum＝0.133323。

根据上述方法，同样可以得到Mx取其他值时的residual_sum。根据公式

Mx存在一个取值Mx′，使得residual_sum最小。本文把Mx＝Mx′时的推断结果，作为最终的推断结果。

此外，Mx依次取0.01，0.02，...，0.50进行上述计算方法，绘制成图像如图3所示。

在图3中，黑点为图像的纵坐标的最低点。那么，Mx取值0.3时，即Mx′＝30％，residual_sum最小，最小为0.300111。那么Mx＝0.3时的推断结果为最终的结果，Mx′＝30％为表11示例数据对应的混合样本中含量最少的供体在混合样本中的比例。最终分析结果如表19所示，表19中最后一行为所有基因座的最小residual之和。

表19

相比前面的对比例，本发明的一个实施例中，引入粒子群算法加速计算：

首先初始化(t＝0)粒子群算法参数：

粒子数n＝3，粒子初始速度v(0)取值范围为[-0.01，0.01]，惯性权值ω＝0.005，学习因子c₁＝c₂＝2，迭代次数为10次(，这些参数值是经过多次实验后得出的最佳取值，此为经验值)。

粒子初始位置X(0)的取值范围为(即Mx的取值范围)[0.0001，0.5]，(这个参数是根据具体的应用问题选取的)。

r₁，r₂为每次迭代生成的0到1之间的随机数。

目标函数f(x)为公式

在这里，粒子群算法的目标是求上述公式的最小值。

在[0.0001，0.5]区间内随机生成n＝3个粒子的位置X_i(0)(i＝1，2，3)(即Mx的三个取值)0.31088753，0.04369095，0.07871866。根据上述方法，分别计算出3个粒子的目标函数值，即residual_sum值，分别为0.03234149，0.73296422，0.53499878。此时3个粒子个体最优位置X^pb分别为0.31088753，0.04369095，0.07871866，最优目标函数值f(X^pb)(即residual_sum值)分别为0.03234149，0.73296422，0.53499878。此时，0.03234149是三个粒子residual_sum值中的最小值，所以粒子群全局最优位置X^gb为0.31088753，粒子群全局最优目标函数值f(X^gb)(即Mx＝0.31088753时，residual_sum值)为0.03234149。

在[-0.01，0.01]区间内随机生成n＝3个粒子的初始速vi(0)(i＝1，2，3)2.15070013e-03，7.01666399e-03，-7.67238216e-05。

初始化(t＝0)阶段参数列表如表20所示。

表20

第1次迭代：

在0到1之间分别生成3对随机数r₁，r₂，分别为0.38859983、0.90158338，0.09569076、0.04808691，0.26583436、0.80850666。

将第1对随机数0.38859983、0.90158338即(r₁＝0.38859983，r₂＝0.90158338)和上次迭代后的参数值(即表23中的ω，c₁，c₂，v₁(0)，X₁(0)，

X^gb)代入公式

和X_i(t+1)＝X_i(t)+v_i(t+1)计算可得粒子1在第1次迭代后的速度v₁(1)和位置X₁(1)。

v₁(0)＝1.07535007e-05

X₁(1)＝0.31089828

将第2对随机数0.09569076、0.04808691即(r₁＝0.09569076，r₂＝0.04808691)和上次迭代后的参数值(即表23中的ω，c₁，c₂，v₂(0)，X₂(0)，

X^gb)代入公式计算可得粒子2在第1次迭代后的速度v₂(1)和位置X₂(1)。

v₂(1)＝2.57323983e-02

X₂(1)＝0.06942334

将第3对随机数0.26583436、0.80850666即(r₁＝0.26583436，r₂＝0.80850666)和上次迭代后的参数值(即表23中的ω，c₁，c₂，v₃(0)，X₃(0)，

Xgb)代入公式计算可得粒子3在第1次迭代后的速度v₃(1)和位置X₃(1)。

v₃(1)＝3.75419763e-01

X₃(1)＝0.45413843

将X₁(1)、X₂(1)、X₃(1)按照附录2中的计算方法计算目标函数f(X₁(1))、f(X₂(1))、f(X₃(1))，即Mx＝0.31089828、0.06942334、0.45413843时，residual_sum值。

f(X₁(1))＝0.03234514

f(X₂(1))＝0.5834667

f(X₃(1))＝0.20163704

本发明要找的是目标函数的最小值，即residual_sum最小值。由表20知，

则第1次迭代后，

的值0.31088753依然是粒子1的个体最优位置

则第1次迭代后，X₂(1)的值为粒子2的个体最优位置

即

则第1次迭代后，X₃(1)的值为粒子3的个体最优位置

即

因为

则第1次迭代后，X^gb的值0.31088753依然是粒子群的全局最优位置X^gb(再一步是将

比较大小，若f(X^gb)数值最小，则X^gb不变；若存在

数值最小，则将X^gb更新为

的值)。

经过第1次迭代后粒子群算法参数值如表21所示。

表21

用上述同样的计算过程计算第2次迭代后的参数值，直到第10次迭代结束。经过10次迭代后粒子群算法参数值如表22所示。

表22

10次迭代的目标函数值residual_sum变化趋势图如图4所示。

本发明将表22中X^gb的值视为使residual_sum最小的Mx的值，即Mx’＝X^gb＝0.30042744，将Mx′＝0.30042744时的分析结果作为最终的分析结果。

可以发现，0.30042744与方法一中计算的Mx′＝0.3很接近，但是不用粒子群算法则需要进行至少50次residual_sum的计算，使用粒子群算法后只需计算33次，减少近四成的计算量。

最后，需要说明的是，本领域的普通技术人员在本说明书的启示下和在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。

Claims (9)

1.一种基于粒子群算法的改进最小残差法分析混合STR图谱的方法，其特征在于：

所述方法用于法医DNA鉴定的STR图谱；

所述方法包括如下步骤：

S100：输入来自两个供体的混合STR图谱，其包括基因座、等位基因和峰面积；

S200：在(0，0.5)区间内随机生成3个粒子的初始位置，并在[-0.01，0.01]区间内随机生成所述3个粒子的3个初始速度，以此执行粒子群算法，其中，所述随机生成的3个粒子的初始位置还作为混合比例Mx的三个初始取值；

S300：根据粒子群算法计算出3个粒子各自对应的、所有基因座的残差的最小值之和residual_sum；

S400：根据3个粒子所对应的各个residual_sum的最小值，计算3个粒子各自的个体最优位置，并以此作为混合比例Mx的3个个体极值；

S500：根据3个粒子所对应的各个residual_sum的最小值，进一步获得其中最小的一个residual_sum，以此作为整个粒子群体的全局最优位置，并将该全局最优位置作为混合比例Mx的全局极值；

S600：3个粒子中的每一个粒子，通过跟踪个体极值和全局极值来更新自己在解空间中的位置与速度，最终找到全局最优；其中，若粒子群算法的迭代次数已达到阈值，则获得residual_sum最小值、混合比例预测值Mx和建立在二者之上的基因座的拆分结果；

S700：若粒子群算法的迭代次数未达到阈值，则更新所述3个粒子的速度和位置，迭代执行步骤S300至S600。

2.如权利要求1所述的方法，其中，优选的，

混合比例Mx的取值范围为(0，50％)。

3.如权利要求1所述的方法，其中，步骤S300中，

其中，

表示混合STR图谱中第i个基因座在基因型组合为com时的残差值，n是基因座i的等位基因数量，

分别表示等位基因j的峰面积观察值与期望值；

基因座i的拆分结果combination_i为：

基因座i的残差的最小值residual_i为：

4.如权利要求3所述的方法，其中，对每个基因座的等位基因峰面积降序排序后对基因座进行归一化之后的结果作为

采用简化后的等位基因模型获取

5.如权利要求4所述的方法，其中，所述简化后的等位基因模型中重新定义了Mx为混合样本中含量最少的供体在混合样本中的比例，并将每个基因座中的等位基因按照其峰面积大小降序排列。

6.如权利要求1所述的方法，其中，

粒子i每次迭代更新的位置均记录在集合P_i中，经过t次迭代后P_i为

P_i＝{X_i(o)，X_i(1)，X_i(2)，…，X_i(t)}

那么粒子i的个体最优位置X_i ^pb为

其中arg best就是使f(x)达到最优值时的自变量X的取值，目标函数f(x)为公式

7.如权利要求1所述的方法，其中，

粒子群中的每个粒子每次迭代更新的位置均记录在集合G中，经过t次迭代后G为

其中，n为粒子的数量，那么粒子群全局最优位置X^gb为

其中arg best就是使f(x)达到最优值时的自变量X的取值，目标函数f(x)为公式

8.如权利要求1所述的方法，其中，

速度更新公式为：

粒子位置更新公式为：

X_i(t+1)＝X_i(t)+v_i(t+1)

其中，粒子i的取值范围是1到n，n为粒子数量，t为当前迭代次数，t＝0表示初始化阶段，即第1次迭代前，ω为惯性权值，r₁，r₂为0到1之间的随机数，X_i(t)为粒子i在第t次迭代时的位置，X_i(0)表示粒子i初始化位置，

为粒子i当前的个体最优位置，X^gb为粒子群全局最优位置，v_i(t)为粒子i在第t次迭代时的速度，v_i(0)为粒子i初始化速度，c₁，c₂是学习因子。

9.如权利要求1所述的方法，其中，

如果达到迭代次数阈值，则获得residual_sum最小值、混合比例Mx的预测值Mx′和拆分结果具体为：

在Mx的取值范围中，寻找使得整个图谱的所有基因座的最小残差值之和最小时的数值，将该数值赋值给Mx′，再根据Mx′分别计算每个基因座的残差的最小值residual_i，取残差值为residual_i时的基因组合作为拆分结果，

其中，n为STR图谱中基因座的数量，residual_sum为整个STR图谱中所有基因座的残差最小值之和，同样，residual_sum越小，当前Mx′越接近真实的混合比例，混合比例为Mx′时的拆分结果正确的概率越大，

此时基因座i的拆分结果combination_i为：

Images