具体实施方式

为了使本发明的发明目的、技术方案和有益技术效果更加清晰，以下结合实施例对本发明进行进一步详细说明，熟悉此技术的人士可由本说明书所揭露的内容容易地了解本申请发明的其他优点及功效。

本发明第一方面提供一种高通量测序突变检测方法，包括：

S1)获取基因样本中各候选突变位点的特征信息，所述特征信息来源于基因样本的高通量测序数据；

S2)基于各候选突变位点的特征信息和目标值s，将各候选突变位点进行层次聚类分析，其中，目标值s由各候选突变位点与其所在的分类内部其他候选突变位点的距离d的平均值a和各候选突变位点与最近的一个分类中的候选突变位点的距离d的平均值b计算获得；

S3)根据层次聚类分析结果，基于候选突变位点的突变丰度，获取各类候选突变位点的背景值pbg；

S4)根据层次聚类分析结果，基于候选突变位点的深度及其所在的分类的背景值pbg，获取各类候选突变位点的最低检测下限；

S5)根据各候选突变位点所在类别的背景值pbg及各候选突变位点的最低检测下限，确定各候选突变位点的检测阈值。

本发明所提供的高通量测序突变检测方法，可以包括：获取基因样本中各候选突变位点的特征信息，所述特征信息来源于基因样本的高通量测序数据。通常来说，本领域技术人员可选择合适的方法，根据基因样本的高通量测序数据，确定所需针对的各候选突变位点，并进一步获取基因样本中各候选突变位点的特征信息。例如，还可以包括：将基因样本的高通量测序数据与人类参考基因组数据进行比对，以鉴别基因样本中的候选突变位点。基因样本的高通量测序数据可以是Fastq文件等，其通常可以从高通量测序的下机数据(例如，BCL格式的数据等)转换获得(例如，通过BCL2fastq等软件)，将基因样本的高通量测序数据与人类参考基因组数据进行比对(例如，通过BWA算法等进行对比)，比对结果可以进一步被转换为BAM文件(例如，通过samtools等软件进行转换)，根据比对结果即可以进一步鉴别基因样本中的候选突变位点(例如，可以通过VarDict算法等进行鉴别)。再例如，特征信息通常可以从适当的文件中提取(例如，BAM文件等)，特征信息具体可以包括深度、变异深度、正链参考碱基深度、负链参考碱基深度、正链变异碱基深度、负链变异碱基深度、基因型、突变丰度、链偏差、读取片段上的位置、读取片段上的位置的标准差、平均碱基质量分数、碱基质量分数标准差、比对质量、高质量片段比例、高质量片段突变丰度、是否为微卫星位点、微卫星位点单元长度、片段上的错配总数、5`端的序列、3`端的序列、突变类型、重复比例等中的一种或多种的组合，这些特征的计算方式可以参照Lai Z,Markovets A,Ahdesmaki M,Chapman B,Hofmann O,McEwen R,Johnson J,Dougherty B,Barrett JC,DryJR.VarDict:a novel and versatile variant caller for next-generationsequencing in cancer research.Nucleic Acids Res.2016Jun 20；44(11):e108.doi:10.1093/nar/gkw227.Epub 2016Apr 7.PMID:27060149；PMCID:PMC4914105.等。各特征信息的中英文对应名称及特征描述可以如表1所示。

表1

上述高通量测序突变检测方法中，所获得的各候选突变位点的特征信息，可以以矩阵的方式给出，例如，可以构成一个所有候选突变位点的特征矩阵M(Mi,j)，其中，Mi,j为第i个候选突变位点的第j个特征信息的具体数值。

本发明所提供的高通量测序突变检测方法，还可以包括：基于各候选突变位点的特征信息和目标值s，将各候选突变位点进行层次聚类分析(Hierarchical clustering)，其中，目标值s由各候选突变位点与其所在的分类内部其他候选突变位点的距离d的平均值a和各候选突变位点与最近的一个分类中的候选突变位点的距离d的平均值b计算获得。在层次聚类分析中，聚类层高阈值通常可以为h，通常可以选取h采用自适应的方式进行优化，h的取值通常满足min(h)<＝h<＝max(h)，h的取值范围通常可以根据待优化的集(例如，上述的特征矩阵M(Mi,j))进行确定，上述算法可以来自于层次聚类分析，可以通过软件R的hclust函数得到，优化的目标值则如上所述可以为目标值s，当h在一定的取值范围内变化时，不同的层次聚类分析结果均会对应不同的目标值s。通常来说，目标值s越小则说明得到的聚类内部差异越小，说明聚类效果更好。例如，所提供的层次聚类分析结果为目标值s最小的情况下的层次聚类分析结果。再例如，目标值s的计算方法可以如下所示：

其中，a为各候选突变位点与其所在的分类内部其他候选突变位点的距离d的平均值；

b为各候选突变位点与最近的一个分类中的候选突变位点的距离d的平均值；

上述公式中，两个位点之间的距离d通常为两个位点的丰度的差值，更具体可以为两个位点的丰度的差值的绝对值。

本发明所提供的高通量测序突变检测方法，还可以包括：根据层次聚类分析结果，基于候选突变位点的突变丰度，获取各类候选突变位点的背景值pbg。在获取了层次聚类分析结果以后，可以根据各候选突变位点的突变丰度以及他们的分类结果，基于同一类的候选突变位点的突变丰度，获取各类候选突变位点的背景值pbg。例如，各类候选突变位点的背景值pbg可以为该类中各候选突变位点的突变丰度的中位数。

本发明所提供的高通量测序突变检测方法，还可以包括：根据层次聚类分析结果，基于候选突变位点的深度及其所在的分类的背景值pbg，获取各类候选突变位点的最低检测下限。在获取了层次聚类分析结果以后，可以根据他们的分类结果，基于各候选突变位点的深度及其所在的分类的背景值pbg，获取各类候选突变位点的最低检测下限。例如，各候选突变位点的最低检测下限的计算方法可以为：

f＝ln(1-p)/-n

其中，f为候选突变位点的最低检测下限；

p为候选突变位点所在的分类的背景值pbg；

n为候选突变位点的深度。

本发明所提供的高通量测序突变检测方法，还可以包括：根据各候选突变位点所在类别的背景值pbg及各候选突变位点的最低检测下限，确定各候选突变位点的检测阈值。通常来说，各候选突变位点其所在类别的背景值pbg以及各候选突变位点的最低检测下限两者相比较大的数值，可以作为各候选突变位点的检测阈值，其原因在于，最低检测下限为该位点深度下能达到的理论最低值，背景值决定了背景噪音的大小，最低检测下限小于背景值的话以最低检测下限为阈值，否则最低只能达到背景值为下限。

本发明所提供的高通量测序突变检测方法，还可以包括：根据各候选突变位点的检测阈值，获取各候选突变位点的突变检测结果。通常来说，各候选突变位点的检测阈值可以对应于各候选突变位点特征信息中的突变丰度，并可以根据两者的比较结果，获取各候选突变位点的突变检测结果。例如，当候选突变位点的突变丰度≥候选突变位点的检测阈值，则可以认为基因样本中候选突变位点的突变为阳性。再例如，当候选突变位点的突变丰度＜候选突变位点的检测阈值，则可以认为基因样本中候选突变位点的突变为阴性。

本发明第二方面提供一种计算机可读存储介质，其上存储有计算机程序，所述程序被处理器执行时实现本发明第一方面所提供的高通量测序突变检测方法。

本发明第三方面提供一种设备，包括：处理器和存储器，所述存储器用于存储计算机程序，所述处理器用于执行所述存储器存储的计算机程序，以使所述设备执行本发明第一方面所提供的高通量测序突变检测方法。

本发明第四方面提供一种装置，所述装置包括：

特征信息获取模块，用于获取基因样本中各候选突变位点的特征信息，所述特征信息来源于基因样本的高通量测序数据；

层次聚类分析模块，用于基于各候选突变位点的特征信息和目标值s，将各候选突变位点进行层次聚类分析，其中，目标值s由各候选突变位点与其所在的分类内部其他候选突变位点的距离d的平均值a和各候选突变位点与最近的一个分类中的候选突变位点的距离d的平均值b计算获得；

背景值pbg计算模块，用于根据层次聚类分析结果，基于候选突变位点的突变丰度，获取各类候选突变位点的背景值pbg；

最低检测下限计算模块，用于根据层次聚类分析结果，基于候选突变位点的深度及其所在的分类的背景值pbg，获取各类候选突变位点的最低检测下限；

检测阈值计算模块，用于根据各候选突变位点所在类别的背景值pbg及各候选突变位点的最低检测下限，确定各候选突变位点的检测阈值；

优选的，还包括候选突变位点鉴别模块，用于将基因样本的高通量测序数据与人类参考基因组数据进行比对，以鉴别基因样本中的候选突变位点；

优选的，还包括突变检测结果计算模块，用于根据各候选突变位点的检测阈值，获取各候选突变位点的突变检测结果。

本发明中，上述装置中各模块的运行原理可以参照本发明第一方面所提供的高通量测序突变检测方法，在此不做赘述。

本申请所提供的高通量测序突变检测方法，可以通过自适应优化的算法，针对每个样本特有的数据特点智能的设定分析的阈值，在保证阳性位点检出的前提下尽可能的去除噪音造成的假阳性，可以更加精准的从样本(例如，肿瘤样本)中检测突变，从而具有良好的产业化前景。

下面通过实施例对本申请予以进一步说明，但并不因此而限制本申请的范围。

实施例1

以肿瘤样本检测数据下机开始的整个分析流程为例，具体如下：

1)测序下机数据通过BCL2fastq从测序BCL数据中分离样本数据并转换成fastq文件，使用参数bcl2fastq--barcode-mismatches 1-o./demultiplex--ignore-missing-bcls--no-lane-splitting；

2)Fastq文件通过BWA算法比对到人类参考基因组并利用samtools软件转换为BAM文件，使用参数bwa mem-t 16-R"@RG\tID:DNA\tLB:DNA\tSM:S2100019497-Plasma\tPL:ILLUMINA"-Mhuman_g1k_v37_decoy.fasta；

3)使用VarDict算法鉴别出样本的所有候选突变，使用参数VarDict-b bam-p-GREF-c1-S 2-E 3-g 5；

4)对于该样本中任意的一个候选突变，使用VarDict提取表一中的所有特征，构建特征矩阵；

5)使用R,hclust函数基于特征矩阵进行层次聚类，以0.01作为h变化的步长，计算每个h下的s值，并根据s，确定最佳的分类方式；

6)定义每个分类内的变异的丰度值中位数为该分类的背景值，并结合突变深度定义AF阈值，阈值定义为所在类别的背景值pbg以及各候选突变位点的最低检测下限两者相比较大的数值；

7)筛选候选突变列表，标记高于该阈值的突变为真实突变。

8)使用已知突变位点的标准样品(例如Horizon HD780)进行标准建库实验和测序，得到的测序数据经过上述处理，比对结果如表1所示：

表1

可见，本发明提供的筛选算法可以在检出真阳性位点的前提下明显降低假阳性位点的检出。

综上所述，本发明有效克服了现有技术中的种种缺点而具高度产业利用价值。

上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。