CN106909806B - 定点检测变异的方法和装置 - Google Patents

Abstract

本发明公开了一种定点检出变异的方法，包括：基于变异的已知信息，确定变异的指定位点和包含该变异的参考序列；获取待测样本的核酸的测序数据，测序数据包括多个读段；提取测序数据中包含指定位点的读段，获得指定读段；以指定读段上的指定位点为中心，往两端方向各延伸N个bp，获得指定片段，4≤N≤10；将指定片段与包含变异的参考序列进行比对，获得支持读段，支持读段为与包含变异的参考序列匹配的指定片段所在的读段；统计支持读段的量，基于支持读段的量判断所述变异是否存在。该方法基于关注读段中是否存在发生变异后应当具有的序列特征来进行变异定点检测，能够规避变异位点附近比对质量下降的问题，能够快速精确的检出变异。

Description

定点检测变异的方法和装置

技术领域

本发明涉及生物信息领域，具体的，本发明涉及定点检测变异的方法和装置，更具体的，本发明涉及一种定点检出变异的方法、一种定点检出变异的装置、一种检测融合基因突变的方法和一种检测融合基因突变的装置。

背景技术

癌症由遗传基因的改变导致，不同癌症、不同患者具有不同类型的基因变异，找到癌症患者的基因突变类型是个体化的治疗的基础，同时能够帮助我们更清晰的认识癌症的机理。

目前，临床上主要通过armsPCR方法来检测SNV、INDEL，通过FISH的方法来检测基因融合，这两种实验方法价格高，探针是针对特定突变设计的，难增加新的突变检测位点。

随着基因组学和生物信息学的不断发展，NGS高通量方法逐渐在这个领域内得到应用。利用高通量方法同时对患者的癌症组织和正常血细胞对照进行测序，首先在癌症组织中检测变异，然后去掉在对照中存在的germline变异(生殖细胞变异)，从而得到最终的somatic变异(体细胞变异)。在这种情况下，检测结果中会包含大量的临床意义未明的变异，这类变异对临床医生并没有有效的指导作用；检测过程中同时需要癌症组织和血细胞进行测序，增加了工作量；更重要的是INDEL附近的碱基的比对质量会下降，例如对EGFRc.2238_2248>GC这类肺癌中存在的复杂INDEL(complex INDEL)变异，缺失(deletion)后插入的GC碱基可能会比对到不同的位置，传统的变异检测方法对这种变异的检测很困难。

发明内容

依据本发明的一方面提供一种定点检出变异的方法，该方法包括：基于所述变异的已知信息，确定所述变异的指定位点和包含所述变异的参考序列；获取待测样本的核酸的测序数据，所述测序数据包括多个读段；提取所述测序数据中包含所述指定位点的读段，获得指定读段；以所述指定读段上的指定位点为中心，往两端方向各延伸N个bp，获得指定片段，4≤N≤10；将所述指定片段与所述包含所述变异的参考序列进行比对，获得支持读段，所述支持读段为与所述参考序列匹配的指定片段所在的读段；统计所述支持读段的量，基于所述支持读段的量判断所述变异是否存在。

依据本发明的另一方面提供一种计算机可读存储介质，用于存储供计算机执行的第一程序，本领域普通技术人员可以理解，在执行该第一程序时，通过指令相关硬件可完成上述定点检出变异的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的再一方面提供一种定点检出变异的装置，该装置包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；处理器，用于执行第一计算机可执行程序，所述第一计算机可执行程序的执行包括完成上述本发明一方面的定点检出变异的方法；存储单元，与所述数据输入单元、数据输出单元和处理器相连，用于存储数据，其中包括所述第一计算机可执行程序。

上述本发明一方面的方法、计算机可读存储介质和/或装置，基于关注读段中是否存在发生变异后应当具有的序列特征来进行定点变异检测，能够规避变异位点附近比对质量下降、变异位点周边比对存在干扰等问题，能够快速精确的检出变异。

依据本发明的一方面提供一种检测融合基因突变的方法，该方法包括：获取待测样本的测序结果，所述测序结果包括多个读段；提取所述测序结果中的割裂读段，所述割裂读段为同一读段的两部分分别匹配到参考序列两个不同位置的读段；分析匹配到所述参考序列上相同位置的割裂读段的数量，确定候选断点；定义所述参考序列上候选断点相应位置为第一融合基因位置，截取匹配到所述第一融合基因位置的割裂读段的不匹配所述第一融合基因位置的部分，以获得第一割裂片段，将所述第一割裂片段进行组装，获得第一一致性序列；将所述第一一致性序列与所述参考序列进行比对，定义所述第一一致性序列与所述参考序列匹配的位置为第二融合基因位置；截取匹配到所述第二融合基因位置的割裂读段的不匹配所述第二融合基因位置的部分，获得第二割裂片段，将所述第二割裂片段进行组装，获得第二一致性序列；将所述第二一致性序列与所述参考序列进行比对，若所述第二一致性序列与所述参考序列匹配的位置为所述第一融合基因位置，确定存在所述融合基因突变。

依据本发明的另一方面提供一种计算机可读存储介质，用于存储供计算机执行的第二程序，本领域普通技术人员可以理解，在执行该第二程序时，通过指令相关硬件可完成上述检测融合基因突变的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的再一方面提供一种检测融合基因突变的装置，该装置包括：数据输入模块，用于输入数据；数据输出模块，用于输出数据；处理器，用于执行第二计算机可执行程序，所述第二计算机可执行程序的执行包括完成上述本发明一方面的检测融合基因突变的方法；存储模块，与所述数据输入模块、数据输出模块和处理器相连，用于存储数据，其中包括所述第二计算机可执行程序。

利用上述本发明一方面的方法、计算机可读存储介质和/或装置，能够准确高效的检测融合基因突变。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解，其中：

图1显示本发明的一个实施例中的定点检出变异的方法的流程。

图2显示本发明的一个实施例中的定点检出变异的装置的示意图。

图3显示本发明的一个实施例中的基于不同的测序深度、利用模型公式进行计算绘制的ROC曲线。

图4显示本发明的一个实施例中的基于不同的等位基因频率、利用模型公式进行计算绘制的ROC曲线。

图5显示本发明的一个实施例中的基于BGISEQ-100测序平台的单样本测序数据确定变异检出的流程。

图6显示本发明的一个实施例中的变异检出部分结果。

图7显示本发明的一个实施例中自动生成的样本检测报告的示意图。

图8显示本发明的一个实施例中基于读段比对的SNV位点的检出图。

具体实施方式

参见图1，根据本发明的实施例提供的一种定点检出变异的方法，该方法包括以下步骤：

S10确定变异的指定位点和包含该变异的参考序列。

基于所述变异的已知信息，例如变异在参考基因组上的位置、类型、等位基因突变频率等，确定所述变异的指定位点和包含所述变异的参考序列，包括确定变异存在时应当出现的序列、序列的起始位置和序列的终止位置等。

所述变异选自SNP和INDEL中的至少一种。所称变异的指定位点指，存在该变异时具有的特征序列的至少一部分，指定位点可以是单核苷酸，也可以是多个核苷酸。所称的参考序列指预先确定的序列，可以是预先获得的待测样本所属生物类别的任意参考模板，例如，若待测样本来源的为人类个体，参考序列可选择NCBI数据库提供的HG19，进一步地，也可以预先配置包含更多参考序列的资源库，例如依据待测样本来源个体的状态、地域等因素选择或是测定组装出更接近的序列作为参考序列。所称包含变异的参考序列为存在该变异的参考序列，例如参考基因组存在该变异后变成的序列。

S12获取待测样本的核酸的测序数据。

获取待测样本的核酸的测序数据，所述测序数据包括多个读段。

所称的测序数据通过对核酸序列进行测序文库制备、上机测序获得。根据本发明的实施例，获取所述测序数据，包括：获取待测样本中的核酸，制备所述核酸的测序文库，对所述测序文库进行测序。测序文库的制备方法根据所选择的测序方法的要求进行，测序方法依据所选的测序平台的不同，可选择但不限于Illumina公司的Hisq2000/2500测序平台、Life Technologies公司的Ion Torrent平台和单分子测序平台，测序方式可以选择单端测序，也可以选择双末端测序，获得的下机数据是测读出来的片段，称为读段(reads)。

需要说明的是，上述S10和S12之间无先后顺序的限制，可先进行S10再进行S12，也可先进行S12再进行S10。

S14提取所述测序数据中包含所述指定位点的读段。

提取所述测序数据中包含所述指定位点的读段，获得指定读段。所称的指定读段也包含存在变异时应当具有的特征序列的至少一部分。

根据本发明的一个实施例，所述提取测序数据中包含所述指定位点的读段，获得指定读段，包括：将所述测序数据与包含所述变异的参考序列进行比对，获得比对结果，将比对结果中的比对到所述参考序列上对应的指定位点位置的读段为所述指定读段。比对可以利用已知比对软件进行，例如SOAP、BWA和TeraMap等。

根据本发明的一个实施例，在获得比对结果后，对比对结果中的reads进行去重，去除重复的reads，例如去除由于测序文库构建过程中的扩增而带来的重复片段，能够减小后续处理依据的数据量，利于基于比对结果进行快速定点检测检测。

S16以指定读段上的指定位点为中心，往两端方向各延伸N个bp，获得指定片段。

以所述指定读段上的指定位点为中心，往两端方向各延伸N个bp，获得指定片段，4≤N≤10。发明人经过大量分析验证，确定N的数值范围。确定的该延伸长度的范围，使获得的指定片段能够用于后续高效筛选确定出可靠的特定读段，以用于定点变异检出。N若小于4，会使后续获得的比对结果复杂度增加，增加后续分析难度；而N若大于10，会使后续获得的比对结果中的特定读段的数量大大减少，不利于后续基于统计准确判定变异是否存在。根据本发明的一个较佳实施例，使N＝5，即使得获得的指定片段达11bp左右，利于后续快速确定出可靠的、数目足够的特定读段以用于变异判定。

S18将所述指定片段与所述包含所述变异的参考序列进行比对，获得支持读段。

将所述指定片段与所述包含所述变异的参考序列进行比对，获得支持读段，所述支持读段为与所述参考序列匹配的指定片段所在的读段。

根据本发明的较佳实施例，获得指定读段之后，对所述指定读段进行过滤，其中包括过滤掉指定位点位于读段的末端N bp内的指定读段。如此，除去相对不可靠或者说难以确定是否可靠的数据，利于后续步骤的高效快速进行。

所称的匹配意同比对上。具体比对时，可以利用已知比对软件进行，例如SOAP、BWA和TeraMap等，本发明对此不作限制。在比对过程中，根据比对参数的设置，一对或一条reads最多允许有n个碱基错配(mismatch)，例如设置n为1或2，若reads中有超过n个碱基发生错配，则视为该对Reads无法比对到参考序列，或者，若错配的n个碱基全部位于reads对中的一个reads，则视为该reads对中的该reads无法比对到参考序列。

根据本发明的一个较佳实施例，所称的匹配为完全匹配，即指定片段与含变异的参考序列零错配，包含这些指定片段的读段为支持读段。即支持读段为支持变异的读段，为包含该变异发生时应当具有的特征序列的读段。

S20基于所述支持读段的量判断所述变异是否存在。

统计所述支持读段的量，基于所述支持读段的量判断所述变异是否存在。所称的支持读段的量，包括支持读段的数目、数目所占的比例、其上特定碱基的测序深度、碱基测序错误率等

根据本发明的一个实施例，所称的测序数据中的读段的长度不相同，例如测序数据是利用Life Technologies公司的Ion Torrent系列中的Proton测序平台进行测序获得的。

发明人发现，肺癌变异检测试剂盒针对特定位点的特定变异进行检测，现有的变异检出程序均是对整个区域所有位点进行循环遍历，为了保证大范围检出的准确性而将检出条件设置很高，但对特定位点的特定变异检出，检出精度可以提高。发明人还发现，Proton测序或者BGISEQ-100测序中，由于测序文库构建中，插入片段两端的接头(P接头和A接头)的不对称性，会使得测到的正链reads和负链reads必定来自于不同的模板，不会是同一个模板的PCR产物；故对于基因组上特定位置的特定变异，例如EGFR L858R、KRASG12C等，限定链偏向(strand bias)的限制可以很大程度上保证变异检测的真实性。所称的正链reads和负链reads是相对的，互为反向互补。

发明人据上述发现以及为解决以上问题而建立参考值模型，确定阈值(cutoff 值)，并证明了在此模型和cutoff值的情况下变异检测具有高的灵敏度(sensitivity)和特 异性(specificity)。所称模型基于以下两点假设而建立：(1)对于任一位点，假设参考基因 组对应的碱基为r∈{A,T,C,G}，(2)对于任一位点，假设覆盖该位点的所有读段的对应碱基 为b_i，碱基质量值为q_i，则对应的碱基错误率为i＝1,2,...,d d表示该位点对应 的测序深度。所称模型可表示为其中M₀表示所述位点的变异 不存在，该位点与参考基因组碱基不同是由于系统误差导致的；表示所述位点的变异 真实存在，所述变异为r突变为m，f为等位基因突变频率，既不为r、也不为m的碱基是系统误 差造成的；L(M₀)表示所述位点的测序数据的分布情况符合M₀的概率，表示该所述位点的测序数 据的分布情况符合的概率，

根据本发明的一个实施例，所述基于支持读段的数量判断所述变异是否存在，包括将所述支持读段的量代入到所述模型中，包括将所述支持读段的量代入到上述式I以确定L(M₀)，将所述支持读段的量代入到上述式III以确定计算获得参考值LOD(m,f)；将所述参考值与阈值比较，当所述参考值大于所述阈值，判定所述待测样本存在所述变异。

所称阈值的取值范围为0至10。利用所称模型确定阈值，可以通过设置置信度，例如通常设置为95％或99％，确定该置信度对应的参考值为阈值。当利用上述步骤检测某一待测样本计算得的参考值大于所述阈值时，表明定点检出的变异95％或99％可信。

根据本发明的实施例，利用ROC分析进行评估，确定阈值。ROC曲线(receiveroperating characteristic curve，接收者操作特征曲线)，是一种二元分类模型，即输出结果只有两种类别的模型。考虑一个二分问题，即将实例分成正类(positive)或负类(negative)，对一个二分问题来说，会出现四种情况：如果一个实例是正类并且也被预测成正类，即为真正类(True positive，TP)，如果实例是负类被预测成正类，称之为假正类(False positive，FP)，相应地，如果实例是负类被预测成负类，称之为真负类(Truenegative，TN)，正类被预测成负类则为假负类(false negative，FN)。TP：正确肯定的数目；FN：漏报，没有正确找到的匹配的数目；FP：误报，给出的匹配是不正确的；TN：正确拒绝的非匹配对数。为了形象化这一变化，引入ROC，ROC曲线可以用于评价一个分类器，即评价引入不同候选阈值的变异判定模型。AUC(Area Under roc Curve)为ROC曲线下方的面积，AUC值介于0.5到1.0之间，AUC越大，分类器分类效果越好。

发明人根据以上方法欲检测样本中是否存在某个/些已知变异，例如EGFR EX19某一区域是否存在15～18bp的缺失(deletion)，针对于之类已知的变异的基因型，该方法提高了检出的敏感性和特异性。若根据该方法检测结果判定不存在目标已知变异基因型，则可以利用敏感性和特异性相对较低的已知的变异检测方法/软件补充检测，判断在该区域是否有未知基因型存在。

根据本发明的实施例提供的一种计算机可读存储介质，用于存储供计算机执行的第一程序，本领域普通技术人员可以理解，在执行该第一程序时，通过指令相关硬件可完成上述定点检出变异的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

参见图2，根据本发明的实施例提供的一种定点检出变异的装置100，该装置100包括：数据输入单元110，用于输入数据；数据输出单元120，用于输出数据；处理器130，用于执行第一计算机可执行程序，所述第一计算机可执行程序的执行包括完成上述本发明一方面或者任一实施例的定点检出变异的方法；存储单元140，与所述数据输入单元、数据输出单元和处理器相连，用于存储数据，其中包括所述第一计算机可执行程序。

上述本发明任一实施例中的方法、计算机可读存储介质和/或装置，基于关注读段中是否存在发生变异后应当具有的序列特征来进行定点变异检测，能够规避变异位点附近比对质量下降、变异位点周边比对存在干扰等问题，能够快速精确的检出变异。

根据本发明的实施例提供的一种检测融合基因突变的方法，该方法包括：获取待测样本的测序结果，所述测序结果包括多个读段；提取所述测序结果中的割裂读段(soft-clippedreads)，所述割裂读段为同一读段被切割为两段、两段分别能匹配到参考序列两个不同位置的读段；分析匹配到所述参考序列上相同位置的割裂读段的数量，确定候选断点；定义所述参考序列上候选断点相应位置为第一融合基因位置，截取匹配到所述第一融合基因位置的割裂读段的不匹配所述第一融合基因位置的部分，以获得第一割裂片段，将所述第一割裂片段进行组装，获得第一一致性序列；将所述第一一致性序列与所述参考序列进行比对，定义所述第一一致性序列与所述参考序列匹配的位置为第二融合基因位置；截取匹配到所述第二融合基因位置的割裂读段的不匹配所述第二融合基因位置的部分，获得第二割裂片段，将所述第二割裂片段进行组装，获得第二一致性序列；将所述第二一致性序列与所述参考序列进行比对，若所述第二一致性序列与所述参考序列匹配的位置为所述第一融合基因位置，确定存在所述融合基因突变。所称的断点，指两个基因发生融合突变的位置。

根据本发明的一个实施例，所述测序结果包含的数据量不少于30x，即测序深度不小于30x；所述分析匹配到所述参考序列上相同位置的割裂读段的数量，确定候选断点，包括：确定匹配到所述参考序列上相同位置的割裂读段的数量不小于10条的位置为所述候选断点。

根据本发明的实施例提供的一种计算机可读存储介质，用于存储供计算机执行的第二程序，本领域普通技术人员可以理解，在执行该第二程序时，通过指令相关硬件可完成上述检测融合基因突变的方法的全部或部分步骤。所称存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

根据本发明的实施例提供的一种检测融合基因突变的装置，该装置包括：数据输入模块，用于输入数据；数据输出模块，用于输出数据；处理器，用于执行第二计算机可执行程序，所述第二计算机可执行程序的执行包括完成上述本发明一方面的检测融合基因突变的方法；存储模块，与所述数据输入模块、数据输出模块和处理器相连，用于存储数据，其中包括所述第二计算机可执行程序。

利用上述本发明任一实施例的检测融合基因突变的方法、计算机可读存储介质和/或装置，能够准确高效的检测融合基因突变。

根据本发明的实施例，发明人结合传统实验方法和高通量测序方法，提出了基于分析高通量测序数据的方法对重要突变位点，包括SNV、INDEL和融合突变(FUSION)进行快速精确检测的方法。该方法概括地说，包括使用寡核苷酸探针捕获技术或PCR多重扩增的方式来获取基因组上的目标序列，对目标序列产物进行高通量测序，从中识别DNA样品中的碱基序列及变异信息。根据示例，在对SNV、INDEL的检测中，针对试剂盒已知突变位点的特性，根据PCR方法的检测位点、COSMIC数据库中记载的致病变异，推算出拥有该变异时测序reads应当具有的序列，然后在测序数据中对这种序列进行检测。在这种情况下，对变异的检测不再关注其具体的比对位置和比对形式，而是关注在测序得到的reads中是否存在发生这种变化后应当具有的序列特征，从而规避了INDEL特别是复杂INDEL(complex INDEL)附近比对质量下降等的情况。利用上述提供的定点检出变异和检测基因融合突变的方法研究某一疾病样本，例如研究肺癌样本时，无需对照样本，基于对单个样本进行突变检测，就能一次性获得该样本的多个及多种有意义的突变信息，利于肺癌的个体化治疗。

以下结合附图和具体实施例对本发明的突变检测方法和/或装置进行详细的描述。下面示例，仅用于解释本发明，而不能理解为对本发明的限制。

需要说明的是，在本文中所使用的术语“第一”、“第二”等仅用于方便描述目的，而不能理解为指示或暗示相对重要性，也不能理解为之间有先后顺序关系。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

除另有交待，以下实施例中涉及的未特别交待的试剂、序列(接头、标签和引物)、软件及仪器，都是常规市售产品或者开源的，比如购自Life Technologies等。

实施例一

(一)参考值模型的构建

1，参考值模型构建的假设基础

1.1，对于任一位点，假设参考基因组对应的碱基为r∈{A,T,C,G}；

1.2，对于任一位点，假设覆盖该位点的所有reads的对应碱基为b_i，碱基质量值为q_i，则对应的碱基错误率为i＝1,2,...,d，d表示该位点对应的测序深度。

2，模型的建立

对于每一个位点的数据分布情况分为两种模型来解释：

模型M₀：这个位点不存在变异，跟参考基因组不同的那些碱基都是系统误差导致的；

模型这个位点的突变r→m是真实存在的，并且等位基因突变频率为f，对于那些既不为r，也不为m的碱基当做系统误差处理。

该位点的数据分布情况能够当作模型M₀来处理的概率为：

其中：

该位点的数据分布能够当做模型来处理的概率为

其中：

到此，变异检出的问题就转换为判断位点的数据分布情况更偏向于哪个模型，也即对两个概率L(M₀)和进行比较，于是建立如下的变异检出模型。

一般情况下，与L(M₀)的差异都是数量级上的差异，因此的值会很大，所以会对其采取取对数的操作。

其中，为参考值，θ即为对应的cut off值。

(二)模型用于变异检测的检测灵敏度

灵敏度(sensitivity)＝f(等位基因allelic fraction,测序深度depth,碱基测序错误率local sequencing error rate,确定的cut off)。

在此模型下，sensitivity表现为所有的使得成立的带有突变的reads数出现的概率的总和。

针对θ∈(0,0.1,0.2,...,10)里的每一个θ值，计算满足的最小k值，(k为带有突变的reads数)，即

然后通过二项分布的概率计算公式：

此处f(1-e)+(1-f)e为reads带有突变的概率。

计算相应的概率，即为灵敏性。

(三)模型用于变异检测的检测特异性

特异性(Specificity)＝f(local sequencing error rate,depth,evidence cutoff)。

在此模型下，Specificity表现为所有的使得成立的系统错误造成的带有突变的reads数出现的概率的总和

针对θ∈(0,0.1,0.2,...,10)里的每一个θ值，计算满足的最小k值，(k为不带突变的reads数)，即

然后通过二项分布的概率计算公式：此处1-e为reads不带有突变的概率。

计算相应的概率，即为特异性。

(四)ROC曲线(ROC curve)和cut off的确定

分别对测序深度和等位基因频率进行限定，根据上述公式进行计算，绘制ROC曲线图，结果如图3和图4所示。

从图3和图4中可以看出在大于200x，突变频率大于0.03的情况下选择cut offθ＝2可以满足要求。随着频率和测序深度的继续增大，AUC也在继续增大，故cut offθ＝2对测序深度大于200x，突变频率大于0.03的变异检出均成立。

实施例二

在获得测序下机数据后，以获得BGISEQ-100平台下机数据为例，进行变异检测一般包括如下部分：

1、变异已知信息处理以及测序数据预处理

1.1将要检出的变异类型转换成检出程序识别的格式，生成待测变异list列表。

1.2将下机数据与参考基因组比对。对BGISEQ-100有效测序数据使用tmap工具比对到参考基因组上，得到精确的比对结果。其中tmap工具源自：https://github.com/iontorrent/TS/tree/master/Analysis/TMAP。

排序。使用samtools sort对利用tmap比对后的结果(bam文件)进行排序：按照染色体编号和所在染色体上的位置按照从小到大的顺序进行排序。

去掉比对结果的PCR重复片段。对排序后的结果(bam文件)使用BamDuplicates工具去除PCR重复片段。其中，BamDuplicates工具的著作权源自Ion Torrent Systems,Inc.。

建索引。对去掉PCR重复片段之后的bam利用samtools index建立相应的索引。

1.3对处理好的bam文件进行QC质控，合格的文件将进行后续的步骤。

2、已知变异位点的检测

使用自主脚本lungSnvIndel.pl来检测去重之后的bam文件中那些已知的位点是否存在相应的变异，包括snv和indel。具体的，调用该脚本包括实现以下：

数据输入，从list文件中提取待测位点信息，读入bam文件中指定位点的测序reads。

过滤，依次过滤掉具有以下特征的reads：MapQ<30(比对质量小于30)、CIGAR中包含MIDS外的标记、指定位点位于reads末端5bp。

变异验证，以指定位点为中心，5`和3`方向各延伸5个bp，共11个bp与list中的记 录做比对，若相同则为支持reads。对所有覆盖reads循环操作，统计支持reads数量，代入以 下参考值模型进行验证判断，若成 立，则判定该变异存在。

结果输出。

3、未知INDEL类型的检测

如果上面步骤2中没有检测到INDEL的存在，则可以使用tvc来对上述已知INDEL的区域进行检测，给出其他的INDEL基因型。Tvc是lifeTechnologies公司针对proton数据开发的snv或indel检测程序。

使用TVC工具(http://ioncommunity.lifetechnologies.com/community/products/torrent-variant-caller)，使用其默认参数json文件检测肿瘤相关的SNV。

使用TVC工具，调整其参数json文件，参数data_quality_stringency由8.5改为6，参数filter_unusual_predictions由0.25改为0.3，检测肿瘤组织indel。

4、对融合基因进行检测

使用程序seekSV来对指定的融合基因进行检测，seekSV为华大基因(BGI)自主开发的融合基因检测软件。运行该程序包括实现以下：

4.1遍历bam文件，提取出有soft clip标记的测序reads，确定参考基因组每个位置覆盖的soft clip reads(割裂读段)数量。

4.2对每个位置覆盖的soft clip reads数量进行分析，确定breakpoint(断点)候选位点。

4.3假定breakpoint候选位点处reads比对的基因是融合基因中一对基因中的一个融合partner 1，切掉breakpoint候选位点的切割读段中非匹配到该候选位点的reads片段进行组装，形成一致性序列。

4.4将新生成的一致性序列在参考基因组上进行比对，若能比对上，假定新比对位置是breakpoint候选位点的另一个融合partner 2，对原比对到此处的soft clip reads重复4.3中的比对操作，若能比对到partner 1，则确认存在基因融合变异。

4.5结果输出。

5、可选择进行部分

变异注释。使用自主脚本Annotation.pl为检出的变异添加注释信息，同时添加患者信息。

生成报告。读入患者信息，利用已知的注释信息生成html版报告，同时自动化生成pdf版本的报告，并绘制检出位点的reads图。

图3显示基于BGISEQ-100测序平台的单样本测序数据确定变异检出的流程。

需要说明的：

在第1部分中的1.1，首先收集了不同试剂盒所检测的变异类型，同时收录了COSMIC数据库中记载的临床意义已明的变异类型，然后通过特定的算法，推定这些变异发生后测序reads应当具有的序列。举个例子，对于变异EGFR c.2235_2249del15，根据变异的描述，结合参考基因组hg19，可知正常测序reads在2235-2249间的碱基序列，当变异发生后，2235-2249间的碱基缺失，5`和3`端的碱基序列将直接连在一起，即生成新的目标序列，也即是发生变异后应当出现的序列，如下表1所示。所称的特定算法模拟以上变化过程，并向5`和3`两个方向各延伸5bp。

表1

表1中的“发生变异后应当出现的序列”、“序列起始位置”和“序列结束位置”即为发生变异后测序reads应当具有的序列特征。

在上述第1部分中的1.2中，针对BGISEQ-100测序数据的特点，对测序结果进行标准的处理，包括比对，去重等。

在上述第1部分中的1.3，对标准处理后数据进行QC质控，只有合格的数据才能参与到后续的变异分析中。

在上述第2-4部分中，通过第1部分中建立的列表，在测序数据中检测列表中的变异类型是否存在。若未检测到已知的INDEL基因型，则利用TVC在这个区域内寻找是否有新的INDEL突变基因型存在。若存在新的基因型且被判定为致病，则这个新的变异基因型也会加入到已知变异检测列表中。同时融合基因也会在这一步中被检测。

在上述第5部分中，对前边检测出来的变异进行注释并格式化输出。

采用上述的已知致病变异检出方法，有效的整合了实验方法和高通量测序方法的优点。相比实验方法，本发明经济成本更底，对检测的内容更加灵活，能够方便的增加新的检测基因型；相比传统的高通量方法，本方法只关注特定位点的确切变异，检测速度更快，灵敏度更高。由于对SNV和INDEL采用了新的检测策略，有效的解决了在INDEL区域比对质量下降对变异检出的影响，同时能够在相同的比对质量下更好的检测到complex INDEL变异，并针对BGISEQ-100的数据做了专门的优化。由于所要检测的变异已知，因此只需要分析单样本。此外，流程中还整合了融合基因检测。在流程检出后，根据收录的用药信息给予注释，直接生成pdf版本的临床检测报告。因此本发明的有益效果为在更低的经济投入，针对临床意义已明的致癌变异，利用BGISEQ-100高通量测序方法快速检测已知的致病变异(包含complex和fusion)，并直接给出pdf版本临床报告。

实施例三

该示例利用一名女性左上肺腺癌患者的FFPE组织样本，对其进行目标区域捕获以及BGISEQ-100平台测序，对测序下机的有效数据通过tmap比对、samtools sort排序、BamDuplicates去重、samtools index建索引、已知位点的变异检测、变异注释、生成报告等步骤，最终得到该患者的已知位点的变异检测报告。

将上述的变异检测方法的各部分流程都整合到软件Otype中，软件的运行环境为Linux操作系统，具体操作步骤如下：

在Linux操作系统计算机终端中输入如下的命令行：

perlOtype.pl–lsample.list–o outdir–O run.sh，会生成相应的运行脚本。

sh run.sh运行脚本。

Otype的命令行参数含义具体见表2的参数说明。

表2参数说明

结果分析：

1、分析QC质控信息，判断数据能否用于变异检出。数据统计信息如表3所示，数据质量满足后续分析要求。

表3

2、统计病人在每个已知位点的变异情况

如图6所示，如果最后一列为“KEEP”表示在特定的位点存在相应的变异，如果最后一列为“REJECT”表示在特定的位点不存在相应的变异。

如图6中的第32行，当EGFR L858R变异发生时，该样本测序结果比对到参考基因组上后，在chr7:55259510-55259521间应当具有TGGGCGGGCCA的序列，对覆盖该区域的258条reads进行过滤，剩余193条，检索目标序列，得到31条匹配reads，其中14条为正向链。提取质量值进行LOD计算，lod>2，通过检验，确认该变异存在。

上图中30行，当EGFR G719S变异发生时，该样本测序结果比对到参考基因组上后，在chr7:55241702-55241713间应当具有TCCTGAGCTCC的序列，对覆盖该区域的262条reads进行过滤，剩余218条，检索目标序列，得到1条正向链匹配reads，。提取质量值进行LOD计算，lod<2，不通过检验，确认该变异不存在。

3、自动化生成pdf的报告。包括已知位点的变异情况以及相应的一些靶向药物的信息，具体如以下图7所示。结合reads图，例如结合图8对检出变异进行确认。图8示例比对结果中，参考基因组chr7上的T突变为G碱基，与变异检测结果一致。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (9)

1.一种定点检出变异的方法，其特征在于，包括：

基于所述变异的已知信息，确定所述变异的指定位点和包含所述变异的参考序列；

获取待测样本的核酸的测序数据，所述测序数据包括多个读段；

提取所述测序数据中包含所述指定位点的读段，获得指定读段；

以所述指定读段上的指定位点为中心，往两端方向各延伸N个bp，获得指定片段，4≤N≤10；

将所述指定片段与所述包含所述变异的参考序列进行比对，获得支持读段，所述支持读段为与所述参考序列匹配的指定片段所在的读段；

统计所述支持读段的量，基于所述支持读段的量判断所述变异是否存在；

所述基于支持读段的量判断所述变异是否存在，包括：

将所述支持读段的量代入到以下模型，获得参考值LOD(m,f)，

所述模型基于以下两点假设而建立：

(1)对于任一位点，假设参考基因组对应的碱基为r∈{A,T,C,G}，

(2)对于任一位点，假设覆盖该位点的所有读段的对应碱基为b_i，碱基质量值为q_i，则对应的碱基错误率为d表示该位点对应的测序深度，其中

M₀表示所述位点的变异不存在，该位点与参考基因组碱基不同是由于系统误差导致的，

表示所述位点的变异真实存在，所述变异为r突变为m，f为等位基因突变频率，既不为r、也不为m的碱基是系统误差造成的，

L(M₀)表示所述位点的测序数据的分布情况符合M₀的概率，将所述支持读段的量代入式I以确定L(M₀)，

表示该所述位点的测序数据的分布情况符合的概率，将所述支持读段的量代入式III以确定

将所述参考值与阈值比较，当所述参考值大于所述阈值，判定所述待测样本存在所述变异。

2.权利要求1所述的方法，其特征在于，所述提取测序数据中包含所述指定位点的读段，获得指定读段，包括：

将所述测序数据与包含所述变异的参考序列进行比对，比对到所述参考序列上对应的指定位点位置的读段为所述指定读段。

3.权利要求2所述的方法，其特征在于，获得指定读段之后，对所述指定读段进行过滤，其中包括过滤掉指定位点位于读段的末端N bp内的指定读段。

4.权利要求1所述的方法，其特征在于，N＝5。

5.权利要求1-4任一所述的方法，其特征在于，所述变异选自SNP和INDEL中的至少一种。

6.一种定点检出变异的装置，其特征在于，包括：

数据输入单元，用于输入数据；

数据输出单元，用于输出数据；

处理器，用于执行第一计算机可执行程序，所述第一计算机可执行程序的执行包括完成权利要求1-5任一所述的方法；

存储单元，与所述数据输入单元、数据输出单元和处理器相连，用于存储数据，其中包括所述第一计算机可执行程序。

7.一种检测融合基因突变的方法，其特征在于，包括：

获取待测样本的测序结果，所述测序结果包括多个读段；

提取所述测序结果中的割裂读段，所述割裂读段为同一读段的两部分分别匹配到参考序列两个不同位置的读段；

分析匹配到所述参考序列上相同位置的割裂读段的数量，确定候选断点；

截取匹配到第一融合基因位置的割裂读段的不匹配所述第一融合基因位置的部分，以获得第一割裂片段，将所述第一割裂片段进行组装，获得第一一致性序列，定义所述参考序列上候选断点相应位置为第一融合基因位置；

将所述第一一致性序列与所述参考序列进行比对，定义所述第一一致性序列与所述参考序列匹配的位置为第二融合基因位置；

截取匹配到所述第二融合基因位置的割裂读段的不匹配所述第二融合基因位置的部分，获得第二割裂片段，将所述第二割裂片段进行组装，获得第二一致性序列；

将所述第二一致性序列与所述参考序列进行比对，若所述第二一致性序列与所述参考序列匹配的位置为所述第一融合基因位置，确定存在所述融合基因突变。

8.权利要求7所述的方法，其特征在于，所述测序结果包含的数据量不少于30x，所述分析匹配到所述参考序列上相同位置的割裂读段的数量，确定候选断点，包括：

确定匹配到所述参考序列上相同位置的割裂读段的数量不小于10条的位置为所述候选断点。

9.一种检测融合基因突变的装置，其特征在于，包括：

数据输入模块，用于输入数据；

数据输出模块，用于输出数据；

处理器，用于执行第二计算机可执行程序，所述第二计算机可执行程序的执行包括完成权利要求7或8所述的方法；

存储模块，与所述数据输入模块、数据输出模块和处理器相连，用于存储数据，其中包括所述第二计算机可执行程序。