一种校正线粒体基因组测序突变的方法、装置和存储介质与流程

1.本技术涉及线粒体基因组突变检测技术领域，特别是涉及一种校正线粒体基因组测序突变的方法、装置和存储介质。


背景技术：

2.通过对人类血液或其他组织样本，采用探针杂交捕获技术进行线粒体基因组dna捕获和富集，然后进行高通量测序，即可得到线粒体基因组序列。通过对测序数据进行质控，去除重复序列，使用比对软件进行序列比对后，使用突变检测软件对线粒体基因组进行突变检测，即可得到线粒体基因组的突变数据，即线粒体基因组测序突变数据。但是，研究显示，这些突变数据有相当一部分是由于测序错误、实验方法不当或比对算法的系统误差导致；也就是说，这些突变有相当部分并非人类线粒体基因组中真实存在的突变；这部分突变会给下游的分析加重负担，并且干扰线粒体基因组测序突变数据检测结果。因此需要对这些突变进行校正和过滤，才能减少下游分析的工作量，提高分析结果的准确性。
3.目前点突变校正和过滤的方法主要使用机器学习训练模型或者使用硬指标来过滤。其中，硬指标过滤，例如突变质量值、突变的链分布偏差等。而主要应用机器学习模型的方法是vqsr(variant quality score recalibration突变质量值重校正)。
4.vqsr的核心原理是利用gmm(gaussian mixture model高斯混合模型)模型构造一个区分真变异和假变异的分类器，在构造该模型的时候使用和已知突变数据集一致的位点，并分配相应的可信度权重来进行训练。这些已知且被严格验证的变异被当作真实的变异，来训练区分真变异的gmm，接着对输入的数据进行打分，然后把评分最低的突变作为假变异的集合，来构造区分假变异的gmm。最后用构造好的两个gmm模型，即真变异gmm和假变异gmm，同时对变异进行打分，判定输入的突变属于真变异的分类还是假变异的分类。
5.为了保证分类结果可靠，vqsr在进行模型训练的时候有最低可用位点数目的要求，例如真和假变异可供训练的数目须超过5000个，否则无法训练出可用的模型。一般人类的全基因组测序能满足最低可用位点数目要求，但对于单个样本的人类外显子组测序、线粒体基因组测序或其他小型基因组测序，则无法满足以上要求。所以单个样本的线粒体基因组测序无法使用vqsr进行突变校正和过滤，即使采取多样本进行突变校正和过滤的方法，所需的样本数量也是极大的，这使得vqsr突变校正和过滤的效率极大受限。
6.另外，也有研究显示，可以使用突变质量值或突变的链分布偏差等指标进行过滤；但是，现有的基于这些硬指标进行过滤的校正方法普遍存在准确性较差、假阳性多等不足。
7.因此，如何更有效和准确的去除由于实验错误和测序错误导致的假阳性，得到更可靠的突变合集，仍然是线粒体基因组测序突变数据分析的重要难题。


技术实现要素：

8.本技术的目的是提供一种新的校正线粒体基因组测序突变的方法、装置和存储介质。
9.为了实现上述目的，本技术采用了以下技术方案：
10.本技术的第一方面公开了一种校正线粒体基因组测序突变的方法，其包括，将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。
11.其中，经过验证的已知突变集合，可以是数据库中的突变集合，例如，线粒体相关疾病数据库mseqdr、临床疾病和表型相关突变数据库clinvar，也可以是其他经过严格验证的已知突变的集合。
12.需要说明的是，本技术的校正方法，先通过已知突变白名单过滤部分突变，仅仅对没有出现在已知突变白名单中的突变进行进一步的分析和过滤，减小了后续分析的数据量和成本；并且，本技术综合各方面因素，采用11个过滤指标，去除由于实验错误和测序错误导致的假阳性，得到最终真实的突变合集。并且，本技术的校正方法，突变的过滤不需要考虑位点数量的多少，不需要同时对多样本进行过滤，单个样本即可进行突变过滤和校正；本技术的突变过滤指标更多，且构建的已知突变白名单保证了不会遗漏有意义的突变，从而使结果更可靠。可以理解，本技术的校正方法，不仅仅能够用于人类线粒体基因组突变数据的过滤校正，也可应用于白鼠、斑马鱼等其他物种的线粒体基因组或果蝇基因组等小型基因组的突变数据校正。此外，本技术的校正方法虽然能够很好的适用于单样本的过滤和校正，自然也能够用于多样本的突变过滤和校正，尤其是数量较少的情况下的多样本校正。
13.还需要说明的是，对于人类全基因组或全外显子组测序突变数据，白名单的构建作用不显著，因为突变数量过大，白名单无法包含海量的突变；但是，对于线粒体基因组这种小型的基因组，白名单突变能保证不会遗漏绝大多数已发现的可能致病突变，使得整体结果更准确可信；并且，还能够减少下游分析的工作量，提高分析结果的准确性。因此，本技术的校正方法，更适合用于线粒体基因组等小型的基因组的测序突变校正；而现有的其他校正方法更适用于较大型的基因组。
14.本技术的一种实现方式中，根据1)突变频率对突变进行过滤包括，计算所有突变的突变频率，对于同一突变位点，只保留同一突变位点中突变频率最高的突变。
15.本技术的一种实现方式中，根据2)突变所在读段的位置的平均值对突变进行过滤包括，统计同一突变位点在不同读段中的位置，计算这些位置的平均值，过滤去除突变频率小于突变频率阈值，且位置的平均值小于位置阈值的突变。
16.本技术的一种实现方式中，根据3)突变是否都位于同一位置对突变进行过滤包括，对于同一突变位点，过滤去除突变位点在不同的读段中的位置都相同的突变。
17.本技术的一种实现方式中，根据4)突变的平均碱基质量值对突变进行过滤包括，计算突变的平均碱基质量值，过滤去除平均碱基质量值低于碱基质量值阈值的突变。
18.本技术的一种实现方式中，根据5)读段比对质量值的均方根对突变进行过滤包括，计算突变所在的所有读段的对比质量值的均方根，过滤去除对比质量值均方根小于均方根阈值的突变。
19.本技术的一种实现方式中，根据6)突变的单位深度质量值对突变进行过滤包括，计算突变的单位深度质量值，过滤去除单位深度质量值小于单位深度质量值阈值的突变。
20.本技术的一种实现方式中，根据7)读段平均比对错配数量对突变进行过滤包括，计算突变所在读段的平均比对错配数量，过滤去除平均比对错配数量大于比对错配数量阈值的突变。
21.本技术的一种实现方式中，根据8)信噪比对突变进行过滤包括，统计突变所在的读段中，碱基质量值大于碱基质量值阈值的读段数量，碱基质量值小于碱基质量值阈值的读段数量，将碱基质量值大于碱基质量值阈值的读段数量除以碱基质量值小于碱基质量值阈值的读段数量的商作为信噪比，过滤去除信噪比值小于信噪比阈值的突变。
22.本技术的一种实现方式中，根据9)杂合位点读段的碱基质量值秩和检验对突变进行过滤包括，对突变碱基所在读段的碱基质量值和参考碱基所在读段的碱基质量值作秩和检验，过滤去除秩和小于质量值秩和阈值的突变。
23.本技术的一种实现方式中，根据10)杂合位点突变位置的秩和检验对突变进行过滤包括，对突变碱基所在读段中的位置和参考碱基所在读段中的位置作秩和检验，过滤去除秩和小于位置秩和阈值的突变。
24.本技术的一种实现方式中，突变频率的计算公式为，
25.af＝alt/覆盖该位点的读段总数，
26.或者，
27.af＝alt/(ref+所有alt)
28.其中，af为突变频率，alt为同一位点突变碱基的数量，ref同一位点参考碱基的数量，所有alt是指在同一位置有多个突变的情况下所有突变的alt。
29.本技术的一种实现方式中，突变频率阈值为0.35，位置阈值为8。
30.本技术的一种实现方式中，碱基质量值阈值为22.5。
31.本技术的一种实现方式中，均方根阈值为40。
32.本技术的一种实现方式中，突变的单位深度质量值的计算方式为，将突变的质量值除以所有含有突变的样本的深度之和的商作为单位深度质量值。
33.本技术的一种实现方式中，单位深度质量值阈值为2。
34.本技术的一种实现方式中，比对错配数量阈值为5.25。
35.本技术的一种实现方式中，根据8)信噪比对突变进行过滤时，如果碱基质量值小于碱基质量值阈值的读段数量为0，则将其赋值为0.5。
36.本技术的一种实现方式中，信噪比阈值为1.5。
37.本技术的一种实现方式中，质量值秩和阈值为-12.5。
38.本技术的一种实现方式中，位置秩和阈值为-8。
39.需要说明的是，以上使用的具体的过滤阈值可根据实验条件和测序平台进行调整，比如突变平均碱基质量值过滤阈值可上调或下调。
40.本技术的一种实现方式中，根据11)突变的链分布偏差对突变进行过滤包括，只使
用平均碱基质量值大于22.5的读段，计算突变的频率，标记为hiaf；分别计算在正链的参考碱基数量，标记为srf，在负链的参考碱基数量，标记为srr，在正链的突变碱基数量，标记为saf，在负链的突变碱基数量，标记为sar；对于参考碱基，当srf和srr任意一个为0，且srf+srr《＝12时，参考碱基偏差值refbias为0；计算srf/(srf+srr)，srr/(srf+srr)，当以上值都大于0.05，且srf和srr都大于2，参考碱基偏差值refbias为2，否则为1；对于突变碱基，当saf和sar任意一个为0，且saf+sar《＝12时，突变碱基偏差值altbias等于0；计算saf/(saf+sar)，sar/(saf+sar)，当以上值都大于0.05，且saf和sar都大于2，突变碱基偏差值altbias为2，否则为1；对srf、srr、saf和sar进行fisher精确检验，得到p值和优势比值odd ratio；当hiaf小于0.25，refbias为2，且altbias为1，p值小于等于0.01，odd ratio大于5或等于0，且突变长度小于100bp时，过滤去除该突变。
41.本技术的一种实现方式中，由线粒体相关疾病数据库mseqdr和临床疾病和表型相关突变数据库clinvar的线粒体突变数据，保留致病、可能致病或意义不明的突变，进行合并和去除重复，构成已知突变白名单。
42.需要说明的是，本技术的已知突变白名单由线粒体相关疾病数据库mseqdr和临床疾病和表型相关突变数据库clinvar的线粒体突变数据，经过合并和去除重复获得；可以理解，原则上，也可以将其他经过严格验证的突变加入到本技术的已知突变白名单中，不仅限于mseqdr数据库和clinvar数据库。因此，本技术的已知突变白名单可以根据现有技术不断调整。
43.本技术的第二方面公开了一种校正线粒体基因组测序突变的装置，包括白名单比对模块、突变过滤模块、校正后数据获取模块；其中，白名单比对模块，包括用于将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；突变过滤模块，包括用于对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；校正后数据获取模块，包括用于合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。
44.需要说明的是，本技术校正线粒体基因组测序突变的装置，实际上就是通过各模块分别实现本技术的校正线粒体基因组测序突变的方法；因此，各模块的具体限定可以参考本技术的校正线粒体基因组测序突变的方法，例如，11个过滤指标具体如何进行过滤、各过滤指标中的具体阈值等都可以参考本技术的校正方法，在此不累述。
45.本技术的第三方面公开了一种校正线粒体基因组测序突变的装置，该装置包括存储器和处理器；存储器，包括用于存储程序；处理器，包括用于通过执行该存储器存储的程序以实现本技术的校正线粒体基因组测序突变的方法。
46.本技术的第四方面公开了一种计算机可读存储介质，该存储介质中存储有程序，该程序能够被处理器执行以实现本技术的校正线粒体基因组测序突变的方法。
47.由于采用以上技术方案，本技术的有益效果在于：
48.本技术的线粒体基因组测序突变校正方法，先通过已知突变白名单保留可能致病的突变，避免遗漏已知的可能致病突变，再通过11个过滤指标去除了因实验方法不当、测序误差或比对软件的系统误差导致的假阳性突变，提高了最终获得的突变集合的可信度和真实性。本技术的校正方法，在减少突变解读工作量的同时，提高了分析结果的准确性；并且，本技术的校正方法不受样本数量的限制，能够对单样本进行校正。
附图说明
49.图1是本技术实施例中线粒体基因组测序突变校正方法的流程框图；
50.图2是本技术实施例中线粒体基因组测序突变校正装置的结构框图；
51.图3是本技术实施例中vcf文件示意图；
52.图4是本技术实施例中bam文件示意图。
具体实施方式
53.下面通过具体实施方式结合附图对本技术作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本技术能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他装置、材料、方法所替代。在某些情况下，本技术相关的一些操作并没有在说明书中显示或者描述，是为了避免本技术的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。
54.现有的线粒体基因组测序突变校正方法，例如，vqsr无法对单样本进行校正和过滤，并且，样本数量较少的情况下，也难以采用vqsr进行校正和过滤。此外，虽然有研究采用突变质量值或突变的链分布偏差进行过滤，但是，总体来说，现有的过滤校正方法都存在准确性较差、假阳性多等问题。
55.基于以上研究和认识，本技术创造性的提出了一种新的校正线粒体基因组测序突变的方法，如图1所示，包括白名单比对步骤11、突变过滤步骤12、校正后数据获取步骤13。
56.白名单比对步骤11，包括将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集。
57.突变过滤步骤12，包括对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集。
58.校正后数据获取步骤13，包括合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变。本技术中，已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。
59.本技术的关键在于：
60.1)获取人类线粒体基因组已知突变数据，对数据进行格式转换和去除重复，构建已知突变白名单。
61.2)保留输入样本数据出现在白名单中的突变。对于未出现在白名单中的突变，使用突变所在读段的位置的平均值、突变的链分布偏差、突变的平均碱基质量值等共11个过滤指标，去除由于实验错误和测序错误导致的假阳性，得到最终真实的突变合集。
62.本技术的一种实现方式中，本技术校正线粒体基因组测序突变的方法，具体包括：
63.1、获取mseqdr(线粒体相关疾病数据库)和clinvar(临床疾病和表型相关突变数据库)的线粒体突变数据，保留致病、可能致病或意义不明的突变，进行合并和去除重复，构成已知突变白名单。
64.2、若样本突变出现在已知突变白名单中，将不会进行以下的过滤步骤，否则，对没有出现在已知突变白名单中的突变进行以下所有过滤指标的过滤：
65.1)同一位点只保留突变频率最高的突变：记同一位点突变碱基的数量为alt，同一位点参考碱基的数量为ref，则突变的频率af＝alt/覆盖该位点的读段总数＝alt/(ref+所有alt)，只保留同一位点中突变频率最高的突变。
66.其中，所有alt是引物，alt或突变在同一位置可能有多个，这种情况下，需要将同一位置的多个突变的alt相加进行计算。
67.2)突变所在读段的位置的平均值：统计同一突变位点在不同读段中的位置，计算这些位置的平均值，过滤去除突变频率小于突变频率阈值，且位置的平均值小于位置阈值的突变。
68.本技术的一种实现方式中，突变频率阈值为0.35，位置阈值为8。即过滤去除突变频率af《0.35，且突变位置发生在读段的前7个bp的位置。需要说明的是，位置阈值为8是本技术根据经验和统计研究获得，一般读段的长度都是百余bp，8这个数值较小，表示这个突变出现在读段的前面，研究显示，当af较小，如af《0.35，且突变出现在前面，则表明这个突变很可能是由于测序错误引起，所以将其过滤去除。
69.3)突变是否都位于同一位置：统计突变在所有读段的位置，若突变在所有读段都位于同一位置，则该突变被过滤去除。
70.需要说明的是，本技术研究认为，如果突变在所有读段中的位置都相同，那意味着这个突变是由于测序错误引起的，因此本技术创造性的过滤去除突变位点在不同的读段中的位置都相同的突变。
71.4)突变的平均碱基质量值：计算突变的平均碱基质量值，若突变的平均碱基质量值低于碱基质量值阈值，则该突变被过滤去除。
72.本技术的一种实现方式中，碱基质量值阈值为22.5，即过滤去除平均碱基质量值低于22.5的突变。需要说明的是，22.5这个值也是本技术根据经验和统计研究获得；研究显示，如果平均碱基质量值过小，如小于22.5，说明那些位置的测序质量较差，则这个突变很可能为假阳性，所以将其过滤去除。
73.5)读段比对质量值的均方根：计算突变所在的所有读段的对比质量值的均方根，若该值小于均方根阈值，则该突变被过滤去除。
74.本技术的一种实现方式中，均方根阈值为40，即过滤去除对比质量值的均方根小于40的突变。需要说明的是，本技术的均方根阈值参考gatk的vqsr的rmsmq值分布，参考该
分布，本技术创造性的提出均方根阈值为40时，能有效的把假突变过滤。
75.6)突变的单位深度质量值：突变的质量值除以所有含有突变的样本的深度之和的商，即单位深度质量值，若该值小于单位深度质量值阈值，则该突变被过滤去除。
76.本技术的一种实现方式中，单位深度质量值阈值为2，即过滤去除单位深度质量值小于2的突变。本技术的单位深度质量值阈值是参考gatk的vqsr的qd值分布，参考该分布，本技术创造性的提出单位深度质量值阈值为2时，能有效的过滤去除假突变。
77.7)读段平均比对错配数量：当突变所在的读段比对错配数量平均值大于比对错配数量阈值时，该突变很可能是由于对比错误而出现，因此将该突变过滤去除。
78.本技术的一种实现方式中，比对错配数量阈值为5.25，即过滤去除比对错配数量平均值大于5.25的突变。需要说明的是，比对错配数量阈值是本技术根据经验和统计研究获得，研究显示，错配数量过多，很可能说明这些读段比对的位置是错误的，即这些snv产生于序列比对，而非真实存在；因此，本技术研究认为过滤去除比对错配数量平均值大于5.25的突变，能够有效去除由此导致的假阳性。
79.8)信噪比：突变所在的读段中，碱基质量值大于碱基质量值阈值的读段数量，除以碱基质量值小于碱基质量值阈值的读段数量(若该值为0，则赋予0.5的值)，即为信噪比，若该比值小于信噪比阈值，则该突变被过滤去除。
80.本技术的一种实现方式中，碱基质量值阈值为22.5，信噪比阈值为1.5。即将碱基质量值大于碱基质量值阈值的读段数量除以碱基质量值小于碱基质量值阈值的读段数量的商作为信噪比，过滤去除信噪比值小于1.5的突变。需要说明的是，本技术的信噪比阈值也是由经验和统计研究获得，研究显示，信噪比越高，意味着该突变是由测序错误引起的概率越低；反之信噪比越低则突变由测序错误引起的概率越高；因此，本技术研究认为，过滤去除信噪比值小于1.5的突变，能够有效的去除测序错误引起突变。
81.9)杂合位点读段的碱基质量值秩和检验：对突变碱基所在读段的碱基质量值和参考碱基所在读段的碱基质量值作秩和检验，若秩和小于质量值秩和阈值，则该突变被过滤去除。
82.本技术的一种实现方式中，质量值秩和阈值为-12.5，即过滤去除秩和小于-12.5的突变。需要说明的是，本技术的质量值秩和阈值参考gatk的vqsr的mqranksum分布，参考该分布，本技术创造性的提出质量值秩和阈值为-12.5时，能有效的把假突变过滤。
83.10)杂合位点突变位置的秩和检验：对突变碱基所在读段中的位置和参考碱基所在读段中的位置作秩和检验，若秩和小于位置秩和阈值，则过滤去除该突变。
84.本技术的一种实现方式中，位置秩和阈值为-8，即过滤去除秩和小于-8的突变。需要说明的是，该位置秩和阈值参考gatk的vqsr的readposranksum分布，参考该分布，本技术创造性的提出位置秩和阈值为-8时，能有效的过滤去除假突变。
85.11)突变的链分布偏差：
86.只使用平均碱基质量值大于22.5的读段，计算出突变的频率，标记为hiaf。hiaf的计算公式与本技术的af相同，不同之处在于，前面的使用所有读段，而hiaf只使用平均碱基质量值大于22.5的读段。
87.分别计算在正链的参考碱基数量srf、在负链的参考碱基数量srr、在正链的突变碱基数量saf和在负链的突变碱基数量sar。
88.对于参考碱基，当srf和srr任意一个为0且srf+srr《＝12时，参考碱基偏差值refbias为0。计算srf/(srf+srr)，srr/(srf+srr)，当以上值都大于0.05且srf和srr都大于2，参考碱基偏差值refbias为2，否则为1。
89.对于突变碱基，当saf和sar任意一个为0且saf+sar《＝12时，突变碱基偏差值altbias等于0。计算saf/(saf+sar)，sar/(saf+sar)，当以上值都大于0.05且saf和sar都大于2，突变碱基偏差值altbias为2，否则为1。
90.对srf、srr、saf和sar进行fisher精确检验，得到p值和odd ratio(优势比值)。
91.当hiaf小于0.25，refbias为2且altbias为1，p值小于等于0.01，odd ratio大于5(或odd ratio等于0)，且突变长度小于100bp时，说明突变碱基极端分布在正链或负链，则过滤去除该突变。
92.需要说明的是，根据突变的链分布偏差对突变进行过滤中，所有采用的具体数值都是本技术由经验和统计研究获得；其中，p值《0.01是统计学上表示差异显著的阈值；优势比odd ratio》5或等于0，同样表示两组对象差异显著。
93.本领域技术人员可以理解，上述方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述方法中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述方法中全部或部分功能。
94.因此，基于本技术的校正线粒体基因组测序突变的方法，本技术提出了一种校正线粒体基因组测序突变的装置，如图2所示，包括白名单比对模块21、突变过滤模块22和校正后数据获取模块23。
95.白名单比对模块21，包括用于将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集。
96.突变过滤模块22，包括用于对没有出现在已知突变白名单中的突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集。
97.校正后数据获取模块23，包括用于合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变。
98.本技术的另一种实现方式中还提供了一种校正线粒体基因组测序突变的装置，该装置包括存储器和处理器；存储器，包括用于存储程序；处理器，包括用于通过执行存储器存储的程序以实现以下方法：将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；对
没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。
99.本技术另一种实现方式中还提供一种计算机可读存储介质，该存储介质中包括程序，该程序能够被处理器执行以实现如下方法：将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。
100.下面通过具体试验对本技术作进一步详细说明。以下试验仅对本技术进行进一步说明，不应理解为对本技术的限制。
101.实施例
102.首先，获取mseqdr(线粒体相关疾病数据库)和clinvar(临床疾病和表型相关突变数据库)的线粒体突变数据，保留致病、可能致病或意义不明的突变，进行合并和去除重复，构成已知突变白名单。
103.具体的，通过编写python脚本，获取mseqdr数据库的全部数据，保留chr列为mt且effect列为+/+(即致病突变)的行，将这些数据放入临时文件1。编写python脚本将临时文件1的数据进行格式转换，从而获取每个突变的起始，终止，参考碱基，突变碱基等信息，得到临时文件2。从clinvar的ftp网址下载variant_summary.txt文件，编写python脚本进行格式转换，只保留chromosome列为mt，clinicalsignificance列包含pathogenic或uncertain字段，不区分大小写，且reviewstatus列不包含no assertion字段的行，从而获取临时文件3。将临时文件2和临时文件3合并，进行去重，生成最终的已知突变白名单。
104.本例获得的已知突变白名单总共包含1473个突变，都是已知确定的人类线粒体基因组的致病、可能致病或意义未明的突变。
105.将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集。
106.对于没有出现在已知突变白名单中的突变，进行如下过滤：
107.1)同一位点只保留突变频率最高的突变：记同一位点突变碱基的数量为alt，同一位点参考碱基的数量为ref，则突变的频率af＝alt/覆盖该位点的读段总数＝alt/(ref+所有alt)，只保留同一位点中突变频率最高的突变。
108.其中，所有alt是引物，alt或突变在同一位置可能有多个，这种情况下，需要将同
一位置的多个突变的alt相加进行计算。
109.2)突变所在读段的位置的平均值：统计同一突变位点在不同读段中的位置，计算这些位置的平均值，过滤去除突变频率小于0.35，且位置的平均值小于8的突变。
110.3)突变是否都位于同一位置：统计突变在所有读段的位置，若突变在所有读段都位于同一位置，则该突变被过滤去除。
111.4)突变的平均碱基质量值：计算突变的平均碱基质量值，若突变的平均碱基质量值低于22.5，则该突变被过滤去除。
112.5)读段比对质量值的均方根：计算突变所在的所有读段的对比质量值的均方根，若该值小于40，则该突变被过滤去除。
113.6)突变的单位深度质量值：突变的质量值除以所有含有突变的样本的深度之和的商，即单位深度质量值，若该值小于2，则该突变被过滤去除。
114.7)读段平均比对错配数量：当突变所在的读段比对错配数量平均值大于5.25时，该突变很可能是由于对比错误而出现，因此将该突变过滤去除。
115.8)信噪比：突变所在的读段中，碱基质量值大于22.5的读段数量，除以碱基质量值小于22.5的读段数量(若该值为0，则赋予0.5的值)，即为信噪比，若该比值小于1.5，则该突变被过滤去除。
116.9)杂合位点读段的碱基质量值秩和检验：对突变碱基所在读段的碱基质量值和参考碱基所在读段的碱基质量值作秩和检验，若秩和小于-12.5，则该突变被过滤去除。
117.10)杂合位点突变位置的秩和检验：对突变碱基所在读段中的位置和参考碱基所在读段中的位置作秩和检验，若秩和小于-8，则过滤去除该突变。
118.11)突变的链分布偏差：
119.只使用平均碱基质量值大于22.5的读段，计算出突变的频率，标记为hiaf。hiaf的计算公式与本技术的af相同，不同之处在于，前面的使用所有读段，而hiaf只使用平均碱基质量值大于22.5的读段。
120.分别计算在正链的参考碱基数量srf、在负链的参考碱基数量srr、在正链的突变碱基数量saf和在负链的突变碱基数量sar。
121.对于参考碱基，当srf和srr任意一个为0且srf+srr《＝12时，参考碱基偏差值refbias为0。计算srf/(srf+srr)，srr/(srf+srr)，当以上值都大于0.05且srf和srr都大于2，参考碱基偏差值refbias为2，否则为1。
122.对于突变碱基，当saf和sar任意一个为0且saf+sar《＝12时，突变碱基偏差值altbias等于0。计算saf/(saf+sar)，sar/(saf+sar)，当以上值都大于0.05且saf和sar都大于2，突变碱基偏差值altbias为2，否则为1。
123.对srf、srr、saf和sar进行fisher精确检验，得到p值和odd ratio(优势比值)。
124.当hiaf小于0.25，refbias为2且altbias为1，p值小于等于0.01，odd ratio大于5(或odd ratio等于0)，且突变长度小于100bp时，说明突变碱基极端分布在正链或负链，则过滤去除该突变。
125.经过以上11个指标的过滤，剩余的突变即第二突变集。合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变。
126.以上过程，输入样本的vcf文件以及sam或bam文件，vcf文件可由突变检测软件比
如freebayes，gatk，vardict等生成，sam或bam文件可由序列比对文件比如bwa，bowtie2等生成，执行过滤程序，通过白名单和过滤指标的过滤，得到过滤后的vcf文件，即校正后的线粒体基因组测序突变。vcf文件一般如图3所示，sam文件一般如图4所示。
127.本例具体选取了10个线粒体基因组二代测序样本，统计过滤前后突变数量、阳性突变的突变数量、精确率和虚警率的变化，结果如表1所示。其中，精确率是指，真阳性突变在过滤前/后的突变的比例；虚警率是指假阳性突变在过滤前/后的突变的比例。
128.表1线粒体基因组二代测序突变校正前后的结果统计
[0129][0130]
表1的结果显示，过滤后平均每个样本减少23.1个突变，精确率平均提升27.2％，虚警率平均降低26.3％，显著地减少了下游需要分析的数据量，提高了真实阳性突变的比例，降低了虚假突变的比例。虽然平均每个样本有1.1个突变会被过滤去除，但经过验证，被遗漏的突变都是人群多态性的突变，即没有致病性的突变。比如as75058样本的nc_012920.1:m.16179c》t突变。因此，本例的线粒体基因组测序突变校正方法，能够有效的去除因实验方法不当、测序误差或比对软件的系统误差导致的假阳性突变，得到可信度高的突变集合，在减少突变解读工作量的同时，提高了分析结果的准确性。并且，本例的校正方法能够对单样本或数量较少的样本进行校正。
[0131]
以上内容是结合具体的实施方式对本技术所作的进一步详细说明，不能认定本技术的具体实施只局限于这些说明。对于本技术所属技术领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干简单推演或替换。

技术特征：
1.一种校正线粒体基因组测序突变的方法，其特征在于：包括以下步骤，将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；所述已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。2.根据权利要求1所述的方法，其特征在于：根据1)突变频率对突变进行过滤包括，计算所有突变的突变频率，对于同一突变位点，只保留同一突变位点中突变频率最高的突变；优选的，根据2)突变所在读段的位置的平均值对突变进行过滤包括，统计同一突变位点在不同读段中的位置，计算这些位置的平均值，过滤去除突变频率小于突变频率阈值，且位置的平均值小于位置阈值的突变；优选的，根据3)突变是否都位于同一位置对突变进行过滤包括，对于同一突变位点，过滤去除突变位点在不同的读段中的位置都相同的突变；优选的，根据4)突变的平均碱基质量值对突变进行过滤包括，计算突变的平均碱基质量值，过滤去除平均碱基质量值低于碱基质量值阈值的突变；优选的，根据5)读段比对质量值的均方根对突变进行过滤包括，计算突变所在的所有读段的对比质量值的均方根，过滤去除对比质量值均方根小于均方根阈值的突变；优选的，根据6)突变的单位深度质量值对突变进行过滤包括，计算突变的单位深度质量值，过滤去除单位深度质量值小于单位深度质量值阈值的突变；优选的，根据7)读段平均比对错配数量对突变进行过滤包括，计算突变所在读段的平均比对错配数量，过滤去除平均比对错配数量大于比对错配数量阈值的突变；优选的，根据8)信噪比对突变进行过滤包括，统计突变所在的读段中，碱基质量值大于碱基质量值阈值的读段数量，碱基质量值小于碱基质量值阈值的读段数量，将碱基质量值大于碱基质量值阈值的读段数量除以碱基质量值小于碱基质量值阈值的读段数量的商作为信噪比，过滤去除信噪比值小于信噪比阈值的突变；优选的，根据9)杂合位点读段的碱基质量值秩和检验对突变进行过滤包括，对突变碱基所在读段的碱基质量值和参考碱基所在读段的碱基质量值作秩和检验，过滤去除秩和小于质量值秩和阈值的突变；优选的，根据10)杂合位点突变位置的秩和检验对突变进行过滤包括，对突变碱基所在读段中的位置和参考碱基所在读段中的位置作秩和检验，过滤去除秩和小于位置秩和阈值的突变。3.根据权利要求2所述的方法，其特征在于：所述突变频率的计算公式为，af＝alt/覆盖该位点的读段总数，或者，af＝alt/(ref+所有alt)
其中，af为突变频率，alt为同一位点突变碱基的数量，ref同一位点参考碱基的数量，所有alt是指在同一位置有多个突变的情况下所有突变的alt。4.根据权利要求2所述的方法，其特征在于：所述突变频率阈值为0.35，所述位置阈值为8；优选的，所述碱基质量值阈值为22.5；优选的，所述均方根阈值为40；优选的，所述突变的单位深度质量值的计算方式为，将突变的质量值除以所有含有突变的样本的深度之和的商作为单位深度质量值；优选的，所述单位深度质量值阈值为2；优选的，所述比对错配数量阈值为5.25；优选的，根据8)信噪比对突变进行过滤时，如果碱基质量值小于碱基质量值阈值的读段数量为0，则将其赋值为0.5；优选的，所述信噪比阈值为1.5；优选的，所述质量值秩和阈值为-12.5；优选的，所述位置秩和阈值为-8。5.根据权利要求1-4任一项所述的方法，其特征在于：根据11)突变的链分布偏差对突变进行过滤包括，只使用平均碱基质量值大于22.5的读段，计算突变的频率，标记为hiaf；分别计算在正链的参考碱基数量，标记为srf，在负链的参考碱基数量，标记为srr，在正链的突变碱基数量，标记为saf，在负链的突变碱基数量，标记为sar；对于参考碱基，当srf和srr任意一个为0，且srf+srr<＝12时，参考碱基偏差值refbias为0；计算srf/(srf+srr)，srr/(srf+srr)，当以上值都大于0.05，且srf和srr都大于2，参考碱基偏差值refbias为2，否则为1；对于突变碱基，当saf和sar任意一个为0，且saf+sar<＝12时，突变碱基偏差值altbias等于0；计算saf/(saf+sar)，sar/(saf+sar)，当以上值都大于0.05，且saf和sar都大于2，突变碱基偏差值altbias为2，否则为1；对srf、srr、saf和sar进行fisher精确检验，得到p值和优势比值odd ratio；当hiaf小于0.25，refbias为2，且altbias为1，p值小于等于0.01，odd ratio大于5或等于0，且突变长度小于100bp时，过滤去除该突变。6.根据权利要求1-4任一项所述的方法，其特征在于：由线粒体相关疾病数据库mseqdr和临床疾病和表型相关突变数据库clinvar的线粒体突变数据，保留致病、可能致病或意义不明的突变，进行合并和去除重复，构成所述已知突变白名单。7.一种校正线粒体基因组测序突变的装置，其特征在于：包括白名单比对模块、突变过滤模块、校正后数据获取模块；所述白名单比对模块，包括用于将待校正的线粒体基因组测序突变与已知突变白名单进行比较，保留出现在已知突变白名单中的待校正的线粒体基因组测序突变，作为第一突变集；所述突变过滤模块，包括用于对没有出现在已知突变白名单中的待校正的线粒体基因组测序突变，根据1)突变频率、2)突变所在读段的位置的平均值、3)突变是否都位于同一位
置、4)突变的平均碱基质量值、5)读段比对质量值的均方根、6)突变的单位深度质量值、7)读段平均比对错配数量、8)信噪比、9)杂合位点读段的碱基质量值秩和检验、10)杂合位点突变位置的秩和检验以及11)突变的链分布偏差对突变进行过滤，将过滤后剩余的突变作为第二突变集；所述校正后数据获取模块，包括用于合并第一突变集和第二突变集，即获得校正后的线粒体基因组测序突变；所述已知突变白名单为经过验证的已知为致病、可能致病或意义不明的突变集合。8.根据权利要求7所述的装置，其特征在于：根据1)突变频率对突变进行过滤包括，计算所有突变的突变频率，对于同一突变位点，只保留同一突变位点中突变频率最高的突变；优选的，根据2)突变所在读段的位置的平均值对突变进行过滤包括，统计同一突变位点在不同读段中的位置，计算这些位置的平均值，过滤去除突变频率小于突变频率阈值，且位置的平均值小于位置阈值的突变；优选的，根据3)突变是否都位于同一位置对突变进行过滤包括，对于同一突变位点，过滤去除突变位点在不同的读段中的位置都相同的突变；优选的，根据4)突变的平均碱基质量值对突变进行过滤包括，计算突变的平均碱基质量值，过滤去除平均碱基质量值低于碱基质量值阈值的突变；优选的，根据5)读段比对质量值的均方根对突变进行过滤包括，计算突变所在的所有读段的对比质量值的均方根，过滤去除对比质量值均方根小于均方根阈值的突变；优选的，根据6)突变的单位深度质量值对突变进行过滤包括，计算突变的单位深度质量值，过滤去除单位深度质量值小于单位深度质量值阈值的突变；优选的，根据7)读段平均比对错配数量对突变进行过滤包括，计算突变所在读段的平均比对错配数量，过滤去除平均比对错配数量大于比对错配数量阈值的突变；优选的，根据8)信噪比对突变进行过滤包括，统计突变所在的读段中，碱基质量值大于碱基质量值阈值的读段数量，碱基质量值小于碱基质量值阈值的读段数量，将碱基质量值大于碱基质量值阈值的读段数量除以碱基质量值小于碱基质量值阈值的读段数量的商作为信噪比，过滤去除信噪比值小于信噪比阈值的突变；优选的，根据9)杂合位点读段的碱基质量值秩和检验对突变进行过滤包括，对突变碱基所在读段的碱基质量值和参考碱基所在读段的碱基质量值作秩和检验，过滤去除秩和小于质量值秩和阈值的突变；优选的，根据10)杂合位点突变位置的秩和检验对突变进行过滤包括，对突变碱基所在读段中的位置和参考碱基所在读段中的位置作秩和检验，过滤去除秩和小于位置秩和阈值的突变；优选的，所述突变频率的计算公式为，af＝alt/覆盖该位点的读段总数，或者，af＝alt/(ref+所有alt)其中，af为突变频率，alt为同一位点突变碱基的数量，ref同一位点参考碱基的数量，所有alt是指在同一位置有多个突变的情况下所有突变的alt；优选的，所述突变频率阈值为0.35，所述位置阈值为8；
优选的，所述碱基质量值阈值为22.5；优选的，所述均方根阈值为40；优选的，所述突变的单位深度质量值的计算方式为，将突变的质量值除以所有含有突变的样本的深度之和的商作为单位深度质量值；优选的，所述单位深度质量值阈值为2；优选的，所述比对错配数量阈值为5.25；优选的，根据8)信噪比对突变进行过滤时，如果碱基质量值小于碱基质量值阈值的读段数量为0，则将其赋值为0.5；优选的，所述信噪比阈值为1.5；优选的，所述质量值秩和阈值为-12.5；优选的，所述位置秩和阈值为-8；优选的，根据11)突变的链分布偏差对突变进行过滤包括，只使用平均碱基质量值大于22.5的读段，计算突变的频率，标记为hiaf；分别计算在正链的参考碱基数量，标记为srf，在负链的参考碱基数量，标记为srr，在正链的突变碱基数量，标记为saf，在负链的突变碱基数量，标记为sar；对于参考碱基，当srf和srr任意一个为0，且srf+srr<＝12时，参考碱基偏差值refbias为0；计算srf/(srf+srr)，srr/(srf+srr)，当以上值都大于0.05，且srf和srr都大于2，参考碱基偏差值refbias为2，否则为1；对于突变碱基，当saf和sar任意一个为0，且saf+sar<＝12时，突变碱基偏差值altbias等于0；计算saf/(saf+sar)，sar/(saf+sar)，当以上值都大于0.05，且saf和sar都大于2，突变碱基偏差值altbias为2，否则为1；对srf、srr、saf和sar进行fisher精确检验，得到p值和优势比值odd ratio；当hiaf小于0.25，refbias为2，且altbias为1，p值小于等于0.01，odd ratio大于5或等于0，且突变长度小于100bp时，过滤去除该突变；优选的，由线粒体相关疾病数据库mseqdr和临床疾病和表型相关突变数据库clinvar的线粒体突变数据，保留致病、可能致病或意义不明的突变，进行合并和去除重复，构成所述已知突变白名单。9.一种校正线粒体基因组测序突变的装置，其特征在于：所述装置包括存储器和处理器；所述存储器，包括用于存储程序；所述处理器，包括用于通过执行所述存储器存储的程序以实现权利要求1-6任一项所述校正线粒体基因组测序突变的方法。10.一种计算机可读存储介质，其特征在于：所述存储介质中存储有程序，所述程序能够被处理器执行以实现权利要求1-6任一项所述校正线粒体基因组测序突变的方法。

技术总结
本申请公开了一种校正线粒体基因组测序突变的方法、装置和存储介质。本申请的线粒体基因组测序突变校正方法包括，先通过已知突变白名单保留可能致病的突变，再使用突变所在读段的位置的平均值、突变的链分布偏差、突变的平均碱基质量值等11个过滤指标，去除由于实验错误和测序错误导致的假阳性，从而得到可信度高的突变集合。本申请的校正方法在减少突变解读工作量的同时，提高了分析结果的准确性。并且，本申请的校正方法不受样本数量的限制，能够对单样本或数量较少的样本进行校正。够对单样本或数量较少的样本进行校正。够对单样本或数量较少的样本进行校正。


技术研发人员：黄凯 窦浩宇 刘永初 燕攀 刘阳 李阳
受保护的技术使用者：深圳雅济科技有限公司
技术研发日：2023.07.17
技术公布日：2023/10/5