CN107480470B - 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 - Google Patents
基于贝叶斯与泊松分布检验的已知变异检出方法和装置 Download PDFInfo
- Publication number
- CN107480470B CN107480470B CN201610407552.3A CN201610407552A CN107480470B CN 107480470 B CN107480470 B CN 107480470B CN 201610407552 A CN201610407552 A CN 201610407552A CN 107480470 B CN107480470 B CN 107480470B
- Authority
- CN
- China
- Prior art keywords
- model
- variation
- probability
- site
- 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
Links
- 238000012360 testing method Methods 0.000 title claims abstract description 46
- 238000001514 detection method Methods 0.000 title claims abstract description 42
- 238000012163 sequencing technique Methods 0.000 claims abstract description 102
- 238000000034 method Methods 0.000 claims abstract description 31
- 230000035772 mutation Effects 0.000 claims description 49
- 230000009897 systematic effect Effects 0.000 claims description 15
- 230000001186 cumulative effect Effects 0.000 claims description 11
- 108700028369 Alleles Proteins 0.000 claims description 7
- 238000012217 deletion Methods 0.000 claims description 7
- 230000037430 deletion Effects 0.000 claims description 7
- 230000036438 mutation frequency Effects 0.000 claims description 7
- 108090000623 proteins and genes Proteins 0.000 claims description 6
- 238000003780 insertion Methods 0.000 claims description 4
- 230000037431 insertion Effects 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 4
- 239000002773 nucleotide Substances 0.000 claims description 3
- 125000003729 nucleotide group Chemical group 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 abstract description 11
- 238000005516 engineering process Methods 0.000 description 10
- 238000012165 high-throughput sequencing Methods 0.000 description 7
- 239000000523 sample Substances 0.000 description 6
- 238000004364 calculation method Methods 0.000 description 5
- 238000003908 quality control method Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 230000002759 chromosomal effect Effects 0.000 description 3
- 210000000349 chromosome Anatomy 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 206010028980 Neoplasm Diseases 0.000 description 2
- 201000011510 cancer Diseases 0.000 description 2
- JJWKPURADFRFRB-UHFFFAOYSA-N carbonyl sulfide Chemical compound O=C=S JJWKPURADFRFRB-UHFFFAOYSA-N 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003752 polymerase chain reaction Methods 0.000 description 2
- 210000001519 tissue Anatomy 0.000 description 2
- 208000010507 Adenocarcinoma of Lung Diseases 0.000 description 1
- 102100028630 Cytoskeleton-associated protein 2 Human genes 0.000 description 1
- 101150039808 Egfr gene Proteins 0.000 description 1
- 206010064571 Gene mutation Diseases 0.000 description 1
- 101000766848 Homo sapiens Cytoskeleton-associated protein 2 Proteins 0.000 description 1
- 108091034117 Oligonucleotide Proteins 0.000 description 1
- 108020005187 Oligonucleotide Probes Proteins 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- WOERBKLLTSWFBY-UHFFFAOYSA-M dihydrogen phosphate;tetramethylazanium Chemical compound C[N+](C)(C)C.OP(O)([O-])=O WOERBKLLTSWFBY-UHFFFAOYSA-M 0.000 description 1
- 102000052116 epidermal growth factor receptor activity proteins Human genes 0.000 description 1
- 108700015053 epidermal growth factor receptor activity proteins Proteins 0.000 description 1
- 108700021358 erbB-1 Genes Proteins 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 238000009396 hybridization Methods 0.000 description 1
- 150000002500 ions Chemical class 0.000 description 1
- 201000005249 lung adenocarcinoma Diseases 0.000 description 1
- YOHYSYJDKVYCJI-UHFFFAOYSA-N n-[3-[[6-[3-(trifluoromethyl)anilino]pyrimidin-4-yl]amino]phenyl]cyclopropanecarboxamide Chemical compound FC(F)(F)C1=CC=CC(NC=2N=CN=C(NC=3C=C(NC(=O)C4CC4)C=CC=3)C=2)=C1 YOHYSYJDKVYCJI-UHFFFAOYSA-N 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 239000002751 oligonucleotide probe Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000001717 pathogenic effect Effects 0.000 description 1
- 238000002203 pretreatment Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000007480 sanger sequencing Methods 0.000 description 1
- 230000000392 somatic effect Effects 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
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
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种基于贝叶斯与泊松分布检验的已知变异检出方法和装置。该方法包括：提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；将所述假设已知变异位点存在时推算出的测序读长序列与所述测序读长碱基序列进行比对检测；进行贝叶斯检验模型和泊松分布检验模型的判断。本发明的方法能够高特异性、高敏感性地实现SNV或InDel变异检出。
Description
技术领域
本发明涉及核苷酸序列的变异位点检测技术领域，尤其涉及一种基于贝叶斯与泊松分布检验的已知变异检出方法和装置。
背景技术
在科研及临床转化领域的技术中，对于体细胞基因突变的检测方法主要有PCR-Sanger法、PCR-Mass法、ARMS-PCR法和高通量捕获测序法等。这几种方法的样本均来源于癌症患者手术时所切除的病灶组织样本。PCR-Sanger测序法的灵敏度低，仅能检测突变率大于20％的高频突变，且难以胜任多基因多位点的检测；PCR-Mass法不能检测位置突变，也不能检测插入删除变异(InDel)，灵敏度在5％左右，该技术局限性较大；ARMS-PCR法特异性好，灵敏度高，能检测1％-5％的低频突变，然而该方法只能检测特定的已知突变，且不能同时检测多个基因多个位点。
BGISEQ-100平台的数据特点，导致传统的变异检出软件(如SOAPsnp、VarScan、GATK等)的检出效果不理想。BGISEQ-100平台数据在InDel附近的碱基比对质量会下降，尤其对EGFR基因c.2238_2248>GC这类常见的复杂变异，删除后插入的GC碱基可能会比对到不同的位置，无法正确比对到参考基因组，从而导致传统变异检出软件的检测效果不理想。
发明内容
本发明提供一种基于贝叶斯与泊松分布检验的已知变异检出方法和装置，能够高特异性、高敏感性地实现SNV或InDel变异检出。
根据本发明的第一方面，本发明提供一种基于贝叶斯与泊松分布检验的已知变异检出方法，包括：提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；将上述假设已知变异位点存在时推算出的测序读长序列与上述测序读长碱基序列进行比对检测，找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列；针对上述变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与上述参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由上述参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差，判断上述模型
的概率与模型M0的概率之比值与第一阈值的关系；针对上述变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，且上述泊松分布累计概率值大于上述第二阈值，判断该位点为强阳性变异；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，或上述泊松分布累计概率值大于上述第二阈值，判断该位点为弱阳性变异；若上述模型
的概率与模型M0的概率之比值小于上述第一阈值，且上述泊松分布累计概率值小于等于上述第二阈值，判断该位点为阴性无变异。
根据本发明的第二方面，本发明提供一种基于贝叶斯与泊松分布检验的已知变异检出装置，包括：数据输入单元，用于提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；比对检测单元，用于将上述假设已知变异位点存在时推算出的测序读长序列与上述测序读长碱基序列进行比对检测，找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列；模型存储单元，用于存储贝叶斯检验模型和泊松分布检验模型，其中，针对上述变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与上述参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由上述参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差；针对上述变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n；变异判断单元，用于判断上述模型
的概率与模型M0的概率之比值与第一阈值的关系，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，且上述泊松分布累计概率值大于上述第二阈值，判断该位点为强阳性变异；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，或上述泊松分布累计概率值大于上述第二阈值，判断该位点为弱阳性变异；若上述模型
的概率与模型M0的概率之比值小于上述第一阈值，且上述泊松分布累计概率值小于等于上述第二阈值，判断该位点为阴性无变异；数据输出单元，用于输出上述变异判断单元所判断的变异数据结果。
本发明的基于贝叶斯与泊松分布检验的已知变异检出方法，不关注变异具体的比对位置和比对形式，而是关注在测序得到的读长中是否存在该变异发生时的序列特征，从而规避了InDel附近比对质量下降的情况。同时，对于检出的变异结果，采用泊松分布检验及贝叶斯检验模型对检出结果进行监测，以降低假阳性。
附图说明
图1是依据本发明的一种实施方式的检出方法的流程示意图；
图2是依据本发明的一种实施方式的检出装置的结构框图。
具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。
结合传统实验技术和高通量测序技术的特点，为了更快速、准确地同时检测多基因多位点，本发明提出了利用高通量测序技术产生的测序数据，对重要的单核苷酸变异(SNV)和/或插入删除变异(InDel)变异位点进行快速、精确检测的方法。针对BGISEQ-100测序平台数据特点，使用寡核苷酸探针捕获技术或PCR多重扩增的方式来获取基因组上的目标序列，对目标序列产物进行高通量测序，从中识别DNA样品中的碱基序列和变异信息。在对SNV、InDel的检测中，本发明针对已知变异位点的特性，根据不同的检测变异位点和COSMIC(Catalogue of Somatic Mutations in Cancer)数据库中记载的致病变异位点，推算出该变异位点存在时测序读长(即reads，也称为“读段”)的碱基序列，然后在BGISEQ-100数据中对这种序列进行检测。在这种情况下，复杂的InDel变异的检测，不再关注其具体的比对位置和比对形式，而是关注在测序得到的读长中是否存在该变异发生时的序列特征，从而规避了InDel附近比对质量下降的情况。同时，对于检出的变异结果，采用泊松分布检验及贝叶斯检验模型对检出结果进行监测，以降低假阳性。
依据本发明的一种实施方式，提供一种基于贝叶斯与泊松分布检验的已知变异检出方法，参考图1，包括如下步骤：
S101：提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算
出的测序读长序列。
测序读长碱基序列，即利用高通量测序技术产生的测序序列，尤其是利用BGISEQ-100测序平台产生的测序数据。本发明的优选实施例中，特别采用的是高通量捕获测序法，该方法利用设计的芯片探针，对样本DNA中感兴趣的目标基因区域进行杂交捕获和高通量测序，对测序结果进行生物信息学分析，从而得到变异位点的检出信息。该方法较其它方法具有灵敏度高、通量高、可同时检测几十甚至上百个基因和位点的优势。
本发明中的样本DNA，可以来源于任何待检测个体基因组DNA，尤其是来源于人体的基因组DNA。相应地，参考基因组序列可以是与待检测个体基因组DNA对应的物种的基因组参考序列，尤其是人类基因组序列。例如，参考基因组序列可以选择为NCBI数据库中版本37.3(hg19；NCBI Build 37.3)的人类基因组参考序列。
本发明使用的测序读长碱基序列可以是测序下机数据。然而，为了提高测序下机数据使用的有效性，避免不必要的冗余和对无效数据不必要的分析，本发明使用的测序读长碱基序列优选是对测序下机数据进行预处理和质控以后得到的数据。这样的预处理优选包括：与参考基因组比对，排序，去重，建索引。
对于与参考基因组比对，例如，对BGISEQ-100测序平台的有效测序数据使用Tmap工具比对到参考基因组(例如hg19)上，得到精确的比对结果。其中Tmap工具源自：https://github.com/iontorrent/TS/tree/master/Analysis/TMAP。
对于排序，例如使用samtools sort对Tmap工具比对后的结果(bam文件)进行排序：根据染色体编号和所在染色体上的位置按从小到大的顺序进行排序。
对于去重，即去掉比对结果的重复片段，例如对排序后的结果(bam文件)使用BamDuplicates工具去除重复片段。其中BamDuplicates工具源自Ion Torrent Systems,Inc。
对于建索引，例如对去重后的bam文件利用samtools index建立相应的索引。
对于质控，例如对建索引后的bam文件进行QC质控，合格的文件将进行后续的变异检出步骤。
本发明的方法主要是针对已知变异进行检测，即检测待检测个体基因组DNA中是否存在已知变异。如果待检测个体基因组DNA中存在已知变异，则会在测序读长碱基序列中存在该变异发生时的序列特征，即变异特征。例如，假设已知变异是缺失(或称“删除”)了一个特定碱基，那么序列特征即表现为相对于参考基因组序列缺失了该特定碱基；同样地，假设已知变异是插入了一个特定碱基，那么序列特征即表现为相对于参考基因组序列插入了该特定碱基。又例如，假设需要检测EGFR_L858R，首先在COSMIC中得知其具体的变异为c.2573T>G，将其转换为染色体上的位点为chr7_55259515_T>G，根据染色体位点和参考基因组hg19的序列，得到变异位点的前后5bp的序列，则可知变异后的位点加上前后5bp的序列为TGGGCGGGCCA。因此，根据已知变异即可推算出假设已知变异位点存在时的测序读长序列。问题就转化为，在测序读长碱基序列中寻找假设已知变异位点存在时推算出的测序读长序列。
S102：将所述假设已知变异位点存在时推算出的测序读长序列与所述测序读长碱
基序列进行比对检测。
在将推算出的测序读长序列比对到测序读长碱基序列时，可使用各种现有的比对软件，包括但不限于Tmap，BWA(Burrows-Wheeler Aligner)，SOAP(Short OligonucleotideAnalysis Package)，samtools等，本实施方式对此不作限定。本领域技术人员也可以根据具体需要编程这样的软件用于实现将假设已知变异位点存在时推算出的测序读长序列与测序读长碱基序列进行比对检测。
通过将假设已知变异位点存在时推算出的测序读长序列与测序读长碱基序列进行比对检测，即可找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列。
例如，对于QC质控之后的bam文件，使用比对软件将假设已知变异位点存在时推算出的测序读长序列比对到bam文件中，检测已知变异位点是否存在，即检测在bam文件中是否存在该变异发生时的序列特征。然后，对于检测出的已知变异位点，根据已知变异位点，在bam文件中找到能覆盖到该变异位点的所有测序读长序列，并且统计该变异位点上的总读长序列覆盖深度(即该位点的测序深度)、发生变异的正负链读长数(包括发生变异的正链读长数和发生变异的负链读长数)，以及覆盖该位点的所有读长在该位点的碱基质量值。接下来即可根据上述统计结果进行贝叶斯检验模型和泊松分布检验模型的判断，以确切地确定出每一变异位点的具体变异情况。
S103：进行贝叶斯检验模型和泊松分布检验模型的判断。
对于贝叶斯检验模型和泊松分布检验模型分别作如下假设和模型的建立，以及根据所建立的模型进行变异类型的判断。
针对找到的变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差，判断所述模型
的概率与模型M0的概率之比值与第一阈值的关系。
其中，模型
的概率与模型M0的概率之比值，可以是指模型
的概率与模型M0的概率直接相除得到的比值。然而，考虑到模型
的概率与模型M0的概率直接相除得到的比值可能很大，例如大于100，因此模型
的概率与模型M0的概率之比值，也可以是指对模型
的概率与模型M0的概率直接相除得到的比值进行相关计算(例如取对数)得到的数值，例如，对模型
的概率与模型M0的概率直接相除得到的比值取以10为底的Log值(即Log10)，如果模型
的概率与模型M0的概率直接相除得到的比值为100，那么取以10为底数的对数以后得到的数值即为2。相应地，相对于模型
的概率与模型M0的概率直接相除得到的比值，在取对数前后，第一阈值是不同的。例如，如果以模型
的概率与模型M0的概率直接相除得到的比值作为比较对象，第一阈值可以设置为100；如果以模型
的概率与模型M0的概率直接相除得到的比值取以10为底的Log值，第一阈值可以设置为2。
针对找到的变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系。
根据上述贝叶斯检验模型和泊松分布检验模型，若模型
的概率与模型M0的概率之比值大于等于第一阈值，且泊松分布累计概率值大于第二阈值，判断该位点为强阳性变异；若模型
的概率与模型M0的概率之比值大于等于第一阈值，或泊松分布累计概率值大于第二阈值，判断该位点为弱阳性变异；若模型
的概率与模型M0的概率之比值小于第一阈值，且泊松分布累计概率值小于等于第二阈值，判断该位点为阴性无变异。
本发明中，第一阈值和第二阈值均可以通过对大量测序数据和最终确定的变异类型进行经验性地总结和分析得到。
一般而言，第一阈值的设置需要考虑假阳性和假阴性出现的可能性大小。如果第一阈值的设置偏大，就意味着对于某一特定位点而言，模型
的概率需要远大于模型M0的概率，才能判断该特定位点是阳性变异位点，可能导致将弱阳性变异判断为阴性无变异，即可能导致假阴性出现。相反，如果第一阈值的设置偏小，就意味着对于某一特定位点而言，模型
的概率只要略大于模型M0的概率，就能判断该特定位点是阳性变异位点，可能导致将阴性无变异判断为弱阳性变异，即可能导致假阳性出现。
同理，第二阈值的设置也需要考虑假阳性和假阴性出现的可能性大小。如果第二阈值的设置偏大，则可能导致将弱阳性变异判断为阴性无变异，即可能导致假阴性出现。相反，如果第二阈值的设置偏小，可能导致将阴性无变异判断为弱阳性变异，即可能导致假阳性出现。
在本发明的优选的实施例中，模型
的概率与模型M0的概率之比值是对模型
的概率与模型M0的概率直接相除得到的比值取以10为底数的对数得到的数据。
具体地，在本发明的优选的实施例中，按照如下假设和模型建立贝叶斯检验模型。
(1)模型的假设：
(1.1)对于任一位点，假设参考基因组对应的碱基为r∈{A,T,C,G}；
(1.2)对于任一位点，假设覆盖该位点的所有读长的对应碱基为bi，碱基质量值为qi，则对应的碱基错误率为
(d表示该位点对应的测序深度)。
(2)模型的建立：
对于每一个位点的数据分布情况分为两种模型来解释：模型M0和模型
模型M0：这个位点不存在变异，与参考基因组不同的那些碱基都是系统误差导致的；
模型
这个位点的变异r→m是真实存在的，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差处理。
该位点的数据分布情况能够当作模型M0来处理的概率为：
该位点的数据分布能够当做模型
来处理的概率为：
其中，
到此，变异检出的问题就转换为判断位点的数据分布情况更偏向于哪个模型，也即对两个概率L(M0)和
进行比较，于是建立如下的变异检出模型：
一般情况下，
与L(M0)的差异都是数量级上的差异，因此
的值会很大，所以会对其取对数(以10为底的对数)的操作。
其中，
为参考值，θ即为对应的cut off值(第一阈值)。
本发明的上述优选的实施例中，θ的取值为2，能够尽可能地避免假阳性和假阴性结果的出现。
在本发明的一个实施方式中，给出一种θ的确定过程，具体如下：
步骤一：
针对θ∈(0,0.1,0.2,...,10)里的每一个θ值，计算满足
的最小k值，其中k为带有突变的reads数，即
然后通过二项分布的概率计算公式计算相应的灵敏性：
其中，f(1-e)+(1-f)e为reads带有突变的概率。
步骤二：
针对θ∈(0,0.1,0.2,...,10)里的每一个θ值，计算满足
的最小k值，其中k为不带突变的reads数，即
然后通过二项分布的概率计算公式计算相应的特异性：
步骤三：
这样可以得到每个θ对应的灵敏性和特异性数值，作灵敏性和特异性的ROC曲线，选择最佳的θ值(权衡最佳的灵敏性和特异性)即是2。然后在实际运用中，通过大量的回顾性样本测试阈值取2是符合要求的。
在本发明的优选的实施例中，按照如下假设和模型建立泊松分布检验模型。
(1)当测序深度一定时已知变异位点发生测序错误的读长条数(λ)服从泊松分布；
(2)其对应的概率分布列为：
其中，k代表测序深度；
(3)建立假设：假设带有已知变异特征的读长是由测序错误导致的且读长条数为n，那么n应该服从参数为λ的泊松分布，即累计概率值
应该小于等于第二阈值；
(4)如果计算出来的累积概率值大于第二阈值，则拒绝原假设，即带有已知变异特征的读长不是由测序错误导致的，而是真实变异的存在。
本发明的上述优选的实施例中，第二阈值取值为0.95，能够尽可能地避免假阳性和假阴性结果的出现。就是说，如果计算出来的累积概率值小于等于0.95，那么带有已知变异特征的读长是由测序错误导致，判断该位点为阴性无变异；如果计算出来的累积概率值大于0.95，那么带有已知变异特征的读长不是由测序错误导致的，而是真实变异的存在，判断该位点为阳性变异(强阳性变异或弱阳性变异)。
本发明的已知变异检出方法，有效地整合传统实验技术和高通量测序技术的优点。相比传统实验技术，本发明的经济成本更低，对检测的内容更加灵活，方便增加新的检测突变；相比传统的高通量测序技术，本发明只关注特定位点的确切变异，检测速度更快，灵敏度更高。由于对SNV和InDel采用了新的检测策略，有效地解决了InDel区域比对质量下降对变异检出的影响，同时能够在相同的比对质量下更好地检测到复杂的InDel变异，并针对BGISEQ-100的数据做了专门的优化。本发明的已知变异检出方法能够提高数据的利用率，降低数据分析的时间。
本领域普通技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。
请参考图2，依据本发明的另一方面还提供一种基于贝叶斯与泊松分布检验的已知变异检出装置，包括：数据输入单元201，用于提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；比对检测单元202，用于将上述假设已知变异位点存在时推算出的测序读长序列与上述测序读长碱基序列进行比对检测，找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列；模型存储单元203，用于存储贝叶斯检验模型和泊松分布检验模型，其中，针对上述变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与上述参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由上述参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差；针对上述变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n；变异判断单元204，用于判断上述模型
的概率与模型M0的概率之比值与第一阈值的关系，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，且上述泊松分布累计概率值大于上述第二阈值，判断该位点为强阳性变异；若上述模型
的概率与模型M0的概率之比值大于等于上述第一阈值，或上述泊松分布累计概率值大于上述第二阈值，判断该位点为弱阳性变异；若上述模型
的概率与模型M0的概率之比值小于上述第一阈值，且上述泊松分布累计概率值小于等于上述第二阈值，判断该位点为阴性无变异；数据输出单元205，用于输出上述变异判断单元所判断的变异数据结果。
以下结合具体实施例对依据本发明的已知变异检出方法及其运行结果进行详细的描述。此处的实例仅仅是用来解释本发明，并不用于限定本发明。
实施例1
本实施例的检出方法所使用的具体参数设置为：
样本DNA：女性左上肺腺癌患者的FFPE组织样本；
样本处理以及测序下机数据：目标区域捕获以及BGISEQ-100平台测序；
测序下机数据预处理和质控：测序下机的有效数据通过tmap比对参考基因组、samtools sort排序、BamDuplicates去重、samtools index建索引、QC质控，合格的文件进行后续的变异检出步骤；
参考基因组：NCBI数据库中版本37.3(hg19；NCBI Build 37.3)的人类基因组参考序列；
假设已知变异位点存在时推算出的测序读长序列：按照表1中野生型和变异列数据和参考基因组序列推算；
贝叶斯检验模型：
(L(M0)和
参考上文中描述)；
第一阈值：取值为2；
泊松分布检验模型：累计概率值
第二阈值：取值为0.95。
统计样本在每个已知位点的变异情况，结果如表2所示(注意：表2中行顺序与表1中行顺序对应，例如表2中第2行对应表1中第2行，其余类似)。最后一列为“SCX”表示在特定的位点存在相应的变异，为强阳性结果；最后一列为“WCX”表示在特定的位点可能存在相应的变异，最好进行第三方验证，为弱阳性结果；最后一列为“NEG”表示在特定的位点不可能存在相应的变异，为阴性结果。
表1
表2
以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。
Claims (10)
1.一种基于贝叶斯与泊松分布检验的已知变异检出方法，其特征在于，所述方法包括：
提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；
将所述假设已知变异位点存在时推算出的测序读长序列与所述测序读长碱基序列进行比对检测，找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列；
针对所述变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与所述参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由所述参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差，判断所述模型
的概率与模型M0的概率之比值与第一阈值的关系；
针对所述变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系；
若所述模型
的概率与模型M0的概率之比值大于等于所述第一阈值，且所述泊松分布累计概率值大于所述第二阈值，判断该位点为强阳性变异；若所述模型
的概率与模型M0的概率之比值大于等于所述第一阈值且所述泊松分布累计概率值小于等于所述第二阈值，或所述模型
的概率与模型M0的概率之比值小于所述第一阈值且所述泊松分布累计概率值大于所述第二阈值，判断该位点为弱阳性变异；若所述模型
的概率与模型M0的概率之比值小于所述第一阈值，且所述泊松分布累计概率值小于等于所述第二阈值，判断该位点为阴性无变异。
2.根据权利要求1所述的方法，其特征在于，所述模型
的概率与模型M0的概率之比值为LOD(m,f)，其满足如下公式(1)：
其中，L(M0)和
分别表示模型M0和模型
的概率；
对于任一位点，假设参考基因组对应的碱基为r∈{A,T,C,G}；假设覆盖该位点的所有读长的对应碱基为bi，碱基质量值为qi，对应的碱基错误率为
d表示该位点对应的测序深度；
L(M0)和
分别满足如下公式(2)～(3)和(4)～(5)：
其中,
其中，
3.根据权利要求2所述的方法，其特征在于，所述第一阈值为2。
4.根据权利要求1所述的方法，其特征在于，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，服从泊松分布，其对应的概率分布列为如下公式：
其中，k代表测序深度；
所述泊松分布累计概率值为
5.根据权利要求4所述的方法，其特征在于，所述第二阈值为0.95。
6.根据权利要求1-5任一项所述的方法，其特征在于，所述变异包括单核苷酸变异和/或插入删除变异。
7.根据权利要求1-5任一项所述的方法，其特征在于，所述测序读长碱基序列是人基因序列，所述参考基因组序列是人类基因组hg19序列。
8.根据权利要求1-5任一项所述的方法，其特征在于，所述方法还包括：对测序下机数据进行预处理和质控以得到所述测序读长碱基序列。
9.根据权利要求8所述的方法，其特征在于，所述预处理包括：与参考基因组比对，排序，去重，建索引。
10.一种基于贝叶斯与泊松分布检验的已知变异检出装置，其特征在于，所述装置包括：
数据输入单元，用于提供测序读长碱基序列、参考基因组序列和假设已知变异位点存在时推算出的测序读长序列；
比对检测单元，用于将所述假设已知变异位点存在时推算出的测序读长序列与所述测序读长碱基序列进行比对检测，找到每一位点变异发生时的变异特征并找到能覆盖到该位点的所有测序读长碱基序列；
模型存储单元，用于存储贝叶斯检验模型和泊松分布检验模型，其中，针对所述变异特征对应的每一位点，在贝叶斯检验模型下，假设模型M0代表该位点不存在变异，与所述参考基因组序列不同的碱基是系统误差，假设模型
代表该位点由所述参考基因组碱基r变异为m真实存在，并且等位基因突变频率为f，对于既不为r也不为m的碱基当作系统误差；针对所述变异特征对应的每一位点，在泊松分布检验模型下，假设当测序深度一定时已知变异位点发生测序错误的读长条数为λ，假设带有已知变异特征的读长是由测序错误导致的且读长条数为n；
变异判断单元，用于判断所述模型
的概率与模型M0的概率之比值与第一阈值的关系，判断n服从参数为λ的泊松分布累计概率值与第二阈值的关系；若所述模型
的概率与模型M0的概率之比值大于等于所述第一阈值，且所述泊松分布累计概率值大于所述第二阈值，判断该位点为强阳性变异；若所述模型
的概率与模型M0的概率之比值大于等于所述第一阈值且所述泊松分布累计概率值小于等于所述第二阈值，或所述模型
的概率与模型M0的概率之比值小于所述第一阈值且所述泊松分布累计概率值大于所述第二阈值，判断该位点为弱阳性变异；若所述模型
的概率与模型M0的概率之比值小于所述第一阈值，且所述泊松分布累计概率值小于等于所述第二阈值，判断该位点为阴性无变异；
数据输出单元，用于输出所述变异判断单元所判断的变异数据结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610407552.3A CN107480470B (zh) | 2016-06-08 | 2016-06-08 | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610407552.3A CN107480470B (zh) | 2016-06-08 | 2016-06-08 | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107480470A CN107480470A (zh) | 2017-12-15 |
| CN107480470B true CN107480470B (zh) | 2020-08-11 |
Family
ID=60593595
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610407552.3A Active CN107480470B (zh) | 2016-06-08 | 2016-06-08 | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107480470B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108763872B (zh) * | 2018-04-25 | 2019-12-06 | 华中科技大学 | 一种分析预测癌症突变影响lir模体功能的方法 |
| CN111919257B (zh) * | 2018-07-27 | 2021-05-28 | 思勤有限公司 | 降低测序数据中的噪声的方法和系统及其实施和应用 |
| CN109411023B (zh) * | 2018-09-30 | 2022-03-18 | 华中农业大学 | 一种基于贝叶斯网络推理的基因间交互关系挖掘方法 |
| CN109658983B (zh) * | 2018-12-20 | 2019-11-19 | 深圳市海普洛斯生物科技有限公司 | 一种识别和消除核酸变异检测中假阳性的方法和装置 |
| CN109727638B (zh) * | 2018-12-27 | 2021-08-17 | 北京优迅医学检验实验室有限公司 | 测序深度的矫正方法及装置 |
| CN109637586B (zh) * | 2018-12-27 | 2020-11-17 | 北京优迅医学检验实验室有限公司 | 测序深度的矫正方法及装置 |
| CN109637585B (zh) * | 2018-12-27 | 2020-11-17 | 北京优迅医学检验实验室有限公司 | 测序深度的矫正方法及装置 |
| CN109785899B (zh) * | 2019-02-18 | 2020-01-07 | 东莞博奥木华基因科技有限公司 | 一种基因型校正的装置和方法 |
| CN114743597A (zh) * | 2022-03-30 | 2022-07-12 | 深圳华大医学检验实验室 | 一种基于碱基序列分析物种的方法及装置 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ES2741966T3 (es) * | 2011-12-31 | 2020-02-12 | Bgi Genomics Co Ltd | Método para detectar una variación genética |
| US20150154352A1 (en) * | 2012-06-21 | 2015-06-04 | Gigagen, Inc. | System and Methods for Genetic Analysis of Mixed Cell Populations |
| SI3011051T1 (sl) * | 2013-06-21 | 2019-05-31 | Sequenom, Inc. | Postopek za neinvazivno oceno genetskih variacij |
-
2016
- 2016-06-08 CN CN201610407552.3A patent/CN107480470B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN107480470A (zh) | 2017-12-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107480470B (zh) | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 | |
| CN106909806B (zh) | 定点检测变异的方法和装置 | |
| CN108573125B (zh) | 一种基因组拷贝数变异的检测方法及包含该方法的装置 | |
| CN104302781B (zh) | 一种检测染色体结构异常的方法及装置 | |
| JP6314091B2 (ja) | Ｄｎａ配列のデータ分析 | |
| CN110692101B (zh) | 用于比对靶向的核酸测序数据的方法 | |
| CN111755068B (zh) | 基于测序数据识别肿瘤纯度和绝对拷贝数的方法及装置 | |
| CN110060733B (zh) | 基于单样本的二代测序肿瘤体细胞变异检测装置 | |
| CN108304694B (zh) | 基于二代测序数据分析基因突变的方法 | |
| CN110621785A (zh) | 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置 | |
| AU2016355983A1 (en) | Methods for detecting copy-number variations in next-generation sequencing | |
| CN117106870B (zh) | 胎儿浓度的确定方法及装置 | |
| WO2019132010A1 (ja) | 塩基配列における塩基種を推定する方法、装置及びプログラム | |
| CN114067908B (zh) | 一种评估单样本同源重组缺陷的方法、装置和存储介质 | |
| JP2021503128A (ja) | アラインされていないシーケンシングデータの高速品質管理のためのｋ−ｍｅｒの使用 | |
| WO2019213810A1 (zh) | 检测染色体非整倍性的方法、装置及系统 | |
| CA3096353C (en) | Determination of frequency distribution of nucleotide sequence variants | |
| US10937523B2 (en) | Methods, systems and computer readable storage media for generating accurate nucleotide sequences | |
| Zachariasen et al. | Identification of representative species-specific genes for abundance measurements | |
| CN114613434A (zh) | 基于群体样本深度信息检测基因拷贝数变异的方法及系统 | |
| CN110021342B (zh) | 用于加速变异位点的识别的方法及系统 | |
| CN113380324B (zh) | 一种T细胞受体序列motif组合识别检测方法、存储介质及设备 | |
| CN114067909B (zh) | 一种矫正同源重组缺陷评分的方法、装置和存储介质 | |
| US20170226588A1 (en) | Systems and methods for dna amplification with post-sequencing data filtering and cell isolation | |
| CN114242170B (zh) | 一种同源重组修复缺陷的评估方法、装置和存储介质 |
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 |