CN108427864B - 一种拷贝数变异的检测方法、装置以及计算机可读介质 - Google Patents

一种拷贝数变异的检测方法、装置以及计算机可读介质 Download PDF

Info

Publication number
CN108427864B
CN108427864B CN201810151291.2A CN201810151291A CN108427864B CN 108427864 B CN108427864 B CN 108427864B CN 201810151291 A CN201810151291 A CN 201810151291A CN 108427864 B CN108427864 B CN 108427864B
Authority
CN
China
Prior art keywords
sample
copy number
tested
average
sequencing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810151291.2A
Other languages
English (en)
Other versions
CN108427864A (zh
Inventor
邵阳
汪笑男
吴雪
常志力
刘思思
那成龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Shihe gene Biotechnology Co., Ltd
Original Assignee
南京世和基因生物技术有限公司
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 南京世和基因生物技术有限公司 filed Critical 南京世和基因生物技术有限公司
Priority to CN201810151291.2A priority Critical patent/CN108427864B/zh
Priority to PCT/CN2018/090086 priority patent/WO2019157791A1/zh
Publication of CN108427864A publication Critical patent/CN108427864A/zh
Application granted granted Critical
Publication of CN108427864B publication Critical patent/CN108427864B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B45/00—ICT specially adapted for bioinformatics-related data visualisation, e.g. displaying of maps or networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/10—Ploidy or copy number detection
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations

Abstract

本发明涉及一种基因片段拷贝数变异的检测方法、装置以及计算机可读介质，特别是涉及了一种基于高通量测序（NGS）数据的检测基因片段拷贝数变异的技术，属于生物信息学技术领域。将本发明提供的基因片段拷贝数变异的检测方法应用于多个检测区域构成的一组基因（panel）的NGS测序过程中时，通过采用T分布拟合获得平均测序深度，进而计算出拷贝数基线及分析出存在异常的基因片段。该分析方法有效地消除各个检测区域因NGS试验方法导致的测序深度的波动性对检测结果的影响，避免拷贝数变异引起的其它区域的测序深度增大/减小而引起的拷贝数值波动性大、及在不同panel中检测结果不一致的问题。

Description

一种拷贝数变异的检测方法、装置以及计算机可读介质
技术领域
本发明涉及一种拷贝数变异的检测方法、装置以及计算机可读介质，特别是涉及了一种基于高通量测序数据的检测拷贝数变异的技术，属于生物信息学技术领域。
背景技术
人类基因组上存在着大量的变异。根据发生变异的碱基数目，基因组上的遗传变异分为单核苷酸变异(Single Nucleotide Variants，SNVs)和结构变异(StructuralVariations，SVs)。拷贝数变异(Copy Number Variations，CNVs)是结构变异的一种形式，是指与参照基因组相比，大小在50bp到几Mb的DNA片段的缺失、插入、复制和复杂多位点变异。近年来的研究表明，基因组片段的CNVs通过改变基因剂量或染色体构象影响基因的表达，导致生物体病变和影响疾病的发展，在表型多态性和进化研究中占据越来越重要的地位。目前在全基因组范围内寻找CNVs主要基于的两种技术，分别是基因芯片技术(DNAchip)及新一代测序技术(Next Generation Sequencing，NGS)。
基于NGS的检测方法是近年来快速发展的CNVs检测新方法，由基因组上某段区域的测序片段数目表征该段区域基因的含量，从而确定各基因含量异常的区域。这种方法的优势在于，在测序深度足够的条件下利用NGS测序数据能够获得比芯片检测更准确的断点位置和检测分辨率，实现在全基因组范围内更细致的CNVs检测。下一代测序技术又可以分为基于读深(read-depth based)、基于读对(paired-read based)、基于序列组装(sequence assembly based)和基于分裂读取(split-read based)四种方法；从所用的核心检测技术可分为基于概率统计模型和基于机器学习方法两类。统计类方法主要是根据读深信号的统计特征来检测拷贝数变异区域。统计类方法的一个假设前提是认为测序过程是均匀的，即沿着染色体的窗口的读深信号服从某种分布，比如泊松分布，而且读深信号与拷贝数数目之间是呈线性关系的。所以连续窗口的读深信号的增加或减少就预示着拷贝数的增加或减少，也即预示着拷贝数变异区域。
现有的进行高通量测序过程中，通常会将多个基因或多个敏感位点区域构成panel进行检测，但是当一些基因上发生了拷贝数变异之后，其它的一些区域上的高通量测序的测序深度会受到影响而出现测序深度的增大或减小的现象，会导致在计算panel上的平均测序深度时的不稳定性，进而导致计算结果在不同的panel上的结果不一致的问题。
发明内容
本发明要解决的是对于由多个检测区域构成的panel检测过程中特定基因的拷贝数计算结果不稳定的问题。
本发明的第一个方面，提供了用于对一个或多个检测区域/位点/基因的拷贝数进行计算的方法，该方法是将区域/位点/基因作为一个检测区域进行检测，检测区域可以与其它的位点/区域共同构成panel后通过NGS方法进行同时测序分析，技术方案是：
一种拷贝数变异的检测方法，包括如下步骤：
S1，将待检测区域和其它区域组成panel，对待测样本、对照样本和正常样本采用所述的panel进行NGS测序；
S2，根据步骤S1的结果，计算出每份正常样本和待测样本上的待检测区域的拷贝数；
S3，将正常样本中所述的待检测区域的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测区域的平均拷贝数基线；
S4，将待测样本中所述的待检测区域的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测区域的拷贝数的计算方法是：根据样本与对照样本的NGS测序结果，比对至参考基因组序列，获得唯一比对至所述的检测区域的读段(reads)，并计算出每个样本中每个检测区域上的测序深度；对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；根据平均测序深度计算出每份样本上的待检测区域的拷贝数。
在一个实施例中，panel中检测区域的数量为3～50000个。
在一个实施例中，拷贝数的计算是公式是：拷贝数＝(样本的待检测区域的测序深度/样本的平均测序深度)/(对照样本的待检测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
在一个实施例中，所述的基线数值是由平均值(mean)±标准偏差(sd)所构成。
在一个实施例中，所述的待检测区域包括但不限于以下基因中的任意一个片断：AKT1、AKT2、AKT3、CD274(PD-L1)、DDR2、EGFR、ERBB2、FGFR1、FGFR3、FLT4、HGF、MET、MTOR、MYC、PDCD1LG2(PD-L2)、PDGFRA、PDGFRB、SOX2、TP53、VEGFA、BRAF、FLT1、HRAS、KDR、KRAS、MAP2K1、MAP2K2、NRAS、PIK3CA、RB1、TOP1、VEGFA、BCL2、BCL6、CCND1、CCND3、CDK6、CEBPA、CEBPB、CEBPD、HOXA10、IL3、IRF4、KMT2A、LYL1、MUC1、MYC、NOTCH1、SETBP1、TAL1、ZAP70。
所述的拷贝数变异的检测方法是用于非治疗与诊断目的。
本发明的第二个方面，提供了用于检测基因的拷贝数变异的方法，这种方法主要是应用于含有多个待检测区域(例如外显子区域)的基因的拷贝数的检测，本方法是将该基因中的多个待检测区域分别进行拷贝数的计算，并将这几个待检测区域拷贝数的平均值作为该基因的拷贝数，并进行比较分析，技术方案是：
一种拷贝数变异的检测方法，包括如下步骤：
S1，从待测基因序列中选出至少两个待检测区域，并与其它区域组成panel，对待测样本、对照样本和正常样本采用所述的panel进行NGS测序；
S2，根据步骤S1的结果，计算出每份待测样本和正常样本上的待检测基因的拷贝数；
S3，将正常样本中待检测基因的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测基因的平均拷贝数基线；
S4，将待测样本中所述的待检测基因的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测基因的拷贝数是指每个待检测区域的拷贝数的平均值；
待检测区域的拷贝数的计算方法是：根据样本与对照样本的NGS测序结果，比对至参考基因组序列，获得唯一比对至所述的检测区域的读段(reads)，并计算出每个样本中每个检测区域上的测序深度；对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；根据平均测序深度计算出每份样本上的待检测区域的拷贝数。
在一个实施例中，panel中检测区域的数量为3～50000个。
在一个实施例中，拷贝数的计算是公式是：拷贝数＝(样本的待检测区域的测序深度/样本的平均测序深度)/(对照样本的待检测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
在一个实施例中，每个待检测区域的拷贝数的平均值是指算术平均值、正态分布平均值卡方分布平均值、F分布平均值或者T分布平均值。
待检测基因包括但不限于：AKT1、AKT2、AKT3、CD274(PD-L1)、DDR2、EGFR、ERBB2、FGFR1、FGFR3、FLT4、HGF、MET、MTOR、MYC、PDCD1LG2(PD-L2)、PDGFRA、PDGFRB、SOX2、TP53、VEGFA、BRAF、FLT1、HRAS、KDR、KRAS、MAP2K1、MAP2K2、NRAS、PIK3CA、RB1、TOP1、VEGFA、BCL2、BCL6、CCND1、CCND3、CDK6、CEBPA、CEBPB、CEBPD、HOXA10、IL3、IRF4、KMT2A、LYL1、MUC1、MYC、NOTCH1、SETBP1、TAL1、ZAP70。
本发明的第三个方面，提供了一种用于对一个或多个检测区域/位点/基因的拷贝数进行计算的方法的检测装置，该检测装置是将区域/位点/基因作为一个检测区域进行检测，检测区域可以与其它的位点/区域共同构成panel后通过NGS方法进行同时测序分析，技术方案是：
一种拷贝数变异的检测装置，包括：
测序数据获取模块，用于对待测样本、对照样本和正常样本采用由待检测区域和其它区域组成的panel进行NGS测序，并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段(reads)；
测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；
拷贝数计算模块，用于计算出每份正常样本和待测样本上的待检测区域的拷贝数；
基线计算模块，用于将正常样本中待检测区域的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测区域的平均拷贝数基线；
分析模块，用于将待测样本中所述的待检测区域的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异。
在一个实施例中，拷贝数是通过如下公式计算得到：拷贝数＝(样本的待检测区域的测序深度/样本的平均测序深度)/(对照样本的待检测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
在一个实施例中，panel中检测区域的数量为3～50000个。
本发明的第四个方面，提供了一种用于基因的拷贝数变异的检测装置，该检测装置是应用于含有多个待检测区域(例如外显子区域)的基因的拷贝数的检测，技术方案是：
一种拷贝数变异的检测装置，包括：
测序数据获取模块，用于对待测样本、对照样本和正常样本采用panel进行NGS测序，并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段(reads)；所述的panel是由至少两个待检测区域和其它区域组成，所述的待检测区域选自待测基因序列；
测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；
拷贝数计算模块，用于计算出每份正常样本和待测样本上的待检测基因的拷贝数；
基线计算模块，用于将正常样本中所述的待检测基因的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测基因的平均拷贝数基线；
分析模块，用于将待测样本中所述的待检测基因的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测基因的拷贝数是指每个待检测区域的拷贝数的平均值。
在一个实施例中，拷贝数是通过如下公式计算得到：拷贝数＝(样本的待检测区域的测序深度/样本的平均测序深度)/(对照样本的待检测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
在一个实施例中，panel中检测区域的数量为3～50000个。
本发明的第五个方面，提供了：
一种计算机可读介质，记载有可以运行上述拷贝数变异的检测方法的程序。
有益效果
本发明提供的拷贝数变异的检测方法应用于多个检测区域构成的panel的NGS测序过程中时，可以有效地避免拷贝数变异引起的其它区域的测序深度增大/减小而引起的拷贝数值波动性大、在不同panel中检测结果不一致的问题；通过采用了T分布拟合获得平均测序深度可以有效地消除各个检测区域测序深度的波动性对计算结果的影响，通过计算出拷贝数基线可以有效地分析出存在异常的检测区域。
附图说明
图1是本发明提供的检测方法的流程图。
图2是本发明提供的检测装置的结构图。
图3是本发明提供的检测方法的MYCN基因检测结果对比图。
图4是本发明提供的检测方法的MET基因检测结果对比图。
图5是本发明提供的检测方法的CDKN2A基因检测结果对比图。
具体实施方式
下面通过具体实施方式对本发明作进一步详细说明。但本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。
本文使用的词语“包括”、“包含”、“具有”或其任何其他变体意欲涵盖非排它性的包括。例如，包括列出要素的工艺、方法、物品或设备不必受限于那些要素，而是可以包括其他没有明确列出或属于这种工艺、方法、物品或设备固有的要素。除非上下文明确规定，否则单数形式“一个/种”和“所述(该)”包括复数个讨论对象。
本发明所提出的拷贝数变异检测方法，是基于读深方法进行判定各个检测区域的拷贝数变异，进行检测的过程中，各个检测的区域共同构成panel进行高通量测序，这些区域可以是外显子，也可以是其它的关注的基因。检测的原理为：如果染色体某一区域发生了拷贝数变异，则高通量测序时该区域的序列片段分布将发生变化，即，拷贝数缺失-序列密度将变小，拷贝数扩增-序列密度将变大。本发明提供的检测方法仅仅是用于用于通过测序结果判定是否存在着拷贝数变异现象，并非是用于非治疗或诊断目的。
在实际检测过程中，为了发挥高通量测序的通量大、效率高的优点，通常是将许多个基因片段、区域共同组成检测panel，同时对这些区域进行检测，在一些通常的panel中，检测区域的数量可以是几十个，也可以是几百、上千个。有可能有一次检测的过程中，较为关心的只是这个panel中的一部分基因或位点，这部分也就是属于待检测的基因或位点区域。本发明中所述的“待检测区域”是指进行拷贝数变异检测时需要检测的区域。本发明中所述的“其它区域”是指在中与待检测区域一起加入至panel的那些区域，这些区域可能是在其它的条件下较为重要、需要检测，也可以是在本次检测中的后续步骤中再进行拷贝数变异的检测，“其它区域”的加入是为了提高高通量测序的效率而进行加入。根据上述的定义，本发明说明书中当提及“每个检测区域”、“panel中的检测区域”时，是指“待检测区域”和“其它区域”，因为这两个区域都属于在一次高通量测序中都进行检测的区域，并且“其它区域”中的测序深度等数据也是要应用于检测的计算过程。
在本文中所使用的术语“高通量测序”、“下一代测序”等，指的是第二代高通量测序技术及之后发展的更高通量的测序方法。下一代测序平台包括但不限于Illumina(Miseq、Hiseq2000、Hiseq2500、Hiseq3000、Hiseq4000、HiseqX Ten等)、ABI-Solid和Roche-454测序平台等。随着测序技术的不断发展，本领域技术人员能够理解的是还可以采用其他方法的测序方法和装置进行本检测。根据本发明的具体示例，可以将根据本发明实施例的核酸标签用于Illumina、ABI-Solid和Roche-454测序平台等的至少一种进行测序。下一代测序技术，例如Illumina测序技术具有以下优势：(1)高灵敏度：下一代测序，例如Miseq的测序通量大，目前一次实验流程可以产生最多15G碱基数据，高的数据通量可以在测序序列数一定的情况下，使得每条序列获得更高的测序深度，所以可以检测到含量更低的突变，同时因其测序深度高，突变位点被多次覆盖，其测序结果也更为可靠。(2)高通量，低成本：利用根据本发明实施例的标签序列，通过一次测序可以检测上万份样本，从而大大降低了成本。
本发明中“待测样本”是指需要进行检测，并判定该样本上的一个或者多个基因区域是否存在有拷贝数变异。“正常样本”是指用于计算判定基线的样本，这些样本可以是选择具有一定数量的正常人的血液样本进行统计分析(可以采用几百例正常人样本，可以根据实际情况使样本数量具有统计意义)。在计算拷贝数时，还需要使用“对照样本(control)”，以对待测样本/正常样本的测序深度进行对比，并计算拷贝数，这里的“对照样本”可以是使用与待测样本配对的样本，当待测样本存在着自己的对照时，可以直接将其作为对照样本，如果没有对照样本时，可以采用标准品细胞DNA样本作为对照，例如当待测样本为人类来源时，可使用健康北京人B淋巴细胞DNA样本NA18535等。
如本文所用的，术语“比对”是指将读数或标签与参考序列进行比较并且由此确定该参考序列是否含有该读数序列的过程。如果该参考序列含有该读数，则可以将该读数映射至参考序列，或者在某些实施方案中，映射至参考序列中的某个特定位置。在有些情况下，比对简单地告知读数是或不是特定参考序列的成员(即，该读数是存在还是不存在于参考序列中)。举例来说，将读数与人类染色体13的参考序列进行比对，将告知该读数是否存在于针对染色体13的参考序列中。提供这种信息的工具可以被判定集合成员身份测试器。在有些情况下，比对另外指示参考序列中读数或标签可以映射至其中的位置。举例来说，如果参考序列是全人类基因组序列，那么比对可以指示读数存在于染色体13上，并且可以进一步指示读数是在染色体13的具体股和/或位点上。术语“参考基因组”或“参考序列”是指任何生物体或病毒的任何具体的已知基因组序列(无论是部分的或完整的)，它可以用于对来自受试者的识别的序列进行参比。例如，用于人类受试者以及许多其他生物体的参考基因组可见于美国国家生物技术信息中心(ncbi.nlm.nih.gov)，对于人的样品来说，参照序列可以是人基因组hg18或hg19的序列。目前hg19的相关数据库相对较多且hg19测出来的碱基量比hg18要多，即样品比对率会相对较高，所以优先选择hg19。
术语“读数(read)”是指来自核酸样品的一部分的序列读数。典型地，虽然并不一定是，读数代表样品中的相邻碱基对的短序列。读数可以通过样品部分的碱基对序列(以ATCG表示)以符号代表。它可以存储在存储设备中，且酌情处理，以确定它是否与参考序列匹配或者满足其它标准。读数可以直接地从测序装置中或者间接地从涉及样品的存储的序列信息中获得。在有些情况下，读数是足够长度(例如，至少约30bp)的DNA序列，其可以用于识别更大的序列或区域，例如，其可以与染色体或基因组区域或基因进行比对并且具体地分配给染色体或基因组区域或基因。
每个目标区域上的序列信息是比对结果中包含该位点的测序片段，位点的测序深度信息是比对结果中包含该位点的测序片段数目。
本发明提供的方法进行拷贝数变异的检测可以是以下的一些基因：
ABL1、APC、ARID2、AURKA、BCL2、BLM、BTK、CCND2、CDC73、CDK8、CEBPA、CRKL、CTNNB1、EGFR、EPHB1、ESR1、FANCC、FANCL、FGF23、FGFR2、FLT4、GID4、GPR124、IDH1、IL7R、JAK2、KDM5C、KLHL6、MAP2K4、MED12、MLH1、MSH2、MYCL1、NFE2L2、NPM1、NUP93、PDGFRA、PIK3R1、PRKDC、RAD51、RICTOR、SF3B1、SMO、SPOP、SUFU、TOP1、VHL、ZNF703、AKT1、AR、ASXL1、AURKB、BCL2L2、BRAF、CARD11、CCND3、CDH1、CDKN1B、CHEK1、CRLF2、DAXX、EMSY、ERBB2、EZH2、FANCD2、FBXW7、FGF3、FGFR3、FOXL2、GNA11、GRIN2A、IDH2、INHBA、JAK3、KDM6A、KRAS、MAP3K1、MEF2B、MLL、MSH6、MYCN、NFKBIA、NRAS、PAK3、PDGFRB、PIK3R2、PTCH1、RAF1、RNF43、SMAD2、SOCS1、SRC、TET2、TP53、WISP3、BRCA1、AKT2、ARAF、ATM、AXL、BCL6、CSF1R、CBFB、CCNE1、CDK12、CDKN2A、CHEK2、FGF10、DDR2、EP300、ERBB3、FAM123B、FANCE、IGF1R、FGF4、FGFR4、GATA1、GNA13、GSK3B、MEN1、IRF4、JUN、KDR、LRP1B、MCL1、PALB2、MLL2、MTOR、MYD88、NKX2-1、NTRK1、SMAD4、PDK1、PPP2R1A、PTEN、RARA、RPTOR、BRCA2、SOX10、STAG2、TGFBR2、TSC1、WT1、CTCF、AKT3、ARFRP1、ATR、BAP1、BCOR、FGF14、CBL、CD79A、CDK4、CDKN2B、CIC、IKBKE、DNMT3A、EPHA3、ERBB4、FAM46C、FANCF、MET、FGF6、FLT1、GATA2、GNAQ、HGF、PAX5、IRS2、KAT6A、KEAP1、MAP2K1、MDM2、SMARCA4、MPL、MUTYH、NF1、NOTCH1、NTRK2、BRIP1、PIK3CA、PRDM1、PTPN11、RB1、RUNX1、CTNNA1、SOX2、STAT4、TNFAIP3、TSC2、XPO1、FGF19、ALK、ARID1A、ATRX、BARD1、BCORL1、IKZF1、CCND1、CD79B、CDK6、CDKN2C、CREBBP、MITF、DOT1L、EPHA5、ERG、FANCA、FANCG、PBRM1、FGFR1、FLT3、GATA3、GNAS、HRAS、SMARCB1、JAK1、KDM5A、KIT、MAP2K2、MDM4、MRE11A、MYC、NF2、NOTCH2、NTRK3、PIK3CG、PRKAR1A、RAD50、RET、SETD2、SPEN、STK11、TNFRSF14、TSHR和ZNF217。
也可以是以下一些癌症相关的敏感基因，例如：
肺癌相关的热点基因AKT1、AKT2、AKT3、CD274(PD-L1)、DDR2、EGFR、ERBB2、FGFR1、FGFR3、FLT4、HGF、MET、MTOR、MYC、PDCD1LG2(PD-L2)、PDGFRA、PDGFRB、SOX2、TP53、VEGFA。
肠癌相关的热点基因AKT1、AKT2、AKT3、BRAF、EGFR、ERBB2、FLT1、FLT4、HRAS、KDR、KRAS、MAP2K1、MAP2K2、MET、MTOR、NRAS、PIK3CA、RB1、TOP1、VEGFA。
血癌相关的热点基因BCL2、BCL6、CCND1、CCND3、CDK6、CEBPA、CEBPB、CEBPD、FGFR3、HOXA10、IL3、IRF4、KMT2A、LYL1、MUC1、MYC、NOTCH1、SETBP1、TAL1、ZAP70。
所检测的区域可以使任何感兴趣的区域，不局限于外显子与基因，长度大小也都没有限制。上面的检测区域数量，可以是3～50000个。
例如：如果检测若干个感兴趣的区域时，可以将它们与其它数量不限的区域共同构成panel同时进行NGS测序，测序过程中这些感兴趣的区域都通过其参考基因序列设计探针进行捕获并检测；如果是对一个较大的基因作为感兴趣的区域进行检测时，也可以只检测这个基因中的若干个子区域(例如是外显子区域)，将这些子区域与其它数量不限的区域共同构成panel同时时行NGS测序，测序过程中根据这些子区域的参考序列设计探针对子区域进行捕获和检测，并将子区域的拷贝数的平均值作为该基因的拷贝数(这里的平均值可以是指算术平均值、正态分布平均值卡方分布平均值、F分布平均值或者T分布平均值)。
本发明提供的拷贝数变异检测方法流程如图1所示：
对于步骤101，需要对正常样本、待测样本、对照样本进行NGS测序，通过高通量测序的方法获得相应位点的序列信息可以按照常规的实验方法、教科书、探针设计方法、测序仪使用手册中的描述进行，主要的流程包括：对每个待测样本和正常样本的组织样本或者全血样本进行DNA提取，获取基因组DNA；对DNA片段过大的样本，通过超声破碎，将样本机械力打断至200-350碱基对；对片段化的DNA分子执行末端修复、添加腺嘌呤、文库接头连接等操作；获得的DNA文库与长度为120碱基的单链生物素标记DNA探针分子杂交，再以链霉亲和素包裹的磁珠分离捕获的DNA文库分子；在illumina下一代测序仪上进行测序。测序反应获得的数据通过生物信息学分析。在获得了相应的测序信息后，可以采用常规方法做数据进行预处理，这里的处理主要是对测序所得的每个样本序列分别进行过滤，以去除掉不合格的序列和接头序列，其中，样本包括目标样本(即，变异组织)和对照样本(即，正常组织)；具体地，对高通量测序后的样本序列进行过滤，去除不合格的序列及接头序列，其中，不合格序列可以为下列情况中的至少一种：测序质量低于某一阈值的碱基个数超过整条序列碱基个数的一定比例(例如，50％)和序列中测序结果不确定的碱基(例如，IlluminaGA测序结果中的N)个数超过整条序列碱基个数的一定比例(例如，10％)。其中，高通量测序技术可以为IlluminaGA或者HiSeq测序技术，也可以为现有的其他高通量测序技术，低质量阈值可以由具体测序技术和测序环境确定。在对读段进行了预处理之后，将过滤后的每个样本序列分别比对到参考基因组序列，对比对后的每个样本序列分别进行筛选以得到唯一比对的样本序列，确定每个唯一比对的样本序列相对于参考基因组序列的位置信息，并对位置信息进行排序；具体地：(1)首先可以通过任何一种短序列映射程序(例如，短寡核苷酸分析包(Short Oligo nucleotide Analysis Package,SOAP))将过滤得到的每个样本序列(即，由多个测序片段数据构成的序列)分别比对到参考基因组序列(例如，人类基因组参考序列)得到每个样本序列在参考基因组上的位置情况；(2)然后，对比对结果进行一系列的筛选，例如，去除比对到多个位置的序列(因为这个序列已无法准确唯一的提供比对位置信息)、去除重复出现的序列(因为这些序列可能是由于前期实验引入的误差，如由测序错误引起，为使检测结果更加精准，故去除)，以得到唯一比对的序列结果。
对于步骤102，需要根据以上获得的唯一比对的序列统计出在各个样本各个目标区域上的测序深度，也就是统计在各个目标区域上的读数。由于本发明中是对多个目标区域构成的panel进行整体检测，在一些位置上的拷贝数变异会导致该样本的总体平均测序深度发生变化，因此计算panel上各个区域的平均测序深度时就会出现会随着选取的目标区域的不同而发生计算值不一致的情况，也就导致了拷贝数计算结果的不稳定性，并且随着目标区域的数量的减少，影响程度呈增大趋势。或者随着拷贝数变异程度增大，影响程度也会增大。对于由若干个待检测区域和其它一些区域共同构成的panel来说，这里是要计算出分别在这些待测区域和其它区域上的测序深度；而对于包含一个基因中的多个子区域的panel来说，这里也是同样地获得这些子区域以及其它一些区域的测序深度。
对于步骤103，需要对同一份样本上的全部的检测区域的测序深度进行统计，目标是为了获得该样本上的平均测序深度，本发明中意外地发现针对步骤102中所描述的问题，采用t分布拟合计算平均值的方法可以非常有效地解决计算结果稳定性问题。对于由若干个待检测区域和其它一些区域共同构成的panel来说，将一份样本中这些全部的检测区域的测序深度进行T分布统计，并拟合计算出分布平均值，就该值作为这份样本的平均测序深度；同样地，对于包含一个基因中的多个子区域的panel来说，也是要将这些全部的子区域和其它的区域的测序深度进行T分布统计，获得平均测序深度；在该步骤中，可以针对每个正常样本、待测样本和对照样本分别统计出它们的平均测序深度。
对于步骤104，通过比值比(odds ratio,OR)的方法计算出每个样本上的每个待测区域的拷贝数。
对于由若干个待检测区域和其它一些区域共同构成的panel来说，该计算方法是通过以下公式得到的：
对于待测样本，拷贝数＝(待测样本的待测区域的测序深度/待测样本的平均测序深度)/(对照样本的待测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
对于正常人群样本，拷贝数＝(正常样本的待测区域的测序深度/正常样本的平均测序深度)/(对照样本的待测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
对于包含一个基因中的多个子区域的panel来说，该计算方法是也是类似地通过以下公式得到的：
对于待测样本，拷贝数＝(待测样本的待测子区域的测序深度/待测样本的平均测序深度)/(对照样本的待测子区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
对于正常人群样本，拷贝数＝(正常样本的待测子区域的测序深度/正常样本的平均测序深度)/(对照样本的待测子区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
使用使用比值比进行计算可以达到以下几个目的，首先可以降低区域的GC含量对其测序深度的影响，这种影响会在计算待测样本和对照样本在同一区域的比值比时，这种影响在除法中抵消了，其次，消除了该区域自身由于自身特性比如不容易捕获而导致的测序深度较其他区域低的影响，通常会可能会被认为是发生了拷贝数降低，但通过比值比的计算后只关注待测样本相对于对照样本的比值比的变化而不是实际的值，因此依然可以准确得到拷贝数变异区域。最后，如果采用对照样本时，可以对不同的对照样本实现比较，因为它们都是相对于同一个标准做的比较。
其中，如果是对包含多个子区域的基因的拷贝数进行检测时，对该基因上的子区域进行检测，然后根据每个子区域的拷贝数做平均值作为该基因上的拷贝数。
对于步骤105，其目的是为了获得用于判定目标区域是否存在拷贝数变异的基线，所使用的方法是将全部的正常样本上的待测个区域(或者是包含子区域的基因)的拷贝数进行T分布拟合，求出平均值，作为拷贝数基线，这里的基线数值可以是采用平均值(mean)±标准偏差(sd)所构成。
对于步骤106，将计算得到的待测样本上的某一个区域(或者是包含子区域的基因)的拷贝数值与该区域的基线值进行比较，大于或小于某一个数值时，则判定存在拷贝数增加/缺失。例如，可以判定当拷贝数<＝0.65或>＝1.6时，存在着缺失/增加，这里的mean±sd可以是Mean±2SD、Mean±2.5SD或Mean±3SD等，采取的平均值加减sd的标准使判定结果具有统计学意义，例如mean±2sd包含了96％的样本，3sd包含99％样本。如果超过这个范围，可以再通过假设检验计算p值，考察其显著性。
基于上述的方法，本发明还提供了用于检测拷贝数变异的装置，该装置中可以多个模块构成，如图2所示：
当用于检测一个待测区域时，装置包括：测序数据获取模块，用于对待测样本和正常样本采用包含多个检测区域的panel进行NGS测序，并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段(reads)；所述的panel是由待检测区域和其它区域共同构成；测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；拷贝数计算模块，用于计算出每份正常样本和待测样本上的待测区域的拷贝数；基线计算模块，用于将正常样本中待检测区域的拷贝数进行T分布拟合，将分布曲线的平均值作为该检测区域的平均拷贝数基线；分析模块，用于将待测样本中所述的待检测区域的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异。在一个实施例中，拷贝数是通过如下公式计算得到：拷贝数＝(样本的待测区域的测序深度/样本的平均测序深度)/(对照样本的待测区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
当用于检测一个包含若干子区域的基因时，装置包括：测序数据获取模块，用于对待测样本和正常样本采用包含多个检测区域的panel进行NGS测序，并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段(reads)；所述的多个检测区域是由待检测基因的待测子区域和其它区域构成；测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；拷贝数计算模块，用于计算出每份正常样本和待测样本上的待检测基因的拷贝数；
基线计算模块，用于将正常样本中所述的待检测基因的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测基因的平均拷贝数基线；分析模块，用于将待测样本中所述的待检测基因的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；其中，待测样本和正常样本的待检测基因的拷贝数是根据每个待测子区域的拷贝数的平均值。在一个实施例中，拷贝数是通过如下公式计算得到：拷贝数＝(样本的待测子区域的测序深度/样本的平均测序深度)/(对照样本的待测子区域的测序深度/对照样本的平均测序深度)×倍体数(Ploidy)。
需要说明的是，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。
显然，本领域的技术人员应该明白，上述的本发明的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明不限制于任何特定的硬件和软件结合。
实施例1 MYCN基因拷贝数的检测
本实施例中检测的目标区域是含有2个外显子区域的MYCN基因的拷贝数进行检测，目标区域的主要信息如下表所示：
表1 MYCN基因的2个外显子
区域名称 在染色体上的起始位置 在染色体上的结束位置
MYCN:NM_005378:e1 chr2：16082186 chr2：16082976
MYCN:NM_005378:e2 chr2：16085614 chr2：16086219
这里的正常人群共349位，他们的MYCN基因OR值组成的z分布的曲线平均值为1.08，sd为0.075，这里的病人样本用的是https://www.horizondiscovery.com公司提供的HD786标准品，它包含MYCN扩增，官方提供的拷贝数OR为4.75，我们利用NGS的方法以及以上描述的算法得到的该基因的拷贝数OR为3.95左右。这里采用的对照样本为NA18535标准品细胞。
进行检测的panel中还包括了其它的一些敏感基因的外显子区域，分别取不同数量的区域与MYCN基因2个外显子共同构成panel，共选择了4684个其它备选区域，例如其中100个其它备选区域的信息如下：
表2 panel中其它检测区域信息
序号 区域名称 序号 区域名称 序号 区域名称
1 MYCN:NM_005378:e1 35 BRAF:NM_004333:e17 68 CCND3:NM_001136017:e5
2 MYCN:NM_005378:e2 36 BRAF:NM_004333:e16 69 CCND3:NM_001136017:e2
3 AKT1:NM_005163:e14 37 BRAF:NM_004333:e3 70 CDK6:NM_001259:e7
4 AKT1:NM_005163:e8 38 FLT4:NM_182925:e27 71 CDK6:NM_001259:e5
5 AKT3:NM_181690:e13 39 FLT4:NM_182925:e17 72 CDK6:NM_001259:e3
6 AKT3:NM_181690:e6 40 FLT4:NM_182925:e14 73 CDK6:NM_001145306:e5
7 AKT3:NM_181690:e2 41 FLT4:NM_182925:e12 74 CEBPA:NM_001285829:e1
8 CEP170:NM_014812:e7 42 KRAS:NM_033360:e5 75 CEBPA:NM_004364:e1
9 BRF1:NM_001242787:e11 43 KRAS:NM_033360:e3 76 CEBPD:NM_005195:e1
10 BRF1:NM_001242786:e17 44 KRAS:NM_004985:e4 77 FGFR3:NM_000142:e2
11 CD274:NM_014143:e2 45 MAP2K1:NM_002755:e3 78 FGFR3:NM_000142:e4
12 CD274:NM_014143:e6 46 MAP2K1:NM_002755:e5 79 FGFR3:NM_000142:e9
13 CD274:NM_001267706:e2 47 MAP2K1:NM_002755:e8 80 FGFR3:NM_000142:e15
14 EGFR:NM_201284:e2 48 NRAS:NM_002524:e 81 FGFR3:NM_001163213:e7
15 EGFR:NM_201284:e7 49 NRAS:NM_002524:e3 82 IRF4:NM_001195286:e4
16 EGFR:NM_201283:e1 50 NRAS:NM_002524:e2 83 IRF4:NM_001195286:e7
17 EGFR:NM_201282:e14 51 PIK3CA:NM_006218:e4 84 IRF4:NM_002460:e4
18 FGFR1:NM_015850:e16 52 PIK3CA:NM_006218:e6 85 IRF4:NM_002460:e8
19 FGFR1:NM_015850:e13 53 PIK3CA:NM_006218:e15 86 LYL1:NM_005583:e4
20 FLT1:NM_001160030:e13 54 top1MT:NM_001258447:e11 87 LYL1:NM_005583:e3
21 FLT1:NM_001160030:e10 55 top1MT:NM_001258447:e8 88 LYL1:NM_005583:e2
22 FLT1:NM_001160030:e5 56 top1MT:NM_001258446:e7 89 MUC1:NM_001204291:e7
23 HGFAC:NM_001528:e2 57 MTOR:NM_004958:e55 90 MUC1:NM_001204291:e5
24 HGFAC:NM_001528:e11 58 MTOR:NM_004958:e52 91 MUC1:NM_001204291:e3
25 PDCD1LG2:NM_025239:e2 59 MTOR:NM_004958:e41 92 MUC1:NM_001204290:e7
26 PDCD1LG2:NM_025239:e3 60 MTOR:NM_004958:e24 93 NOTCH1:NM_017617:e33
27 PDCD1LG2:NM_025239:e5 61 BCL2L15:NM_001010922:e4 94 NOTCH1:NM_017617:e29
28 PDCD1LG2:NM_025239:e6 62 BCL2L15:NM_001010922:e3 95 NOTCH1:NM_017617:e26
29 PDCD1LG2:NM_025239:e7 63 BCL2L15:NM_001010922:e2 96 NOTCH1:NM_017617:e21
30 PDGFRA:NM_001347827:e3 64 CCND1:NM_053056:e1 97 NOTCH1:NM_017617:e13
31 PDGFRA:NM_001347827:e4 65 CCND1:NM_053056:e2 98 SETBP1:NM_015559:e2
32 PDGFRA:NM_001347827:e17 66 CCND1:NM_053056:e4 99 SETBP1:NM_015559:e4
33 TP53INP1:NM_033285:e4 67 CCND1:NM_053056:e5 100 TAL1:NM_001287347:e4
34 TP53INP1:NM_033285:e3
在进行计算时，将MYCN基因的2个外显子分别和10个、50个、100个、300个、400个、……4000个、4684个其它的热点基因的外显子区域共同构成panel同时进行NGS测序，具体的panel中的区域数量如表3所示，其中采用10个、50个、100个其它区域时，是采用的表2中序号为1～10、1～50、1～100的外显子区域，其它的区域在此不一一描述。
采用的检测步骤如下：
1、采用Illumina高通量测序技术对待测序列进行测序，接收测序序列后，对测序序列进行过滤，去除不合格的序列，将样本接头序列从序列片段中去除，其中，不合格序列包括：测序质量值低于5的碱基个数超过整条序列碱基个数的50％或序列中测序结果中N的个数超过整条序列碱基个数的10％；
2、采用短寡核苷酸分析包(SOAP)映射程序将高通量测序技术得到的癌症样本和癌旁样本测序片段比对到人类参考基因组序列上，筛选掉比对结果中多重比对的序列、去除重复出现的测序序列、一端测序数据中碱基序列一样的序列(仅保留一份)，以降低结果中的假阳性，最后，根据需求提取比对结果中所需要处理的染色体编号及位置信息，染色体位置。
3、根据待测样本、正常样本与对照样本的各个检测区域的读段(reads)计算出各个区域上的测序深度；对于每一个单独的样本，将样本上各个检测区域的测序深度进行T分布拟合(例如，当需要检测100个区域时，将某一份样本上这100个区域的测序深度进行T分布统计)，并计算出T分布曲线的平均值，将这个平均值作为该单独样本的平均测序深度；按照上述的方法，再对其它的样本进行平均测序深度的计算，此时，就获得了待测样本集、正常样本集与对照样本集中的每一个样本的平均测序深度。
4、通过比值比的方法计算出每一个样本上的每一个检测区域的拷贝数，计算方法是：
对于待测样本，某一个特定检测区域的拷贝数＝(待测样本上检测区域上的测序深度/待测样本上的平均测序深度)/(对照样本上某区域上的测序深度/对照样本上的平均测序深度)×倍体数(Ploidy)。
对于正常人群样本，拷贝数＝(正常样本上检测区域上的测序深度/正常样本上的平均测序深度)/(对照样本上某区域上的测序深度/对照样本上的平均测序深度)×倍体数(Ploidy)。
可以获得了每个待测样本和正常样本上的各个区域上的拷贝数。
本实施例中，MYCN基因具有2个外显子，对这两个外显子区域计算出拷贝数之后，将两者的算术平均值作为基因的拷贝数。
5、把正常样本集中每个样本上的MYCN基因的拷贝数进行T分布统计，并计算出分布曲线上的平均值，作为MYCN基因的拷贝数基线值，并以mean±3sd数值范围作为判定区间；
6、将待测样本的MYCN基因的拷贝数与步骤5得到的MYCN基因拷贝数基线进行比较，判定是否出现拷贝数异常。
作为对照，步骤3中计算每个样本的平均测序深度中，采用的是将一个样本中的全部检测区域的测序深度作算术平均值，作为该样本的平均测序深度。
计算结果如表4所示，区域数量与拷贝数值之间进行作图后如图3所示。
表4 MYCN基因拷贝数
从表中可以看出，采用本发明提供的T分布计算得到的MYCN基因的2个外显子区域的拷贝数数值与panel中的其它检测区域数量之间没有直接关联，panel中其它区域的数量多少对检测出的MYCN基因拷贝数并不产生影响，该数值始终在3.94-4.27之间，变化很小，并且与标准品的官方提供的拷贝数为4.75接近；而与之截然不同的是，采用算术平均值计算得到的拷贝数值随着panel中的其它检测区域数量的变化会发生明显的变化，从图3中就可以看出，只要在panel中的检测区域数量非常大时，这种方法计算得到的结果才能与标准品的数值保持一致，而在小panel的情况下，由于基因拷贝数对样本测序深度的影响更大，使计算结果与真实值之间出现了明显偏差。
另外，还采用了2个MYCN外显子分别与10个不同的其它区域构成panel来进行拷贝数的检测，分别采用表2中序号为1～10、21～30、41～50、61～70、81～90号的外显子区域与MYCN组成panel(分为称为1～5组)，依照同样的方法对MYCN拷贝数进行检测，结果如表5所示：
表5 MYCN基因拷贝数
从表中可以看出，当在MYCN基因与其它10个检测区域构成panel进行检测时，在panel中的其它区域的不同，会对采用平均值方法检测得到的结果产生明显不同，因此，这种方法的检测结果的可信度较低；而采用T分布拟合方法计算得到的结果与panel中其它的区域并没有明显的关联，得到的检测结果始终保持稳定，可信度高。
实施例2 MET基因拷贝数的检测
本实施例中检测的目标区域是含有20个外显子区域的MET基因的拷贝数进行检测，目标区域的主要信息如下表所示：
表6
区域名称 在染色体上的起始位置 在染色体上的结束位置
MET:NM_000245:e2 chr7：116339138 chr7：116340338
MET:NM_000245:e3 chr7：116371721 chr7：116371913
MET:NM_000245:e4 chr7：116380003 chr7：116380138
MET:NM_000245:e5 chr7：116380905 chr7：116381079
MET:NM_000245:e6 chr7：116395408 chr7：116395569
MET:NM_000245:e7 chr7：116397490 chr7：116397593
MET:NM_000245:e8 chr7：116397691 chr7：116397828
MET:NM_000245:e9 chr7：116398512 chr7：116398674
MET:NM_000245:e10 chr7：116399444 chr7：116399544
MET:NM_000245:e11 chr7：116403103 chr7：116403322
MET:NM_000245:e12 chr7：116409698 chr7：116409845
MET:NM_000245:e13 chr7：116411551 chr7：116411708
MET:NM_000245:e14 chr7：116411902 chr7：116412043
MET:NM_000245:e15 chr7：116414934 chr7：116415165
MET:NM_000245:e16 chr7：116417442 chr7：116417523
MET:NM_000245:e17 chr7：116418829 chr7：116419011
MET:NM_000245:e18 chr7：116422041 chr7：116422151
MET:NM_000245:e19 chr7：116423357 chr7：116423523
MET:NM_000245:e20 chr7：116435708 chr7：116435845
MET:NM_000245:e21 chr7：116435940 chr7：116436178
这里的正常人群共349位，他们的MET基因OR值组成的z分布的曲线平均值为1.01，sd为0.039，这里的病人样本用的是https://www.horizondiscovery.com公司提供的HD786标准品，它包含MET扩增，官方提供的拷贝数OR为2.25，我们利用NGS的方法以及以上描述的算法得到的该基因的OR为1.72。这里采用的对照样本为NA18535标准品细胞。
进行检测的panel中还包括了其它的一些敏感基因的外显子区域，分别取不同数量的区域与MET基因20个外显子共同构成panel，共选择了4666个其它备选区域，例如其中100个其它备选区域的信息如表2所示。
在进行计算时，将MET基因20个外显子分别和10个、50个、100个、300个、400个、……4000个、4666个其它的热点基因的外显子区域共同构成panel同时进行NGS测序，具体的panel中的区域数量如表3所示，其中采用10个、50个、100个其它区域时，是采用的表2中序号为1～10、1～50、1～100的外显子区域，其它的区域在此不一一描述。
采用的检测步骤同实施例1。将20个MET外显子区域的拷贝数的算术平均值作为MET基因拷贝数值，同时也在步骤3中采用算术平均值计算每个样本的平均测序深度作为对照。
计算结果如表7所示，区域数量与拷贝数值之间进行作图后如图4所示。
表7 MET基因拷贝数
panel组成 平均值方法 t分布方法
20个MET外显子+4666个其他区域 1.72 1.72
20个MET外显子+4000个其他区域 1.98 1.71
20个MET外显子+3000个其他区域 2.62 1.72
20个MET外显子+2000个其他区域 3.87 1.71
20个MET外显子+1000个其他区域 7.33 1.71
20个MET外显子+900个其他区域 8.14 1.7
20个MET外显子+800个其他区域 9.14 1.72
20个MET外显子+700个其他区域 10.46 1.73
20个MET外显子+600个其他区域 11.73 1.72
20个MET外显子+500个其他区域 13.74 1.71
20个MET外显子+400个其他区域 16.77 1.72
20个MET外显子+300个其他区域 21.32 1.75
20个MET外显子+200个其他区域 27.84 1.74
20个MET外显子+100个其他区域 32.21 1.72
20个MET外显子+50个其他区域 43.49 1.7
20个MET外显子+10个其他区域 57.41 1.73
这里的计算结果也与实施例1类似，本实施例中，MET基因具有20个外显子，对这两个外显子区域计算出拷贝数之后，将两者的算术平均值作为基因的拷贝数。采用本发明的方法可以避免在panel中的其它检测区域的数量的影响，使每次检测结果的数值与对标准品的实际值一致。
还采用了MET外显子分别与10个不同的其它区域构成panel来进行拷贝数的检测，分别采用表2中序号为1～10、21～30、41～50、61～70、81～90号的外显子区域与MYCN组成panel(分为称为1～5组)，依照同样的方法对MYCN拷贝数进行检测，结果如表8所示：
表8 MET基因拷贝数
panel组成 平均值方法 t分布方法
第1组 57.41 1.73
第2组 43.21 1.72
第3组 49.74 1.70
第4组 61.21 1.72
第5组 40.55 1.71
从表中可以看出，当在MET基因与其它10个检测区域构成panel进行检测时，在panel中的其它区域的不同，会对采用平均值方法检测得到的结果产生明显不同，因此，这种方法的检测结果的可信度较低；而采用T分布拟合方法计算得到的结果与panel中其它的区域并没有明显的关联，得到的检测结果始终保持稳定，可信度高。
实施例3 CDKN2A基因拷贝数的检测
本实施例中检测的目标区域是含有5个外显子区域的CDKN2A基因的拷贝数进行检测，目标区域的主要信息如下表所示：
表9
区域名称 在染色体上的起始位置 在染色体上的结束位置
CDKN2A:NM_000077:e3 chr9：21968227 chr9：21968241
CDKN2A:NM_001195132:e3 chr9：21968723 chr9：21968770
CDKN2A:NM_000077:e2 chr9：21970900 chr9：21971207
CDKN2A:NM_000077:e1 chr9：21974676 chr9：21974826
CDKN2A:NM_000077:utr3 chr9：21974827 chr9：21975132
这里的正常人群共349位，他们的CDKN2A基因OR值组成的z分布的曲线平均值为1.01，sd为0.054，这里的病人样本用的是https://www.horizondiscovery.com公司提供的HD786标准品，官方没有提供该基因的拷贝数变异，但我们利用NGS的方法以及以上描述的算法得到的该基因的OR为0.62左右。
进行检测的panel中还包括了其它的一些敏感基因的外显子区域，分别取不同数量的区域与CDKN2A基因5个外显子共同构成panel，共选择了4681个其它备选区域，例如其中100个其它备选区域的信息如表2所示。
在进行计算时，将CDKN2A基因20个外显子分别和25个、50个、100个、300个、400个、……4000个、4681个其它的热点基因的外显子区域共同构成panel同时进行NGS测序，具体的panel中的区域数量如表3所示，其中采用25个、50个、100个其它区域时，是采用的表2中序号为1～25、1～50、1～100的外显子区域，其它的区域在此不一一描述。
采用的检测步骤同实施例1。将5个CDKN2A外显子区域的拷贝数的算术平均值作为CDKN2A基因拷贝数值，同时也在步骤3中采用算术平均值计算每个样本的平均测序深度作为对照。
计算结果如表10所示，区域数量与拷贝数值之间进行作图后如图5所示。
表10 CDKN2A基因拷贝数
这里的计算结果也与实施例1类似，本实施例中，CDKN2A基因具有5个外显子，对这5个外显子区域计算出拷贝数之后，将5个拷贝数的算术平均值作为基因的拷贝数。采用本发明的方法可以避免在panel中的其它检测区域的数量的影响，使每次检测结果的数值与对标准品的实际值一致。
还采用了CDKN2A外显子分别与25个不同的其它区域构成panel来进行拷贝数的检测，分别采用表2中序号为1～25、26～50、51～75、76～100号的外显子区域与CDKN2A组成panel(分为称为1～4组)，依照同样的方法对CDKN2A拷贝数进行检测，结果如表11所示：
表11 CDKN2A基因拷贝数
panel组成 平均值方法 t分布方法
第1组 96.3 0.66
第2组 75.4 0.64
第3组 67.7 0.63
第4组 103.3 0.66
从表中可以看出，当在CDKN2A基因与其它25个检测区域构成panel进行检测时，在panel中的其它区域的不同，会对采用平均值方法检测得到的结果产生明显不同，因此，这种方法的检测结果的可信度较低；而采用T分布拟合方法计算得到的结果与panel中其它的区域并没有明显的关联，得到的检测结果始终保持稳定，可信度高。
实施例4食管癌患者样本中PIK3CA基因拷贝数变异的检测
以PIK3CA基因Variation 91720位点作为待检测区域，与实施例3中其它4681基因区域组成panel进行高通量测序，采用本发明的方法对Variation 91720位点的拷贝数进行计算。
同时，作为对照，采用基于TaqMan探针的实时荧光定量PCR检测Variation 91720位点的拷贝数，并采用Copy Caller v2.0软件对拷贝数进行计算。
病例组共33例，来自于中国南方地区的汉族人群，正常人群是以上实施例中349例。病例组年龄中位数56岁，平均55.1岁，病理类型以ESCC为主，占84.8％，其次为腺癌，占12.1％，其它病理类型共占3％。
判定标准：基因的拷贝数大于或小于基线值的mean±3sd被判定是基因存在着拷贝数扩增或缺失；对照样本在检测区域上的平均测序深度达标，大于5x，覆盖比例大于该区域长度的70％。
以上33例病例样本的Variation 91720位点的拷贝数检测结果如表12所示：
表12食管癌患者样本中PIK3CA基因拷贝数
从上表可以看出，本发明的检测结果可以应用于病例的拷贝数变异的检测，检测结果与实时荧光PCR对该基因单独进行检测的结果接近。

Claims (10)

1.一种拷贝数变异的检测方法，其特征在于，包括如下步骤：
S1，采用待检测区域和其它区域组成基因组合，对待测样本、对照样本和正常样本采用所述的基因组合进行NGS测序；
S2，根据步骤S1的结果，计算出每份正常样本和待测样本上的待检测区域的拷贝数；
S3，将正常样本中所述的待检测区域的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测区域的平均拷贝数基线；
S4，将待测样本中所述的待检测区域的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测区域的拷贝数的计算方法是：根据样本与对照样本的NGS测序结果，比对至参考基因组序列，获得唯一比对至所述的检测区域的读段，并计算出每个样本中每个检测区域上的测序深度；对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；根据平均测序深度计算出每份样本上的待检测区域的拷贝数。
2.根据权利要求1所述的拷贝数变异的检测方法，其特征在于，基因组合中检测区域的数量为3～50000个；拷贝数的计算是公式是：拷贝数=（样本的待检测区域的测序深度/样本的平均测序深度）/（对照样本的待检测区域的测序深度/对照样本的平均测序深度）×倍体数；所述的基线数值是由平均值±标准偏差所构成。
3.根据权利要求1所述的拷贝数变异的检测方法，其特征在于，所述的待检测区域选自以下基因中的任意一个片断：AKT1、AKT2、AKT3、CD274、DDR2、EGFR、ERBB2、FGFR1、FGFR3、FLT4、HGF、MET、MTOR、MYC、PDCD1LG2、PDGFRA、PDGFRB、SOX2、TP53、VEGFA、BRAF、FLT1、HRAS、KDR、KRAS、MAP2K1、MAP2K2、NRAS、PIK3CA、RB1、TOP1、VEGFA、BCL2、BCL6、CCND1、CCND3、CDK6、CEBPA、CEBPB、CEBPD、HOXA10、IL3、IRF4、KMT2A、LYL1、MUC1、MYC、NOTCH1、SETBP1、TAL1、ZAP70；所述的拷贝数变异的检测方法是用于非治疗与诊断目的。
4.一种拷贝数变异的检测方法，其特征在于，包括如下步骤：
S1，从待测基因序列中选出至少两个待检测区域，并与其它区域组成基因组合，对待测样本和正常样本采用所述的基因组合进行NGS测序；
S2，根据步骤S1的结果，计算出每份正常样本和待测样本上的待检测基因的拷贝数；
S3，将正常样本中待检测基因的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测基因的平均拷贝数基线；
S4，将待测样本中所述的待检测基因的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测基因的拷贝数是指每个待检测区域的拷贝数的平均值；
待检测区域的拷贝数的计算方法是：根据样本与对照样本的NGS测序结果，比对至参考基因组序列，获得唯一比对至所述的检测区域的读段，并计算出每个样本中每个检测区域上的测序深度；对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；根据平均测序深度计算出每份样本上的待检测区域的拷贝数。
5.根据权利要求4所述的拷贝数变异的检测方法，其特征在于，基因组合中检测区域的数量为3～50000个；拷贝数的计算是公式是：拷贝数=（样本的待检测区域的测序深度/样本的平均测序深度）/（对照样本的待检测区域的测序深度/对照样本的平均测序深度）×倍体数；每个待检测区域的拷贝数的平均值是指算术平均值、正态分布平均值卡方分布平均值、F分布平均值或者T分布平均值；所述的拷贝数变异的检测方法是用于非治疗与诊断目的。
6.一种拷贝数变异的检测装置，其特征在于，包括：
测序数据获取模块，用于对待测样本和正常样本采用由待检测区域和其它区域组成的基因组合进行NGS测序，并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段；
测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；
拷贝数计算模块，用于计算出每份正常样本和待测样本上的待检测区域的拷贝数；
基线计算模块，用于将正常样本中待检测区域的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测区域的平均拷贝数基线；
分析模块，用于将待测样本中所述的待检测区域的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异。
7.根据权利要求6所述的拷贝数变异的检测装置，其特征在于，拷贝数是通过如下公式计算得到：拷贝数=（样本的待检测区域的测序深度/样本的平均测序深度）/（对照样本的待检测区域的测序深度/对照样本的平均测序深度）×倍体数。
8.一种拷贝数变异的检测装置，其特征在于，包括：
测序数据获取模块，用于对待测样本、对照样本和正常样本采用基因组合进行NGS测序，所述的基因组合是由至少两个待检测区域和其它区域组成，所述的待检测区域选自待测基因序列；并将样本的测序下机数据与参考基因组序列对比，获得唯一比对至所述的检测区域的读段；
测序深度计算模块，用于根据获得的唯一比对至所述的检测区域的读段计算出每个样本中每个检测区域上的测序深度，并对同一份样本中的每个检测区域的测序深度进行T分布拟合，将分布曲线的平均值作为该样本上的平均测序深度；
拷贝数计算模块，用于计算出每份正常样本和待测样本上的待检测基因的拷贝数；
基线计算模块，用于将正常样本中所述的待检测基因的拷贝数进行T分布拟合，将分布曲线的平均值作为该待检测基因的平均拷贝数基线；
分析模块，用于将待测样本中所述的待检测基因的拷贝数与平均拷贝数基线进行比较，判定是否存在拷贝数变异；
其中，待测样本和正常样本的待检测基因的拷贝数是指每个待检测区域的拷贝数的平均值。
9.根据权利要求8所述的拷贝数变异的检测装置，其特征在于，拷贝数是通过如下公式计算得到：拷贝数=（样本的待检测区域的测序深度/样本的平均测序深度）/（对照样本的待检测区域的测序深度/对照样本的平均测序深度）×倍体数。
10.一种计算机可读介质，记载有可以运行权利要求1～5任一项所述的拷贝数变异的检测方法的程序。
CN201810151291.2A 2018-02-14 2018-02-14 一种拷贝数变异的检测方法、装置以及计算机可读介质 Active CN108427864B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201810151291.2A CN108427864B (zh) 2018-02-14 2018-02-14 一种拷贝数变异的检测方法、装置以及计算机可读介质
PCT/CN2018/090086 WO2019157791A1 (zh) 2018-02-14 2018-06-06 一种拷贝数变异的检测方法、装置以及计算机可读介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810151291.2A CN108427864B (zh) 2018-02-14 2018-02-14 一种拷贝数变异的检测方法、装置以及计算机可读介质

Publications (2)

Publication Number Publication Date
CN108427864A CN108427864A (zh) 2018-08-21
CN108427864B true CN108427864B (zh) 2019-01-29

Family

ID=63157045

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810151291.2A Active CN108427864B (zh) 2018-02-14 2018-02-14 一种拷贝数变异的检测方法、装置以及计算机可读介质

Country Status (2)

Country Link
CN (1) CN108427864B (zh)
WO (1) WO2019157791A1 (zh)

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109637585B (zh) * 2018-12-27 2020-11-17 北京优迅医学检验实验室有限公司 测序深度的矫正方法及装置
CN110246543B (zh) * 2019-06-21 2021-02-26 元码基因科技（北京）股份有限公司 基于二代测序技术利用单样本检测拷贝数变异的方法和计算机系统
CN112151112A (zh) * 2019-06-27 2020-12-29 天津中科智虹生物科技有限公司 一种遗传基因检测的方法和装置
CN110241191A (zh) * 2019-06-28 2019-09-17 中国人民解放军第四军医大学 一种基于NGS同时检测mtDNA拷贝数和突变的方法
CN114502744B (zh) * 2019-12-11 2023-06-23 深圳华大基因股份有限公司 一种基于血液循环肿瘤dna的拷贝数变异检测方法和装置
CN110993022B (zh) * 2019-12-20 2023-09-05 北京优迅医学检验实验室有限公司 检测拷贝数扩增的方法和装置及建立检测拷贝数扩增的动态基线的方法和装置
CN112885406B (zh) * 2020-04-16 2023-01-31 深圳裕策生物科技有限公司 检测hla杂合性缺失的方法及系统
CN111508559B (zh) * 2020-04-21 2021-08-13 北京橡鑫生物科技有限公司 检测目标区域cnv的方法及装置
CN111863124B (zh) * 2020-06-06 2024-01-30 聊城大学 一种拷贝数变异检测方法、系统、存储介质、计算机设备
CN112349346A (zh) * 2020-10-27 2021-02-09 广州燃石医学检验所有限公司 检测基因组区域中的结构变异的方法
CN112435711B (zh) * 2020-11-11 2022-04-01 赛福解码(北京)基因科技有限公司 一种改善小panel数据中大型cnv检测效果的方法
CN112768000B (zh) * 2021-01-25 2021-07-20 深圳吉因加医学检验实验室 一种预测met基因拷贝数变化类型的方法及装置
WO2023030233A1 (zh) * 2021-08-30 2023-03-09 广州燃石医学检验所有限公司 一种拷贝数变异的检测方法及其应用
CN114582427B (zh) * 2022-03-22 2023-04-07 成都基因汇科技有限公司 一种渐渗区段鉴定方法及计算机可读存储介质
CN114758720B (zh) * 2022-06-14 2022-09-02 北京贝瑞和康生物技术有限公司 用于检测拷贝数变异的方法、设备和介质
CN115376609B (zh) * 2022-10-24 2023-03-10 广州燃石医学检验所有限公司 一种判别met基因拷贝数扩增类型的方法及装置
CN116386718B (zh) * 2023-05-30 2023-08-01 北京华宇亿康生物工程技术有限公司 检测拷贝数变异的方法、设备和介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106755444A (zh) * 2016-12-31 2017-05-31 东北农业大学 一种大豆基因拷贝数变异分析方法
CN107435070A (zh) * 2012-04-12 2017-12-05 维里纳塔健康公司 拷贝数变异的检测和分类

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101795124B1 (ko) * 2012-04-05 2017-12-01 비지아이 다이어그노시스 씨오., 엘티디. 복제 수 변이를 검측하기 위한 방법 및 시스템
CN105722994B (zh) * 2013-06-17 2020-12-18 维里纳塔健康公司 用于确定性染色体中的拷贝数变异的方法
CN104133914B (zh) * 2014-08-12 2017-03-08 厦门万基生物科技有限公司 一种消除高通量测序引入的gc偏差及对染色体拷贝数变异的检测方法
CN105574361B (zh) * 2015-11-05 2018-11-02 上海序康医疗科技有限公司 一种检测基因组拷贝数变异的方法
CN105349678A (zh) * 2015-12-03 2016-02-24 上海美吉生物医药科技有限公司 一种染色体拷贝数变异的检测方法
CN105760712B (zh) * 2016-03-01 2019-03-26 西安电子科技大学 一种基于新一代测序的拷贝数变异检测方法
CN107423534B (zh) * 2016-05-24 2021-08-06 郝柯 基因组拷贝数变异的检测方法和系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107435070A (zh) * 2012-04-12 2017-12-05 维里纳塔健康公司 拷贝数变异的检测和分类
CN106755444A (zh) * 2016-12-31 2017-05-31 东北农业大学 一种大豆基因拷贝数变异分析方法

Also Published As

Publication number Publication date
CN108427864A (zh) 2018-08-21
WO2019157791A1 (zh) 2019-08-22

Similar Documents

Publication Publication Date Title
CN108427864B (zh) 一种拷贝数变异的检测方法、装置以及计算机可读介质
JP7145907B2 (ja) 疾患細胞不均一性を示す疾患の検出および処置、ならびに通信試験結果のためのシステムおよび方法
CN106047998B (zh) 一种肺癌基因的检测方法及应用
AU2023248050A1 (en) Diagnostic methods
CN109880910A (zh) 一种肿瘤突变负荷的检测位点组合、检测方法、检测试剂盒及系统
JP2019531700A5 (zh)
Zhang et al. Detecting false expression signals in high-density oligonucleotide arrays by an in silico approach
CN108603228A (zh) 通过分析无细胞dna确定肿瘤基因拷贝数的方法
JP2017506875A5 (zh)
CN109427412A (zh) 用于检测肿瘤突变负荷的序列组合和其设计方法
CN105925665A (zh) 试剂盒、建库方法以及检测目标区域变异的方法及系统
CN105555970A (zh) 同时进行单体型分析和染色体非整倍性检测的方法和系统
CN110023509A (zh) 基因型分型测定中的非独特条形码
WO2019046804A1 (en) IDENTIFICATION OF FALSE POSITIVE VARIANTS USING A MODEL OF IMPORTANCE
CN105950709A (zh) 试剂盒、建库方法以及检测目标区域变异的方法及系统
US20230057154A1 (en) Somatic variant cooccurrence with abnormally methylated fragments
KR20220125708A (ko) 차세대 염기서열분석 기반 표적유전자 rna 염기서열 분석 패널 및 분석알고리즘
CN115491423A (zh) 一种用于b细胞淋巴瘤mrd监测的基因组合、试剂盒与应用
TW201920682A (zh) 多型之檢測方法
CN112837749B (zh) 一种癌症筛查用基因芯片探针的优选方法
CN109136371B (zh) 一种放疗疗效和毒性反应相关基因组合、检测探针库以及检测试剂盒
CN113981080A (zh) 晚期三阴性乳腺癌铂类治疗敏感性的预测指标生成分析方法
KR20240049800A (ko) 비정상적으로 메틸화된 단편을 갖는 체세포 변이 동시 발생
WO2023164713A1 (en) Probe sets for a liquid biopsy assay
Mariani et al. Microarray Techniques and Data in Asthma/Chronic Obstructuve Pulmonary Disease

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20180821

Assignee: Nanjing and the medical equipment Co., Ltd.

Assignor: GENESEEQ TECHNOLOGY Inc.

Contract record no.: 2019320000140

Denomination of invention: Detection methods and devices of copy number variations (CNVs) and computer-readable medium

Granted publication date: 20190129

License type: Common License

Record date: 20190514

EE01 Entry into force of recordation of patent licensing contract
CP01 Change in the name or title of a patent holder

Address after: 210032 No. 3-1 Xinjinhu Road, Nanjing High-tech Development Zone, Jiangsu Province

Patentee after: Nanjing Shihe gene Biotechnology Co., Ltd

Address before: 210032 No. 3-1 Xinjinhu Road, Nanjing High-tech Development Zone, Jiangsu Province

Patentee before: GENESEEQ TECHNOLOGY Inc.

CP01 Change in the name or title of a patent holder