一种碱基测序质量评估方法

技术领域

本发明涉及一种碱基测序质量评估方法及系统，属于基因测序领域。

背景技术

核酸测序技术可以探明遗传物质的序列，被广泛应用于临床肿瘤分型、微生物鉴定和遗传病诊断等领域。当今主流的核酸测序技术除了产出被测核酸样品的序列之外，还会给所测得的每个碱基以一个质量值，用来评估其所测的准确性。这个质量值一般以Phred的形式表示：

q＝-10log

式中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

简并多聚物长度(Degenerate

简并多聚物即由简并碱基构成的多聚物，例如字母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

本申请中，纠错码测序具备如下特征：

该测序方法需要多回测序，每回测序得到的信息不完整，而多回测序得到的总的信息是冗余的；利用多回测序的信息冗余来检测和校正潜在的测序错误，得到高准确度的序列。例如的，以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倍大、或至少约10

变异

是指与核酸参考序列不同的核酸序列。典型的变异包括但不限于单核苷酸变异(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

可见，总质量值最高只达到了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

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。