本发明涉及基因检测技术领域，尤其涉及一种基因高通量测序数据突变检测方法。
背景技术：
在过去的肿瘤基因突变检测的临床与科研应用中我们通常只关注肿瘤组织中高丰度基因突变的情况。低丰度突变由于突变核酸含量低，在测序覆盖度较低的情况下极易出现漏检或与假阳性的情况。但在某些应用场景下，例如通过液态活检的方式检测血液中的低丰度肿瘤突变核酸，需要准确的检测出低丰度的突变。通过高通量测序靶向捕获或扩增技术结合高深度的测序，可以提高重要肿瘤突变位点的测序覆盖度，提高检测的灵敏度。但是由于高通量测序中天然存在的噪音，仅从实验角度仍然很难将真实突变和噪音点区分开来，必须通过算法建立降噪和突变检测的模型解决上述问题。
现有方案中采用健康人的测序数据作为背景值，通过正态分布拟合确定每个位点的背景噪音阈值，从而区分真阳性点和噪音。但是此种方案存在以下几个问题：1.高通量测序实验及数据产生存在批次效应，通过健康人的测序数据建立的背景模型能够去除测序系统本身存在的系统误差，但是对于每次实验中随机产生的实验误差无法有效去除；2.建立健康人群的背景数据需要测定大量位点的大样本量的数据，所需成本较高，对于背景数据库中暂未覆盖的位点无法起到降低噪音的作用。
技术实现要素：
针对上述现有技术中的不足，本发明提供一种基因高通量测序数据突变检测方法，采用虚拟分子标签与背景数据库结合的方法降低噪音，提高检测的特异性和敏感性，可在不增加实验成本的前提下能有效的降低实验中的随机误差，结合背景数据库对系统误差的校正，可以达到准确鉴定低丰度突变的目的。
为了实现上述目的，本发明提供一种基因高通量测序数据突变检测方法，包括步骤：
s1：获取一基因样本的高通量测序数据；
s2：生成所述基因样本的高通量测序数据的各基因序列的位置信息标签；
s3：根据所述位置信息标签将各所述基因序列分组并计算获得一突变总量；
s4：将所述突变总量代入一背景模型输出突变检测结果。
优选地，所述s2步骤进一步包括步骤：
s21：通过一序列对比算法将所述各基因序列对比到一参考基因组，形成各所述基因序列的比对信息；
s22：将所述比对信息存储于一sam/bam格式文件中；
s23：根据所述sam/bam格式文件判断各所述基因序列的序列来源的模板链ti，1≤i≤n，n为所述基因序列个数；
s24:根据所述序列来源的模板链ti和所述sam/bam格式文件生成各所述基因序列的位置信息标签。
优选地，所述s23步骤进一步包括步骤：
自所述sam/bam格式文件中提取每一条所述基因序列的一第一比对起始位置pi、一同片段对比序列的一第二对比起始位置qi、正负链信息si和所述基因序列的序列号ri；
当所述基因序列的序列号ri等于所述sam/bam格式文件的read1位置的数值且所述正负链信息si等于所述sam/bam格式文件的foward位置的数值时，或所述基因序列的序列号ri不等于所述sam/bam格式文件的read1位置的数值且所述正负链信息si不等于所述sam/bam格式文件的foward位置的数值时，所述序列来源的模板链ti为正；
当所述基因序列的序列号ri等于所述sam/bam格式文件的read1位置的数值且所述正负链信息si不等于所述sam/bam格式文件的foward位置的数值时，或所述基因序列的序列号ri不等于所述sam/bam格式文件的read1位置的数值且所述正负链信息si等于所述sam/bam格式文件的foward位置的数值时，所述序列来源的模板链ti为负。
优选地，所述位置信息标签表示为(pi,qi,ti)。
优选地，所述s3步骤进一步包括步骤：
s31：将所述位置信息标签一致的所述基因序列分至同一基因组；
s32：统计所述基因组中各所述基因序列与所述参考基因组一目标基因位置gi对应的一当前基因位置为突变型基因型且碱基质量q>30的所述基因序列的突变数vj，j为大于等于1的自然数；
如vj>0，记录所述当前基因位置碱基质量q>30的所述基因序列个数nj；
如vj<f*nj，则vj＝0，其中f为预设的最低碱基一致性比例值；
s33：重复步骤s32获得各所述目标基因位置的突变数vj，并根据所述突变数计算一突变总数
当
当
优选地，所述s4步骤进一步包括步骤：
s41：建立一背景模型，所述背景模型的公式为：
其中，pgi为累积分布频率，γ为第一拟合参数，δ为第二拟合参数，ε为第三拟合参数，λ为第四拟合参数；
根据多个样本数据的拟合获得所述第一拟合参数、所述第二拟合参数、所述第三拟合参数和所述第四拟合参数；
s42：将所述突变总量代入所述背景模型，计算所述累积分布频率；
s43：当所述累积分布频率数值大于0.95时，判定与当前所述位置信息标签对应的一基因位点为阳性位点。
优选地，所述样本数据的个数大于等于1000。
本发明由于采用了以上技术方案，使其具有以下有益效果：
1、在不增加实验步骤和成本的前提下去除高通量测序中的随机噪音。
2、通过对去除随机噪音后的健康人测序数据进行建模，建立了用来判别阳性突变位点的计算模型。
最终，可以在不改变现有实验体系的前提下明显提高低丰度变异检出的敏感性和特异性。
附图说明
图1为本发明实施例的基因高通量测序数据突变检测方法的流程图。
具体实施方式
下面根据附图1，给出本发明的较佳实施例，并予以详细描述，使能更好地理解本发明的功能、特点。
请参阅图1，本发明实施例的一种基因高通量测序数据突变检测方法，包括步骤：
s1：获取一基因样本的高通量测序数据。
s2：生成基因样本的高通量测序数据的各基因序列的位置信息标签。
其中，s2步骤进一步包括步骤：
s21：通过一序列对比算法将各基因序列对比到一参考基因组，形成各基因序列的比对信息；序列对比算法可采用现有的任意序列对比算法，对其不做具体限制；比对信息包括第一比对起始位置信息、第二对比起始位置信息、碱基质量信息、正负链信息和基因序列的序列号信息等；
s22：将比对信息存储于一sam/bam格式文件中；
s23：根据sam/bam格式文件判断各基因序列的序列来源的模板链ti，1≤i≤n，n为基因序列个数；
s24:根据序列来源的模板链ti和sam/bam格式文件生成各基因序列的位置信息标签。
其中，s23步骤进一步包括步骤：
自sam/bam格式文件中提取每一条基因序列的一第一比对起始位置pi、一同片段对比序列的一第二对比起始位置qi、正负链信息si和基因序列的序列号ri；序列来源的模板链ti的逻辑关系式可表示为：
当基因序列的序列号ri等于sam/bam格式文件的read1位置的数值且正负链信息si等于sam/bam格式文件的foward位置的数值时，或基因序列的序列号ri不等于sam/bam格式文件的read1位置的数值且正负链信息si不等于sam/bam格式文件的foward位置的数值时，序列来源的模板链ti为正；
当基因序列的序列号ri等于sam/bam格式文件的read1位置的数值且正负链信息si不等于sam/bam格式文件的foward位置的数值时，或基因序列的序列号ri不等于sam/bam格式文件的read1位置的数值且正负链信息si等于sam/bam格式文件的foward位置的数值时，序列来源的模板链ti为负。
本实施例中，位置信息标签表示为(pi,qi,ti)。该三元组能够唯一标识来自于统一模板核酸的所有序列，并且能够区分模板的正义链和反义链。
s3：根据位置信息标签将各基因序列分组并计算获得一突变总量。
其中，s3步骤进一步包括步骤：
s31：将位置信息标签一致的基因序列分至同一基因组；
s32：统计基因组中各基因序列与参考基因组一目标基因位置gi对应的一当前基因位置为突变型基因型且碱基质量q>30的基因序列的突变数vj，j为大于等于1的自然数；
如vj>0，记录当前基因位置碱基质量q>30的基因序列个数nj；
如vj<f*nj，则vj＝0，其中f为预设的最低碱基一致性比例值；
s33：重复步骤s32获得各目标基因位置的突变数vj，并根据突变数计算一突变总数
当
当
s4：将突变总量代入一背景模型输出突变检测结果。
其中，s4步骤进一步包括步骤：
s41：建立一背景模型，背景模型的公式为：
其中，pgi为累积分布频率，γ为第一拟合参数，δ为第二拟合参数，ε为第三拟合参数，λ为第四拟合参数；
根据1000多个样本数据的拟合获得第一拟合参数、第二拟合参数、第三拟合参数和第四拟合参数；
s42：将突变总量代入背景模型，计算累积分布频率；
s43：当累积分布频率数值大于0.95时，判定与当前位置信息标签对应的一基因位点为阳性位点。
本发明实施例的一种基因高通量测序数据突变检测方法，具有以下有益效果：
1、在不增加实验步骤和成本的前提下去除高通量测序中的随机噪音。
2、通过对去除随机噪音后的健康人测序数据进行建模，建立了用来判别阳性突变位点的计算模型。
最终，可以在不改变现有实验体系的前提下明显提高低丰度变异检出的敏感性和特异性。
以上结合附图实施例对本发明进行了详细说明，本领域中普通技术人员可根据上述说明对本发明做出种种变化例。因而，实施例中的某些细节不应构成对本发明的限定，本发明将以所附权利要求书界定的范围作为本发明的保护范围。