CN115762636A - 一种动态突变拷贝数的方法和系统 - Google Patents

Abstract

本发明提供一种动态突变拷贝数的方法和系统。所述方法包括提取下机数据中的荧光信号矩阵，并进行基线对齐，先检测分子量内标信号峰，根据内标大小使用不同阈值过滤合适的峰，通过距离分类序列匹配分子量内标，并逐步校正匹配后的片段大小，来识别分子量内标。再检测待测样本信号峰，根据分子量内标计算目标峰的片段大小，并根据基线计算片段的拷贝数，最后筛选可能的等位基因，提示样本的基因型。由上，本发明构建的方法和系统可以快速检测并筛选可能的动态突变等位基因峰，直接输出拷贝数，并且适用大小不同的分子量内标，能并行分析多个下机结果，实现流程的高效和可控。

Description

一种动态突变拷贝数的方法和系统

技术领域

本发明涉及罕见病基因检测领域，具体涉及一种从罕见病患者基因的片段毛细管电泳结果中快速计算动态突变拷贝数的方法和系统。

背景技术

动态突变是DNA中的核苷酸重复序列的拷贝数发生扩增而产生的突变，在人类遗传病中主要为三个核苷酸的重复。动态突变的拷贝数会伴随着世代的传递而增加，在达到一定的倍数后，就会导致疾病的发生，发病率和疾病的严重性也会逐代加重。由动态突变导致的罕见病有脊髓小脑共济失调（以下简称SCA）、脆性X综合征、亨廷顿舞蹈病等，如SCA的亚型脊髓小脑共济失调1型（简称SCA1），就是由ATXN1基因的Chr6:16327867-16327953（Hg19）区域的CAG三核苷酸拷贝数超出正常值导致，正常范围是6至39个拷贝（即该区域有6至39个CAG片段重复），致病范围为40至88个拷贝（即该区域有40至88个CAG片段重复），它是一种常染色体显性遗传病，会导致周围神经损害等症状，严重地影响患者和后代的生活质量。确诊动态突变导致的疾病，即准确检测动态突变拷贝数，对于疾病正确的诊断治疗和指导优生优育都有重要的意义。

目前检测动态突变主要采用基于毛细管电泳的片段分析法，即先通过对动态突变核心区域附近设计引物进行PCR，该PCR引物上带有特定的荧光标记，获得带荧光的片段后，和带有另一种荧光的分子量内标（Size Standard）经过毛细管电泳分离，通过与分子量内标比较来确定DNA片段的相对分子量大小和基因型。目前主要是通过界面式软件进行数据分析，如付费GeneMapper或GeneMarker或开源的OSIRIS等，这些软件能够处理下机的原始数据，并计算待测样本包含的片段大小，这些片段在结果中以峰的形式体现，实验人员需要通过峰图人工判断和选择可能的等位基因峰，并通过预先构建的阴性样本基线来计算突变拷贝数。

这种分析流程存在一定的局限性。首先，动态突变的等位基因峰在PCR过程中，本身就会因为DNA聚合酶的特性，导致在其前后3个碱基内出现非等位基因峰（如延伸过程中DNA聚合酶滑动导致两端出现Stutter峰，或者DNA聚合酶会在3’末端添加一个A碱基形成黏性末端导致Plus-A峰）。再加上如果待测峰中噪声很多，可能影响分析人员的判断，导致得到错误的拷贝数结果。其次，如果动态突变拷贝数较多，超出了分子量内标的检测上限，它的信号峰值会远远低于分子量内标检测范围内的峰值，如果分析人员没有放大坐标轴查看，可能会忽略掉这种等位基因峰，导致漏检了致病的等位基因，得到假阴性的结果。此外，界面式的分析流程较难进行流程的搭建和自动化，如果需要批量分析大量的样本，每个样本又分析多个文件的时候，工作量非常大，结果也无法质控。现有主流是界面式的分析软件，开源命令行式工具文献报道极少，报道的软件Fragman经测试，结果和GeneMapper或GeneMarker等行业标准软件差异较大，并且只能适配小的分子量内标（检测上限为500bp），然而动态突变由于拷贝数可能超过500bp，需要使用大分子量内标（检测上限为1200bp）。因此，需要一款能够自动化地快速计算罕见病动态突变拷贝数的方法和系统，并且能同时适用于不同大小的分子量内标，方便搭建自动化流程，提高检测效率和结果的准确性，保证流程可控。

发明内容

针对现有技术的不足，本发明提供一种动态突变拷贝数的方法和系统，能够从毛细管电泳下机结果中快速检测并筛选可能的动态突变等位基因峰，直接输出拷贝数，并且适用大小不同的分子量内标，能并行分析多个下机结果，实现流程的高效和可控。

为实现上述目的，本发明提供一种动态突变拷贝数的方法和系统，包括如下步骤。

S1：提取下机数据中的荧光信号矩阵，并进行基线对齐。其特征在于，从原始测序文件中按照荧光的顺序提取所有的荧光信号矩阵。进一步地，将各个信道的荧光信号减去前20%的信号均值进行信号的归一化。

S2：检测分子量内标信号峰，根据内标大小使用不同阈值过滤合适的峰。其特征在于，先对信号进行平滑，并通过一阶导数找到峰值顶点，再拟合高斯峰，计算所有峰拟合的R方、峰高、峰的起始终止位置以及半峰全宽。进一步地，根据峰的特征进行筛选，针对半峰全宽中值大于10的分子量内标，输出峰高在50至10000 RFU之间、半峰全宽在3至21帧之间，且R方不小于0.8的结果；针对半峰全宽中值小于10的分子量内标，输出峰高在50至10000 RFU之间、半峰全宽不超过10帧，且R方不小于0.7的结果。过滤连续峰和上拉峰后，得到峰信号矩阵，由峰高、峰值位置、峰起始、峰终止四列组成，峰信号矩阵按照峰值位置从小到大排序，每个峰就代表一个片段。

S3：通过距离分类序列匹配分子量内标，并逐步校正匹配后的片段大小。其特征在于，先将分子量内标峰信号矩阵标准化为距离分类序列，并与片段大小的距离分类矩阵进行序列比对。对于比对成功的距离分类值，将该值对应的片段匹配相应的分子量内标，得到分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小五列组成，遍历该矩阵进行片段大小的校正，包括如下步骤。

S301：使用单位距离归一化信号峰，其特征在于，首先计算片段的峰信号差分均值，除以分子量内标片段的差分均值，得到单位距离。将信号差分除以单位距离，使其归一化到与分子量内标片段大小相同的数量级。

S302：将峰信号转换为距离分类序列。其特征在于，计算分子量内标差分的最小值、中值和最大值，分别用A/C/T三个字母表示，将标准化信号和分子量内标片段大小的差分数列与这三个值进行比较并取最接近的值所代表的字母，得到标准化信号和分子量内标的距离分类序列。

S303：比对距离分类序列，将比对成功的片段对应分子量内标片段大小。其特征在于，将两个距离分类序列进行比对，并计算比对成功的序列对应的片段和分子量内标。如果2个不同的片段对应到相同的分子量内标，选择峰信号较高的片段，确保分子量内标对应到唯一的片段，得到分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小（size）五列组成。

S304：逐步校正分子量内标峰信号矩阵。其特征在于，遍历每一个峰信号，一个信号代表一个片段，首先定义该片段左边2个和右边1个已知size的片段构成直线1；定义该片段左边1个和右边2个已知size的片段构成直线2；如果该片段缺少直线1或者直线2，则只取一条直线；如果两条直线都缺少，则取最近的连续三个已知size的片段作为直线。根据电泳中片段迁移相同距离所需的时间取决于片段大小，来计算未知片段的size，如果利用2条直线则求两者的均值。进一步地，求离该值最接近的分子量内标片段，计算两者的偏差，并将偏差最小的分子量内标片段size值作为该片段的size。对于size偏差小于4的片段，进一步比较它与左右两边已知片段的size是否一致。如果该片段size介于左右两边最近的已知size的片段之间，则认为该片段属于分子量内标；如果片段size等于左边最近的已知片段，则选偏差较小的片段作为该size的分子量内标；如果size小于左边最近的已知片段，则该片段不属于分子量内标；如果片段size等于右边最近的片段，则也认为该片段属于分子量内标，在遍历下一个片段的时候再进行校正；如果片段size大于右边最近的片段，则认为该片段不属于分子量内标。经过校正，去掉不属于分子量内标的峰，得到校正后的分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、预期片段大小五列组成。

S4：检测待测样本信号峰，根据内标大小使用不同阈值过滤合适的峰。其特征在于，提取待测样本信道的信号，先对信号进行平滑，并通过一阶导数找到峰值顶点，再拟合高斯峰，计算所有峰拟合的R方、峰高、峰的起始终止位置以及半峰全宽。进一步地，根据峰的特征进行筛选，针对半峰全宽中值大于10的分子量内标，输出半峰全宽在3至21帧之间，且R方不小于0.8的结果；针对半峰全宽中值小于10的分子量内标，输出半峰全宽不超过10帧，且R方不小于0.7的结果。得到峰信号矩阵，由峰高、峰值位置、峰起始、峰终止四列组成，峰信号矩阵按照峰值位置从小到大排序。

S5：根据分子量内标计算目标峰的片段大小，并根据基线计算片段的拷贝数。其特征在于，对于分子量内标上下限100帧内的片段，采用附近的分子量内标拟合的直线进行片段大小的计算，在此之外用所有的分子量内标拟合的线性模型进行片段大小的计算。进一步地，根据基线计算片段峰的拷贝数。基线数据由已知拷贝数的样本在相同实验条件下得到的size构成。

S6：根据峰的特征筛选可能的等位基因，提示样本的基因型。其特征在于，根据峰高、峰的位置、片段大小等特征筛选可能的等位基因峰，包括如下步骤。

S601：先筛选片段不小于300，且峰高在

至

之间的峰，按照峰高降序排列，以第一个峰的30%作为峰高阈值1，以第五个峰高度的5%作为峰高阈值2。进一步地，侯选峰中最高的峰为等位基因峰1，以该峰为标准，遍历该峰左边（即片段小于该峰）和右边（片段大于该峰）的候选峰。先遍历等位基因峰1左边的峰，若第二高的峰高于阈值1，距离等位基因峰1大于2.5个碱基，则标记为等位基因2；若第二高峰不符合条件，接继续往下遍历，以此类推，直到遍历完所有峰或是找到等位基因2为止。

S602：遍历等位基因峰1右边的峰，如果小于500个碱基采用阈值1，大于500个碱基的片段则采用阈值2；若存在峰高大于阈值，且距离等位基因峰1大于2.5个碱基的峰，将其标记为等位基因峰2。

S603：如果左右两边都存在等位基因峰2，则优先选择右边的峰。得到待测样本峰的检测结果矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小、拷贝数和等位基因标记7列组成。

本发明的有益效果是，与现有技术相比，本方法能快速计算罕见病动态突变拷贝数，并判定可能的等位基因峰，不仅可以适配不同大小的分子量内标，还适用于多个下机结果并行分析。此方法减少人工判读对结果的干扰，可快速分析同一个样本的多个文件，同时判断等位基因峰的标准统一而且可控，避免人工判读引入的误差，极大提升动态突变数据分析的效率，并保证分析流程可以追溯，提高了准确率。

附图说明

图1是本发明实施例提供的一种动态突变拷贝数的方法和系统的流程图。

图2是本发明实施例提供的距离分类序列匹配分子量内标的流程图。

图3是本发明实施例1中信号峰拟合结果示意图。

图4是本发明实施例1中距离分类序列匹配分子量内标的算法示意图。

图5是本发明实施例1中分子量内标检测结果拟合直线示意图。

图6是本发明实施例1中待测样本片段大小检测算法示意图。

图7是本发明实施例1中待测样本动态突变检测结果表格示意图。

图8是本发明实施例1中待测样本检出等位基因和分子量内标的示意图。

图9是本发明实施例2中信号峰拟合结果示意图。

图10是本发明实施例2中Fragman软件识别分子量内标示意图。

图11是本发明实施例2中本方法识别分子量内标示意图。

图12是本发明实施例2中待测样本检出等位基因和分子量内标的示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的说明。以下实施例和附图用于说明本发明，但不用来限制本发明的范围。

实施例一

本实施例以脊髓小脑共济失调患者Sample01为例，该样本进行了SCA1/2/3/6/7五种亚型的动态突变的毛细管电泳片段分析，待测样本染蓝色荧光标记；分子量内标使用GeneScan™ 500 LIZ™ （简称GS500），包含 16 个单链标记片段：35、50、75、100、139、150、160、200、250、300、340、350、400、450、490 和 500 个核苷酸，带橙色荧光标记。通过识别文件名可获得全部5个下机数据fsa文件，每个文件都利用本发明所述方法进行处理，即提取下机数据中的荧光信号矩阵，先检测分子量内标的信号峰并匹配片段大小，再检测待测样本信号峰，通过分子量内标计算目标峰对应的片段大小，并根据基线计算片段拷贝数，最后筛选可能的等位基因峰，输出检测结果表格，多种亚型结果可以合并输出。本实施例提供的实现过程主要以检出致病拷贝数的Sample01.SCA7.fsa文件为例，其他文件处理过程同该文件，下文不赘述。

如图1，本发明提供一种动态突变拷贝数的方法和系统包括以下步骤。

S1：提取下机数据中的荧光信号矩阵，并进行基线对齐。其特征在于，从原始测序文件中按照荧光的顺序提取所有的荧光信号矩阵。进一步地，将各个信道的荧光信号减去前20%的信号均值进行信号的归一化。

本实施例中，读取Sample01.SCA7.fsa二进制文件的DATA字段，该字段存储蓝绿黄红橙五种荧光染料的信号，标签分别是DATA.1、DATA.2、DATA.3、DATA.4、DATA.105。按照顺序输出为矩阵，矩阵的每一列分别是一种信道的信号。其中DATA.1是待测样本的信道，DATA.105是分子量内标的信道。

本实施例中，对所有信号进行基线对齐，对于每一列信号，提取前20%求均值，这部分信号不包含片段峰，代表了不同信道间的初始差异，通过减去该值，保证信道能对齐到相同水平。计算公式如下：

其中

代表某一信道上第i个原始信号，j代表该信道前 20%上的第j个信号，N代表该信道所有信号的总数，

代表基线对齐后的信号。

S2：检测分子量内标信号峰，根据内标大小使用不同阈值过滤合适的峰。其特征在于，先对信号进行平滑，并通过一阶导数找到峰值顶点，再拟合高斯峰，计算所有峰拟合的R方、峰高、峰的起始终止位置以及半峰全宽。进一步地，根据峰的特征进行筛选，针对半峰全宽中值大于10的分子量内标，输出峰高在50至10000 RFU之间、半峰全宽在3至21帧之间，且R方不小于0.8的结果；针对半峰全宽中值小于10的分子量内标，输出峰高在50至10000 RFU之间、半峰全宽不超过10帧，且R方不小于0.7的结果。过滤连续峰和上拉峰后，得到峰信号矩阵，由峰高、峰值位置、峰起始、峰终止四列组成，峰信号矩阵按照峰值位置从小到大排序，每个峰就代表一个片段。

本实施例中，选择信号矩阵的第五列，即分子量内标的信号，先用Savitzky-Golay滤波器对信号进行平滑，步长为21。对平滑后数据计算一阶差分的符号差分。遍历符号差分等于-2的信号，选择15、11、21、31四种窗口大小，以该信号为中心提取4种窗口内的所有数据分别进行高斯峰的拟合，计算模型的R方、峰高、峰起始终止位置以及半峰全宽，输出峰高最接近原始数据且R方最高的窗口结果。用一阶差分可以去除驼峰信号，并且拟合得到的峰不受分裂峰型的影响。如图3，第6537帧的符号差分等于-2，为可能的峰信号，以6573为中心，取4种窗口内的数据点拟合高斯峰，空心圆为原始数据点，曲线为拟合结果，拟合结果见横坐标，窗口为11的结果峰高为4136.1 RFU最接近原始峰高4063 RFU，且R方最高，输出该窗口的拟合结果。

本实施例中，使用的分子量内标检测上限是500个碱基，属于小型的分子量内标，即半峰全宽中值小于10帧，过滤标准为峰高在50至10000RFU之间、半峰全宽不超过10帧，且R方不小于0.7。进一步地，过滤连续峰，计算峰的距离，如果发现信号距离小于3.5帧的连续峰，则只保留其中最高的峰。进一步地，过滤上拉峰（Pull-up Peak），由于待测样本超过8000RFU时，超出了机器检测阈值，会导致荧光信号卷积出现错误，造成没有内标的区域出现上拉峰，上拉峰的特征是和待测样本峰具有相同的峰位置，且高度一般不超过待测样本峰的20%。如图4-a，筛选得到峰信号矩阵，由峰高height、峰值位置position、峰起始start、峰终止end四列组成，峰信号矩阵按照position从小到大排序。

S3：通过距离分类序列匹配分子量内标，并逐步校正匹配后的片段大小。其特征在于，先将分子量内标峰信号矩阵标准化为距离分类序列，并与片段大小的距离分类矩阵进行序列比对。对于比对成功的距离分类值，将该值对应的片段匹配相应的分子量内标，得到分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小五列组成，遍历该矩阵进行片段大小的校正，包括如下步骤。

S301：使用单位距离归一化信号峰，其特征在于，首先计算片段的峰信号差分均值，除以分子量内标片段的差分均值，得到单位距离。将信号差分除以单位距离，使其归一化到与分子量内标片段大小相同的数量级。

本实施例中，如图4，计算4-a 的position列的差分均值为293，并除以分子量内标片段的差分均值，分子量内标片段size为35、50、75、100、139、150、160、200、250、300、340、350、400、450、490 和 500，其差分均值为31，得到单位距离为9.45。将position的差分除以单位距离，与分子量内标的距离归一化到相同的数量级，如图4-b-1。

S302：将峰信号转换为距离分类序列。其特征在于，计算分子量内标差分的最小值、中值和最大值，分别用A/C/T三个字母表示，将标准化信号和分子量内标片段大小的差分数列与这三个值进行比较并取最接近的值所代表的字母，得到标准化信号和分子量内标的距离分类序列。

本实施例中，如图4-b-1和4-b-2，先计算分子量内标差分的最小值、中值和最大值分表为10、39和50，分别用A/C/T三个字母表示，将片段位置和分子量内标片段大小的差分数列与这三个值进行比较，并取最接近的值所代表的字母，得到标准化信号和分子量内标的距离分类序列。

S303：比对距离分类序列，将比对成功的片段对应分子量内标片段大小。其特征在于，将两个距离分类序列进行比对，并计算比对成功的序列对应的片段和分子量内标。如果2个不同的片段对应到相同的分子量内标，选择峰信号较高的片段，确保分子量内标对应到唯一的片段，得到分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小（size）五列组成。

本实施例中，如图4-c，标准化后的信号和分子量内标的距离分类序列比对后，计算比对成功的序列对应片段和分子量内标。如第12个序列比对成功，其对应的是信号峰4-a第12和13行，对应分子量内标第8和9个片段，即位置3013对应的片段大小为200，3570对应250。没有对应分子量内标的片段size先使用NA进行替代，结果如图4-a。

S304：逐步校正分子量内标峰信号矩阵。其特征在于，遍历每一个峰信号，一个信号代表一个片段，首先定义该片段左边2个和右边1个已知size的片段构成直线1；定义该片段左边1个和右边2个已知size的片段构成直线2；如果该片段缺少直线1或者直线2，则只取一条直线；如果两条直线都缺少，则取最近的连续三个已知size的片段作为直线。根据电泳中片段迁移相同距离所需的时间取决于片段大小，来计算未知片段的size，如果利用2条直线则求两者的均值。进一步地，求离该值最接近的分子量内标片段，计算两者的偏差，并将偏差最小的分子量内标片段size值作为该片段的size。对于size偏差小于4的片段，进一步比较它与左右两边已知片段的size是否一致。如果该片段size介于左右两边最近的已知size的片段之间，则认为该片段属于分子量内标；如果片段size等于左边最近的已知片段，则选偏差较小的片段作为该size的分子量内标；如果size小于左边最近的已知片段，则该片段不属于分子量内标；如果片段size等于右边最近的片段，则也认为该片段属于分子量内标，在遍历下一个片段的时候再进行校正；如果片段size大于右边最近的片段，则认为该片段不属于分子量内标。经过校正，去掉不属于分子量内标的峰，得到校正后的分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、预期片段大小五列组成。

本实施例中，电泳中片段迁移相同距离所需的时间取决于片段大小，符合如下计算公式：

其中y和x为未知片段的size和迁移时间，c、

和

为常数。将三个已知片段的position和size带入公式（position是峰值位置代表的帧数，即电泳迁移时间），能求出常数c、

和

值，利用公式能求出待测片段的大小。如图4-a，如第1行977处信号，它缺少直线1和2，只能取离它最近的且有片段大小的点即第5、6、7行构成的直线来计算size为4，由于它离最近的分子量内标30相差超过4个碱基，因此认为它不属于分子量内标，将size标注为NA；如第16行4770帧信号，直线1为14、15、17行，求得size为350.8，直线2为15、17、18行，求得size为349.5，计算两条直线的均值为350.1作为该位置对应的片段size。由于它离最近分子量内标350的距离偏差小于4，且其size介于15行和17行片段size之间，因此认为他属于分子量内标，且size校正为350。最终经过校正，去掉size为NA的信号，得到分子量内标峰信号矩阵，由峰高、峰值位置、峰起始、峰终止、预期片段大小五列组成。如图5，实心圆代表检出的分子量内标对应的位置和片段size，校正后的峰的迁移距离和片段size符合线性关系。

S4：检测待测样本信号峰，根据内标大小使用不同阈值过滤合适的峰。其特征在于，提取待测样本信道的信号，先对信号进行平滑，并通过一阶导数找到峰值顶点，再拟合高斯峰，计算所有峰拟合的R方、峰高、峰的起始终止位置以及半峰全宽。进一步地，根据峰的特征进行筛选，针对半峰全宽中值大于10的分子量内标，输出半峰全宽在3至21帧之间，且R方不小于0.8的结果；针对半峰全宽中值小于10的分子量内标，输出半峰全宽不超过10帧，且R方不小于0.7的结果。得到峰信号矩阵，由峰高、峰值位置、峰起始、峰终止四列组成，峰信号矩阵按照峰值位置从小到大排序。

本实施例中，选择S1得到的信号矩阵第1列，即待测样本信号，同S2，先用Savitzky-Golay滤波器对信号进行平滑，步长为21。对平滑后数据计算一阶差分的符号差分。遍历符号差分等于-2的信号，选择15、11、21、31四种窗口大小，以该信号为中心提取4种窗口内的信号拟合高斯峰，计算模型的R方、峰高、峰起始终止位置以及半峰全宽，输出峰高最接近原始数据且R方最高的窗口结果。进一步地，根据峰的特征进行筛选，针对半峰全宽中值小于10的分子量内标，输出半峰全宽不超过10帧，且R方不小于0.7的结果，不根据信号大小和信号间距进行过滤。得到峰信号矩阵，由峰高、峰值位置、峰起始、峰终止四列组成，峰信号矩阵按照峰值位置从小到大排序，如图7的前4列展示了片段大于300的峰的检出结果。

S5：根据分子量内标计算目标峰的片段大小，并根据基线计算片段的拷贝数。其特征在于，对于分子量内标上下限100帧内的片段，采用附近的分子量内标拟合的直线进行片段大小的计算，在此之外用所有的分子量内标拟合的线性模型进行片段大小的计算。进一步地，根据基线计算片段峰的拷贝数。基线数据由已知拷贝数的样本在相同实验条件下得到的size构成。

本实例中，根据分子量内标计算目标峰的片段大小，如图6，实心圆代表分子量内标对应的位置和片段大小，空心圆代表检出的片段，此处示例了4936帧处的待测样本峰，其片段大小由2条直线的均值计算得来，第一条直线由左边的340和350内标以及400的内标构成，第二条直线由左边的350和右边的400、450内标构成，计算方法同S304，得到结果size是363.4。7036帧处的待测样本峰，由于超出了分子量内标的检测上限，并且离500的内标对应的6537帧超过100帧，则采用所有的分子量内标拟合的直线来计算size为545.0。

本实例中，根据基线计算片段峰的拷贝数，基线数据由2个已知拷贝数的纯合样本在相同实验条件下测得的片段大小构成，分别是361.96和367.48，拷贝数分别是10和12。通过以下公式换算得到待测峰的拷贝数。

其中S为片段大小，C为拷贝数，N为重复序列的碱基个数，下标b为基线，t为待测样本。

本实例中，重复序列为3个碱基，因此N=3，可求得待测峰的拷贝数，第4936帧的片段将size=363.4代入公式，得到拷贝数为11，而7036帧片段代入size=545得到拷贝数为71。如图7的前6列，得到待测样本峰的检测结果矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小和拷贝数6列组成。

S6：根据峰的特征筛选可能的等位基因，提示样本的基因型。其特征在于，根据峰高、峰的位置、片段大小等特征筛选可能的等位基因峰，包括如下步骤。

S601：先筛选片段不小于300，且峰高在

至

之间的峰，按照峰高降序排列，以第一个峰的30%作为峰高阈值1，以第五个峰高度的5%作为峰高阈值2。进一步地，侯选峰中最高的峰为等位基因峰1，以该峰为标准，遍历该峰左边（即片段小于该峰）和右边（片段大于该峰）的候选峰。先遍历等位基因峰1左边的峰，若第二高的峰高于阈值1，距离等位基因峰1大于2.5个碱基，则标记为等位基因2；若第二高峰不符合条件，接继续往下遍历，以此类推，直到遍历完所有峰或是找到等位基因2为止。

S602：遍历等位基因峰1右边的峰，如果小于500个碱基采用阈值1，大于500个碱基的片段则采用阈值2；若存在峰高大于阈值，且距离等位基因峰1大于2.5个碱基的峰，将其标记为等位基因峰2。

S603：如果左右两边都存在等位基因峰2，则优先选择右边的峰。得到待测样本峰的检测结果矩阵，由峰高、峰值位置、峰起始、峰终止、片段大小、拷贝数和等位基因标记7列组成。

本实例中，先筛选片段不小于300，且峰高在

至

之间的峰，按照峰高降序排列，如图7，为候选等位基因峰。最高峰为第11行4936帧，峰高33404.91，片段大小为363.4，将该峰标记为等位基因峰1。先在小于4936帧的范围内寻找等位基因峰2，第二高峰,为第10行4936帧，峰高未32855.83，片段大小为363.3，但因为该帧距离等位基因峰1小于2.5个碱基，因此可能是stutter峰，而不是真正的等位基因峰；再寻找第二高峰，为第9行的4905帧，峰高为11668.72，片段大小为360.9，高于最高峰的30%（即峰高阈值1），且距离等位基因峰1大于2.5个碱基，暂时标记该峰为等位基因峰2。

进一步地，在大于4936帧范围寻找第二高峰，为第14行7036帧，峰高6530.46，片段大小为545，超过500，因此峰高阈值不采用峰高阈值1，而是第五个高峰7733.94的5%即386.70作为阈值，该峰高于峰高阈值2，且距离等位基因峰1超过2.5个碱基，暂时标记该峰为等位基因峰2。

进一步地，在两边都存在等位基因峰2的情况下，优先提示大于等位基因峰1的结果。如图7，allele列标记selected，表示该峰被作为可能的等位基因峰，该样本可能的基因型为杂合。图8为待测样本和分子量内标的片段分析结果示意图，灰色线代表分子量内标，黑色线代表待测样本，并标注了峰对应的片段大小。上方的图展示了检出片段和分子量内标的整体关系，下方的2个图则放大视野展示了等位基因峰的检出情况。

实施例二

本实施例以脊髓小脑共济失调患者Sample02为例，该样本进行了SCA2亚型动态突变的毛细管电泳片段分析，待测样本染蓝色荧光标记；分子量内标使用GeneScan™ 1200LIZ™ （简称GS1200），包含68个单链标记片段：20、30、40、60、80、100、114、120、140、160、180、200、214、220、240、250、260、280、300、314、320、340、360、380、400、414、420、440、460、480、500、514、520、540、560、580、600、614、620、640、660、680、700、714、720、740、760、780、800、820、840、850、860、880、900、920、940、960、980、1000、1020、1040、1060、1080、1100、1120、1160和1200个核苷酸，带橙色荧光标记。本实施例主要以Sample02.SCA2.fsa文件为例展示本方法对大分子量内标的检测效果。

S1：提取下机数据中的荧光信号矩阵，并进行基线对齐。

S2：检测分子量内标信号峰，根据内标大小使用不同阈值过滤合适的峰。

本实施例中，检测分子量内标信号峰方法与实施例1相同，如图9，展示了第26265帧处4个窗口的拟合结果，窗口大小为31的结果峰高为198.8 RFU，最接近原始峰高198RFU且R方最高，输出该结果。大分子量内标的峰信号具有更多的数据点，较大的窗口拟合结果相对更好。

本实施例中，由于样本采用了大分子量内标，拟合后的半峰全宽大于10，因此采用适合大分子量内标的过滤阈值，即峰高在50至10000 RFU之间、半峰全宽在3至21帧之间，且R方不小于0.8。经过滤连续峰和上拉峰后得到峰信号矩阵。

S3：通过距离分类序列匹配分子量内标，并逐步校正匹配后的片段大小。

本实例中，通过S2的检测峰的方法，可以有效地消除驼峰，保证S3匹配分子量内标的准确。如图10和11分别展示了Fragman软件和本方法的分子量内标检测结果，横坐标表示片段size，纵坐标为信号高度，黑线代表每一个信号，实心圆代表识别的分子量内标峰。可以看出在600-640片段大小的附近有一个驼峰，Fragman无法准确识别附近的分子量内标，而本方法能够消除驼峰的影响。同时Fragman对大分子量内标无法正确识别，导致影响结果准确性，如20至60之间漏检了一个峰，导致将80的片段识别成60；如600和614间分子量内标，漏检了600的峰，错误地识别614为600。

S4：检测待测样本信号峰，根据内标大小使用不同阈值过滤合适的峰。

S5：根据分子量内标计算目标峰的片段大小，并根据基线计算片段的拷贝数。

S6：根据峰的特征筛选可能的等位基因峰，提示样本的基因型。

S3至S6与实例1相同，片段分析结果如图12，同样，灰色线代表分子量内标，黑色线代表待测样本，并标注了峰对应的片段大小。上方的图展示了检出片段和分子量内标的整体关系，下方的图则放大视野展示了等位基因峰的检出情况。

上述是本发明的优选实施方式，通过本发明所述方法能够快速计算罕见病动态突变的拷贝数，本实施例1中示例文件计算在5秒内完成，本实施例样本全部文件在1分钟内可完成所有文件分析，大大提高了检测效率和准确性。

本发明所公开的实施例和上述说明，用于使本领域专业技术人员能够实现或使用本发明，熟悉本领域的专业技术人员在不违背本发明精神的前提下对本发明做的等同变形或替换的内容，均包含在本申请权利要求所限定的范围内。

Claims (2)

1.一种动态突变拷贝数的方法和系统，其特征在于，包括：

S1：提取下机数据中的荧光信号矩阵，并进行基线对齐；

S2：检测分子量内标信号峰，根据内标大小使用不同阈值过滤合适的峰；

S3：通过距离分类序列匹配分子量内标，并逐步校正匹配后的片段大小；

S4：检测待测样本信号峰，根据内标大小使用不同阈值过滤合适的峰；

S5：根据分子量内标计算目标峰的片段大小，并根据基线计算片段的拷贝数；

S6：根据峰的特征筛选可能的等位基因，提示样本的基因型。

2.根据权利要求1所述的方法，其特征在于，所述S3包括：

S301：使用单位距离归一化信号峰，其特征在于，首先计算片段的峰信号差分均值，除以分子量内标片段的差分均值，得到单位距离；将信号差分除以单位距离，使其归一化到与分子量内标片段大小相同的数量级；

S302：将峰信号转换为距离分类序列，其特征在于，计算分子量内标差分的最小值、中值和最大值，分别用A/C/T三个字母表示，将标准化信号和分子量内标片段大小的差分数列与这三个值进行比较并取最接近的值所代表的字母，得到标准化信号和分子量内标的距离分类序列；

S303：比对距离分类序列，将比对成功的片段对应分子量内标片段大小，其特征在于，将两个距离分类序列进行比对，并计算比对成功的序列对应的片段和分子量内标；如果2个不同的片段对应到相同的分子量内标，选择峰信号较高的片段，确保分子量内标对应到唯一的片段；

S304：逐步校正分子量内标峰信号矩阵，其特征在于，遍历每一个峰信号，一个信号代表一个片段，首先定义该片段左边2个和右边1个已知size的片段构成直线1；定义该片段左边1个和右边2个已知size的片段构成直线2；如果该片段缺少直线1或者直线2，则只取一条直线；如果两条直线都缺少，则取最近的连续三个已知size的片段作为直线；根据电泳中片段迁移相同距离所需的时间取决于片段大小，来计算未知片段的size，如果利用2条直线则求两者的均值；进一步地，求离该值最接近的分子量内标片段，计算两者的偏差，并将偏差最小的分子量内标片段size值作为该片段的size；对于size偏差小于4的片段，进一步比较它与左右两边已知片段的size是否一致；如果该片段size介于左右两边最近的已知size的片段之间，则认为该片段属于分子量内标；如果片段size等于左边最近的已知片段，则选偏差较小的片段作为该size的分子量内标；如果size小于左边最近的已知片段，则该片段不属于分子量内标；如果片段size等于右边最近的片段，则也认为该片段属于分子量内标，在遍历下一个片段的时候再进行校正；如果片段size大于右边最近的片段，则认为该片段不属于分子量内标；经过校正，去掉不属于分子量内标的峰，得到校正后的分子量内标峰信号矩阵。

Images