CN117079720B - 高通量测序数据的处理方法和装置 - Google Patents
高通量测序数据的处理方法和装置 Download PDFInfo
- Publication number
- CN117079720B CN117079720B CN202311335212.0A CN202311335212A CN117079720B CN 117079720 B CN117079720 B CN 117079720B CN 202311335212 A CN202311335212 A CN 202311335212A CN 117079720 B CN117079720 B CN 117079720B
- Authority
- CN
- China
- Prior art keywords
- mutation
- list
- information
- file
- merging
- 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
- 238000012165 high-throughput sequencing Methods 0.000 title claims abstract description 45
- 238000003672 processing method Methods 0.000 title claims abstract description 17
- 230000035772 mutation Effects 0.000 claims abstract description 444
- 238000000605 extraction Methods 0.000 claims abstract description 55
- 238000000034 method Methods 0.000 claims abstract description 43
- 235000019506 cigar Nutrition 0.000 claims description 87
- 108020004705 Codon Proteins 0.000 claims description 37
- 238000003780 insertion Methods 0.000 claims description 37
- 238000012937 correction Methods 0.000 claims description 36
- 230000037431 insertion Effects 0.000 claims description 36
- 238000001914 filtration Methods 0.000 claims description 24
- 238000012545 processing Methods 0.000 claims description 23
- 238000012217 deletion Methods 0.000 claims description 22
- 230000037430 deletion Effects 0.000 claims description 21
- 210000004602 germ cell Anatomy 0.000 claims description 18
- 238000012216 screening Methods 0.000 claims description 13
- 150000001413 amino acids Chemical class 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 4
- 241000703228 Heterocodon Species 0.000 claims description 3
- 230000011218 segmentation Effects 0.000 claims 1
- 238000012163 sequencing technique Methods 0.000 abstract description 25
- 238000001514 detection method Methods 0.000 abstract description 20
- 239000003814 drug Substances 0.000 abstract description 13
- 229940079593 drug Drugs 0.000 abstract description 10
- 238000012790 confirmation Methods 0.000 abstract description 7
- 239000000523 sample Substances 0.000 description 13
- 238000010586 diagram Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 7
- 230000003321 amplification Effects 0.000 description 5
- 239000013068 control sample Substances 0.000 description 5
- 210000000265 leukocyte Anatomy 0.000 description 5
- 238000003199 nucleic acid amplification method Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 206010028980 Neoplasm Diseases 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 239000012634 fragment Substances 0.000 description 4
- 239000002773 nucleotide Substances 0.000 description 4
- 125000003729 nucleotide group Chemical group 0.000 description 4
- 230000037429 base substitution Effects 0.000 description 3
- 238000004590 computer program Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 210000004881 tumor cell Anatomy 0.000 description 3
- WSFSSNUMVMOOMR-UHFFFAOYSA-N Formaldehyde Chemical compound O=C WSFSSNUMVMOOMR-UHFFFAOYSA-N 0.000 description 2
- 108091081062 Repeated sequence (DNA) Proteins 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 210000000349 chromosome Anatomy 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 230000000295 complement effect Effects 0.000 description 2
- 238000001647 drug administration Methods 0.000 description 2
- 102000052116 epidermal growth factor receptor activity proteins Human genes 0.000 description 2
- 108700015053 epidermal growth factor receptor activity proteins Proteins 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 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 2
- 108090000623 proteins and genes Proteins 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- QIVUCLWGARAQIO-OLIXTKCUSA-N (3s)-n-[(3s,5s,6r)-6-methyl-2-oxo-1-(2,2,2-trifluoroethyl)-5-(2,3,6-trifluorophenyl)piperidin-3-yl]-2-oxospiro[1h-pyrrolo[2,3-b]pyridine-3,6'-5,7-dihydrocyclopenta[b]pyridine]-3'-carboxamide Chemical compound C1([C@H]2[C@H](N(C(=O)[C@@H](NC(=O)C=3C=C4C[C@]5(CC4=NC=3)C3=CC=CN=C3NC5=O)C2)CC(F)(F)F)C)=C(F)C=CC(F)=C1F QIVUCLWGARAQIO-OLIXTKCUSA-N 0.000 description 1
- 102100021147 DNA mismatch repair protein Msh6 Human genes 0.000 description 1
- 101000968658 Homo sapiens DNA mismatch repair protein Msh6 Proteins 0.000 description 1
- 108091081021 Sense strand Proteins 0.000 description 1
- 108010078814 Tumor Suppressor Protein p53 Proteins 0.000 description 1
- 102000015098 Tumor Suppressor Protein p53 Human genes 0.000 description 1
- 230000000692 anti-sense effect Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000014509 gene expression Effects 0.000 description 1
- 238000010295 mobile communication Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000012188 paraffin wax Substances 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000013518 transcription Methods 0.000 description 1
- 230000035897 transcription Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
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
- G16B30/10—Sequence alignment; Homology search
-
- 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
- G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
- G16B15/30—Drug targeting using structural data; Docking or binding prediction
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Chemical & Material Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Analytical Chemistry (AREA)
- Medicinal Chemistry (AREA)
- Pharmacology & Pharmacy (AREA)
- Crystallography & Structural Chemistry (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种高通量测序数据处理方法和装置。该处理方法包括：对高通量测序数据进行比对，获得比对变异bam文件；对比对变异bam文件中的序列信息进行提取，得到序列提取文件；合并序列提取文件中的变异位点，得到变异位点合并文件；根据变异位点合并文件确定最终变异位点，输出变异结果文件。通过在变异检测确认之前，对比对获得的比对文件中存在变异信息的序列进行提取，并对其中存在重叠或临近的变异位点的位置信息进行合并，获得变异位点合并文件，利用该变异位点合并文件进行后续的变异位点确认，变异结果更准确，更接近真实突变，因而能够更准确地用于指导临床用药。
Description
技术领域
本发明涉及高通量测序数据处理领域，具体而言，涉及一种高通量测序数据处理的方法和装置。
背景技术
目前针对杂交捕获技术开发的第二代高通量测序数据突变检测的软件与方法有很多，但是准确性仍有提升空间，对于复杂区域（是指多种变异类型的区域）无法做到准确识别真实突变形式。且针对肿瘤临检FFPE(***固定石蜡包埋)类型样本，数据质量在无法保证的情况下更容易出现假阴性，假阳性等漏检或多检现象。在已有的近似方案中，北京雅康博生物科技有限公司针对多重扩增测序法为了改善假阳性位点多的缺陷，提供了一种高通量测序数据处理方法：包括获取二级测序序列，二级测序序列为高通量测序数据中能够被目的片段扩增引物识别，且去除对应的扩增引物后的测序序列；比对二级测序序列与参考基因组序列，得到初级变异结果；利用已知突变数据中的突变数据修正初级变异结果，得到处理结果。
现有技术存在如下几方面的缺陷：
1）针对杂交捕获技术开发的第二代高通量测序数据突变检测软件与方法多数仅能检出简单突变，涉及复杂突变等情况会出现报出检出位点不准确的情况，不精确的位点信息在涉及用药位点时会导致无法准确指导用药。
2）针对杂交捕获技术开发的第二代高通量测序数据突变检测软件与方法对indel(***/缺失)突变的检测也存在遗漏风险。
3）受到测序数据质量以及肿瘤细胞含量等影响，对于样本极为重要的突变位点也存在被参数阈值过滤，从而出现假阴性结果的风险。
由此可见，在可能存在复杂突变（即同时包括多种变异类型的突变）时，如何使检测得到的变异结果具有相对更高的准确性，目前尚无能够有效的解决方案。
发明内容
本发明的主要目的在于提供一种高通量测序数据处理方法和装置，以解决现有技术中变异分析结果准确性低的问题。
为了实现上述目的，根据本发明的一个方面，提供了一种高通量测序数据的处理方法，包括如下步骤：S101，对高通量测序数据进行比对，获得比对变异bam文件；S102，对比对变异bam文件中的序列信息进行提取，得到序列提取文件；S103，合并序列提取文件中的变异位点，得到变异位点合并文件；S104，根据变异位点合并文件确定最终变异位点，输出变异结果文件。
进一步地，S101包括：将高通量测序数据与参考基因组序列进行初次比对，获得初始比对bam文件；将初始比对bam文件与参考基因组进行局部重比对，获得比对变异bam文件。
进一步地，S102包括：从比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值；3）flag值；4）CIGAR标签；5）MD标签；6）起始位点；7）比对质量；根据信息提取文件对比对变异bam文件中的序列进行过滤，剔除如下至少之一的序列，获得序列提取文件：1）比对质量小于比对阈值的序列；2）在参考基因组上的比对位置不唯一的序列；3）重复序列。
进一步地，S103包括：根据信息提取文件中的CIGAR标签拆分CIGAR标签中的突变信息，输出带有第一突变信息的CIGAR变异列表，CIGAR变异列表中包含匹配M、缺失D、***I信息；根据信息提取文件中的MD标签拆分MD标签中的突变信息，输出带有第二突变信息的MD变异列表，MD变异列表中包含匹配M、缺失D和替换R信息；循环比对CIGAR变异列表和MD变异列表中的变异信息，对CIGAR变异列表和MD变异列表的变异及位置信息执行合并操作，得到合并列表；对合并列进行校正，得到变异位点合并文件。
进一步地，对CIGAR变异列表和MD变异列表的变异及位置信息执行合并操作，得到合并列表包括：如果CIGAR变异列表中存在软剪切信息，则去除软剪切信息后，输出调整匹配信息为变异起始位点，并输出至合并列表；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置相同，将将CIGAR变异列表的匹配M信息输出至合并列表；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至合并列表；如果CIGAR变异列表中的缺失D的位置同MD变异列表中的缺失D的位置相同，将CIGAR变异列表中的缺失D的信息输出至合并列表；如果CIGAR变异列表中的***I的位置同MD变异列表中的匹配M的位置重合，以***I的起始位置信息为分割，将在***I的位置***MD变异列表中的匹配M的位置中作为修正，输出至合并列表，并记录***长度；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的替换R的位置重合，以替换R的起始位置信息为分割，将替换R的位置***CIGAR变异列表中的匹配M的位置中作为修正，输出至合并列表。
进一步地，对合并序列进行校正，得到变异位点合并文件，包括：对合并列表中的缺失D信息进行校正；当合并列表中前一个突变信息为***I并且后一个突变信息为替换R时，根据序列信息对提取的CIGAR变异列表和MD变异列表进行整体校正；对合并列表中的匹配M的位置和长度进行验证，得到校正后的合并列表；根据校正后的合并列表，当突变信息符合如下至少之一的条件时，执行合并操作，获得变异位点合并文件：1）两个替换R的突变信息相距4bp以内时；2）***I的突变信息同其他任意突变信息相距20bp内时；3）缺失D的突变信息同其他任意突变信息相距20bp内时。
进一步地，S104包括：判断变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点；当存在情形1）时，则将分属不同密码子的SNV位点取消合并，删除不同密码子合并的位点，保留取消合并的不同密码子上单个SNV位点信息，形成第一更新后的变异位点合并文件；当存在情形2）时，则从第一更新后的变异位点合并文件中去除胚系变异位点，得到第二更新后的变异位点合并文件；利用第二更新变异位点合并文件进行突变重注释，获得重注释变异位点；对重注释变异位点进行位点筛选，获得最终变异位点，并输出变异结果文件；其中，位点筛选包括：1）去除比对到参考基因组的正链与负链的序列的比例大于等于90%的序列支持的位点；2）除去低于过滤阈值的位点；3）从低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变；变异结果文件中包括热点突变和非热点突变，热点突变包括待验证的低频突变，待验证的低频突变指低于所述过滤阈值的热点突变。
进一步地，判断变异位点合并文件是否存在情形1）包括：对变异位点合并文件中的突变信息进行注释，得到变异注释文件；根据变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子。
为了实现上述目的，根据本发明的第二个方面，提供了一种高通量测序数据的处理装置，包括：比对单元，被设置为对高通量测序数据进行比对，获得比对变异bam文件；序列提取单元，被设置为对比对变异bam文件中的序列信息进行提取，得到序列提取文件；变异位点合并单元，被设置为合并序列提取文件中的变异位点，得到变异位点合并文件；变异确定单元，被设置为根据变异位点合并文件确定最终变异位点，输出变异结果文件。
进一步地，比对单元包括：第一比对模块，被设置为将高通量测序数据与参考基因组序列进行初次比对，获得初始比对bam文件；第二比对模块，被设置将初始比对bam文件与参考基因组进行局部重比对，获得比对变异bam文件。
进一步地，序列提取单元包括：提取模块，被设置为从比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值；3）flag值；4）CIGAR标签；5）MD标签；6）起始位点；7）比对质量；筛选单元，被设置为根据信息提取文件对比对变异bam文件中的序列进行筛选，剔除如下至少之一的序列，获得序列提取文件：1）比对质量小于比对阈值的序列；2）在参考基因组上的比对位置不唯一的序列；3）重复序列。
进一步地，变异位点合并单元包括：CIGAR变异输出模块，被设置为根据信息提取文件中的CIGAR标签拆分CIGAR标签中的突变信息，输出带有第一突变信息的CIGAR变异列表，CIGAR变异列表中包含匹配M、缺失D、***I信息；MD变异输出模块，被设置为根据信息提取文件中的MD标签拆分MD标签中的突变信息，输出带有第二突变信息的MD变异列表，MD变异列表中包含匹配M、缺失D和替换R信息；合并模块，被设置为循环比对CIGAR变异列表和MD变异列表中的变异信息，对CIGAR变异列表和MD变异列表的变异及位置信息执行合并操作，得到合并列表；校正模块，被设置为对合并列进行校正，得到变异位点合并文件。
进一步地，合并模块包括：第一判断输出子模块，被设置为在CIGAR变异列表中存在软剪切信息时，去除软剪切信息后，调整匹配信息为变异起始位点，并输出至合并列表；第二判断输出子模块，被设置为在CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置相同，将将CIGAR变异列表的匹配M信息输出至合并列表；第三判断输出子模块，被设置为在CIGAR变异列表的匹配M的位置同MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至合并列表；第四判断输出子模块，被设置为在CIGAR变异列表中的缺失D的位置同MD变异列表中的缺失D的位置相同时，将CIGAR变异列表中的缺失D的信息输出至合并列表；第五判断输出子模块，被设置为在CIGAR变异列表中的***I的位置同MD变异列表中的匹配M的位置重合，以***I的起始位置信息为分割，将在***I的位置***MD变异列表中的匹配M的位置中作为修正，输出至合并列表，并记录***长度；第六判断输出模块，被设置为在CIGAR变异列表中的匹配M的位置同MD变异列表中的替换R的位置重合，以替换R的起始位置信息为分割，将替换R的位置***CIGAR变异列表中的匹配M的位置中作为修正，输出至合并列表。
进一步地，校正模块包括：第一校正子模块，被设置为对合并列表中的缺失D信息进行校正；第二校正子模块，被设置为当合并列表中前一个突变信息为***I并且后一个突变信息为替换R时，根据对序列信息提取的CIGAR变异列表和MD变异列表整体校正；第三校正子模块，被设置为对合并列表中的匹配M的位置和长度进行验证，得到校正后的合并列表；合并子模块，被设置为根据校正后的合并列表，当突变信息符合如下至少之一的条件时，执行合并操作，获得变异位点合并文件：1）两个替换R的突变信息相距4bp以内时；2）***I的突变信息同其他任意突变信息相距20bp内时；3）缺失D的突变信息同其他任意突变信息相距20bp内时。
进一步地，变异确定单元包括：情形判断模块，被设置为判断变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点；第一更新模块，被设置为在存在情形1）时，则将分属不同密码子的SNV位点取消合并，删除不同密码子合并的位点，保留取消合并的不同密码子上单个SNV位点信息，形成第一更新后的变异位点合并文件；第二更新模块，被设置为在存在情形2）时，则从第一更新后的变异位点合并文件中去除胚系变异位点，得到第二更新后的变异位点合并文件；重注释模块，被设置为利用第二更新变异位点合并文件进行突变重注释，获得重注释变异位点；位点筛选及输出模块，被设置为对重注释变异位点进行位点筛选，获得最终变异位点，并将输出变异结果文件；其中，位点筛选包括：1）去除比对到参考基因组的正链与负链的序列的比例大于等于90%的序列支持的位点；2）除去低于过滤阈值的位点；3）从低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变；变异结果文件中包括热点突变和非热点突变，热点突变包括待验证的低频突变，待验证的低频突变指低于所述过滤阈值的热点突变。
进一步地，情形判断模块包括：情形1）判断子模块和情形2）判断子模块，其中，情形1）判断子模块包括：注释元件，被设置为对变异位点合并文件中的突变信息进行注释，得到变异注释文件；异密码子判定模块，被设置为根据变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子。
根据本发明的另一方面，提供了一种计算机可读存储介质，计算机可读存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行上述任一种高通量测序数据的处理方法。
根据本发明的另一方面，提供了一种处理器，处理器用于运行程序，其中，程序运行时执行上述任一种高通量测序数据的处理方法。
应用本发明的技术方案，通过在变异检测确认之前，对比对获得的比对文件中存在变异信息的序列进行提取，并对其中存在重叠或临近的变异位点的位置信息进行合并，获得变异位点合并文件，合并后的变异位点避免了将同一密码子内的变异分属于单个突变位点，而难以准确指导临床用药的情况。利用该变异位点合并文件进行后续的变异位点确认，变异结果更准确，更接近真实突变，因而能够更准确地用于指导临床用药。
附图说明
构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1示出了根据本发明的实施例1的高通量测序数据处理方法的流程示意图；
图2示出了根据本发明的实施例3的高通量测序数据处理装置的结构示意图。
图3示出了根据本发明的一种具体实施例中测序数据与参考基因组序列比对后得到的BAM文件中的读段（reads）信息示意图。
图4示出了上述图3所示的实施例中CIGAR标签和MD标签两个标签的变异信息示意图。
图5示出了上述图3所示的实施例中两个标签的变异信息的合并过程示意图。
图6示出了上述图5的合并结果示意图。
图7示出了上述图6的合并结果的展示示意图。
图8示出了根据本发明的一种具体实施例中对于分属不同密码子的SNV位点取消合并的具体操作示意图。
图9示出了根据本发明的一种具体实施例中对重注释中3’端右移进行说明的示意图。
图10至图16分别示出了根据本发明的实施例2中各示例的结果图。
图17示出了本发明的一个优选的实施例中所提供的一种高通量测序数据处理方法的终端的硬件结构框图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合实施例来详细说明本发明。
术语解释：
bam文件说明：bam文件是一种二进制的压缩文件，需要通过特定的软件来进行查看，bam文件通常由12个字段组成。bam格式分为header section（头部分，注释信息，以@开头，可有可无）和alignment section（比对结果）两个部分。alignment section由11个字段组成：
第一列：序列的名字，也就是reads的名称；
第二列：FLAG，每一个read的比对情况可以用十进制数字（或者十六进制数字）表示，如果比对情况有多个，将多个比对情况所代表的十进制数字加和就是这一行的FLAG。以下为十进制数字代表含义：
1：该read是成对的paired reads中的一个。
2：paired reads中每个都正确比对到参考序列上。
4：该read没比对到参考序列上。
8：与该read成对的matepair read没有比对到参考序列上。
16：该read其反向互补序列能够比对到参考序列。
32：与该read成对的matepair read其反向互补序列能够比对到参考序列。
64：在paired reads中，该read是与参考序列比对的第一条。
128：在paired reads中，该read是与参考序列比对的第二条。
256：该read是次优的比对结果。
512：该read没有通过质量控制。
1024：由于PCR或测序错误产生的重复reads。
2048：补充匹配的read。
第三列：参考序列的名字；比对上的参考基因组序列的名字，通常为染色体名称。
第四列：在参考序列上的位置：read比对到的参考序列最左侧的位置坐标，也是CIGAR中第一个比对标识“M（Match，匹配）”对应的最左侧碱基在参考序列的位置。
第五列：比对质量（mapping quality），计算方法为比对错误率的-10*log10的值，一般是四舍五入的整数值，如果是255，说明该比对值无效。
第六列：代表比对结果的CIGAR字符串，也叫CIGAR标签，表示比对的具体情况。CIGAR标识符表示read中每个碱基的比对情况，主要有以下标识符：
M: alignment match，read上的碱基与参考序列完全匹配，碱基一一对应，同时囊括正确匹配碱基与错误匹配碱基。
I: insertion to the reference，read上的碱基相对于参考序列存在***碱基。
D: deletion from the reference，read上的碱基相对于参考序列存在删除碱基。
S: soft clipping，软剪切，序列信息在bam中仍然完整保留，不会对序列进行修改，表示软剪切的部分序列无法与read的其他部分序列在参考序列中比对至同一位置。
H: hard clipping，硬剪切，read的开头或者结尾部分没有比对到参考序列上，在bam文件中不会保留未比对的部分序列。
举例：如15S37M1D2M1I，这段字符的意思是开头为15个碱基无法与后续的40个碱基比对至同一位置，前15个碱基为软剪切，在该比对位置无可匹配的参考序列，37个碱基与参考基因组完全匹配，1个参考序列上的删除(缺少一个参考基因组碱基)，2个匹配，1个参考序列上的***(比参考基因组多一个碱基)。
第七列：mate read，该read的mate read比对上的参考序列的名字。
如果和该read所在行的第三列参考序列名字一样，则用“=”表示，说明这对read比对到了同一条参考序列上；
如果mate read没有比对可比对的参考序列，第七列则用“*”表示；
如果这对read没有比对到同一条参考序列，那么这一列则是mate read所能比对的参考序列名称。
第八列：mate 序列在参考序列上的位置；
第九列：测序片段长度，可以简单理解为文库片段长度；
第十列：read的序列信息；
第十一列：read序列对应的ASCII码格式的碱基质量值；
第十二列以及后续其他列：为可选的区域 header section，包括MC、MD、NM、AS、XS等标签。其中，MD标签中，“^”代表删除（Deletion）变异；“A/T/C/G”代表单核苷酸变异，即碱基替换；***数字代表匹配参考基因组碱基数量。MD标签中不存在***（Insertion）变异。
如背景技术所提到的，现有的高通量测序数据进行肿瘤变异检测分析的方法存在无法准确识别复杂突变（即同时包含***、缺失或替换等变异类型的突变），而需要人工确认来提高检测准确性的问题，为改善这一状况，在本申请中发明人经分析发现，现有的高通量测序数据处理方法通常情况下无法对相邻突变进行准确合并报出，而是单独报出各个突变位点，这使得致报出的位点无法准确指导用药。基于这一发现，发明人提出了对相邻突变进行准确合并的改进思路，并在改进思路下，提出了本申请的一系列保护方案。
实施例1
在本申请第一种典型的实施例中，提供了一种高通量测序数据的处理方法，该处理方法如图1所示，包括如下步骤：
S101，对高通量测序数据进行比对，获得比对变异bam文件；
S102，对比对变异bam文件中的序列信息进行提取，得到序列提取文件；
S103，合并序列提取文件中的变异位点，得到变异位点合并文件；
S104，根据变异位点合并文件确定最终变异位点，输出变异结果文件。
该处理方法实质上是一种自动化筛选变异位点的处理方法，该方法通过在变异检测确认之前，对比对获得的比对文件中存在变异信息的序列进行提取，并对其中存在重叠或临近的变异位点的位置信息进行合并，获得变异位点合并文件，合并后的变异位点避免了将同一密码子内的变异分属于单个突变位点，难以准确表示真实突变形式，进而无法指导临床用药的情况以及存在假阴性和假阳性位点较多的情况。利用本申请中的该变异位点合并文件进行后续的变异位点确认，变异结果更准确，更接近真实突变，因而能够更准确地用于指导临床用药。
上述S101是高通量测序数据与参考基因组比对，获取比对变异bam文件的步骤，具体的比对方式采用本领域的常规比对操作即可。比如，先进行全基因组比对，再进行局部重比对。在一种优选的实施例中，上述S101包括：将高通量测序数据与参考基因组序列进行初次比对，获得初始比对bam文件；将初始比对bam文件与参考基因组进行局部重比对，获得比对变异bam文件。局部重比对后经重新排序（常规操作，后续所有数据的处理均需要bam排序建立索引后才能使用软件快速定位到指定位置）生成索引文件。
需要说明的是，本申请的高通量测序数据可适配于不同类型高通量测序数据，且不区分测序探针是否为panel还是其他WES等不同探针。具体地，包括但不限于针对已知突变而设计的panel产出的高通量测序数据，也适用于基于未知突变无法准确分型且不适用于检测范围更大的其他类型探针产出的高通量测序数据。
上述步骤S102的目的为后续可能会出现的需要合并变异位点的情形进行前期准备，通过对比对变异bam文件中的序列信息进行提取，根据提取变异信息进行合并，能够使合并的准确性相对较高，进而有利于后续最终变异位点报出的准确性。具体提取与变异有关的序列信息可以从bam文件中合理选择得到。
在一种优选的实施例中，上述步骤S102包括：从比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值（第11列）；3）flag值（第2列）；4）CIGAR标签（第6列）；5）MD标签（第12列）；6）起始位点（第4列）；7）比对质量（第5列）；根据信息提取文件对比对变异bam文件中的序列进行筛选，剔除如下至少之一的序列，获得序列提取文件：1）比对质量小于比对阈值的序列，优选阈值小于60；2）在参考基因组上的比对位置不唯一的序列（可根据第2列的FLAG的数字来确定）；3）重复序列（根据第2列的FLAG值确定）（具体地，可以根据flag信息，使用pysam程序包中的is_reverse功能来实现）。
上述S103的主要目的是对相邻突变位点进行合并，以避免最终分析报出的变异位点为单独位点，而难以确定真实的突变形式，进而难以指导临床用药的情况。因而，在进行上述S102的变异序列信息提取之后，根据提取所得文件中的变异信息指标，可选择合适的方式来确定具体需合并的变异序列列表。
在一种优选的实施例中，S103包括：根据信息提取文件中的CIGAR标签拆分CIGAR标签中的突变信息，CIGAR变异列表中包含Match(匹配) ，Delete(缺失) 和Insert(***)信息，输出带有第一突变信息的CIGAR变异列表，根据信息提取文件中的MD标签拆分MD标签中的突变信息，MD变异列表中包含Match(匹配)，Delete(缺失)和替换（Replace）信息；输出带有第二突变信息的MD变异列表，循环比对CIGAR变异列表和MD变异列表中的变异信息，对CIGAR变异列表和MD变异列表中的变异及位置信息执行合并操作，得到合并列表，对合并列进行校正，得到变异位点合并文件。
上述优选实施例根据前述提取文件中提取的CIGAR标签和MD标签，其中CIGAR标签包涵了M-匹配，I-***，D-缺失，S-单核苷酸变异（又记为SV）等与参考基因组比对上（比对上也叫匹配）及比对不上（也叫不匹配，意味着存在变异）的变异类型及具体碱基数量。MD标签中，匹配M、缺失D和替换R信息。具体地，“^”代表删除（Deletion）变异，“A/T/C/G”代表单核苷酸变异，***数字代表匹配参考基因组碱基数量。
具体地，结合以下示例对上述根据两种标签对变异信息进行拆分和合并的过程进行详细展示和说明。图3所示的是测序数据与参考基因组序列比对后得到的BAM文件中的读段（reads）信息。其中CIGAR标签所标记的变异列表中包含匹配M、缺失D和***I等变异信息。MD变异列表中包含匹配M、缺失D和替换R等变异信息（其具体展现形式与CIGAR变异列表有所不同）。在图3中，两个标签列表分别如下：
CIGAR：67M4D8M4I71M 表示67个匹配，4个缺失，8个匹配，4个***及71个匹配。
MD：67^CCTC4T0C73
67表示67个匹配，^CCTC 表示4个缺失，4表示4个匹配，T表示1个碱基替换，0表示0个匹配，C表示1个碱基替换，73表示73个匹配。
因此，上述两个标签的变异信息图示例见图4（该段序列比对的是基因组反义链，而描述突变信息通常使用的都是正义链碱基。因此显示的是互补配对碱基）。合并过程见图5。合并结果见图6，结果展示见图7。
为了进一步使合并后的结果更准确，在上述合并操作的过程中，根据不同的情况，合并操作也存在差异。在一种优选的实施例中，对CIGAR变异列表和MD变异列表执行变异位点合并操作，得到合并列表包括：如果所述CIGAR变异列表中存在软剪切信息，则去除所述软剪切信息后，输出调整匹配信息为变异起始位点，并输出至合并列表；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置相同，将CIGAR变异列表的匹配M信息输出至合并列表；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至合并列表；如果CIGAR变异列表中的缺失D的位置同MD变异列表中的缺失D的位置相同，将CIGAR变异列表中的缺失D的信息输出至合并列表；如果CIGAR变异列表中的***I的位置同MD变异列表中的匹配M的位置重合，以***I的起始位置信息为分割，将在***I的位置***MD变异列表中的匹配M的位置中作为修正，输出至合并列表，并记录***长度；如果CIGAR变异列表中的匹配M的位置同MD变异列表中的替换R的位置重合，以替换R的起始位置信息为分割，将替换R的位置***CIGAR变异列表中的匹配M的位置中作为修正，输出至合并列表（参见图5-7）。
为了进一步提高合并后的序列的准确性，本申请的处理方法还包括对合并序列进行校正的步骤。在一种优选的实施例中，对合并序列进行校正，得到变异位点合并文件包括：对合并列表中缺失D信息校正；当合并列表中前一个突变信息为***I并且后一个突变信息为替换R时，根据对序列信息对提取的CIGAR变异列表和MD变异列表整体校正（目的是防止异常情况导致信息失准）；对合并列表中的匹配M的位置和长度进行验证，得到校正后的合并列表；根据校正后的合并列表，当突变信息符合如下至少之一的条件时，执行合并操作，获得变异位点合并文件：1）两个Replace(碱基替换)突变信息相距4bp以内时；2）Inserttion(***)突变信息同其他任意突变信息相距20bp内时；3）Delete(缺失)突变信息同其他任意突变信息相距20bp内时。
在一种优选的实施例中，上述S104包括：判断变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点（具体参照本领域常规的胚系变异的鉴定方法进行判定即可，比如，如果为配对样本流程，读取白细胞对照样本的bam文件，根据提取的白细胞对照样本序列信息提取白细胞对照样本的CIGAR标签和MD标签信息，根据CIGAR变异列表和MD变异列表信息执行变异列表的合并和重定位，合格的（符合质量阈值判断的归为合格）突变信息列表输出至白细胞对照样本的变异统计表中，根据白细胞对照样本的变异统计表，删除上一步中变异统计表中的胚系位点（除了胚系位点，也可能是胚系中也同时存在的伪影或者测序错误等情况）；当存在情形1）时，则将分属不同密码子的SNV位点取消合并，删除不同密码子合并的位点，保留取消合并的不同密码子上单个SNV位点信息，形成第一更新后的变异位点合并文件；当存在情形2）时，则从第一更新后的变异位点合并文件中去除胚系变异位点，得到第二更新后的变异位点合并文件；当存在情形1）和情形2），利用第二更新变异位点合并文件进行突变注释，获得最终变异位点，并将输出变异位点结果文件进行突变注释，获得重注释变异位点；对重注释变异位点进行位点筛选（包括1）去除比对到参考基因组的正链的序列与比对到参考基因组的负链的序列的比例大于等于90%的位点，也即，覆盖某位点的双端测序的读段reads中90%比对到了参考基因组的其中一个方向的链上，该位点需要去除；2）去除低于过滤阈值的位点；3）从低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变），获得最终变异位点，并输出变异结果文件，其中，变异结果文件中包括热点突变和非热点突变，热点突变包括待验证的低频突变位点，待验证低频的突变位点指低于过滤阈值的突变位点，从结果中输出用于后续进行验证。
需要说明的是，在上述对重注释变异位点进行位点筛选（即变异信息确认，保留预期内突变信息）过程中，所涉及的过滤阈值的参数可根据具体测序数据的不同而有所不同（比如，样本类型、样本来源、测序深度等等）。具体地，可根据不同样本类型以及不同需求自定义不同阈值标准。在实际的过滤阈值的参数的设定时，根据不同需求设定reads数与突变丰度的过滤阈值（在本申请的一些优选实施例中，默认的过滤阈值为7条reads支持+1%丰度），同时针对可指导用药的热点突变（即待验证的低频突变）设定更低阈值（比如为3条或4条reads支持+1%丰度），从而避免热点突变遗漏。
对于分属不同密码子的SNV位点取消合并的具体操作示例如图8所示，TP53存在两个相邻的SNV，分别从CC突变成了AT，之前基于相邻将其合并，在该步骤中，基于是否属于不同的密码子内，再将其拆分开来，最后以单独的两个SNV报出结果。具体见下表：
表1：
需要说明的是，热点突变在不同的研究或场景中存在差异，具体可以根据实际需求进行自定义。本申请优先指可以直接指导A级用药的位点。非热点突变指B/C级用药位点。同样地，在不同的场景下，上述报出阈值的具体取值也不同。
上述重注释的过程中，需要用到基因的具体信息，例如，转录本、外显子区域及CDS区域的具体信息，需要同时判断比对序列位置是否为重复序列，如果是重复序列，需要根据基因转录方向，确认重复的***/缺失需要移动的方向。同时根据转录本以及CDS区域具体坐标信息，确认移动后的***/缺失突变注释。
在一种优选的实施例中，判断变异位点合并文件是否存在情形1）包括：对变异位点合并文件中的突变信息进行注释，得到变异注释文件；根据变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子（参见图8）。
该优选实施例通过先进行一次注释，注释后即可确定不同的变异位点属于同一个密码子内，还是分属不同的密码子，这样能更进一步使得最终变异注释的结果更准确，也更能准确指示变异的类型，进而能更准确地用于指导临床用药。
其中的注释通过调用已有的注释模块来进行注释即可。具体包括对突变信息进行基础注释、dup突变信息注释的格式化以及基于突变位点根据基因表达方向进行3’端右移等校正操作，最后输出注释后信息。
下面结合MSH6位点对上述重注释中3’端右移进行示例性说明，如图9及下表所示。
表2：
实施例2
本实施例通过测试真实2067例样本，其中携带复杂突变样本376例，携带低频用药位点突变样本245例。经人工校对，复杂突变准确合并报出比例达100%，低频用药位点突变防遗漏提示比例100%。
对部分复杂突变自动合并符合结论的测序数据的测试结果示例如下：
示例-1 ：snv与indel自动合并注释（见图10及下表）
表3：
示例-2：相同密码子内snv自动合并注释（见图11及下表）
表4：
。
示例-3：多个snv，indel自动合并注释（见图12及下表）
表5：
。
需要说明的是，本实施例中所有示例的表中，每一列的信息均分别为：1. 样本示例；2. 突变位点信息；3. 染色体(chr)；4. 突变起始位置(start)；5. 突变终止位置(end)；6. 野生型(ref)；7. 突变型(alt) ；8. 注释信息(annotation)。下面各表不再解释。
测试结论：自动合并功能完善，重新注释位点符合HGVS国际标准。
示例-4：对部分低频突变测试并进行3D PCR验证结果如下：
（1）KRAS低频，仅4条reads支持该突变，经3D验证为真，见图13及下表。
表6：
（2）KRAS低频，仅4条reads支持该突变，经3D验证为假，见图14及下表。
表7：
（3）EGFR低频，仅3条reads支持该突变，经3D验证为真，见图15及下表。
表8：
（4）EGFR低频，有6条reads支持该突变，经3D验证为假，见图16及下表。
表9：
/>
测试结论：由于肿瘤细胞含量低以及数据质量差等因素影响，导致部分可指导用药位点存在支持读段（reads）少，突变丰度低的情况，仅仅通过阈值进行筛选均可能存在假阳性与假阴性的情况，为临床检测样本负责，数据分析应全面提示可能存在的有意义低频位点，经其他方法进行验证，防止出现假阳性与假阴性的情况。
需要说明的是，本申请的高通量测序数据的处理方法能够自动合并同读段（reads）上相同密码子内点突变，自动合并带有indel(***/缺失)突变的复杂突变，重注释自动进行3’端右移，区分热点与非热点突变的不同位点质量评估，提示低频热点突变。
现有技术受限于检测范围内的已知突变且仅适用于多重扩增测序法，针对未知突变无法赋予准确分型且不适用于检测范围更大的其他类型panel。本申请的方法可适配于不同类型高通量测序数据，且不区分测序探针是否为panel还是其他WES等不同探针。对于数据中存在的复杂突变都能做到自行合并，以符合HGVS的国际标准报出准确突变信息，同步提示低频用药位点，防止假阴性遗漏。
因此，本申请的方法可进一步扩展至所有使用杂交捕获测序技术的SNV(点突变)/InDel(小片段***缺失)变异检测，调参测试后也可适用于多重扩增测序法数据，适用于需要准确提供突变信息的临检方向。无需人工重新校对真实突变形式，报出突变符合HGVS等国际标准。
需要说明的是，对于前述的各方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明并不受所描述的动作顺序的限制，因为依据本发明，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明所必须的。
通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的检测仪器等硬件设备的方式来实现。基于这样的理解，本申请的技术方案中数据处理的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本申请各个实施例或者实施例的某些部分的方法。
本申请可用于众多通用或专用的计算***环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器***、基于微处理器的***、置顶盒、可编程的消费电子设备、网络PC、小型计算机、大型计算机、包括以上任何***或设备的分布式计算环境等等。
本申请所提供的方法可以在终端、计算机终端或者类似的运算装置中执行。以运行在终端上为例，图17是本发明实施例的一种方法的终端的硬件结构框图。如图17所示，终端可以包括一个或多个（图17中仅示出一个）处理器102（处理器102可以包括但不限于微处理器MCU或可编程逻辑器件FPGA等的处理装置）和用于存储数据的存储器104，可选地，上述终端还可以包括用于通信功能的传输设备106以及输入输出设备108。本领域普通技术人员可以理解，图17所示的结构仅为示意，其并不对上述终端的结构造成限定。例如，终端还可包括比图17中所示更多或者更少的组件，或者具有与图4所示不同的配置。
存储器104可用于存储计算机程序，例如，应用软件的软件程序以及模块，如本发明实施例中的读段拼接、分簇、一致性处理等方法对应的计算机程序，处理器102通过运行存储在存储器104内的计算机程序，从而执行各种功能应用以及数据处理，即实现上述的方法。存储器104可包括高速随机存储器，还可包括非易失性存储器，如一个或者多个磁性存储装置、闪存、或者其他非易失性固态存储器。在一些实例中，存储器104可进一步包括相对于处理器102远程设置的存储器，这些远程存储器可以通过网络连接至终端。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。
传输设备106用于经由一个网络接收或者发送数据。上述的网络具体实例可包括终端的通信供应商提供的无线网络。在一个实例中，传输设备106包括一个网络适配器（Network Interface Controller，简称为NIC），其可通过基站与其他网络设备相连从而可与互联网进行通讯。在一个实例中，传输设备106可以为射频（Radio Frequency，简称为RF）模块，其用于通过无线方式与互联网进行通讯。
显然，本领域的技术人员应该明白，上述的本申请的部分模块或步骤可以在通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本申请不限制于任何特定的硬件和软件结合。
实施例3
本实施例提供了一种高通量测序数据处理装置，如图2所示，该装置包括：比对单元10、序列提取单元20、变异位点合并单元30及变异确定单元40，其中，
比对单元10被设置为对高通量测序数据进行比对，获得比对变异bam文件；
序列提取单元20被设置为对比对变异bam文件中的序列信息进行提取，得到序列提取文件；
变异位点合并单元30被设置为合并序列提取文件中的变异位点，得到变异位点合并文件；
变异确定单元40被设置为根据变异位点合并文件确定最终变异位点，输出变异结果文件。
进一步地，比对单元包括：第一比对模块，被设置为将高通量测序数据与参考基因组序列进行初次比对，获得初始比对bam文件；第二比对模块，被设置将初始比对bam文件与参考基因组进行局部重比对，获得比对变异bam文件。
进一步地，序列提取单元包括：提取模块，被设置为从比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值；3）flag值；4）CIGAR标签；5）MD标签；6）起始位点；7）比对质量；筛选单元，被设置为根据信息提取文件对比对变异bam文件中的序列进行筛选，剔除如下至少之一的序列，获得序列提取文件：1）比对质量小于比对阈值的序列，优选阈值小于60；2）在参考基因组上的比对位置不唯一的序列；3）重复序列。
进一步地，变异位点合并单元包括：CIGAR变异输出模块，被设置为根据信息提取文件中的CIGAR标签拆分CIGAR标签中的突变信息，“D”代表Deletion缺失；“I”代表Insertion***；“S”代表软剪切（softclip）序列；“M”代表匹配参考基因组，输出带有第一突变信息的CIGAR变异列表；MD变异输出模块，被设置为根据信息提取文件中的MD标签拆分MD标签中的突变信息，其中，“^”代表Deletion缺失；“A/T/C/G”代表单核苷酸变异（即替换Replace）；***数字代表匹配参考基因组碱基数量，输出带有第二突变信息的MD变异列表；合并模块，被设置为循环比对CIGAR变异列表和MD变异列表中的变异信息，对CIGAR变异列表和MD变异列表的变异及位置信息执行变异位点合并操作，得到合并列表；校正模块，被设置为对合并列进行校正，得到变异位点合并文件。
进一步地，合并模块包括：第一判断输出子模块，被设置为在CIGAR变异列表中存在软剪切信息时，去除软剪切信息后，调整匹配信息为变异起始位点，并输出至合并列表；第二判断输出子模块，被设置为在CIGAR变异列表中的匹配M的位置同MD变异列表中的匹配M的位置相同，将将CIGAR变异列表的匹配M信息输出至合并列表；第三判断输出子模块，被设置为在CIGAR变异列表的匹配M的位置同MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至合并列表；第四判断输出子模块，被设置为在CIGAR变异列表中的缺失D的位置同MD变异列表中的缺失D的位置相同时，将CIGAR变异列表中的缺失D的信息输出至合并列表；第五判断输出子模块，被设置为在CIGAR变异列表中的***I的位置同MD变异列表中的匹配M的位置重合，以***I的起始位置信息为分割，将在***I的位置***MD变异列表中的匹配M的位置中作为修正，输出至合并列表，并记录***长度；第六判断输出模块，被设置为在CIGAR变异列表中的匹配M的位置同MD变异列表中的替换R的位置重合，以替换R的起始位置信息为分割，将替换R的位置***CIGAR变异列表中的匹配M的位置中作为修正，输出至合并列表。
进一步地，校正模块包括：第一校正子模块，被设置为采用循环验证的方法对合并列表中的缺失Delete信息校正；第二校正子模块，被设置为当合并列表中前一个为Insert(***)并且后一个为Replace(碱基替换)时，根据对序列信息提取的CIGAR变异列表和MD变异列表整体校正（防止异常情况导致信息失准）；第三校正子模块，被设置为对合并列表中的Match(匹配)信息位置和长度进行验证，得到校正后的合并列表；合并子模块，被设置为根据校正后的合并列表，当突变信息符合如下至少之一的条件时，执行变异信息合并操作，获得变异位点合并文件：1）两个Replace(碱基替换)突变信息相距4bp以内时；2）Insert(***)突变信息同其他任意突变信息相距20bp内时；3）Delete(缺失)突变信息同其他任意突变信息相距20bp内时。
进一步地，变异确定单元包括：情形判断模块，被设置为判断变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点；第一更新模块，被设置为在存在情形1）时，则将分属不同密码子的SNV位点拆分开，并根据拆分后开的突变信息与变异位点合并文件中的信息进行比对，若存在匹配的，则进行加和；若存在不匹配的，则将不匹配的并入变异位点合并文件中，形成第一更新后的变异位点合并文件；第二更新模块，被设置为在存在情形2）时，则从第一更新后的变异位点合并文件中去除胚系变异位点，得到第二更新后的变异位点合并文件；重注释输出模块，被设置为利用第二更新后的变异位点合并文件进行突变重注释，获得重注释变异位点；位点筛选及输出模块，被设置为对述重注释变异位点进行位点筛选（包括1）去除比对到参考基因组的正链的序列与比对到参考基因组的负链的序列的比例大于等于90%的位点，也即，覆盖某位点的双端测序的读段reads中90%比对到了参考基因组的其中一个方向的链上，该位点需要去除；2）去除低于过滤阈值的位点；3）从低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变），获得最终变异位点，并将输出变异结果文件；其中，变异结果文件中包括热点突变和非热点突变最终变异位点，热点突变包括待验证的低频突变位点，待验证的低频突变位点指低于报出阈值的突变位点。
进一步地，情形判断模块包括：情形1）判断子模块和情形2）判断子模块，其中，情形1）判断子模块包括：注释元件，被设置为对变异位点合并文件中的突变信息进行注释，得到变异注释文件；异密码子判定模块，被设置为根据变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子。
实施例4
本实施例提供了一种计算机可读存储介质，计算机可读存储介质包括存储的程序，其中，在程序运行时控制存储介质所在设备执行前述高通量测序数据处理方法。
本实施例还提供了一种处理器，处理器用于运行程序，其中，程序运行时执行前述高通量测序数据处理方法。
从以上的描述中，可以看出，本发明上述的实施例实现了如下技术效果：目前基于第二代高通量测序数据进行肿瘤变异检测分析的方法存在无法准确识别复杂突变的困难，在得到初步的变异检测结果后仍需要人工确认真实突变形式，本发明的方法和装置可准确提供复杂突变的真实突变形式。在肿瘤临床检测中的可指导用药的重要突变位点上，由于测序质量以及肿瘤细胞含量等影响因素都会对变异检测准确率造测成极大影响，因此本发明的方法和装置在矫正复杂突变同时，基于指定参数重新定义突变热点检出下限，减少突变检测假阴性情况的发生。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (12)
1.一种高通量测序数据的处理方法，其特征在于，包括如下步骤：
S101，将高通量测序数据与参考基因组进行比对，获得比对变异bam文件；
S102，对所述比对变异bam文件中的序列信息进行提取，得到序列提取文件；
S103，合并所述序列提取文件中的变异位点，得到变异位点合并文件；
S104，根据所述变异位点合并文件确定最终变异位点，输出变异结果文件；
所述S102包括：
从所述比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值；3）flag值；4）CIGAR标签；5）MD标签；6）起始位点；7）比对质量；
根据所述信息提取文件对所述比对变异bam文件中的序列进行过滤，剔除如下至少之一的序列，获得所述序列提取文件：1）比对质量小于比对阈值的序列；2）在所述参考基因组上的比对位置不唯一的序列；3）重复序列；
所述S103包括：
根据所述信息提取文件中的所述CIGAR标签拆分CIGAR标签中的突变信息，输出带有第一突变信息的CIGAR变异列表，所述CIGAR变异列表中包含匹配M、缺失D、***I信息；
根据所述信息提取文件中的所述MD标签拆分MD标签中的突变信息，输出带有第二突变信息的MD变异列表，所述MD变异列表中包含匹配M、缺失D和替换R信息；
循环比对所述CIGAR变异列表和所述MD变异列表中的变异信息，对所述CIGAR变异列表和所述MD变异列表的变异及位置信息执行合并操作，得到合并列表；
对所述合并列表进行校正，得到所述变异位点合并文件；
所述S104包括：
判断所述变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点；
当存在所述情形1）时，则将分属不同密码子的SNV位点取消合并，删除不同密码子合并的位点，保留取消合并的不同密码子上单个SNV位点信息，形成第一更新后的所述变异位点合并文件；
当存在所述情形2）时，则从第一更新后的所述变异位点合并文件中去除所述胚系变异位点，得到第二更新后的所述变异位点合并文件；
利用所述第二更新后的所述变异位点合并文件进行突变重注释，获得重注释变异位点；
对所述重注释变异位点进行位点筛选，获得所述最终变异位点，并输出所述变异结果文件；
其中，所述位点筛选包括：1）去除比对到所述参考基因组的正链与负链的序列的比例大于等于90%的序列支持的位点；2）除去低于过滤阈值的位点；3）从所述低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变；
所述变异结果文件中包括热点突变和非热点突变，所述热点突变包括待验证的低频突变，所述待验证的低频突变指低于所述过滤阈值的热点突变。
2.根据权利要求1所述的处理方法，其特征在于，所述S101包括：
将所述高通量测序数据与所述参考基因组进行初次比对，获得初始比对bam文件；
将所述初始比对bam文件与所述参考基因组进行局部重比对，获得所述比对变异bam文件。
3.根据权利要求1所述的处理方法，其特征在于，对所述CIGAR变异列表和所述MD变异列表的变异及位置信息执行合并操作，得到合并列表包括：
如果所述CIGAR变异列表中存在软剪切信息，则去除所述软剪切信息后，输出调整匹配信息为变异起始位点，并输出至所述合并列表；
如果所述CIGAR变异列表中的匹配M的位置同所述MD变异列表中的匹配M的位置相同，将所述CIGAR变异列表的匹配M信息输出至所述合并列表；
如果所述CIGAR变异列表中的匹配M的位置同所述MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至所述合并列表；
如果所述CIGAR变异列表中的缺失D的位置同所述MD变异列表中的缺失D的位置相同，将所述CIGAR变异列表中的缺失D的信息输出至所述合并列表；
如果所述CIGAR变异列表中的***I的位置同所述MD变异列表中的匹配M的位置重合，以所述***I的起始位置信息为分割，将在所述***I的位置***所述MD变异列表中的匹配M的位置中作为修正，输出至所述合并列表，并记录***长度；
如果所述CIGAR变异列表中的匹配M的位置同所述MD变异列表中的替换R的位置重合，以所述替换R的起始位置信息为分割，将所述替换R的位置***所述CIGAR变异列表中的匹配M的位置中作为修正，输出至所述合并列表。
4.根据权利要求1所述的处理方法，其特征在于，对所述合并列表进行校正，得到所述变异位点合并文件，包括：
对所述合并列表中的缺失D信息进行校正；
当所述合并列表中前一个突变信息为***I并且后一个突变信息为替换R时，根据所述序列信息对提取的所述CIGAR变异列表和所述MD变异列表进行整体校正；
对所述合并列表中的匹配M的位置和长度进行验证，得到校正后的所述合并列表；
根据校正后的所述合并列表，当突变信息符合如下至少之一的条件时，执行合并操作，获得所述变异位点合并文件：
1）两个替换R的突变信息相距4bp以内时；
2）***I的突变信息同其他任意突变信息相距20bp内时；
3）缺失D的突变信息同其他任意突变信息相距20bp内时。
5.根据权利要求1所述的处理方法，其特征在于，判断所述变异位点合并文件是否存在所述情形1）包括：
对所述变异位点合并文件中的突变信息进行注释，得到变异注释文件；
根据所述变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子。
6.一种高通量测序数据的处理装置，其特征在于，包括：
比对单元，被设置为将高通量测序数据与参考基因组进行比对，获得比对变异bam文件；
序列提取单元，被设置为对所述比对变异bam文件中的序列信息进行提取，得到序列提取文件；
变异位点合并单元，被设置为合并所述序列提取文件中的变异位点，得到变异位点合并文件；
变异确定单元，被设置为根据所述变异位点合并文件确定最终变异位点，输出变异结果文件；
所述序列提取单元包括：
提取模块，被设置为从所述比对变异bam文件中提取如下至少之一的序列信息，得到信息提取文件：1）序列字符；2）碱基质量值；3）flag值；4）CIGAR标签；5）MD标签；6）起始位点；7）比对质量；
筛选单元，被设置为根据所述信息提取文件对所述比对变异bam文件中的序列进行筛选，剔除如下至少之一的序列，获得所述序列提取文件：1）比对质量小于比对阈值的序列；2）在所述参考基因组上的比对位置不唯一的序列；3）重复序列；
所述变异位点合并单元包括：
CIGAR变异输出模块，被设置为根据所述信息提取文件中的所述CIGAR标签拆分CIGAR标签中的突变信息，输出带有第一突变信息的CIGAR变异列表，所述CIGAR变异列表中包含匹配M、缺失D、***I信息；
MD变异输出模块，被设置为根据所述信息提取文件中的所述MD标签拆分MD标签中的突变信息，输出带有第二突变信息的MD变异列表，所述MD变异列表中包含匹配M、缺失D和替换R信息；
合并模块，被设置为循环比对所述CIGAR变异列表和所述MD变异列表中的变异信息，对所述CIGAR变异列表和所述MD变异列表的变异及位置信息执行合并操作，得到合并列表；
校正模块，被设置为对所述合并列表进行校正，得到所述变异位点合并文件；
所述变异确定单元包括：
情形判断模块，被设置为判断所述变异位点合并文件是否存在如下至少之一的情形：情形1）分属不同密码子而被合并的SNV位点；情形2）属于胚系变异位点；
第一更新模块，被设置为在存在所述情形1）时，则将分属不同密码子的SNV位点取消合并，删除不同密码子合并的位点，保留取消合并的不同密码子上单个SNV位点信息，形成第一更新后的所述变异位点合并文件；
第二更新模块，被设置为在存在所述情形2）时，则从第一更新后的所述变异位点合并文件中去除所述胚系变异位点，得到第二更新后的所述变异位点合并文件；
重注释模块，被设置为利用所述第二更新后的所述变异位点合并文件进行突变重注释，获得重注释变异位点；
位点筛选及输出模块，被设置为对所述重注释变异位点进行位点筛选，获得所述最终变异位点，并将输出所述变异结果文件；
其中，所述位点筛选包括：1）去除比对到所述参考基因组的正链与负链的序列的比例大于等于90%的序列支持的位点；2）除去低于过滤阈值的位点；3）从所述低于过滤阈值的位点中查找并保留属于待验证的低频突变的热点突变；
所述变异结果文件中包括热点突变和非热点突变，所述热点突变包括待验证的低频突变，所述待验证的低频突变指低于所述过滤阈值的热点突变。
7.根据权利要求6所述的处理装置，其特征在于，所述比对单元包括：
第一比对模块，被设置为将所述高通量测序数据与所述参考基因组进行初次比对，获得初始比对bam文件；
第二比对模块，被设置将所述初始比对bam文件与所述参考基因组进行局部重比对，获得所述比对变异bam文件。
8.根据权利要求6所述的处理装置，其特征在于，所述合并模块包括：
第一判断输出子模块，被设置为在所述CIGAR变异列表中存在软剪切信息时，去除所述软剪切信息后，调整匹配信息为变异起始位点，并输出至所述合并列表；
第二判断输出子模块，被设置为在所述CIGAR变异列表中的匹配M的位置同所述MD变异列表中的匹配M的位置相同，将所述CIGAR变异列表的匹配M信息输出至所述合并列表；
第三判断输出子模块，被设置为在所述CIGAR变异列表的匹配M的位置同所述MD变异列表中的匹配M的位置不同但重合，将其中短的匹配M信息输出至所述合并列表；
第四判断输出子模块，被设置为在所述CIGAR变异列表中的缺失D的位置同所述MD变异列表中的缺失D的位置相同时，将所述CIGAR变异列表中的缺失D的信息输出至所述合并列表；
第五判断输出子模块，被设置为在所述CIGAR变异列表中的***I的位置同所述MD变异列表中的匹配M的位置重合，以所述***I的起始位置信息为分割，将在所述***I的位置***所述MD变异列表中的匹配M的位置中作为修正，输出至所述合并列表，并记录***长度；
第六判断输出模块，被设置为在所述CIGAR变异列表中的匹配M的位置同所述MD变异列表中的替换R的位置重合，以所述替换R的起始位置信息为分割，将所述替换R的位置***所述CIGAR变异列表中的匹配M的位置中作为修正，输出至所述合并列表。
9.根据权利要求6所述的处理装置，其特征在于，所述校正模块包括：
第一校正子模块，被设置为对所述合并列表中的缺失D信息进行校正；
第二校正子模块，被设置为当所述合并列表中前一个突变信息为***I并且后一个突变信息为替换R时，根据对所述序列信息提取的CIGAR变异列表和MD变异列表整体校正；
第三校正子模块，被设置为对所述合并列表中的匹配M的位置和长度进行验证，得到校正后的所述合并列表；
合并子模块，被设置为根据校正后的所述合并列表，当突变信息符合如下至少之一的条件时，执行合并操作，获得所述变异位点合并文件：
1）两个替换R的突变信息相距4bp以内时；
2）***I的突变信息同其他任意突变信息相距20bp内时；
3）缺失D的突变信息同其他任意突变信息相距20bp内时。
10.根据权利要求6所述的处理装置，其特征在于，所述情形判断模块包括：情形1）判断子模块和情形2）判断子模块，其中，所述情形1）判断子模块包括：
注释元件，被设置为对所述变异位点合并文件中的突变信息进行注释，得到变异注释文件；
异密码子判定模块，被设置为根据所述变异注释文件，以氨基酸位点是否发生变化为判断基准，判断多个SNV位点是否分属不同密码子。
11.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质包括存储的程序，其中，在所述程序运行时控制所述存储介质所在设备执行权利要求1至5中任意一项所述的高通量测序数据的处理方法。
12.一种处理器，其特征在于，所述处理器用于运行程序，其中，所述程序运行时执行权利要求1至5中任意一项所述的高通量测序数据的处理方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311335212.0A CN117079720B (zh) | 2023-10-16 | 2023-10-16 | 高通量测序数据的处理方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311335212.0A CN117079720B (zh) | 2023-10-16 | 2023-10-16 | 高通量测序数据的处理方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117079720A CN117079720A (zh) | 2023-11-17 |
| CN117079720B true CN117079720B (zh) | 2024-01-30 |
Family
ID=88708401
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311335212.0A Active CN117079720B (zh) | 2023-10-16 | 2023-10-16 | 高通量测序数据的处理方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117079720B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117746989B (zh) * | 2024-02-20 | 2024-05-10 | 北京贝瑞和康生物技术有限公司 | 变异描述信息的处理方法、装置及电子设备 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106202991A (zh) * | 2016-06-30 | 2016-12-07 | 厦门艾德生物医药科技股份有限公司 | 一种基因组多重扩增测序产物中突变信息的检测方法 |
| CN108280325A (zh) * | 2017-12-08 | 2018-07-13 | 北京雅康博生物科技有限公司 | 高通量测序数据的处理方法、处理装置、存储介质及处理器 |
| CN109698011A (zh) * | 2018-12-25 | 2019-04-30 | 人和未来生物科技（长沙）有限公司 | 基于短序列比对的Indel区域校正方法及*** |
| CN109767810A (zh) * | 2019-01-10 | 2019-05-17 | 上海思路迪生物医学科技有限公司 | 高通量测序数据分析方法及装置 |
| CN112086131A (zh) * | 2020-08-18 | 2020-12-15 | 西安医学院 | 一种高通量测序中假阳性变异位点的筛选方法 |
| CN112111565A (zh) * | 2019-06-20 | 2020-12-22 | 上海其明信息技术有限公司 | 一种细胞游离dna测序数据的突变分析方法和装置 |
| CN112397142A (zh) * | 2020-10-13 | 2021-02-23 | 山东大学 | 面向多核处理器的基因变异检测方法及*** |
| CN114913918A (zh) * | 2021-02-10 | 2022-08-16 | 中国科学院脑科学与智能技术卓越创新中心 | 一种针对孤独症的高通量测序数据分析方法及装置 |
| CN115896256A (zh) * | 2022-11-25 | 2023-04-04 | 臻悦生物科技江苏有限公司 | 基于二代测序技术的rna***缺失突变的检测方法、装置、设备和存储介质 |
| WO2023185559A1 (zh) * | 2022-03-28 | 2023-10-05 | 深圳吉因加医学检验实验室 | 一种用于结构变异检测的方法、装置和存储介质 |
-
2023
- 2023-10-16 CN CN202311335212.0A patent/CN117079720B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106202991A (zh) * | 2016-06-30 | 2016-12-07 | 厦门艾德生物医药科技股份有限公司 | 一种基因组多重扩增测序产物中突变信息的检测方法 |
| CN108280325A (zh) * | 2017-12-08 | 2018-07-13 | 北京雅康博生物科技有限公司 | 高通量测序数据的处理方法、处理装置、存储介质及处理器 |
| CN109698011A (zh) * | 2018-12-25 | 2019-04-30 | 人和未来生物科技（长沙）有限公司 | 基于短序列比对的Indel区域校正方法及*** |
| CN109767810A (zh) * | 2019-01-10 | 2019-05-17 | 上海思路迪生物医学科技有限公司 | 高通量测序数据分析方法及装置 |
| CN112111565A (zh) * | 2019-06-20 | 2020-12-22 | 上海其明信息技术有限公司 | 一种细胞游离dna测序数据的突变分析方法和装置 |
| CN112086131A (zh) * | 2020-08-18 | 2020-12-15 | 西安医学院 | 一种高通量测序中假阳性变异位点的筛选方法 |
| CN112397142A (zh) * | 2020-10-13 | 2021-02-23 | 山东大学 | 面向多核处理器的基因变异检测方法及*** |
| CN114913918A (zh) * | 2021-02-10 | 2022-08-16 | 中国科学院脑科学与智能技术卓越创新中心 | 一种针对孤独症的高通量测序数据分析方法及装置 |
| WO2023185559A1 (zh) * | 2022-03-28 | 2023-10-05 | 深圳吉因加医学检验实验室 | 一种用于结构变异检测的方法、装置和存储介质 |
| CN115896256A (zh) * | 2022-11-25 | 2023-04-04 | 臻悦生物科技江苏有限公司 | 基于二代测序技术的rna***缺失突变的检测方法、装置、设备和存储介质 |
Non-Patent Citations (1)
| Title |
|---|
| 高通量测序数据分析和临床诊断流程的解读;黎籽秀;刘博;徐凌丽;杨琳;王慧君;周文浩;;中国循证儿科杂志(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117079720A (zh) | 2023-11-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN117079720B (zh) | 高通量测序数据的处理方法和装置 | |
| Allam et al. | Karect: accurate correction of substitution, insertion and deletion errors for next-generation sequencing data | |
| US11538557B2 (en) | System and methods for indel identification using short read sequencing | |
| Ewing et al. | Base-calling of automated sequencer traces using phred. II. Error probabilities | |
| CN107391965A (zh) | 一种基于高通量测序技术的肺癌体细胞突变检测分析方法 | |
| JP6314091B2 (ja) | Ｄｎａ配列のデータ分析 | |
| CN103993069B (zh) | 病毒整合位点捕获测序分析方法 | |
| Nelson et al. | Whole-genome validation of high-information-content fingerprinting | |
| CN108280325B (zh) | 高通量测序数据的处理方法、处理装置、存储介质及处理器 | |
| CN108830044B (zh) | 用于检测癌症样本基因融合的检测方法和装置 | |
| CN108256295B (zh) | 一种用于检测基因融合的装置 | |
| CN110993023A (zh) | 复杂突变的检测方法及检测装置 | |
| CN110310702B (zh) | 一种基因组测序组装结果修复的方法、装置和存储介质 | |
| CN110692101A (zh) | 用于比对靶向的核酸测序数据的方法 | |
| CN111326212A (zh) | 一种结构变异的检测方法 | |
| CN104794371A (zh) | 检测逆转座子***多态性的方法和装置 | |
| CN113724791A (zh) | Cyp21a2基因ngs数据分析的方法、装置及应用 | |
| CN111292809B (zh) | 用于检测rna水平基因融合的方法、电子设备和计算机存储介质 | |
| CN115083521A (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** | |
| CN113205857B (zh) | 基因组性染色体非同源区域的鉴定方法和装置 | |
| CN112860957B (zh) | 一种定值单的核对方法、介质及*** | |
| CN107967411B (zh) | 一种脱靶位点的检测方法、装置及终端设备 | |
| CN113782101A (zh) | 高杂合二倍体序列组装结果去冗余的方法、装置及其应用 | |
| US20150142328A1 (en) | Calculation method for interchromosomal translocation position | |
| CN114566214B (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 |