CN116434843A - 一种碱基测序质量评估方法 - Google Patents
一种碱基测序质量评估方法 Download PDFInfo
- Publication number
- CN116434843A CN116434843A CN202310322408.XA CN202310322408A CN116434843A CN 116434843 A CN116434843 A CN 116434843A CN 202310322408 A CN202310322408 A CN 202310322408A CN 116434843 A CN116434843 A CN 116434843A
- Authority
- CN
- China
- Prior art keywords
- base
- sequencing
- nucleic acid
- acid sequence
- sequence
- 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.)
- Pending
Links
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 204
- 238000001303 quality assessment method Methods 0.000 title abstract description 13
- 150000007523 nucleic acids Chemical group 0.000 claims abstract description 123
- 108091028043 Nucleic acid sequence Proteins 0.000 claims abstract description 74
- 238000012549 training Methods 0.000 claims abstract description 11
- 238000003780 insertion Methods 0.000 claims description 71
- 230000037431 insertion Effects 0.000 claims description 71
- 238000006467 substitution reaction Methods 0.000 claims description 66
- 238000012217 deletion Methods 0.000 claims description 65
- 230000037430 deletion Effects 0.000 claims description 65
- 238000000034 method Methods 0.000 claims description 44
- 102000039446 nucleic acids Human genes 0.000 claims description 41
- 108020004707 nucleic acids Proteins 0.000 claims description 41
- 238000012937 correction Methods 0.000 claims description 26
- 108020004414 DNA Proteins 0.000 claims description 23
- 230000007614 genetic variation Effects 0.000 claims description 19
- 229920000642 polymer Polymers 0.000 claims description 19
- 238000007622 bioinformatic analysis Methods 0.000 claims description 17
- 238000004422 calculation algorithm Methods 0.000 claims description 17
- 239000003550 marker Substances 0.000 claims description 12
- 238000007672 fourth generation sequencing Methods 0.000 claims description 10
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 241000701959 Escherichia virus Lambda Species 0.000 claims description 6
- 241000588724 Escherichia coli Species 0.000 claims description 3
- 240000004808 Saccharomyces cerevisiae Species 0.000 claims description 3
- 235000014680 Saccharomyces cerevisiae Nutrition 0.000 claims description 3
- 238000013528 artificial neural network Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims description 3
- 238000012417 linear regression Methods 0.000 claims description 3
- 238000007477 logistic regression Methods 0.000 claims description 3
- 238000007637 random forest analysis Methods 0.000 claims description 3
- 238000012706 support-vector machine Methods 0.000 claims description 3
- 108090000623 proteins and genes Proteins 0.000 abstract description 5
- 239000002585 base Substances 0.000 description 204
- 239000000523 sample Substances 0.000 description 33
- 102000053602 DNA Human genes 0.000 description 20
- 230000035772 mutation Effects 0.000 description 17
- 238000012216 screening Methods 0.000 description 14
- 238000001514 detection method Methods 0.000 description 9
- 238000009792 diffusion process Methods 0.000 description 8
- 239000002773 nucleotide Substances 0.000 description 8
- 125000003729 nucleotide group Chemical group 0.000 description 8
- 230000008569 process Effects 0.000 description 8
- 230000002441 reversible effect Effects 0.000 description 8
- 206010028980 Neoplasm Diseases 0.000 description 6
- 210000004027 cell Anatomy 0.000 description 6
- 150000002500 ions Chemical class 0.000 description 6
- 238000012175 pyrosequencing Methods 0.000 description 6
- 239000000758 substrate Substances 0.000 description 6
- 238000010276 construction Methods 0.000 description 5
- 239000004065 semiconductor Substances 0.000 description 5
- 206010064571 Gene mutation Diseases 0.000 description 4
- 108091034117 Oligonucleotide Proteins 0.000 description 4
- 108091093037 Peptide nucleic acid Proteins 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 4
- 238000003766 bioinformatics method Methods 0.000 description 4
- 125000004122 cyclic group Chemical group 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000006116 polymerization reaction Methods 0.000 description 4
- 229920002477 rna polymer Polymers 0.000 description 4
- 241000894007 species Species 0.000 description 4
- 206010069754 Acquired gene mutation Diseases 0.000 description 3
- 230000005778 DNA damage Effects 0.000 description 3
- 231100000277 DNA damage Toxicity 0.000 description 3
- 201000011510 cancer Diseases 0.000 description 3
- 238000013145 classification model Methods 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 3
- 238000003909 pattern recognition Methods 0.000 description 3
- 238000002360 preparation method Methods 0.000 description 3
- 238000013441 quality evaluation Methods 0.000 description 3
- 238000002864 sequence alignment Methods 0.000 description 3
- 230000037439 somatic mutation Effects 0.000 description 3
- 210000001519 tissue Anatomy 0.000 description 3
- -1 xylose Nucleic Acid Chemical class 0.000 description 3
- 241000143060 Americamysis bahia Species 0.000 description 2
- 102100028630 Cytoskeleton-associated protein 2 Human genes 0.000 description 2
- 101000766848 Homo sapiens Cytoskeleton-associated protein 2 Proteins 0.000 description 2
- 108091092878 Microsatellite Proteins 0.000 description 2
- 241000283907 Tragelaphus oryx Species 0.000 description 2
- SRBFZHDQGSBBOR-IOVATXLUSA-N Xylose Natural products O[C@@H]1COC(O)[C@H](O)[C@H]1O SRBFZHDQGSBBOR-IOVATXLUSA-N 0.000 description 2
- PYMYPHUHKUWMLA-UHFFFAOYSA-N arabinose Natural products OCC(O)C(O)C(O)C=O PYMYPHUHKUWMLA-UHFFFAOYSA-N 0.000 description 2
- SRBFZHDQGSBBOR-UHFFFAOYSA-N beta-D-Pyranose-Lyxose Natural products OC1COC(O)C(O)C1O SRBFZHDQGSBBOR-UHFFFAOYSA-N 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 238000002144 chemical decomposition reaction Methods 0.000 description 2
- 239000003153 chemical reaction reagent Substances 0.000 description 2
- 210000000349 chromosome Anatomy 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 239000005546 dideoxynucleotide Substances 0.000 description 2
- WOERBKLLTSWFBY-UHFFFAOYSA-M dihydrogen phosphate;tetramethylazanium Chemical compound C[N+](C)(C)C.OP(O)([O-])=O WOERBKLLTSWFBY-UHFFFAOYSA-M 0.000 description 2
- 230000001973 epigenetic effect Effects 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 229910052739 hydrogen Inorganic materials 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 239000012528 membrane Substances 0.000 description 2
- 239000003147 molecular marker Substances 0.000 description 2
- 238000013188 needle biopsy Methods 0.000 description 2
- 210000000056 organ Anatomy 0.000 description 2
- 102000054765 polymorphisms of proteins Human genes 0.000 description 2
- 102000004169 proteins and genes Human genes 0.000 description 2
- 238000007480 sanger sequencing Methods 0.000 description 2
- 239000000344 soap Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 241000283690 Bos taurus Species 0.000 description 1
- 241000282472 Canis lupus familiaris Species 0.000 description 1
- 230000004543 DNA replication Effects 0.000 description 1
- 241000283086 Equidae Species 0.000 description 1
- 241000282326 Felis catus Species 0.000 description 1
- 208000033640 Hereditary breast cancer Diseases 0.000 description 1
- 208000026350 Inborn Genetic disease Diseases 0.000 description 1
- 241001494479 Pecora Species 0.000 description 1
- 206010036790 Productive cough Diseases 0.000 description 1
- 241000282887 Suidae Species 0.000 description 1
- 241000700605 Viruses Species 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 239000003513 alkali Substances 0.000 description 1
- 210000004381 amniotic fluid Anatomy 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 239000012491 analyte Substances 0.000 description 1
- 210000003567 ascitic fluid Anatomy 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000013060 biological fluid Substances 0.000 description 1
- 238000001574 biopsy Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000005119 centrifugation Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000009089 cytolysis Effects 0.000 description 1
- 238000007418 data mining Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000007865 diluting Methods 0.000 description 1
- 238000010790 dilution Methods 0.000 description 1
- 239000012895 dilution Substances 0.000 description 1
- 238000004821 distillation Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000002073 fluorescence micrograph Methods 0.000 description 1
- 238000013467 fragmentation Methods 0.000 description 1
- 238000006062 fragmentation reaction Methods 0.000 description 1
- 238000004108 freeze drying Methods 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 208000016361 genetic disease Diseases 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 208000025581 hereditary breast carcinoma Diseases 0.000 description 1
- 238000012165 high-throughput sequencing Methods 0.000 description 1
- 239000001257 hydrogen Substances 0.000 description 1
- 230000002779 inactivation Effects 0.000 description 1
- 238000010348 incorporation Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 244000005700 microbiome Species 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 210000004910 pleural fluid Anatomy 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 230000010076 replication Effects 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 210000003802 sputum Anatomy 0.000 description 1
- 208000024794 sputum Diseases 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 210000002700 urine Anatomy 0.000 description 1
Images
Classifications
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
- C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
- C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
- C12Q1/6869—Methods for sequencing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- 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
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/10—Signal processing, e.g. from mass spectrometry [MS] or from PCR
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Chemical & Material Sciences (AREA)
- Organic Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Biology (AREA)
- Wood Science & Technology (AREA)
- Zoology (AREA)
- General Health & Medical Sciences (AREA)
- Biotechnology (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Immunology (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Bioethics (AREA)
- Signal Processing (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Genetics & Genomics (AREA)
- Biochemistry (AREA)
- Microbiology (AREA)
- Analytical Chemistry (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了一种碱基测序质量评估方法和系统，用于基因测序领域。本发明提供的碱基测序质量评估方法，通过提供标准核酸序列，将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记；训练分类器以拟合碱基的测序信号特征与其标记之间的关系；提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；赋予待测核酸序列中的碱基以至少两种质量值等步骤实现。
Description
技术领域
本发明涉及一种碱基测序质量评估方法及系统，属于基因测序领域。
背景技术
核酸测序技术可以探明遗传物质的序列，被广泛应用于临床肿瘤分型、微生物鉴定和遗传病诊断等领域。当今主流的核酸测序技术除了产出被测核酸样品的序列之外，还会给所测得的每个碱基以一个质量值，用来评估其所测的准确性。这个质量值一般以Phred的形式表示：
q＝-10log10(1-a)
式中a为该碱基的准确率，q为Phred值。例如，准确率99％、99.9％、99.99％对应的Phred值分别为20、30、40。
在对核酸测序数据的生物信息学分析中，质量值起到了非常重要的作用。例如，在鉴定基因突变时，若所测序列上的某碱基与参考序列上的对应碱基不同，则当该碱基的质量值较高时，此处会被判定为基因突变；而当该碱基的质量值较低时，该序列会被认为发生了测序错误、不存在基因突变。
插入错误、缺失错误、替换错误为生物信息学中的常见概念。插入错误指的是测得序列相对参考序列，额外出现一个或多个碱基。缺失错误指的是测得序列相对参考序列，缺失了一个或多个碱基。替换错误指的是测得序列相对参考序列，一个或多个碱基被替换。现有技术中对每个碱基只赋予一个质量值，这种做法适用于Illumina和Ion Torrent的测序技术，因为Illumina的测序错误基本全是替换错误，Ion Torrent的测序错误基本全是插入和缺失错误，所以一个质量值足以满足其后续的生物信息学分析需求。
发明内容
在ECC测序中则面临不一样的错误模式：插入、缺失、替换三种错误大约各占三分之一。且这三种错误在ECC测序中的主要产生原因不同：插入与缺失主要是因为ECC在长同源多聚物上的准确度不及短同源多聚物，而替换主要由建库时的DNA损伤和PCR复制错误导致。这就导致拆分这三种不同的错误类型并以质量值的形式予以标明，有助于下游生物信息学分析时区分背景噪音和真实突变，给出更准确的鉴定结果。因此在ECC测序中，如果按照传统方法给每个碱基仅一个总质量值，会导致如下问题：第一，由于高通量测序仪可以产出大量数据，因此无论是测序仪的随机软件还是常规生物信息学流程都包含大量的数据筛选步骤，根据碱基质量值筛选出高质量的序列数据，提高分析准确率。在ECC测序中，总质量值会使得很多序列在长同源多聚物处的质量值偏低，在数据筛选时损失较多数据，影响下游生物信息学分析。区分三种错误类型并分别赋予质量值，有助于提高碱基质量的区分度，减少在数据筛选中的损失。第二，建库中引入的DNA损伤和PCR引入的DNA复制错误难以消除，却和测序仪无关，影响了测序仪对碱基质量值的判断，使得总质量值并不能如实反映测序仪的真实准确率，而仅仅反映了建库的准确性。
在此，本发明公开了一种对碱基赋予多个质量值的方法以克服上述问题。
具体的，本发明提供了一种碱基测序质量评估方法，其特征在于，包括：
1)提供标准核酸序列，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
根据优选的实施方式，根据第4)步中计算所得的测序信号特征，赋予待测核酸序列中的碱基以三种质量值，分别表示碱基的插入错误率，缺失错误率，替换错误率。
根据优选的实施方式，测序方法包括纠错码测序，即ECC测序。
根据优选的实施方式，测序方法包括纳米孔测序。
根据优选的实施方式，碱基的测序信号特征，指的是测序过程中该碱基发生测序化学反应时产生的信号的特征，包括但不限于：碱基种类，碱基在序列上的位置，碱基所处多聚物的长度，碱基在其所处多聚物中的位置，碱基发生测序化学反应的轮数，信号强度，信号强度(及其邻近信号强度)接近整数的程度，测序信号的参数(单位信号、背景信号、超前系数、滞后系数、衰减系数)，碱基对应的失相程度，等等。
根据优选的实施方式，分类器根据碱基的测序信号特征，将碱基分成若干类，统计每一类碱基的准确率。
根据优选的实施方式，步骤3)拟合碱基的测序信号特征与其标记之间的关系，包括将分类器的拟合结果转化为碱基的测序质量值。
根据优选的实施方式，质量值指的是表征测序准确率的一个数值，选自准确率、错误率、Phred值等。例如的，准确率99％、99.9％、99.99％对应的错误率分别为1％、0.1％、0.01％，对应的Phred值分别为20、30、40。
根据优选的实施方式，质量值对数地基于碱基检出误差概率，并且其中所述质量值包括Q10、Q15、Q20、Q25、Q30、Q35、Q40、Q45、Q50、Q55、Q60。
根据优选的实施方式，分类器包括线性回归，多项式回归，逻辑回归，支持向量机，人工神经网络，随机森林，Phred算法，集成学习等。
根据优选的实施方式，核酸样品包括脱氧核糖核酸(DNA)、核糖核酸(RNA)、肽核酸(PNA)、木糖核酸(XNA)、锁式核酸(LNA)等。
根据优选的实施方式，标准核酸样品，指的是来源和序列均已确定、在基因组的几乎所有位点上均高度纯合的核酸样品，包括大肠杆菌DNA，酿酒酵母DNA、λ噬菌体DNA等。
根据优选的实施方式，核酸序列包括序列确定的碱基序列或者简并碱基序列。
根据优选的实施方式，测序方法包括双脱氧核苷酸终止法(Sanger测序法)、化学降解法(Gilbert法)、焦磷酸测序法(pyrosequencing)、半导体测序法(semiconductorsequencing)、循环可逆终止法(cyclic reversible terminator)、荧光发生测序法(fluorogenic sequencing)、纠错码测序法(error-correction code sequencing)、模糊测序法(fuzzy sequencing)、联合探针锚定连接法(combinatorial probe-anchorligation)、联合探针锚定聚合法(combinatorial probe-anchor polymerization)、寡核苷酸连接检测测序法(sequencing by oligonucleotide ligation and detection)、边结合边测序法(sequencing-by-binding)、单分子荧光测序法、单分子实时测序、纳米孔测序法等。
根据优选的实施方式，标准核酸样品和待测核酸样品的测序方法相同，例如的，二者均使用ECC测序法或者均使用前述的任一测序方法。
根据优选的实施方式，质量评估方法还包括：6)对赋予质量值后的核酸序列进行生物信息学分析。
根据优选的实施方式，生物信息学分析包括，根据所赋予的质量值筛选高质量的核酸序列。筛选方法包括但不限于，筛选全部质量值均高于或低于某一阈值的核酸序列，筛选全部质量值的均值均高于或低于某一阈值的核酸序列，筛选核酸序列中质量值均高于或低于某一阈值的区域，筛选核酸序列中质量值的均值均高于或低于某一阈值的区域，等。
根据优选的实施方式，生物信息学分析包括，根据所赋予的质量值，将核酸序列比对到参考序列上。比对是生物信息学中的常规概念，可以利用Smith-Waterman算法、Bowtie、BWA、SOAP、Needleman-Wunch算法、Bowtie2、BLAST、ELAND、TMAP、MAQ、minimap2、SHRiMP等进行。
根据优选的实施方式，生物信息学分析包括，根据比对结果及被比对序列所赋予的质量值，鉴定基因变异。基因变异是生物学中的常规概念，包括但不限于单核苷酸多态性、拷贝数变异、表观遗传学变异、大范围结构变异等。
根据优选的实施方式，生物信息学分析包括，鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入变异时，使用表征插入错误率的质量值来计算变异存在的概率；当鉴定缺失变异时，使用表征缺失错误率的质量值来计算变异存在的概率。
根据优选的实施方式，生物信息学分析包括，鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入缺失变异时，使用表征插入和/或缺失错误率的质量值来计算变异存在的概率。
根据优选的实施方式，鉴定基因变异时，可以利用比对结果的某些特征，来去除潜在的假阳性或假阴性结果。
根据优选的实施方式，生物信息学分析包括，根据所赋予的质量值，将核酸序列组装为较长的核酸序列。
根据优选的实施方式，步骤1)-3)可预先完成，并将训练好的分类器作为配置文件存于系统中，在执行步骤4)-5)时调取即可。
根据优选的实施方式，标准核酸样品和待测核酸样品可带上不同的分子标记，并混合在一起同时测序。测序结束后，先利用分子标记将两种样品拆分出来，完成步骤1)-3)，得到训练好的分类器，再应用在待测核酸样品上。
本发明还提供一种用于评估核酸测序碱基识别的准确度的系统，其包括：
机器可读存储器；和
被配置成执行机器可读指令的处理器，指令在由处理器执行时使系统执行一方法，该方法包括：处理器接收
1)标准核酸序列，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
根据优选的实施方式，根据第4)步中计算所得的测序信号特征，赋予待测核酸序列中的碱基以三种质量值，分别表示该碱基的插入错误率，缺失错误率，替换错误率。
本发明的有益效果
本发明公开的对一个碱基赋予多种质量值的方法，相比于现有技术只赋予一种质量值的方法，具有如下明显优势：
1.由于样品制备过程一般难以引入插入/缺失错误，所以ECC测序的插入/缺失质量值可以很高，轻易达到Q60，有助于精确地鉴定插入/缺失变异。
2.由于长DPL的序列测序并不容易发生替换错误，所以随着DPL的增加，ECC测序的替换质量值下降并不明显，有助于鉴定长DPL处的替换变异。
附图说明
本发明的新颖特征在所附权利要求书中具体阐述。将参考以下详细描述和附图来获得对本发明特征和优势的更好理解，以下详细描述阐述利用本发明原理的说明性实施例，在附图中：
图1示出了测序反应的三种错误类型。
图2说明了测序信号特征的实例。
图3说明了根据一个实施例的碱基测序质量值，包括：插入质量值、缺失质量值、替换质量值和总质量值。
图4说明了根据一个实施例的不同简并多聚物长度对应的4种碱基质量值。
图5说明了序列重比对的示意图。
具体实施方式
除非另外定义，否则本文使用的所有科学和技术术语的含义与本领域普通技术人员通常理解的含义相同。
术语解释
每个
如本文所用，当参考项目的集合使用时，术语“每个”旨在识别集合中的单个项目，但不一定是指集合中的每个项目。如果明确公开或上下文另有明确规定，则可能会出现例外情况。
包括
术语“包括”在本文中旨在为开放式的，不仅包括所列举的元素，而且还涵盖任何附加的元素。
简并碱基
本发明中，按照IUPAC符号命名规则(Nucleic acid notation)，使用下面表1的字母表示简并碱基，例如字母M表示A和/或C。
表1
| 字母 | 所代表的 |
| M | A/C |
| K | G/T |
| R | A/G |
| Y | C/T |
| W | A/T |
| S | C/G |
| B | C/G/T |
| D | A/G/T |
| H | A/C/T |
| V | A/C/G |
简并多聚物长度(Degenerate
Polymer
Length，DPL)
简并多聚物即由简并碱基构成的多聚物，例如字母M表示A和/或C，简并多聚物MMKKK，其长度为5，即DPL为5。
碱基在其所处多聚物中的位置
即多聚物中的碱基与多聚物中最近的一个末端之间的距离，举例说明如下，左边所示为多聚物序列，右边所示为每个碱基在该多聚物中的位置：
A:0
AA:0，0
AAA:0，1，0
AAAA:0，1，1，0
AAAAA:0,1,2,1,0
AAAAAA:0,1,2,2,1,0.
纠错码(Error
Correcting
Code,ECC)测序以及纠错校正
本申请中，纠错码测序具备如下特征：
该测序方法需要多回测序，每回测序得到的信息不完整，而多回测序得到的总的信息是冗余的；利用多回测序的信息冗余来检测和校正潜在的测序错误，得到高准确度的序列。例如的，以2+2测序为例，将测序试剂按对偶碱基分为两两匹配的三组(例如的，分别为MK、RY、WS三组)，并对待测DNA序列进行三回独立测序，继而产生三条简并序列编码，这三条编码可互为校验，后续不但能够通过解码推导出真实碱基序列信息，而且具备对单回测序错误位点的校正能力。此校正过程即为纠错校正。
失相校正
在高通量测序中每个测序单元中包含有数千至数万条不等的DNA分子，它们由同一条DNA模板复制而得，具有相同的序列，将测序信号放大到可被准确检测的水平。然而在测序过程中有些DNA分子并不能保证每轮都充分反应，导致滞后现象，另一些会因为底物掺杂或错配等提前反应，导致超前现象。这些DNA分子的延伸会随着测序进行逐渐不再同步，这一现象被称作测序中的“失相”。在算法上将失相的测序信号根据测序反应的模型进行拟合和重新纠正的过程则称作“失相校正”。
归一化信号
即经过归一化后的信号，原始测序信号经过衰减校正后得到较为准确的单位信号，每个测序位点衰减校正后的测序信号与此位点的单位信号的比值即为此位点在每个测序轮中的归一化信号。
比对
比对(align或alignment)是生物信息学中的常见概念，在生物信息学中，比对经常用于比较不同核酸之间或者不同蛋白质之间的相似性。本发明中的比对指的是将测序得到的碱基序列和参考序列进行比较，从而确定测序所得的碱基序列正确与否。常用的序列比对算法及软件包括但不限于，例如的，Smith-Waterman算法、Bowtie、BWA、SOAP、Needleman-Wunch算法、Bowtie2、BLAST、ELAND、TMAP、MAQ、minimap2、SHRiMP等。
参考序列
参考序列是指可用于参考来自受试者的已鉴定序列的任何生物体的任何特定已知基因组序列，无论是部分的还是完整的。例如，可在美国国家生物技术信息中心(National Center for Biotechnology Information)找到用于人类受试者以及许多其他生物体的参考基因组。“基因组”是指以核酸序列表达的生物体或病毒的完整遗传信息。基因组既包括基因又包括DNA的非编码序列。参考序列可大于与其比对的读段。例如，参考序列可为比对读段的至少约100倍大、或至少约1000倍大、或至少约104倍大、或至少约105倍大、或至少约106倍大、或至少约107倍大。参考序列的其他示例包括其他物种的基因组，以及任何物种的染色体、亚染色体区域(诸如链)等。本发明中，参考序列对应的物种优选的来自包括大肠杆菌、λ噬菌体、酿酒酵母等物种的核酸序列，这些物种的核酸序列已知，且在基因组的几乎所有位点上均高度纯合；λ噬菌体优选的来自New England Biolabs公司生产的λ噬菌体DNA。
变异
是指与核酸参考序列不同的核酸序列。典型的变异包括但不限于单核苷酸变异(SN)、短缺失和插入多态性(Indel)、拷贝数变异(CNV)、表观遗传学变异、微卫星标记或短串联重复序列和结构变异。体细胞变异检出是识别以低频率存在于DNA样本中的变异的工作。体细胞变异检出在癌症治疗的背景下是引人关注的。癌症是由DNA中突变的积聚引起的。来自肿瘤的DNA样本通常是异质的，包括一些正常细胞、癌症进展早期的一些细胞(具有较少突变)和一些晚期细胞(具有较多突变)。由于这种异质性，当对肿瘤(例如，来自FFPE样本)测序时，体细胞突变将通常以低频率出现。例如，可在覆盖给定碱基的读段的仅10％中看到SNV。
位置
术语“位置”是指核苷酸序列内的一个或多个核苷酸的位置或坐标，也指核苷酸序列中的一个或多个碱基对的位置或坐标。
分类器
分类器分类是数据挖掘的非常重要的方法，在机器学习中,分类器的作用是在标记好类别的训练数据基础上判断一个新的观察样本所属的类别。
分类器的构造和实施大体需经过以下几个步骤：
选定样本(包含正样本和负样本)，将所有样本分成训练样本和测试样本两部分；
在训练样本上执行分类器算法，生成分类模型；
在测试样本上执行分类模型，生成预测结果。
在优选的实施方式中，根据预测结果，计算必要的评估指标，评估分类模型的性能。
需要说明的是，本发明中上述指出的术语，并非是独特的第一次发明的术语或者定义。申请人之前的专利或者本领域的基础知识中，均为通用的解释。不需要将该定义作为独特特征。涉及到数学的用语，实际也是数学领域的常见用语。其并未包含特殊含义。
发明详述
插入错误、缺失错误、替换错误为生物信息学中的常见概念。插入错误指的是测得序列相对参考序列，额外出现一个或多个碱基。缺失错误指的是测得序列相对参考序列，缺失了一个或多个碱基。替换错误指的是测得序列相对参考序列，一个或多个碱基被替换。图1展示了这三种测序错误的例子，从左到右，测得序列相对参考序列，分别插入了一个G、缺失了一个T、以及G被替换成T。在ECC测序中，三种错误均有发生，且比例相当，有必要分别对其进行质量评估，本发明即公开了这样一种碱基测序质量评估方法，对待测核酸序列中的每个碱基赋予至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率，以提高碱基读出的准确性，更具针对性地进行后续的生物信息学分析。
具体的，本发明的第一方面公开了一种碱基测序质量评估方法，其特征在于，包括：
1)提供标准核酸序列，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；
3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
在优选的实施方式中，赋予所测核酸序列中的碱基以三种质量值，分别表示该碱基的插入错误率，缺失错误率，替换错误率。
本发明中，核酸包括脱氧核糖核酸(DNA)、核糖核酸(RNA)、肽核酸(PNA)、木糖核酸(XNA)、锁式核酸(LNA)等。核酸样品是指包含核酸或核酸混合物的样本，通常来源于生物流体、细胞、组织、器官或生物体，该核酸或核酸混合物包含待测序和/或定相的至少一种核酸序列。此类样本包括但不限于血液、血液级分、痰/口腔液、羊水、细针活检样本(例如，外科活检、细针活检等)、尿液、腹膜液、胸膜液、组织外植体、器官培养物和任何其他组织或细胞制剂，或其级分或衍生物，或从其分离的级分或衍生物。虽然样本通常取自人类受试者(例如，患者)，但样本可取自具有染色体的任何生物体，包括但不限于牛、马、猪、羊、狗、猫等。样本可按从生物来源获得的原样直接使用，或者经过预处理以改变样本的性质后使用。例如，此类预处理可包括由血液制备血浆、稀释粘性流体等。预处理的方法还可涉及但不限于过滤、沉淀、稀释、蒸馏、混合、离心、冷冻、冻干、浓缩、扩增、核酸片段化、干扰组分的灭活、添加试剂、裂解等。
本发明中，测序方法包括双脱氧核苷酸终止法(Sanger测序法)、化学降解法(Gilbert法)、焦磷酸测序法(pyrosequencing)、半导体测序法(semiconductorsequencing)、循环可逆终止法(cyclic reversible terminator)、荧光发生测序法(fluorogenic sequencing)、纠错码测序法(error-correction code sequencing)、模糊测序法(fuzzy sequencing)、缺失测序(专利CN202210104037.3)、联合探针锚定连接法(combinatorial probe-anchor ligation)、联合探针锚定聚合法(combinatorial probe-anchorpolymerization)、寡核苷酸连接检测测序法(sequencing by oligonucleotideligation and detection)、边结合边测序法(sequencing-by-binding)、单分子荧光测序法、单分子实时测序、纳米孔测序法等。
根据优选的实施方式，测序方法是纠错码(ECC)测序，在ECC测序中，插入错误、缺失错误、替换错误三种类型所占的比例相当，如果按照传统方式，一个碱基只赋予一个总的质量值，那么ECC测序存在两个问题：第一，质量值最高只能达到Q50。这是因为样品制备过程会引入一些替换错误(PCR错误、DNA损伤等)，难以制备高纯度样品；第二，质量值随简并多聚物长度(DPL)的增加而快速下降，这是因为长DPL测不准，容易发生插入和缺失错误。可以看到，ECC测序发生替换和插入/缺失错误的原因是不同的，因此需要对每个碱基赋予多个质量值(至少2种，例如的，2种或3种)，分别表示碱基的替换错误率和插入/缺失错误率，以克服ECC测序中存在的上述问题。
在一些实施方式中，对标准核酸样品进行ECC测序，对待测核酸样品进行ECC测序，两种样品的测序反应可以同时进行，也可以先后进行。
根据优选的实施方式，测序方法是纳米孔测序，在纳米孔测序中，纳米孔被固载在薄膜上，薄膜一侧的离子在浓度差的作用下穿过纳米孔，扩散到另一侧，形成扩散电流。待测DNA穿过纳米孔，阻碍了离子的穿孔，造成扩散电流下降。不同的碱基因其体积和化学性质的不同，对离子穿孔的阻碍大小不同，造成的扩散电流下降也不同，因此纳米孔测序通过检测扩散电流的变化就可以获得待测DNA的序列信息。由于单分子运动的随机性、扩散电流检测的误差等因素，纳米孔测序的错误率非常高，常常在10％以上，且插入、缺失、替换等多种错误模式并存。纳米孔测序对扩散电流进行高频检测，因此通常多个扩散电流值对应一个被测出的碱基，所述多个扩散电流值中不仅包含碱基种类的信息，也包含其潜在的错误类型的信息。因此，在纳米孔测序中，需要对每个碱基赋予多种质量值(至少2种，例如的，2种或3种)，分别表示碱基的替换错误率，插入/缺失错误率，以充分挖掘所检测的扩散电流值所包含的信息，提示可能存在的错误模式，提高下游生物信息学分析的准确率。
本发明中，得到核酸序列的过程也就是碱基识别(base calling)的过程，测序仪在测序运行期间产生原始数据，例如的，这些原始数据可以包括荧光图像数据，可以是单色的、多色的荧光图像，以循环可逆终止法测序(cyclic reversible terminator)、焦磷酸测序法(pyrosequencing)、荧光发生测序法(fluorogenic sequencing)、纠错码测序法(error-correction code sequencing)为代表；这些原始数据还可以包括电压变化，以IonTorrent的半导体测序法为典型代表，碱基并入释放的氢离子导致pH变化，最终反映为与掺入的核苷酸的数量成比例的电压变化；另选的，原始数据还包括以牛津纳米孔测序技术(ONT)为代表的测序方法中的电流信号，依赖于纳米孔感测使用生物传感器来测量当分析物穿过纳米孔或靠近其孔口时电流的中断，同时确定碱基的种类。
在优选的实施方式中，碱基识别得到的核酸序列是确定的碱基序列，即由A,G,C,T表示的序列，或者由A,G,C,U表示的序列。
在一些实施方式中，碱基识别得到的核酸序列是简并碱基序列，即由M,K,R,Y,W,S,B,D,H,V等表示的简并碱基，以MK测序为例，M底物的A和C分别用两种不同的荧光染料标记，K底物的G和T分别用两种不同的荧光染料标记，对于每个测序反应循环，可以得到延伸的两种碱基的数量，此简并碱基序列简称为双色简并碱基序列。
可以理解的，核酸序列还可以是确定的碱基序列和简并碱基序列的组合，例如的，以A、B为底物进行测序，以A为底物的测序反应得到的序列是确定的碱基序列，以B为底物的反应得到的序列是简并碱基序列。
本发明中，碱基的测序信号特征，指的是测序过程中，被测序列上的该碱基发生测序化学反应时产生的信号的特征，图2给出了测序信号特征的实例，包括但不限于：该碱基的种类，即碱基属于A,G,C,T(或U)的哪一种；该碱基在序列上的位置，即碱基在其所在核苷酸序列上的位置位次，例如的，对于单端测序，位置靠前的碱基的测序质量值通常高于位置靠后的碱基；该碱基所处多聚物的长度，即碱基所处的同源多聚物或简并多聚物的碱基的数量，通常的，多聚物长度短，测序质量值高；该碱基在其所处多聚物中的位置，即碱基与其所处的同源多聚物或简并多聚物的最近一个末端的距离；该碱基发生测序化学反应的轮数，即该碱基并入核苷酸链时对应的cycle数，通常的，其对应的cycle数小，质量值高；信号强度，可以是测序仪直接采集到的信号的强度，包括亮度、电压水平或电流水平等，可以是归一化信号，可以是失相校正后的信号；信号强度(及其邻近信号强度)接近整数的程度，即归一化信号或失相校正后信号或纠错校正后的信号与最接近整数之间的差值，通常的，差值小的，准确度更高；测序信号的参数，即单位信号、背景信号、超前系数、滞后系数、衰减系数等；测到该碱基时的失相程度，通常的，失相程度低，准确度更高；等等。
在具体的实施方式中，将标准核酸序列比对到其对应的参考序列，得到比对结果，再根据比对结果将碱基标记为测序正确或测序错误；优选的，从比对结果中进一步筛选出高质量比对的碱基序列，再将高质量比对的碱基序列中的碱基标记为测序正确或测序错误，忽略无法确定的碱基(即无法成功比对到参考序列上的碱基或比对质量较低的碱基)。根据比对结果，将比对结果为“匹配”的碱基标记为“测序正确”，将比对结果为“错配”的碱基标记为“替换错误”，将比对结果为“插入”的碱基标记为“插入错误”，将比对结果为“缺失”的碱基标记为“缺失错误”；可选的，将比对结果为“错配”的碱基标记为“替换错误”，将比对结果为“插入”或“缺失”的碱基标记为“插入或缺失错误”。本发明中所述的高质量比对，需要根据所用的比对软件或算法来具体选择质量值范围；例如的，当使用BWA进行序列比对时，高质量比对的碱基序列指的是，比对质量大于0、或大于等于10、或大于等于20、或大于等于30、或大于等于40、或大于等于50、或大于等于60的碱基序列。
本发明中，分类器是模式识别领域的常规概念，包括线性回归，多项式回归，逻辑回归，支持向量机，人工神经网络，随机森林，Phred算法、集成学习等。随着模式识别领域的发展，近年来有多种新颖的分类器算法提出。使用新颖的分类器算法并不改变本发明的实质。
根据优选的实施方式，分类器可以根据碱基的测序信号特征，将碱基分成若干类，统计每一类碱基的准确率。例如的，可以将位于序列1-20、21-40、41-60、61-80、81-100bp的碱基分别划为一类，或将位于长度为1、2、3、4、5及5以上的多聚物中的碱基分别划为一类。当使用多种测序信号特征时，可进行正交划分，例如位于序列1-20bp内且位于1倍多聚物中的碱基划为一类，位于序列21-40bp内切位于2倍多聚物中的碱基划为另一类，以此类推。
在优选的实施方式中，训练分类器，来拟合碱基的测序信号特征与其标记之间的关系，包括将分类器的拟合结果转化为质量值。存在大量文献报道如何将分类器的预测结果转化为质量值。以著名的softmax算法为例，设某分类器的输出为(a,b)，其中(1,0)表示正确，(0,1)表示错误。由于分类器训练的精度或预测时的计算误差等因素，分类器在预测时的输出并不总恰好是(1,0)或(0,1)，而是(0.9,0.05)或(0.1,0.99)这样较为接近(1,0)或(0,1)的数值。此时softmax算法利用下式将输出(a,b)转化为正确率：
随着模式识别领域的发展，近年来有多种新颖的转化算法提出，例如的，包括Sparse-softmax，log-softmax,Taylor softmax,log-Taylor softmax,soft-marginsoftmax,SM-Taylor softmax等。使用新颖的转化算法并不改变本发明的实质。
质量值是表征测序准确率的一个数值，质量值可以不同的数学方式表达，如准确率、错误率、Phred值等。例如，准确率99％、99.9％、99.99％对应的错误率分别为1％、0.1％、0.01％，对应的Phred值分别为20、30、40。在某些实现中，为了便于记录和存储，会将Phred值加上33后转为ASCII码，例如Phred值20、30、40会分别转为字符’5’、’？’、’I’。质量值表达形式的不同不影响本发明的实质。
在优选的实施方式中，质量值对数地基于碱基检出误差概率，并且其中质量值包括Q10、Q15、Q20、Q25、Q30、Q35、Q40、Q45、Q50、Q55、Q60。
在优选的实施方式中，步骤1)-3)可预先完成，并将训练好的分类器作为配置文件存于系统中，在执行步骤4)-5)时调取即可。
在优选的实施方式中，标准核酸样品和待测核酸样品可带上不同的分子标记，并混合在一起同时测序。测序结束后，先利用分子标记(例如的，index)将两种样品拆分出来，对于标准核酸序列执行完成步骤1)-3)，得到训练好的分类器，再应用在待测核酸样品上，即执行步骤4)-5)。
根据优选的实施方式，质量评估方法还包括序列比对，比对时的罚分依赖于单个碱基的至少两种质量值，不必进行重比对(realignment)步骤，简化生物信息学流程，加快分析速度。具体实施方式为：比对软件不再采用固定的罚分规则，而使罚分依赖于碱基的不同质量值。即，若一个碱基的替换质量值较低，则在判定该碱基存在替换时赋予较低的罚分；若一个碱基的插入质量值较低，则在判定该碱基存在插入时赋予较低的罚分；若一个碱基的缺失质量值较低，则在判定该碱基存在缺失时赋予较低的罚分。
根据优选的实施方式，质量评估方法还包括：6)对赋予质量值后的核酸序列进行生物信息学分析。
在优选的实施方式中，生物信息学分析，可以是根据碱基所赋予的质量值，筛选高质量的核酸序列。筛选方法包括但不限于，筛选全部质量值均高于或低于某一阈值的核酸序列，筛选全部质量值的均值均高于或低于某一阈值的核酸序列，筛选核酸序列中质量值均高于或低于某一阈值的区域，筛选核酸序列中质量值的均值均高于或低于某一阈值的区域，等。所述阈值可基于经验分析而改变，例如的，阈值设置为15，或20，或25，或30。
在优选的实施方式中，生物信息学分析指的是，根据比对结果及被比对序列所赋予的质量值，鉴定基因变异。
根据优选的实施方式，生物信息学分析包括鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入变异时，使用表征插入错误率的质量值来计算变异存在的概率；当鉴定缺失变异时，使用表征缺失错误率的质量值来计算变异存在的概率。
根据优选的实施方式，生物信息学分析包括鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入缺失变异时，使用表征插入和/或缺失错误率的质量值来计算变异存在的概率。
在一些实施方式中，生物信息学分析，可以是在鉴定基因变异时，给出零假设：该位点不存在基因变异。根据质量值和比对结果，计算零假设成立的概率，若该概率大于给定的显著性水平，则接受零假设，否则拒绝零假设、认为该位点存在基因变异。在计算所述零假设成立的概率时：
1.当鉴定插入变异时，仅使用插入质量值；
2.当鉴定缺失变异时，仅使用缺失质量值；
3.当鉴定替换变异时，仅使用替换质量值。
在ECC测序中，当DPL较大时，若只使用一个质量值来评估碱基质量，则可能因为质量值较低而无法检测到较低频的基因突变；而本发明使用多个质量值来评估碱基测序质量，其中，替换质量值受长DPL的影响较小，质量值较高，更容易成功检出较低频的突变。
根据优选的实施方式，鉴定基因变异时，可以利用比对结果的某些特征，来去除潜在的假阳性或假阴性结果。这些都是生物信息学中的常规操作，其增添均不影响本发明的实质。这样的特征包括但不限于：
1.该基因变异集中出现在正向或反向比对的序列上，而在反向或正向比对的序列上较少出现；
2.该基因变异集中出现在序列的两端，而在序列的中央较少出现；
3.当使用双端测序(paired-end sequencing)时，read1测到该位点主要为G变T、而read2测到该位点主要为C变A，或read1测到该位点主要为C变T、而read2测到该位点主要为G变A；
4.该基因变异附近频繁出现其他不同的基因变异。
在某些实施方式中，插入质量值和缺失质量值可以合并为插入缺失质量值。即每个碱基被赋予插入缺失质量值和替换质量值两种质量值，分别表示该碱基发生插入或缺失错误，或者替换错误的概率。前述相关操作均可自然推广到两种质量值的情况，在此不再赘述。
本发明的第二方面提供一种用于评估核酸测序碱基识别的准确度的系统，其包括：
机器可读存储器；和
被配置成执行机器可读指令的处理器，所述指令在由所述处理器执行时使所述系统执行一方法，所述方法包括：所述处理器接收
1)标准核酸序列，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
在优选的实施方式中，根据第4)步中计算所得的测序信号特征，赋予待测核酸序列中的碱基以三种质量值，即：插入质量值、缺失质量值、替换质量值，分别表示该碱基的插入错误率，缺失错误率，替换错误率。
在本发明的碱基测序质量评估方法的具体实施部分中所讨论的特征中的特征同样适用于用于评估核酸测序碱基识别的准确度的系统的具体实施。如上所示，所有其他特征在此处不再重复，并且应被视为以引用方式重复。本领域普通技术人员将理解在这些具体实施中识别的特征可如何容易地与在其他具体实施中识别的基本特征组组合。
实施例1
从New England Biolabs公司购买λ噬菌体的基因组DNA，建库后进行纠错码测序，得到碱基序列(即标准核酸序列)。对每个碱基，计算如下5种预测器的值：
1.该碱基在序列中的位置；
2.该碱基所处同源多聚物的长度；
3.该碱基所处失相校正信号与纠错校正信号间的差异；
4.该碱基所处失相校正信号与其取整信号间的欧氏距离；
5.该碱基所在序列全部失相校正信号与其取整信号间的欧氏距离。
利用BWA-MEM软件将所得DNA序列比对到参考基因组上，忽略未比对序列。根据比对结果，将碱基分别标记为“匹配”、“替换”、“插入”、“缺失”中的一种。利用上述5种预测器，通过Phred算法构建质量评估表。构造时，对碱基标记分别采取以下4种做法：
1.仅将“插入”视为错误，其他标记视为正确；
2.仅将“缺失”视为错误，其他标记视为正确；
3.仅将“替换”视为错误，其他标记视为正确；
4.仅将“匹配”视为正确，其他标记视为错误。
从而得到4张质量评估表，分别为插入质量值、缺失质量值、替换质量值、总质量值。统计这4种质量评估表的区分度，如图3所示：(即质量值在某一阈值以上的碱基的比例)。
具体的，质量值不低于20、30、40、50、60的比例如表2所示：
表2
| 质量值类型 | 插入质量值 | 缺失质量值 | 替换质量值 | 总质量值 |
| AQ20 | 100％ | 91.56％ | 92.56％ | 85.79％ |
| AQ30 | 82.49％ | 74.37％ | 82.63％ | 73.58％ |
| AQ40 | 73.34％ | 63.56％ | 65.85％ | 58.08％ |
| AQ50 | 58.25％ | 59.40％ | 6.58％ | 1.95％ |
| AQ60 | 50.78％ | 49.36％ | 0％ | 0％ |
可见，总质量值最高只达到了1.95％Q50，而插入质量值和缺失质量值最高不仅均达到了Q60，而且比例接近50％。替换质量值虽然最高也是Q50，但其比例大幅提高至6.58％。因此，本发明可以显著提升碱基质量值的区分度，其中插入和缺失质量值可以达到Q60不是显而易见的。
图4从上到下依次为碱基所处同源多聚物长度和其插入质量值、缺失质量值、替换质量值、总质量值之间的关系(小提琴图)，可见，当同源多聚物较长时，插入质量值、缺失质量值、总质量值均较低，而替换质量值受到的影响则较小。在检测长同源多聚物上的替换型突变时，若按传统方法使用总质量值，会因质量值较低而无法检测较低频的突变。而若使用替换质量值，则会因质量值较高而成功检出较低频的突变。这一效果也不是显而易见的。
实施例2
比对是一个试图找到所测序列与参考序列之间最佳碱基对应关系的过程。在一个比对结果中，如果两个碱基之间是匹配的关系，则会有奖励分，而错配、插入和缺失则会造成罚分。所谓最佳碱基对应关系就是综合比对结果中所有奖励分和罚分之后，得分最高的比对结果。比对软件在计算得到比对结果的过程中，通常会有一套固定的罚分规则，这套规则仅和碱基配对关系有关，而和碱基质量值无关。除此之外，大部分比对软件往往还存在另一个问题，即它们通常是独立处理每条序列的比对的，即一条序列的比对结果和其他序列无关。这就造成当出现插入/缺失突变、或连续多个多种突变时，不同序列容易出现不同的比对结果，难以鉴定到底是何种突变。传统方法中，是在比对之后增加一个“重比对”(realignment)的步骤，具体过程参见图5，根据覆盖到同一位点的所有序列的比对情况，局部更正每条序列的比对结果，以增加对插入/缺失、或连续突变的鉴定准确度。
若使用本发明所描述的方法，每个碱基分别赋予多个不同的质量值，则可规避上述问题，不必进行重比对步骤，简化生物信息学流程，加快分析速度。具体实施方式为：比对软件不再采用固定的罚分规则，而使罚分依赖于碱基的不同质量值。即，若一个碱基的替换质量值较低，则在判定该碱基存在替换时赋予较低的罚分；若一个碱基的插入质量值较低，则在判定该碱基存在插入时赋予较低的罚分；若一个碱基的缺失质量值较低，则在判定该碱基存在缺失时赋予较低的罚分。
实施例3
1.对HBL101细胞系反复传代后提取基因组DNA，使用艾吉泰康遗传性乳腺癌Panel试剂盒构建其靶向测序文库。
2.对该文库进行ECC测序，并对测得的每个碱基赋予插入质量值、缺失质量值、替换质量值、总质量值四种质量值，分别得到4个fastq文件。
3.使用fastp软件对这4个fastq文件进行预处理，去除低质量数据，然后使用BWA-MEM软件将序列比对到hg19参考基因组上。
4.使用pysam统计比对结果的覆盖均一度和背景噪音，并用GATK鉴定文库中存在的基因变异。
5.分析结果如表3所示，在使用总质量值、插入质量值、缺失质量值的结果中，由于在长同源多聚物处的质量值相对较低，此处的序列大多被fastp软件去除，导致比对后在基因组上的覆盖均一度较低，只有0.6左右。而在使用替换质量值的结果中，长、短同源多聚物处的质量值持平，在fastp软件处理时均被保留，因此比对后在基因组上的覆盖均一度较高，可达到0.83。覆盖均一度的计算方式为覆盖深度的Lorenz曲线与x轴围成的面积。在GATK的鉴定结果中，分别比较不同变异类型的检出数量，由于共有4种质量值，在使用插入/缺失/替换质量值时，仅检测对应的插入/缺失/替换变异，不检测另外2种变异，使用总质量值时检测上述3种变异类型。鉴定结果如表3所示，使用替换质量值的结果比使用总质量值的结果可以多鉴定出20％的变异。
表3
| 质量值类型 | 覆盖均一度 | 插入变异数 | 缺失变异数 | 替换变异数 |
| 插入质量值 | 0.62 | 8 | 0 | 0 |
| 缺失质量值 | 0.61 | 0 | 5 | 0 |
| 替换质量值 | 0.83 | 0 | 0 | 30 |
| 总质量值 | 0.59 | 6 | 4 | 25 |
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
Claims (15)
1.一种碱基测序质量评估方法，其特征在于，包括：
1)提供标准核酸序列，在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；
3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
2.根据权利要求1所述的方法，其特征在于，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以三种质量值，分别表示碱基的插入错误率，缺失错误率，替换错误率。
3.根据权利要求1或2所述的方法，其特征在于，所述测序方法包括纠错码测序、纳米孔测序。
4.根据权利要求3所述的方法，其特征在于，所述碱基的测序信号特征指的是，测序过程中该碱基发生测序化学反应时产生的信号的特征，包括但不限于：碱基种类，碱基在序列上的位置，碱基所处多聚物的长度，碱基在其所处多聚物中的位置，碱基发生测序化学反应的轮数，信号强度，信号强度(及其邻近信号强度)接近整数的程度，测序信号的参数，碱基对应的失相程度，等。
5.根据权利要求1-4任一项所述的方法，其特征在于，所述分类器根据碱基的测序信号特征，将碱基分成若干类，统计每一类碱基的准确率。
6.根据权利要求5所述的方法，其特征在于，所述拟合碱基的测序信号特征与其标记之间的关系，包括将分类器的拟合结果转化为碱基的测序质量值。
7.根据权利要求6所述的方法，其特征在于，所述分类器包括线性回归，多项式回归，逻辑回归，支持向量机，人工神经网络，随机森林，Phred算法，集成学习等。
8.根据权利要求1-7任一项所述的方法，其特征在于，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；所述标准核酸样品，指的是来源和序列均已确定、在基因组的几乎所有位点上均高度纯合的核酸样品，包括λ
噬菌体DNA，大肠杆菌DNA，酿酒酵母DNA等。
9.根据权利要求1所述的方法，其特征在于，还包括：6)对赋予质量值后的核酸序列进行生物信息学分析。
10.根据权利要求9所述的方法，其特征在于，所述生物信息学分析包括，鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入变异时，使用表征插入错误率的质量值来计算变异存在的概率；当鉴定缺失变异时，使用表征缺失错误率的质量值来计算变异存在的概率。
11.根据权利要求9所述的方法，其特征在于，所述生物信息学分析包括，鉴定基因变异：当鉴定替换变异时，使用表征替换错误率的质量值来计算变异存在的概率；当鉴定插入缺失变异时，使用表征插入和/或缺失错误率的质量值来计算变异存在的概率。
12.根据权利要求1所述的方法，其特征在于，还包括序列比对，比对时的罚分依赖于单个碱基的至少两种质量值。
13.根据权利要求1-12任一项所述的方法，其特征在于，步骤1)-3)预先完成，并将训练好的分类器作为配置文件存于系统中，在执行步骤4)-5)时调取所述配置文件。
14.一种用于评估核酸测序碱基识别的准确度的系统，其包括：
机器可读存储器；和
被配置成执行机器可读指令的处理器，所述指令在由所述处理器执行时使
所述系统执行一方法，所述方法包括：所述处理器接收
1)标准核酸序列，所述标准核酸序列是对标准核酸样品测序得到的碱基序列；
在所述标准核酸序列中，以碱基为基本单元，计算碱基的测序信号特征；
2)将所述标准核酸序列比对到参考序列上，根据比对结果对标准核酸序列中的碱基赋予标记，所述标记选自测序正确、插入错误、缺失错误或替换错误；
3)训练分类器以拟合碱基的所述测序信号特征与其标记之间的关系；
4)提供待测核酸序列，在所述待测核酸序列中，以碱基为基本单元，计算碱基的所述测序信号特征；
5)利用训练好的分类器，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以至少两种质量值，分别表示该碱基的插入和/或缺失错误率，替换错误率。
15.根据权利要求14所述的系统，其特征在于，根据第4)步中计算所得的测序信号特征，赋予所述待测核酸序列中的碱基以三种质量值，分别表示碱基的插入错误率，缺失错误率，替换错误率。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310322408.XA CN116434843A (zh) | 2023-03-29 | 2023-03-29 | 一种碱基测序质量评估方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310322408.XA CN116434843A (zh) | 2023-03-29 | 2023-03-29 | 一种碱基测序质量评估方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116434843A true CN116434843A (zh) | 2023-07-14 |
Family
ID=87091875
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310322408.XA Pending CN116434843A (zh) | 2023-03-29 | 2023-03-29 | 一种碱基测序质量评估方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116434843A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117594130A (zh) * | 2024-01-19 | 2024-02-23 | 北京普译生物科技有限公司 | 纳米孔测序信号评价方法、装置、电子设备和存储介质 |
-
2023
- 2023-03-29 CN CN202310322408.XA patent/CN116434843A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117594130A (zh) * | 2024-01-19 | 2024-02-23 | 北京普译生物科技有限公司 | 纳米孔测序信号评价方法、装置、电子设备和存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| NZ759659A (en) | Deep learning-based variant classifier | |
| CN112020565A (zh) | 用于确保基于测序的测定的有效性的质量控制模板 | |
| AU2023251452A1 (en) | Validation methods and systems for sequence variant calls | |
| US20190206510A1 (en) | Validation methods and systems for sequence variant calls | |
| GB2620315A (en) | Determination of base modifications of nucleic acids | |
| US20210065847A1 (en) | Systems and methods for determining consensus base calls in nucleic acid sequencing | |
| CN110016497B (zh) | 一种检测肿瘤单细胞基因组拷贝数变异的方法 | |
| CN111321209A (zh) | 一种用于循环肿瘤dna测序数据双端矫正的方法 | |
| CN116434843A (zh) | 一种碱基测序质量评估方法 | |
| CN108460248B (zh) | 一种基于Bionano平台检测长串联重复序列的方法 | |
| CN114220487A (zh) | 一种新型9基因risk急性髓系白血病预后模型的构建方法 | |
| CN116246703A (zh) | 一种核酸测序数据的质量评估方法 | |
| CN108728515A (zh) | 一种使用duplex方法检测ctDNA低频突变的文库构建和测序数据的分析方法 | |
| CN114420214A (zh) | 核酸测序数据的质量评估方法和筛选方法 | |
| CN114107525B (zh) | 一种铜绿假单胞菌的mnp标记位点、引物组合物、试剂盒及其应用 | |
| WO2024007971A1 (en) | Analysis of microbial fragments in plasma | |
| CN114790487B (zh) | 一种霍氏鲍特菌的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN114790488B (zh) | 一种金黄色葡萄球菌的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN115029454B (zh) | 一种卡他莫拉菌的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN114790494B (zh) | 一种水痘-带状疱疹病毒的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN114790485B (zh) | 一种不动杆菌属的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN114214435B (zh) | 一种肺炎支原体的mnp标记组合、引物对组合、试剂盒及其应用 | |
| CN114790493B (zh) | 一种单纯疱疹病毒的mnp标记位点、引物组合物、试剂盒及其应用 | |
| CN114836574B (zh) | 一种流行性腮腺炎病毒的mnp标记位点、引物组合物、试剂盒及其应用 | |
| WO2017136606A1 (en) | Apparatus, systems, and methods for dna amplification with post-sequencing data filtering and cell isolation |
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 |