确定读段质量的方法和测序方法

技术领域

本发明涉及数据处理技术领域，具体的，涉及核酸检测数据处理领域，更具体的，涉及确定读段质量的方法及其应用。

背景技术

二代基因测序，也称为高通量测序或者平行大规模测序，广泛使用Q30、Q20等参数来评价测序序列(读段)的质量，该质量分数能够反映测序序列的错误概率以及可靠程度，能用于区分测序错误造成的数据波动还是生物样本本身带来的数据波动[Mbandi,S.K.,etal(2014),A glance at quality score:implication for de novo transcriptomereconstruction of Illumina reads,Frontiers in genetics,5,17.；Garber,M.,et al(2011),Computational methods for transcriptome annotation and quantificationusing RNA-seq,Nature methods,8(6),469-477.；Kelley,D.R.,et al(2010),Quake:quality-aware detection and correction of sequencing errors,Genome biology,11(11),R116.；US8965076B2；US2018051329A1]。

各种测序原理、各类测序平台经常都有相适配或适用的测序序列质量定量评估方案，更可靠和/或更加普适的测序序列的质量评价方法仍是值得关注的问题，特别是针对或者同样适用于非主流测序平台的测序数据的质量评估方法，例如适用于三代测序、单分子测序的测序序列质量的定量评估，或者同时适用于二代测序、三代测序、四代测序、单分子测序、单分子边合成边测序等多种平台的测序序列质量的定量评估方案等。

发明内容

本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。

本申请基于发明人的下列发现和试验测试而作出：

在二代基因测序技术中，通常涉及扩增来放大目标信号，一般使用Q30、Q20等参数评价测序序列的质量来反映测序序列的错误概率和可靠程度，从而助力区分开来源于检测错误的数据波动和来源于生物样本的数据波动[Mbandi,S.K.,et al(2014),A glance atquality score:implication for de novo transcriptome reconstruction ofIllumina reads,Frontiers in genetics,5,17.；Garber,M.,et al(2011),Computational methods for transcriptome annotation and quantification usingRNA-seq,Nature methods,8(6),469-477.；Kelley,D.R.,et al(2010),Quake:quality-aware detection and correction of sequencing errors,Genome biology,11(11),R116.；US8965076B2；US2018051329A1]。然而，将这个质量评估方案、质量分数计算方法应用于定量评估单分子测序所产生的数据，如用来评估真迈生物、太平洋生物、牛津纳米孔等公司的单分子测序平台所产生的测序数据的质量，表现出不适用或者说不能真实地反映出这些测序平台的测序序列的质量，而且，针对单分子测序的测序数据发明人目前还未见有已公开的评估方法。

单分子检测，一般指不涉及克隆如利用扩增来放大目标分子的信号的检测。例如单分子测序，一般指不涉及对待测核酸分子进行扩增的测序。可以理解地，单分子检测的目标信号很弱或者说信噪比很低，很容易由于噪声或杂质信号的干扰而丢失或难被检出。因此，单分子测序所产生的测序数据的错误率通常较高，甚至包含非样本来源的杂质或干扰或噪声序列；这些错误率较高和/或非样本来源的序列会占用大量的数据处理或分析资源例如增加比对计算所需的时间，并且即便利用较佳或较适配的参考序列进行比对，很多情况下根据比对结果也难以区分测序错误带来的数据波动和生物样本/待测样本本身特征造成的数据波动。例如，这些序列中的一部分有可能随机比对到某参考序列如人参考基因组上，在样本情况复杂、涉及不同物种来源的测序中，如检测人源样本中的病原体，这些随机比对上的测序序列可能对不同物种测序序列的拆分以及后续的生信分析造成不良影响或干扰。

鉴于此，本发明的实施方式提供一种确定读段质量的方法，包括：获取来自基于芯片成像检测、对芯片上的模板进行边合成边测序的读段，读段由读取轮数据和轮空轮数据构成，读段的读取轮数据反映对应于图像集的指定位置的模板上存在碱基延伸信号、轮空轮数据反映对应于图像集的指定位置的模板上不存在碱基延伸信号，图像集包括测序运行期间的一轮或多轮测序产生的、包含指定位置及其周围背景信号的多个图像；基于图像集、读取轮数据和轮空轮数据中的至少一类信息，确定预定特征组合中每个特征的值，预定特征组合中的每个特征为与该读段的产生有关的、可量化的特征；基于预定特征组合中的特征的值和预先训练的预测模型，预测读段被分类为指定类别读段的概率；和基于概率，确定读段的测序质量，预先训练的预测模型为能够基于读段的预定特征组合中的特征的值计算该读段属于指定类别读段的概率的分类器，读段的质量得分与该读段被分类为指定类别读段的概率正相关。

本发明的另一实施方式提供一种确定读段质量的方法，所称的读段来自基于芯片成像检测、对芯片上的模板进行边合成边测序而确定的序列，读段由读取轮数据和轮空轮数据构成，每条读段的读取轮数据反映对应于图像集的指定位置的模板上存在碱基延伸信号、轮空轮数据反映对应于图像集的指定位置的模板上不存在碱基延伸信号，图像集包括测序运行期间的一轮或多轮测序产生的、包含指定位置及其周围背景信号的多个图像，该方法包括：(1)针对读段，基于图像、读取轮数据和轮空轮数据中的至少一类信息获取预定特征组合中的每个特征的取值；(2)针对读段，基于特征的取值，确定与取值对应的比对系数和非比对系数；(3)基于比对系数、非比对系数以及预定的比对常数和预定的非比对常数，确定读段的质量得分，其中，预定特征组合包括测序运行期间产生的、与读段有关的多个特征。

本发明的又一实施方式提供一种确定读段质量的方法，所称的读段来自基于芯片成像检测、对芯片上的模板进行边合成边测序而确定的序列，读段由读取轮数据和轮空轮数据构成，每条读段的读取轮数据反映对应于图像集的指定位置的模板上存在碱基延伸信号、轮空轮数据反映对应于图像集的指定位置的模板上不存在碱基延伸信号，图像集包括测序运行期间的一轮或多轮测序产生的、包含指定位置及其周围背景信号的多个图像，该方法包括：(1)针对读段，基于图像、读取轮数据和轮空轮数据中的一类信息获取预定特征的取值；和(2)基于预定特征的取值以及预先构建的特征-质量关系，确定与读段的质量得分，其中，预定特征包括测序运行期间产生的、与读段有关的多个特征，预先构建的特征-质量关系是基于训练读段集确定的，训练读段集包含指定类别读段和非指定类别读段以及这两类读段中的各读段的特征的取值。

在一些实施例中，对于给定的一条读段是否属于指定类别读段，通常存在“是”或“否”两种结论。因此，可以理解地，所称的预先训练的预测模型或者预先构建的特征-质量关系实质为分类器，该分类器可以根据读段的预定特征确定该读段属于指定类别读段的概率的大小，进而来确定该读段的质量得分。也可以理解地，预先训练的预测模型或者预先构建的特征-质量关系为关联了特征与读段属于指定类别读段的概率的量化方案，与所称的预先训练的预测模型或者预先构建的特征-质量关系为关联了特征和读段不属于指定类别读段的概率的量化方案，这两种表述指向相同的技术意思。

本发明的实施方式还提供一种测序方法，包括：采用基于芯片成像检测的边合成边测序技术，对芯片上的模板进行测序，获取图像，并且分析图像以确定测序结果，测序结果包括多个读段；和根据上述任一实施方式的确定读段质量的方法，确定各读段的质量得分。

本发明的实施方式还提供一种计算设备，包括：处理器和存储器；存储器，用于存储计算机程序；处理器，用于执行计算机程序以实现上述任一实施方式的确定读段质量的方法。

本发明的实施方式还提供一种计算机可读存储介质，存储介质包括计算机指令，当指令被计算机执行时，使得计算机实现上述任一实施方式的确定读段质量的方法。

本发明的实施方式还提供一种测序装置，包括上述任一实施方式的计算设备；或者上述任一实施方式的计算机可读存储介质。

采用上述任一实施方式、实施例的方法、设备或装置，能够有效地对利用边合成边测序所得到的测序数据的质量进行准确测定，该方法既适用于二代或高通量测序平台产生的测序数据，也适用于单分子测序平台产生的数据；适用于单色或多色的基于芯片成像检测实现边合成边测序的测序平台所产生的测序数据的质量评估。

上述任一实施方式、实施例的方法根据读段被分类为指定类别读段的概率来确定该读段的质量，可以准确反映该读段的准确性、错误率或可靠程度。如此，根据读段的质量得分评估该次测序运行所产生的测序数据的质量，可以在不需另外的比对、根据比对结果来评估测序数据的质量的情况下，直接根据确定出的质量分值识别出一部分很可能是高错误率或者非样本来源的序列，知晓该次测序运行所产出的数据的质量水平，包括能够是否满足目标应用检测的要求，例如是否满足、其中有多少满足目标检测应用对测序数据质量的要求等；而且，根据质量得分对测序数据进行处理，能够去除掉一部分高错误率序列和/或非样本来源的杂质序列，从而能够提升下游生信分析的速度、降低高错误序列和/或杂质序列对后续生信分析的影响。

在某个示例中，利用下机获得的测序数据进行遗传变异检测，可以在将该些测序数据比对到参考序列获得比对结果之前，根据测序数据中的各读段的质量分数对该些读段进行筛选，如此，能够大幅减少比对和后续分析所需消耗的时间。而且，去除低质量分值或高错误率(大概率不能比对到目标物种参考序列)的读段，也能减少错误序列或非目标物种测序序列对后续分析结果的干扰。例如，在一次测序运行中获得多个样本的混合测序序列，过滤掉质量分数低于预定阈值的读段有利于测序数据准确拆分对应至相应的样本，也利于基于拆分后的测序数据对相应样本进行准确检测。又例如，在待测样本情况复杂、涉及不同物种来源的测序中，如基于测序数据检测人源样本中的病原体，基于带有量化的质量评估分值的读段，去除掉高概率可比对到人参考序列的读段、和/或低概率可比对到目标病原参考序列的读段，可以大大地提升后续序列分析的效率以及病原检测结果的准确性，有效降低人源或低质量病原体来源读段对检测结果造成的不良影响。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本申请实施例的确定读段质量的方法的流程示意图；

图2是根据本申请实施例的确定读段质量的方法的流程示意图；

图3是根据本申请实施例的确定读段质量的方法的流程示意图；

图4是根据本申请实施例，每条序列中读到的碱基在图像中对应模板位置处的平均荧光亮度分布在比对读段和非比对读段之间的分布图；

图5显示了根据本申请一个实施例，人类gDNA混样测试集数据中实际比对读段与实际非比对读段的质量分数值分布；

图6显示了根据本申请一个实施例，人类gDNA混样测试集数据在不同的质量分数过滤阈值之下，对应的比对率及序列过滤和损失比例；

图7显示了根据本申请一个实施例，人类gDNA混样测试集数据中质量分数与对应的唯一比对读段的平均错误率的关联；

图8显示了根据本申请一个实施例，大肠杆菌测序数据中实际比对读段与实际非比对读段的质量分数值分布；

图9显示了根据本申请一个实施例，大肠杆菌测序数据在不同的质量分数过滤阈值之下，对应的比对率及序列过滤和损失比例；

图10显示了根据本申请一个实施例，大肠杆菌测序数据中质量分数与对应的唯一比对读段的平均错误率的关联；

图11显示了根据本申请一个实施例，不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后非比对读段过滤比例的统计结果；

图12显示了根据本申请一个实施例，不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后比对读段损失比例的统计结果；和

图13显示了根据本申请一个实施例，不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后比对率的统计结果。

具体实施方式

下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

在本文中，除非另有说明，以单数形式“一种”、“一个”等限定的方式包括其复数指代物(一个以上)；“一组”或者“多个”指两个或两个以上。

在本文中，除非另有说明，术语“第一”、“第二”、“第三”、“第四”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。

在本文中，除非另有说明，核苷酸指四种天然核苷酸(如dATP、dCTP、dGTP和dTTP或者ATP、CTP、GTP和UTP)或其衍生物，有时也直接以其包含的碱基(A、T/U、C、G)表示。本领域普通技术人员根据上下文记载可以知晓核苷酸或碱基在具体实施方式中的指代。

在本文中，除非另有说明，单链核酸分子或双链核酸分子，包括所称的插入片段、核酸片段、序列、位点、多核苷酸、接头、引物/探针等，均是以5'至3'方向从左到右书写表示的。

在本文中，除非另有说明，“连接”、“固定”等应做广义理解，例如，可以是固定连接，也可以是可逆连接，可以是直接相连，也可以通过中间媒介间接相连，可以是化学键如氢键、磷酸二酯键等。

在本文中，除非另有说明，接头(adapter，也称为衔接子或适配子)、引物或者探针，为序列预设或者说序列已知的寡核苷酸片段，通常而言，接头为单链或双链核酸分子，引物或探针为寡核苷酸单链；对于目前市面上的主流测序平台中，通常是通过处理使来自样本的待测核酸片段(也称为插入片段)的末端带上预设序列(接头)，以及利用能够与接头的至少一部分互补或结合的引物或探针(寡核苷酸链)将待测片段连接或固定到反应器指定位置(如流动池flow cell或芯片的指定表面上固定有该引物或探针)。基于碱基互补原则，接头的至少一部分序列可以用于设计引物/探针，可以作为特定引物/探针的结合位点。

在本文中，所称的“测序”为序列测定，同“核酸测序”或“基因测序”，指核酸序列中碱基次序的测定；包括合成测序(边合成边测序，SBS)和/或连接测序(边连接边测序，SBL)；包括DNA测序和/或RNA测序；包括长片段测序和/或短片段测序，所称的长片段和短片段是相对的，如长于1Kb、2Kb、5Kb或者10Kb的核酸分子可称为长片段，短于1Kb或者800bp的可称为短片段；包括双末端测序、单末端测序和/或配对末端测序等，所称的双末端测序或者配对末端测序可以指同一核酸分子的不完全重叠的任意两段或两个部分的读出。

测序可以通过自动化测序平台进行，根据本申请的实施例，可选择的自动化测序平台包括但不限于Illumina公司的Hiseq、Miseq、Nextseq或Novaseq测序平台、ThermoFisher/Life Technologies公司的Ion Torrent平台、华大基因的BGISEQ和MGISEQ/DNBSEQ平台以及单分子测序平台；测序方式可以选择单端测序、双末端测序或者选择所选用的自动化测序平台可支持的测序方式等。

测序测读出来的序列称为测序序列，也称为读段(reads)，测序序列或读段的长度称为读长。

在某些示例中，利用边合成边测序进行多轮测序来获得测序序列或读段。例如，使待测核酸分子和聚合酶、改造的核苷酸接触并置于适合聚合反应的条件下，可控地使改造的核苷酸掺入到待测核酸分子中或者说可控地实现单碱基延伸，并检测相应的反应信号，基于此信号确定该次反应掺入到待测核酸分子的核苷酸的类型，如此进行多次可控的单碱基延伸和相应信号的检测，以便据反应信号信息检出多次或多轮反应中掺入到待测核酸分子的核苷酸或碱基的类型，以测读出待测核酸分子的一部分序列。

所称的待测核酸分子，也称为待测模板或模板，可以是未经扩增的单分子，也可以是扩增后的包含多个相同多核苷酸分子的分子簇或长链，如主流测序平台采用的经桥式扩增或滚环扩增所形成的克隆簇cluster或DNA纳米球DNB。待测核酸分子可呈现为单链、双链和/或杂交有探针或引物的复合物形态。

相应的反应信号例如可以为荧光信号，也可以转化为采集该些荧光信号所形成的图像数据，由此，处理和分析该些图像数据以检出各次或各轮反应掺入到待测核酸分子的核苷酸，以便确定待测核酸分子的一部分碱基序列。

具体地，在某些示例中，基于表面荧光成像检测实现测序，待测核酸分子连接在固相表面，例如，可改造核苷酸使核苷酸带有或者可结合荧光标记、以及带有可阻止其它核苷酸聚合连接到待测酸分子下一个位置的可切除的抑制基团(该种改造的核苷酸也称为可逆终止子reversible terminator)，并且在每次聚合反应或单碱基延伸反应完成后激发使荧光标记发光，以及采集该些发光信号以获取指定表面位置上发生单碱基延伸反应的待测核酸分子的图像；接着，去除抑制基团和荧光标记等以进行下一次或下一轮聚合反应和信号采集(拍照)，如此重复进行多次或多轮的聚合-拍照-切除，来获得关联于各次单碱基延伸反应连接到待测核酸分子的核苷酸相关的图像集信息。

可以理解地，表面指定位置上的待测核酸分子发生聚合反应发出荧光，一般会在该轮反映所采集到的图像的相应位置上表现为高于背景信号强度的亮点或亮斑(spots)。因此，根据该些图像集包括其中的对应于特定化学特征(发生聚合反应的待测核酸分子)的亮斑的信息，可判断指定位置上的待测核酸分子是否发生聚合反应，结合预设的可区分的荧光发光信号和核苷酸种类的对应关系，可检出发生聚合反应聚合连接到待测核酸分子的核苷酸的类型，由此测定出待测核酸分子的至少一部分序列，获得所称的读段。

需要说明的是，本文所称的核苷酸，包括核糖核酸或脱氧核糖核酸，包括天然核苷酸或其衍生物或其改造物(也称为改造的核苷酸或修饰的核苷酸等)。本文中，有时也以核苷酸所包含的碱基来指代该种核苷酸，本领域技术人员根据常规认识和/或上下文能够清楚知晓。

在本文中，除非另有说明，术语“图像”是指在延伸反应过程中，每一轮延伸反应后，激发荧光信号获得的记录荧光信号的“图像”。这些图像将反映各个位置荧光信号的类型和强度。需要说明的是，按照一次延伸反应的反应混合物中所添加碱基类型的数目，每一轮延伸反应所获得的图像数据将有不同，通常而言，针对不同颜色的荧光信号，可以获取不同的图像。换句话说，本领域技术人员可以理解，在每次延伸反应后，可以进行多次拍照，可以针对不同颜色的荧光采用不同检测通道进行拍照检测。由此，每个测序轮，为了确定延伸碱基的类型，需要收集至少一张图像，甚至多张图像，例如，一个碱基延伸反应可能包含一次图像采集，也可能包含多次图像采集，这些图像的集合作为后续进行测序分析的图像集。

在本文中所使用的术语“位置”是指在测序过程中，为了将荧光信号与靶核酸分子对应，会将固定有靶核酸的位点进行定位不同的坐标，具体的，在边合成边测序的过程中，在每个延伸反应中，所延伸的碱基携带光学可检测标记例如荧光标记，因此，通过激发荧光，可以使得至少一部分核酸分子或克隆簇在图像上表现为亮斑。这些“亮斑”也称为“亮点”(spots或peaks)，在图像上会呈现为发光点，一个发光点占有至少一个像素点。由此，对于一个“给定位置”(也称为“模板位置”)，将与特定的靶核酸分子或者片段匹配。由此，可以通过位置，将所得到的一系列图像集进行分类，以便分别获得各靶核酸分子所对应的测序读段。

在本文中，除非另有说明，术语“预定特征组合”是指本发明的发明人经过大量研究获得的可以用于采用机器学习模型预测测序读段的测序质量的一系列特征，这些特征是可以从相关的图像集、读取轮数据和轮空轮数据获得的。根据本申请的实施例，这些“特征”是可以量化的特征，由此，每一个具体的测序读段都可以产生一系列具体的特征值。根据本申请的实施例，可以采用的特征包括但不限于下列：读长(rl)、平均轮空长度(emp_cyc)、平均连续读取轮数(rep_RG)、平均亮度(int_mean)、重复比例(rep_perc)、平均信噪比(snr_mean)、模板距离(temp_dist)。为了方便理解，下面采用读段_GT__G_ACG_A___AC_T_C_T___T____作为示例，对相关参数的确定方法进行举例如下：

1)读长(rl)：表示每条序列读到的碱基数。经计算，rl＝13；

2)平均轮空长度(emp_cyc)：表示一条测序序列上两个读出的碱基之间轮空轮数的平均值，其中头尾的轮空不算。经计算，emp_cyc＝(0+2+1+0+0+1+3+0+1+1+1+3)/12＝13/12。

例如其中的部分测序序列：_GT__G_ACG_，经计算，emp_cyc＝(0+2+1+0+0)/5＝0.6。

3)平均连续读取轮数(rep_RG)：表示将一条测序序列根据不同荧光基团的颜色拆分为多个部分，分别统计不同荧光基团所对应的连续读取的碱基数，然后将不同荧光基团所对应的连续读取的碱基数计算获得平均值，即为连续读取轮。例如，上述测序中每轮均含有携带红色荧光基团和绿色荧光基团两种碱基，其连续读取轮的计算是通过将上述获得的测序序列拆成红色和绿色两部分，分别统计红色和绿色对应的连续读取的碱基数，然后将红色和绿色的连续读取碱基数计算平均值。

其中，绿色碱基为：_T__C___CTCT_T__；绿色碱基连续亮的个数：[1,1,4,1]；

红色碱基为：G_GAGA_A_______；红色碱基连续亮的个数：[1,4,1]；

经计算，rep_RG＝(1+1+4+1+1+4+1)/7＝13/7。

4)平均亮度(int_mean)：每条测序序列中读出的碱基对应的原始图像中模板处的亮度值。进一步地，计算亮度之前可先用算法优化过图像质量。

5)重复比例(rep_perc)：针对一条测序序列中，在同一轮里不同荧光基团的均被观测到的轮数除以该序列的总读长。例如，上述测序中每轮均含有携带红色荧光基团和绿色荧光基团两种碱基，重复比例为在同一轮里红色碱基和绿色碱基均被观测到的轮数除以该序列的总读长。

其中，绿色碱基为：_T__C___CTCT_T_，红色碱基为：G_GAGA_A_______，红绿同时亮的次数为：1，序列的总读长为13；经计算，rep_perc＝1/13。

6)平均信噪比(snr_mean)：每条测序序列中读到的碱基对应的原始图像中模板位置(或者称为给定位置)处信噪比的平均值。

7)模板距离(temp_dist)：给定位置在对应的红色和绿色图片上的位置距离(单位为像素值)。若某模板仅能在单色图上能找到，无法对应到另一种颜色的图上，则将模板距离标注为预先确定的固定数值(例如0或其他数值)。

在本文中，除非另有说明，术语“比对读段”是指测序反应得到能够与参考序列比对匹配上的读段，即在得到读段后的分析中，将所得到的读段与目标物种的参考基因组进行比对，读段能够比对上，这意味着该读段的质量比较高，可信度高。在测序过程中，尤其是单分子测序过程中，单分子测序系统不依赖样本的扩增，因而能够很好地避免由于扩增而导致的序列偏好性，进而提升下游生信分析的灵敏度。然而无扩增的测序体系意味着有效的测序信号均来源于单个DNA分子，从而单分子测序容易受到噪音及杂质信号的干扰，也容易出现信号丢失的情况。因此，在原始的测序数据中，会有部分错误率较高的序列，甚至是并非来源于样本的杂质序列，这些高错误率或者非样本来源的序列将占用大量测序仪比对计算的时间。通过本发明的方法，能有效地对边合成边测序所产生的读段进行质量评价，且可根据序列的质量评价过滤掉一部分高错误率序列及非样本来源的杂质序列，从而提升下游生信分析的速度、降低高错误序列及杂质序列对后续生信分析的影响，并且具有比较广的适用性。

在本文中，除非另有说明，“预测模型”与“机器学习模型”为一种计算模型或算法，它能够通过对输入数据进行学习和分析，从而自动地进行预测、分类、识别或决策等任务。该模型的学习过程基于统计学原理和数据模式识别，使用训练数据集进行参数调整和模型优化，以提高其预测或推理能力。机器学习模型可以采用各种算法和技术，例如神经网络、支持向量机、决策树、随机森林、深度学习等。这些模型可以通过监督学习、无监督学习或强化学习的方式进行训练和优化。在实际应用中，机器学习模型可用于各种领域，如自然语言处理、图像识别、模式识别、数据挖掘、推荐系统、预测分析等。它在处理大规模数据、自动化决策和智能化系统中具有重要的应用潜力。然而，需要说明的是，在机器学习模型的具体应用中，需要对用于预测的特征和模型等进行深入研究，才能够获得比较令人满意的预测效果，否则会出现各种问题，例如过拟合问题、欠拟合问题、泛化能力较差的问题等。本申请的发明人经过多年的研究结合研究经验，意外发现了一种特征取值，能够结合机器学习模型用于预测读段的质量分数。

在本文中，“模型参数”是指机器学习模型中的可调节参数。通过调整这些参数，改变模型的预测能力和泛化能力，从而得到更好的模型。其中，调整模型参数的过程通常称为模型调优或超参数调优。

本文选择决策树作为机器学习模型，其主要模型参数包括：最大深度(MaxDepth)、最小样本数(Min Samples Split)、特征选择标准(Criterion)等。

在本文中，“训练集(Training Set)”是指用来训练机器学习模型的数据集。在训练阶段，机器学习模型使用训练集中的样本来学习和调整模型的参数，以最小化预测结果与真实标签之间的误差。通过不断地尝试不同的参数和算法，模型逐渐学习到样本之间的规律和特征，从而得到更准确的预测能力。

在本文中，术语“测试集(Test Set)”是用来评估机器学习模型性能的数据集。而且该测试集中的样本在模型的训练阶段没有被用到过，这样可以确保测试集的数据对模型是未知的。在训练完成后，使用测试集来进行模型的性能评估，即将模型对测试集中的样本进行预测，并与真实标签进行比较。通过比较预测结果和真实标签，评估模型在未知数据上的表现，判断其泛化能力和预测准确性。

下面参考附图对根据本申请第一方面的用于确定读段的测序质量的方法进行详细描述。

参考图1，根据本申请的实施例，本申请提出了一种用于确定读段的测序质量的方法，其包括：

S100：获取来自基于芯片成像检测、对芯片上的模板进行边合成边测序的读段。

在该步骤中，获得边合成边测序的读段的方式并不受特别限制。

本领域技术人员能够理解，通常而言，核酸测序包括使核苷酸(包括核苷酸类似物)结合到模板并采集相应的延伸碱基荧光信号的过程。在一些平台中，核苷酸结合到模板和采集相应的荧光信号是非同步或者实时进行的，通常需要多轮测序反应来实现针对模板核酸分子的核苷酸序列或者碱基序列的测定，其中，可以将延伸一个碱基并确定所延伸碱基类型的过程称为一轮测序(cycle)，换句话说，在一轮测序中，完成了针对模板核酸分子上给定位置的碱基类型的测定过程，该方法也被称为边合成边测序(SBS)。根据本申请的实施例，可以将四种核苷酸进行任意组合来实现通过至少一次延伸反应和反应信号的采集来完成所延伸碱基类型的确定。例如，每次延伸反应中添加一种、两种、三种或者四种类型的核苷酸，分别进行延伸反应和反应信号的采集，从而完成延伸反应和延伸碱基类型的确定。

需要说明的是，在一个示例中，所称的一轮测序包括一次碱基延伸反应。例如，可以使四种不同的核苷酸(dATP、dTTP、dGTP和dCTP)与多个待测核酸分子置于相同的聚合反应体系中，使每种核苷酸可被激发出可区分于其他种类核苷酸的信号，从而通过一次重复反应获得的信息确定多个待测分子的一个位置上所掺入或引入的核苷酸的类型。例如，使四种核苷酸分别带有四种不同发光波段的荧光标记，进行四色或四通道单分子测序或高通量测序；又例如，使四种核苷酸分别带有荧光标记F1、荧光标记F2、荧光标记F1及F2以及不带有荧光标记进行双色或双通道高通量测序等，荧光标记F1和F2具有可区分的发光波长，从而通过一次碱基延伸反应获得的信息确定多个待测分子的一个位置上所掺入或引入的核苷酸的类型。

在某些示例中，一轮测序包括使模板与四种核苷酸接触一次进行单碱基延伸以及采集相应的图像以识别碱基延伸信号的过程，经过一轮测序在读段上表现为增加了对应于四种核苷酸的四个数据位。可以理解地，对于单色荧光SBS测序或者多色荧光SBS测序，一轮测序表现在一条读段上的四个数据位必定包含轮空轮数据(位)。

在另一个示例中，所称的一轮测序也可以包括多次碱基延伸反应。例如，可以使四种核苷酸依次与表面上的多个待测核酸分子接触，分别进行碱基延伸和相应的反应信号的采集如拍照，一轮测序包括四次碱基延伸；又例如，使四种核苷酸任意组合与表面上的多个待测核酸分子接触，例如两两组合或者一三组合，两个组合分别进行碱基延伸和相应的反应信号的采集如拍照，一轮测序包括两次碱基延伸。

在某些示例中，一条读段包含读取轮数据和轮空轮数据，所称的读取轮数据反映对应于图像集的指定位置的模板(待测核酸分子)上存在碱基延伸信号、所称的轮空轮数据反映对应于该图像集的该指定位置的模板上不存在碱基延伸信号。所称的存在或不存在，可以是物理意义上的存在或不存在、或者说发生或未发生碱基延伸，也可以指借助本领域常规手段能够或已经识别出、检出或者不能够或难以识别出、检出碱基延伸信号；一种可实现的实施例是，所称的读取轮数据和轮空轮数据为前者，为相对的发生碱基延伸和没有发生碱基延伸；在另一种可实现的实施例是，所称的读取轮数据和轮空轮数据为后者，为相对的检出特定碱基和没有检出特定碱基。

通常地，给定的一条读段一般由四种碱基A、T/U、C和G组成，有时还包含不能确定的碱基N或gap(N或gap选自A、T/U、C和G中的一种)。该条读段包含的碱基的总数为该读段的长度(该读段的读长)。

根据读段的常见表现形式为一条由碱基组成的序列，可以理解地，常见的给定的一条读段只体现所称的读取轮数据，一般不包含所称的轮空轮数据，或者说，一般不会专门体现轮空轮数据。

在一些示例中，为了反映利用边合成边测序测定表面指定位置上的待测核酸分子的一个位置的碱基的过程，包括其与四种碱基或碱基的组合依次或同时接触过才能确定该位置的碱基类型的过程，将读段呈现为由读取轮数据和轮空轮数据组成的形式。

例如，读段可表示为包含多个连续的重复单元，每个重复单元对应待测核酸分子上的一个待测(碱基)位置，每个重复单元包含四个数据位，四个数据位分别对应四种以预定顺序排列的碱基比如ACTG；倘若某个数据位上存在特定的碱基检出，如该数据位显示为相应的特定碱基检出A、C、T或G，则该数据位为所称的读取轮数据，倘若该数据位不存在相应的碱基检出，则该数据位为所称的轮空轮数据。

可以理解地，对于通过边合成边测序的多轮反应，假设在每一轮反应中，待延伸的碱基都能成功发生延伸，测读出的任意一条以常规形式呈现的只包含读取轮数据的读段，都可以将该读段转化成包含读取轮数据和轮空轮数据的读段。

例如，对于给定的一条读段AATGCGCCGT，该读段通过四色荧光单分子测序或者四色荧光高通量测序获得(一种荧光信号对应一种碱基)，按照上面给定顺序排列的四种碱基ACTG(一个重复单元)，可以将该读段转化表示为A___(左侧以下划线表示的三个数据位依次表示未检到的C、T和G)A_____(左侧以下划线表示的五个数据位依次表示未检到的C、T、G、A和C)T____(左侧以下划线表示的四个数据位依次表示未检到的G、A、C和T)G_(左侧以下划线表示的一个数据位表示未检到的A)C_____(左侧以下划线表示的五个数据位依次表示未检到的T、G、A、C和T)G_(左侧以下划线表示的一个数据位表示未检到的A)C___(左侧以下划线表示的三个数据位依次表示未检到的T、G和A)C_____(左侧以下划线表示的五个数据位依次表示未检到的T、G、A、C和T)G__(左侧以下划线表示的两个数据位依次表示未检到的A和C)T_(左侧以下划线表示的一个数据位表示未检到的G)，如此，在不改变该读段对应的物理或化学含义或特征的情况下，使该读段的表现形式转为包含读取轮数据和轮空轮数据。

又例如，对于给定的一条碱基序列仍旧为AATGCGCCGT的读段，该读段通过双色荧光单分子测序获得(一种荧光信号对应一种碱基，A和C带有不同的荧光标记，T和G带有不同的荧光标记，一轮测序或者检出任意模板的一个位置的碱基类型包括模板先接触包含A和C的溶液进行一次碱基延伸、再接触包含T和G的溶液进行再次碱基延伸反应)，按照上面给定顺序排列的四种碱基ACTG(一个重复单元)，可以将该读段转化表示为A___(左侧以下划线表示的三个数据位依次表示未检到的C、T和G)A_(左侧以下划线表示的一个数据位依次表示未检到的C)T____(左侧以下划线表示的四个数据位依次表示未检到的G、A、C和T)G_(左侧以下划线表示的一个数据位表示未检到的A)C_(左侧以下划线表示的一个数据位表示未检到的T)G_(左侧以下划线表示的一个数据位表示未检到的A)C___(左侧以下划线表示的三个数据位依次表示未检到的T、G和A)C_(左侧以下划线表示的一个数据位表示未检到的T)G__(左侧以下划线表示的两个数据位依次表示未检到的A和C)T_(左侧以下划线表示的一个数据位表示未检到的G)，如此，在不改变该读段对应的物理或化学含义或特征的情况下，使该读段的表现形式转为包含读取轮数据和轮空轮数据。

一般情况下，呈现形式包含读取轮数据和轮空轮数据的读段的读长，一般仍是只根据其中的读取轮数据计算。另外，对于本文示例的序列，若无特别说明，以下划线“_”表示轮空轮。

在一些示例中，模板为单分子(非克隆体)，使模板在一次延伸反应中接触两种核苷酸，例如一个延伸反应溶液体系包含A和T，另一个延伸反应溶液体系包含G和C，并且采用两种荧光发光波段来区分开来自混合在一起的两种核苷酸的信号，由此，检测模板上的一个位置的碱基类型可看成包括两次延伸反应，按照上述示例的顺序ATGC为一个重复单元来表示两次延伸反应所测得的序列(读段)，可以理解地，读段将包括读取轮数据和轮空轮数据。

根据本申请的实施例，本申请的技术方案特别适合于“单分子测序”，在本文中所采用的术语“单分子测序”的含义是分辨率为单个或少数几个分子的大小，例如10个、8个、5个、4个或3个以下分子，例如在每次采集的图像上每个像素点(或者称为“亮斑”)的信号是来自于不超过10个分子拷贝的模板的延伸反应，例如来自于不超过8个、5个、4个或3个分子拷贝的模板的延伸反应。由于在单分子测序时，每个像素点或者亮斑的信号来源于较少的核酸模板或者靶核酸分子，不依赖样本的扩增，因而能够很好地避免由于扩增而导致的序列偏好性，进而提升下游生信分析的灵敏度。然而无扩增的测序体系意味着有效的测序信号均来源于单个DNA分子，从而单分子测序容易受到噪音及杂质信号的干扰，也容易出现信号丢失的情况。因此，在原始的测序数据中，会有部分错误率较高的序列，甚至是并非来源于样本的杂质序列，这些高错误率或者非样本来源的序列将占用大量测序仪比对计算的时间。并且即便有参考序列，很多情况下仍然很难区分错误造成的数据波动和生物样本特征造成的数据波动。

在本发明的一个可选实施例中，一轮测序是通过如下(1)-(4)任意一种方式实现：

(1)使所述模板先接触包含四种核苷酸中的两种的溶液进行单碱基延伸以及采集相应的图像、再接触包含四种核苷酸的另外两种的溶液进行单碱基延伸以及采集相应的图像，例如，一般所说的单分子双色测序(两种碱基两两带有相同的荧光标记)；

(2)使所述模板接触包含全部四种核苷酸的溶液进行单碱基延伸以及采集相应的图像，例如，四种碱基的延伸信号可区分的四色测序(每种碱基延伸对应唯一的一种荧光标记发光波段)或双色测序(两种碱基带有两种荧光标记、第三种碱基带有这两种荧光标记、第四种碱基不带荧光标记)或三色测序(三种碱基分别带有三种荧光标记、第四种碱基不带荧光标记)；

(3)使所述模板先接触包含四种核苷酸中的三种的溶液进行单碱基延伸以及采集相应的图像、再接触包含四种核苷酸中的剩下那一种的溶液进行单碱基延伸以及采集相应的图像，例如，四种碱基中的三种带有可区分的发光标记、剩下的那种碱基可带有与前面三种碱基中的任意一种碱基相同的发光标记；

(4)使所述模板依次接触包含四种核苷酸中的一种的溶液进行单碱基延伸以及采集相应的图像，同一般所称的单分子或高通量单色测序(四种碱基均带有相同的荧光标记)。

S200：基于所述图像集、所述读取轮数据和所述轮空轮数据中的至少一类信息，确定预定特征组合中每个特征的值，所述预定特征组合中的每个特征为与该读段的产生有关的、可量化的特征。

根据本申请的实施例，本申请的发明人提出了一种新的获得预定特征组合的方法，得到了可以用于利用预测模型进行测序读段质量判断的特征组合。

简言之，发明人首先在已知的读段所对应的图像集、读取轮数据和轮空轮数据的大量可检测特征选择了多个初步选定特征参数，其选择原则有如下几点：

1、测序软件能够较为容易地在每轮测序过程中获取这些参数的中间值，并在测序结束时快速地计算参数对应的最终数值；

2、各参数之间相对独立；

3、比对序列(本文有时也称比对读段)与非比对序列(本文有时也称非比对读段)在该参数下的分布有较明显的差异。

由此，发明人根据GenoCare测序技术的测序特点以及数据产生的原理，选择了一些潜在可用于区分比对序列和非比对(即无法比对上的)序列的参数作为潜在的特征。例如每条序列红色碱基总数、绿色碱基总数、红色碱基平均轮空轮数、绿色碱基平均轮空轮数、同一轮加碱基红绿碱基同时加上的次数、序列对应碱基的平均亮度、不同轮数亮度的方差、红绿碱基个数比例等等。然后发明人分析了这些参数两两之间的斯皮尔曼排序相关系数，对于相关性高的(斯皮尔曼相关系数高于0.5的)两个参数，仅保留其中的一个。例如发明人发现每条序列的红色碱基总数、绿色碱基总数这两个参数，均与最终序列读长高度正相关，因此最后保留了序列读长这一个参数。又例如发明人发现红色碱基平均轮空长度、绿色碱基平均轮空长度均与红绿整体的平均轮空长度呈较强的正相关关系，因此最终仅保留整体平均轮空长度这一参数。接下来发明人针对各个特征，分别制作了数据库中比对序列和非比对序列的参数值分布图，并判断比对与非比对序列的分布是否有差异。若几乎无差异，则弃用该参数。以图4为例，其显示了发明人使用的参数之一为例画出了分布图，该参数为每条序列读到的碱基在图像中对应模板位置处的平均荧光亮度。从图4可见，对于比对序列和非比对序列，其碱基平均荧光亮度分布有较为明显的差异，因此可以用于后续对模型的训练。

基于发明人进行的大量研究，根据本申请的一些实施例，预定特征组合包括选自下列的至少之一：

(a)读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征，所述读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据的至少一部分确定的；

(b)平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过至少一部分所述读取轮数据和所述轮空轮数据确定的；

(c)平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征，所述平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中的为指定类型碱基的那部分数据确定的，所述指定类型碱基选自A、T或U、C和G中的至少一种；

(d)平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置的信号强度确定的；

(e)重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征，所述重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中、在相同一轮测序读出的多种碱基的那部分数据以及读长确定的，所述读长的值是通过所述读取轮数据的至少一部分确定的；

(f)平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征，所述平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置及其周围背景的信号强度确定的；和

(g)模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征，所述模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置和/或所述轮空轮数据对应的图像的指定位置来确定的。

根据本申请的一些实施例，预定特征组合包括下列特征：选自(b)；(a)、(c)和(e)中的至少之一；以及选自(d)、(f)和(g)中的至少之一。根据本申请的一些实施例，预定特征组合包括(a)～(g)的所有特征。

根据本申请的一些实施例，与读长的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中指定一种碱基的总长度、所述读取轮数据中指定两种碱基的总长度和所述读取轮数据中指定三种碱基的总长度中的至少之一。根据本申请的一些实施例，与平均轮空轮长度的斯皮尔曼相关系数不低于0.5的特征包括：所述读段中指定一种碱基两个读出之间的平均轮空轮长度、所述读段中指定两种碱基两个读出之间的平均轮空轮长度和所述读段中指定三种碱基两个读出之间的平均轮空轮长度中的至少之一。根据本申请的一些实施例，与平均连续读取轮长度的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中指定一种碱基连续读取个数平均值、所述读取轮数据中指定两种碱基连续读取个数平均值和所述读取轮数据中指定三种碱基连续读取个数平均值中的至少之一。根据本申请的一些实施例，与所述平均亮度的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中指定一种碱基对应的图像位置的信号强度的平均值、所述读取轮数据中指定两种碱基对应的图像位置的信号强度的平均值、所述读取轮数据中指定三种碱基对应的图像位置的信号强度的平均值中的至少之一。根据本申请的一些实施例，与所述重复比例的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中在相同一轮测序中读出指定两种碱基所占的比例、所述读取轮数据中在相同一轮测序中读出指定三种碱基所占的比例中的至少之一。根据本申请的一些实施例，与所述平均信噪比的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中指定一种碱基对应的图像位置的信噪比平均值、所述读取轮数据中指定两种碱基对应的图像位置的信噪比平均值、所述读取轮数据中指定三种碱基对应的图像位置的信噪比平均值中的至少之一。根据本申请的一些实施例，与所述模板距离的斯皮尔曼相关系数不低于0.5的特征包括：所述读取轮数据中指定一种碱基对应的多个图像位置之间的间距、所述读取轮数据中指定两种碱基对应的多个图像位置之间的间距、所述读取轮数据中指定三种碱基对应的多个图像位置之间的间距中的至少之一。

根据本申请的实施例，斯皮尔曼相关系数(Spearman correlationcoefficient)，又称为等级相关系数，用于衡量两个变量之间的非线性相关性。斯皮尔曼相关系数的取值范围为-1到1之间，其中-1表示完全的负相关，1表示完全的正相关，0表示没有相关性。斯皮尔曼相关系数的计算基于两个变量的等级之间的差异，并通过对这些差异的平方和进行归一化来计算。简言之，计算斯皮尔曼相关系数的步骤可以包括如下：对两个变量的观测值分别进行排序，得到其等级。计算每对等级之间的差异，并计算差异的平方和，即斯皮尔曼相关系数的分子部分。根据数据的数量n，计算斯皮尔曼相关系数的分母部分，即分母等于n乘以n减去1，再除以2。将分子部分除以分母部分，得到最终的斯皮尔曼相关系数。斯皮尔曼相关系数通常用于非线性关系的变量之间的相关性分析，尤其适用于等级数据、有序分类数据和偏态数据。它比皮尔逊相关系数更鲁棒，能更好地捕捉到变量之间的整体趋势和关联性，而不受离群值和非线性关系的影响。

示例性地，“与读长的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与读长这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与读长的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中指定一种碱基(G、C、T、A中的一种，例如G)的总长度与读长之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中指定一种碱基(G、C、T、A中的一种，例如G)的总长度与读长可作为“与读长的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与平均轮空轮长度的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与平均轮空轮长度这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与平均轮空轮长度的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读段中指定一种碱基(G、C、T、A中的一种，例如G)两个读出之间的平均轮空轮长度与平均轮空轮长度之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读段中指定一种碱基(G、C、T、A中的一种，例如G)两个读出之间的平均轮空轮长度可作为“与平均轮空轮长度的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与平均连续读取轮长度的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与平均连续读取轮长度这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与平均连续读取轮长度的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中指定(G、C、T、A中的一种，例如G)连续读取个数的平均值与平均连续读取轮长度之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中指定(G、C、T、A中的一种，例如G)连续读取个数的平均值可作为“与平均连续读取轮长度的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与平均亮度的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与平均亮度这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与平均亮度的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中指定一种碱基(G、C、T、A中的一种，例如G)所对应的图像位置的信号强度的平均值与平均亮度之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中指定一种碱基(G、C、T、A中的一种，例如G)所对应的图像位置的信号强度的平均值可作为“与平均亮度的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与重复比例的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与重复比例这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与重复比例的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中在相同一轮测序中读出指定两种碱基(G、C、T、A中的两种，例如G和C)所占的比例与重复比例之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中在相同一轮测序中读出指定两种碱基(G、C、T、A中的两种，例如G和C)所占的比例可作为“与重复比例的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与平均信噪比的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与平均信噪比这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与平均信噪比的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中指定一种碱基(G、C、T、A中的两种，例如G)对应的图像位置的信噪比平均值与平均信噪比之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中指定一种碱基(G、C、T、A中的两种，例如G)对应的图像位置的信噪比平均值可作为“与平均信噪比的斯皮尔曼相关系数不低于0.5的特征”。

示例性地，“与模板距离的斯皮尔曼相关系数不低于0.5”是通过计算本申请中某一个特征与模板距离这两个变量之间的差异、并通过对这些差异的平方和进行归一化获得的。当得到的斯皮尔曼相关系数不低于0.5时，其所对应的特征为“与模板距离的斯皮尔曼相关系数不低于0.5的特征”。例如，当计算读取轮数据中指定一种碱基(G、C、T、A中的两种，例如G)对应的多个图像位置之间的间距与模板距离之间的差异、并通过对这些差异的平方和进行归一化处理后，获得的斯皮尔曼相关系数不低于0.5，该读取轮数据中指定一种碱基(G、C、T、A中的两种，例如G)对应的多个图像位置之间的间距可作为“与模板距离的斯皮尔曼相关系数不低于0.5的特征”。

如前所述，特征的选择对于能否成功地应用机器学习模型来预测测序读段质量是非常重要的，本发明的发明人也是通过大量研究提出了根据本申请实施例的预定特征的组合能够有效的应用于测序读段质量的预测，例如通过机器学习模型，诸如贝叶斯模型等。

S300：基于所述预定特征组合中的特征的值和预先训练的预测模型，预测所述读段被分类为指定类别读段的概率。

在获得相关预定特征组合的特征的值后，可以将所得到的特征输入经过预先训练的预测模型(也被称为机器学习模型)中，从而可以输出相关读段能够被归类为比对读段的概率，这个概率可以用来表征该读段的质量，因为这个概率表示该读段在后续分析中能够与参考序列匹配的几率。如果该几率比较低，则该读段的可靠性比较差，为了提高分析的效率和准确性，可以将该读段去除。

根据本申请的实施例，该预测模型选自逻辑回归、贝叶斯函数(或称贝叶斯模型)、决策树、随机森林、支持向量机、神经网络、K近邻算法，其中，优选贝叶斯函数，尤其是朴素贝叶斯函数。

为了方便理解，下面对贝叶斯模型进行简单介绍。

贝叶斯模型是一类基于贝叶斯定理的统计模型，它通过考虑观测数据和先验知识的关系来进行概率推断和预测。贝叶斯模型基于贝叶斯定理，利用已有的先验概率和观测数据的条件概率，更新为后验概率，从而得到对未观测数据的概率预测。贝叶斯定理是概率论中的一个重要定理，描述了在给定观测数据的条件下，更新先验概率为后验概率的过程。贝叶斯定理的数学表达形式如下：

P(A|B)＝(P(B|A)*P(A))/P(B)

其中，P(A|B)表示在给定观测数据B的条件下，事件A发生的概率；P(B|A)表示在事件A发生的条件下，观测数据B出现的条件概率；P(A)和P(B)分别表示事件A和观测数据B的先验概率。

贝叶斯模型的核心思想是将未知的参数或模型看作是随机变量，并使用贝叶斯定理来计算参数或模型的后验概率。贝叶斯模型的优点是能够根据已有的经验和知识对未知参数进行概率推断，并提供了一个灵活的框架来进行预测和决策。此外，贝叶斯模型还能够处理参数不确定性和模型不确定性，并能够灵活地进行模型更新和调整。

根据本申请的实施例，可以采用的贝叶斯模型包括朴素贝叶斯函数、贝叶斯线性回归模型、贝叶斯混合模型等。根据本申请的具体实施例，可以采用朴素贝叶斯函数来进行本申请的预测处理。朴素贝叶斯函数的一个优势是对于小规模的数据预测效果也良好，而且适合增量式训练。对于目前正在不断升级和持续积累数据的测序领域来说，使用朴素贝叶斯算法能够快速为不同版本的测序仪提供有效的质量分数。另一方面，朴素贝叶斯形成的模型在计算质量分数时可以直接使用查表法，运算速度比较快。

具体的，根据本申请的实施例，采用朴素贝叶斯函数，用于预测测序读段(也称为“测序序列”或者“读段”)的质量分数，即质量分数越高，代表对应序列能够被比对上(仅需比对上，无需唯一比对)的概率越高。因此标记数据集数据时，目标值为是否比对。

下面我们对朴素贝叶斯函数的基本公式进行解释：

在下述公式中，y

根据贝叶斯公式：

根据朴素贝叶斯中各参数相互独立的假设，可以认为

由此，可根据如下公式计算P(y

其中，先验概率P(y

需要说明的是，当特征数目较少时，例如仅采用一个特征进行预测时(即n＝1)，只需要对训练数据集的数据进行统计分析就可以得到先验概率P(y

根据本申请的实施例，预测模型是通过下列步骤进行预先训练的：获取训练读段集的比对结果和特征的值，比对结果是基于训练读段集中的读段与参考序列进行比对得到的，训练读段集的比对结果包括指定类别读段和非指定类别读段；将特征的取值和比对结果，输入至预测模型，以便对预测模型进行训练。

在本文中，“训练读段集”是指在对预测模型进行训练时，获得的测序数据集。该训练读段集分为训练集和测试集；其中训练集用于建立预测模型，测试集用于检验最终选择最优的模型的性能。

在本文中，“参考序列”是指理论上训练读段集中的读段需要能够比对上的参考序列，例如与训练读段来自相同物种的基因组序列，也可以已知与训练读段所对应核酸来源具有相同序列的人工序列。

在本发明的一个可选实施例中，该训练包括如下步骤：

1)获取训练读段集：首先，从已知是否属于指定类别的读段中获得一个训练读段集。每个读段都具有比对结果、以及一个或多个特征(具体参见上述预定特征组合中的特征)，并且这些特征是可量化的。

2)数据预处理：对数据进行预处理，包括处理缺失值、处理异常值、归一化或标准化特征等。

3)分割数据集：将数据集分为训练集和测试集。训练集用于构建模型，测试集用于评估模型性能。

4)计算类别先验概率：对于每个分类类别(目标值类别)，如对于指定类别读段和非指定类别读段，计算其在训练集中的先验概率，即该类别的样本占总样本数的比例。

5)计算特征的条件概率：对于每个特征，计算在每个类别下的条件概率。如对于每个特征的取值，计算在指定类别读段和非指定类别读段下的条件概率。朴素贝叶斯算法的"朴素"假设是特征之间相互独立，因此条件概率可以通过分别计算每个特征的概率来得到。

6)建立模型：使用计算得到的先验概率和条件概率建立朴素贝叶斯模型。

7)模型评估：使用测试集来评估模型的性能，若评估通过，则完成模型的训练。通常使用分类准确度、精确度、召回率等指标来评估模型在新数据上的表现。

更具体地，以下对本申请预测模型训练方法进行详细描述。

在本文中，对特征的方式不作限制，生成的特征可以是具有特定物理和/或化学含义的可量化的特征，也可以是不具有任何物理或化学意义的可量化的特征；例如，可以纯粹依据数据表现或呈现出的数字、符号、颜色和/或字符串等一种或多种形式的组合生成的可量化或者说可计算出取值的特征，换句话说，据数据所具有的数学特点和/或存在的规律或规则可以生成的任何可以量化的特征，都可以作为特征。例如，能够或不能够表征数据集中的全部或部分读段所具有的特性的可生成且可量化特征，可作为特征，能够或不能够表征该数据集产生过程中的全部或部分中间数据所具有的特性的可生成且可量化的特征，也可以作为特征。又例如，数据集和/或其产生过程的中间数据可生成的可量化的多个特征的集合、拟合、变形如归一化或标准化处理、任意数学运算等之后确定出的特征，也可以作为特征。再例如，读取轮数据的读长、轮空轮数据的平均轮空轮长度、读取轮数据给定碱基的平均连续读取轮长度、图像集确定给定位置的荧光信号平均亮度、读取轮的重复比例、图像集的给定位置的信号强度、图像集的模板距离等等。

在本发明的一个可选实施例中，通过下列步骤对所述预测模型进行训练：基于比对结果，将训练读段集分类为比对读段和非比对读段；将特征的取值和训练读段集的分类结果，输入至朴素贝叶斯函数中，来计算先验概率和条件概率；将所述先验概率和条件概率输入至朴素贝叶斯函数中，建立所述特征的取值和相应读段属于比对读段或属于非比对读段的概率的关联以便获得所述预先训练的预测模型。

在本申请的一些可选实施例中，将训练读段集的读段与参考序列进行比对，基于比对的结果将训练读段集中的所有读段分为比对读段和非比对读段；例如，当读段能够比对上或者相似性不低于60％(或者不低于70％、不低于80％、不低于90％、不低于95％、不低于99％)，即为比对读段，否则为非比对读段。将本发明上述预定特征组合中每个特征的取值、以及训练读段集的分类结果(即比对读段和非比对读段)作为特征输入至朴素贝叶斯函数，将分类结果(即比对读段和非比对读段)作为标记，对朴素贝叶斯函数的先验概率和条件概率进行训练，然后将得到的先验概率和条件概率输入至朴素贝叶斯函数中，建立特征的取值和相应读段属于比对读段或属于非比对读段的概率的关联，以便得到预先训练的预测模型。

在本发明的一个可选实施例中，先验概率和条件概率包括：训练读段集中比对读段的概率P(y

在本发明的一个可选实施例中，利用高斯函数gaussian(描述正态分布的数学表达式)和洛伦兹函数lorentzian(描述柯西分布的数学表达式，与正态分布相比，柯西分布没有有限的均值和方差，其尾部衰减非常慢。这意味着柯西分布的形状特征是长而尖的，没有明显的峰值)拟合所述训练读段集中特征的取值在比对读段中的频率分布，以计算该特征的取值在比对读段中的概率，以及利用高斯函数和洛伦兹函数拟合所述训练读段集中特征的取值在非比对读段中的频率分布，以计算该特征的取值在非比对读段中的概率。

在本发明的一个可选实施例中，该训练读段集包括多个读段，每个读段均包括多个特征，通过将训练读段集中读段的每个特征的取值输入至高斯函数和洛伦兹函数中，拟合得到该读段的每个特征的取值在比对读段/非比对读段中的频率分布，得到特征的取值在比对读段/非比对读段中的频率，即条件概率。然后将得到的条件概率作为优化后的参数输入至朴素贝叶斯函数中，用于训练预测模型，最终获得预先训练的预测模型。

在本发明的一个可选实施例中，利用该高斯函数和洛伦兹函数得到Voigt函数、Voigt2函数和Voigt2_Delta函数。采用上述Voigt函数、Voigt2函数和Voigt2_Delta函数的至少之一拟合得到特征的取值在比对读段/非比对读段中的频率分布，并计算该特征的取值在比对读段/非比对读段中的概率，作为条件概率。通过采用上述函数对预测模型中条件概率进行优化，可提高对训练模型的拟合能力。

根据本申请的实施例，本申请的发明人发现由于目前测序仪以及配套的试剂盒、测序方案均在不断更新，积累的数据量较少，通常测序数据的噪声较大。如果直接对序列在不同参数下的概率分布进行平滑处理，并用作参数对应的先验概率分布的话，分布中的离群点仍然会对一定范围内的分布图造成较明显的影响，导致预测结果不佳。因此发明人采用将多种函数进行叠加得到复合函数的方式，通过拟合的方式构建特征的取值-比对读段或非比对读段的概率分布图。

在本发明的一个可选实施例中，所述条件概率是通过所述训练读段集构建的特征的取值-比对读段或非比对读段的概率分布确定的，所述特征的取值-比对读段或非比对读段的概率分布是采用Voigt函数、Voigt2函数和Voigt2_Delta函数的至少之一进行拟合的，其中，Vogit函数是包含中心偏移项的单个Vogit函数；Voigt2函数是两个包含中心偏移项的Vogit函数的叠加结果；和Voigt2_Delta函数是Voigt2函数与狄拉克δ函数的叠加结果。其中相关函数的数学表达式如下：

Voigt函数：

其中，Г

Voigt2函数：

voigt2＝b

其中，Г

Voigt2_Delta函数：

voigt_Delta2＝b

+b

其中，Г

在Voigt2_Delta函数中，x表示各特征取值作为自变量，比对读段或非比对读段的概率作为因变量。

图4显示了采用Voigt2_Delta函数拟合得到的特征的取值-比对读段或非比对读段的概率分布图示例(采用平均亮度)，显示拟合函数

对原始数据分布曲线的拟合效果很好。

如前所述，通过分析各个特征各自的特征的取值-比对读段或非比对读段的概率分布图，可以得到先验概率P(y

其中，

P(y

P(y

P(x

P(x

P(y

S400：基于所述概率，确定所述读段的测序质量，所述预先训练的预测模型为能够基于读段的预定特征组合中的特征的值计算该读段属于指定类别读段的概率的分类器，所述读段的质量得分与该读段被分类为指定类别读段的概率正相关。

根据本申请的实施例，在前述步骤中确定读段被归类为比对读段的概率后，可以进一步将相关概率数值转换为通用的测序质量评分。根据本申请的实施例，可以通过P(y

由此，根据本发明的实施例，通过采用本发明实施例的技术方案，本发明的发明人发现能够有效地估算每一条序列(读段)比对成功(例如与参考序列成功匹配或者反应真实的核苷酸序列)的概率及评估序列的错误率，从而在原始测序数据下机后能够根据序列的质量分数过滤掉一部分高错误率序列及非样本来源的杂质序列，从而提升下游生信分析的速度、降低高错误序列及杂质序列对后续生信分析的影响。

另外，根据本申请的实施例，可在进行比对之前筛除掉很大比例的无法比对上的序列以及高错误率序列，大幅减少比对分析和后续生信分析的耗时。根据本申请的实施例，采用本发明的方法可降低高错误率和非比对读段(非比对读段)的输出比例，帮助减少错误序列对测序数据分析结果的影响。例如当测序样本有多种物种来源时(例如以检测病原微生物感染为目的而采集的人类体液样本等)，如果存在大量高错误率、实际应该是非比对读段的测序序列，将很难用比对结果将其区分开。这是由于样本中存在多物种时，非比对读段随机比对上不同序列的概率会增加。因此在比对之前根据质量分数进行一步过滤能够去除掉很大一部分这样的高错误率非比对读段，从而减少这些序列在后续生信分析过程中对分析结果造成的背景噪音影响。

另外，需要说明的是，本领域技术人员能够理解的是，在完成机器学习的训练后，输出了经过迭代优化的先验概率P(y

具体的，参考图2，在本申请的第二方面，本发明提出了一种确定读段质量的方法，根据本申请的实施例，所述读段来自基于芯片成像检测、对芯片上的模板进行边合成边测序而确定的序列，所述读段由读取轮数据和轮空轮数据构成，每条所述读段的读取轮数据反映对应于图像集的指定位置的模板上存在碱基延伸信号、所述轮空轮数据反映对应于所述图像集的所述指定位置的模板上不存在碱基延伸信号，所述图像集包括测序运行期间的一轮或多轮测序产生的、包含所述指定位置及其周围背景信号的多个图像，所述方法包括：

S100A：针对所述读段，基于所述图像、所述读取轮数据和所述轮空轮数据中的至少一类信息获取预定特征组合中的每个特征的取值。

根据本申请的实施例，所述预定特征的取值是基于图像集、读取轮数据和轮空轮数据中的至少一种确定的。

根据本申请的实施例，本申请的发明人获得预定特征组合的方法可参见第一方面所述的方法。

根据本申请的一些可选实施例中，所述预定特征包括如下中的至少之一：

(a)读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征，所述读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据的至少一部分确定的；

(b)平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过至少一部分所述读取轮数据和所述轮空轮数据确定的；

(c)平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征，所述平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中的为指定类型碱基的那部分数据确定的，所述指定类型碱基选自A、T或U、C和G中的至少一种；

(d)平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置的信号强度确定的；

(e)重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征，所述重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中、在相同一轮测序读出的多种碱基的那部分数据以及读长确定的，所述读长的值是通过所述读取轮数据的至少一部分确定的；

(f)平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征，所述平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置及其周围背景的信号强度确定的；和

(g)模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征，所述模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置和/或所述轮空轮数据对应的图像的指定位置来确定的。

根据本申请的一些实施例，预定特征组合包括下列特征：选自(b)；(a)、(c)和(e)中的至少之一；以及选自(d)、(f)和(g)中的至少之一。根据本申请的一些实施例，预定特征组合包括(a)～(g)的所有特征。

由此，根据本发明的实施例，预定特征组合包括所述测序运行期间产生的、与所述读段有关的多个特征。

S200A：针对所述读段，基于所述特征的取值，确定与所述取值对应的比对系数(例如P1m)和非比对系数。

根据本申请的实施例，在该方法中，与取值对应的比对系数和非比对系数是采用训练测序读段集预先确定的，训练读段集具有已知比对结果以及特征的取值。

根据本申请的实施例，在获得特征的取值后，可以将所得到的特征输入经过预先训练的预测模型(例如第一方面所述的方法中的预测模型)中，从而可以输出与取值对应的比对系数(即为P1m)和非比对系数(即为P2m)。

根据本申请的实施例，在该方法中，与取值对应的比对系数和非比对系数是采用朴素贝叶斯函数预先确定的。

根据本申请的一个可选实施例中，比对系数和非比对系数是基于朴素贝叶斯函数的条件概率确定的。根据本申请的实施例，朴素贝叶斯函数的条件概率的获得可参见第一方面所述的方法。

S300A：基于所述比对系数、非比对系数以及预定的比对常数和预定的非比对常数，确定所述读段的质量得分。

根据本申请的实施例，在该方法中，比对常数、非比对常数是采用训练测序读段集预先确定的，训练读段集具有已知比对结果以及特征的取值。

根据本申请的实施例，在获得特征的取值后，可以将所得到的特征输入经过预先训练的预测模型(例如第一方面所述的方法中的预测模型)中，从而可以输出与取值对应的比对常数(先验概率P(y

根据本申请的实施例，在该方法中，比对常数、非比对常数是采用朴素贝叶斯函数预先确定的。

可选的，比对常数和非比对常数是基于朴素贝叶斯函数的先验概率确定的，比对系数和非比对系数是基于朴素贝叶斯函数的条件概率确定的。

根据本申请的实施例，在该方法中，质量是基于下列公式确定的：

其中，

Q表示读段的测序质量；

C

C

P1m表示第m个特征的取值所对应的比对系数；和

P2m表示第m个特征的取值所对应的非比对系数。

[]表示取整函数。

根据本申请的实施例，在该方法中，测序质量可以是基于下列公式确定的：

其中，

Q表示读段的测序质量；

C

C

P1m表示第m个特征的取值所对应的比对系数；和

P2m表示第m个特征的取值所对应的非比对系数；

表示下取整函数。

根据本申请的实施例，在该方法中，进一步包括采用训练读段集确定所述比对常数、所述非比对常数、与所述特征的特定取值对应的所述比对系数和非比对系数。

根据本申请的实施例，在该方法中，比对常数、非比对常数、与取值对应的比对系数和非比对系数是采用训练测序读段集预先确定的，训练读段集具有已知比对结果以及特征的取值。

根据本申请的实施例，在该方法中，比对常数、非比对常数、与取值对应的比对系数和非比对系数是采用朴素贝叶斯函数预先确定的。

可选的，比对常数和非比对常数是基于朴素贝叶斯函数的先验概率确定的，比对系数和非比对系数是基于朴素贝叶斯函数的条件概率确定的。

参考图3，在本申请的第三方面，本申请提出了一种确定读段测序质量的方法，根据本申请的实施例，所述读段来自基于芯片成像检测、对芯片上的模板进行边合成边测序而确定的序列，所述读段由读取轮数据和轮空轮数据构成，每条所述读段的读取轮数据反映对应于图像集的指定位置的模板上存在碱基延伸信号、所述轮空轮数据反映对应于所述图像集的所述指定位置的模板上不存在碱基延伸信号，所述图像集包括测序运行期间的一轮或多轮测序产生的、包含所述指定位置及其周围背景信号的多个图像，所述方法包括：

S100B：针对所述读段，基于所述图像、所述读取轮数据和所述轮空轮数据中的一类信息获取预定特征的取值。

根据本申请的实施例，所述预定特征的取值是基于图像集、读取轮数据和轮空轮数据中的至少一种确定的。

根据本申请的实施例，本申请的发明人获得预定特征组合的方法可参见第一方面所述的方法。

根据本申请的一些可选实施例中，所述预定特征包括如下中的至少之一：

(a)读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征，所述读长或者与所述读长的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据的至少一部分确定的；

(b)平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均轮空长度或者与所述平均轮空长度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过至少一部分所述读取轮数据和所述轮空轮数据确定的；

(c)平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征，所述平均连续读取轮数或者与所述平均连续读取轮数的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中的为指定类型碱基的那部分数据确定的，所述指定类型碱基选自A、T或U、C和G中的至少一种；

(d)平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征，所述平均亮度或者与所述平均亮度的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置的信号强度确定的；

(e)重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征，所述重复比例或者与所述重复比例的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据中、在相同一轮测序读出的多种碱基的那部分数据以及读长确定的，所述读长的值是通过所述读取轮数据的至少一部分确定的；

(f)平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征，所述平均信噪比或者与所述平均信噪比的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置及其周围背景的信号强度确定的；和

(g)模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征，所述模板距离或者与所述模板距离的斯皮尔曼相关系数不低于0.5的任意特征的值是通过所述读取轮数据对应的图像的指定位置和/或所述轮空轮数据对应的图像的指定位置来确定的。

根据本申请的一些实施例，预定特征组合包括下列特征：选自(b)；(a)、(c)和(e)中的至少之一；以及选自(d)、(f)和(g)中的至少之一。根据本申请的一些实施例，预定特征组合包括(a)～(g)的所有特征。

由此，根据本发明的实施例，预定特征组合包括所述测序运行期间产生的、与所述读段有关的多个特征。

S200B：基于所述预定特征的取值以及预先构建的特征-质量关系，确定与所述读段的质量得分。

根据本申请的实施例，针对上述获取的每个特征的取值，预先构建了特征-质量关系(或称特征值-质量关系图)。

根据本申请的一个可选实施例中，特征-质量关系是采用Voigt函数、Voigt2函数和Voigt2_Delta函数的至少之一进行拟合的，其中，Vogit函数、Voigt2函数和Voigt2_Delta函数具体参见第一方面所述的方法。

在本发明的一个可选实施例中，所述预先构建的特征-质量关系是基于训练读段集确定的，所述训练读段集包含指定类别读段和非指定类别读段以及相应类别中的读段的所述特征的取值。

本领域技术人员能够理解，针对第一方面描述的特征和优点，同样适用于第二方面和第三方面，在此不再赘述。需要说明的是，为了提升测序过程中软件对序列质量分数的计算速度，软件在套用朴素贝叶斯模型进行预测时，还可以不直接用上述参数拟合得到的先验分布函数计算对应的概率数值，而是将这些先验分布函数分箱并制表，在测序结束后直接查表得到对应的概率数值。加入质量分数计算需求后，对于标准的16通道自动测序加分析的流程，整体数据分析的总延时可以大大缩短，例如10分钟左右。

在本申请的第四方面，本申请提出了一种测序方法。根据本申请的实施例，该方法包括：采用基于芯片成像检测的边合成边测序技术，对芯片上的模板进行测序，获取图像，并且分析图像以确定测序结果，所述测序结果包括多个读段；和第一方面所述的方法，确定各读段的质量得分。

另外，根据本申请的实施例，在测序后，基于测序质量，将低于预定阈值的测序读段进行过滤。根据本申请的实施例，这里所使用的阈值，是可以根据实际需要，采用对照试验或者平行试验确定的，在此不再赘述。

由此，根据本发明的实施例，通过采用本发明实施例的技术方案，本发明的发明人发现能够有效地估算每一条序列(读段)比对成功(例如与参考序列成功匹配或者反应真实的核苷酸序列)的概率及评估序列的错误率，从而在原始测序数据下机后能够根据序列的质量分数过滤掉一部分高错误率序列及非样本来源的杂质序列，从而提升下游生信分析的速度、降低高错误序列及杂质序列对后续生信分析的影响。

另外，根据本申请的实施例，可在进行比对之前筛除掉很大比例的无法比对上的序列以及高错误率序列，大幅减少比对分析和后续生信分析的耗时。根据本申请的实施例，采用本发明的方法可降低高错误率和非比对读段(非比对读段)的输出比例，帮助减少错误序列对测序数据分析结果的影响。例如当测序样本有多种物种来源时(例如以检测病原微生物感染为目的而采集的人类体液样本等)，如果存在大量高错误率、实际应该是非比对读段的测序序列，将很难用比对结果将其区分开。这是由于样本中存在多物种时，非比对读段随机比对上不同序列的概率会增加。因此在比对之前根据质量分数进行一步过滤能够去除掉很大一部分这样的高错误率非比对读段，从而减少这些序列在后续生信分析过程中对分析结果造成的背景噪音影响。

本发明的发明人在对测序技术进行研究的过程中，本发明的方法尤其适用于单分子双色测序技术。根据本发明的实施例，单分子双色测序技术在测序原理及数据处理方法方面与二代测序技术有很大的差异。简言之，单分子测序系统不依赖样本的扩增，因而能够很好地避免由于扩增而导致的序列偏好性，进而提升下游生信分析的灵敏度。然而无扩增的测序体系同时也意味着有效的测序信号均来源于单个核酸分子或者少数核酸分子，因此单分子测序会更容易受到噪音及杂质信号的干扰，也容易出现信号丢失的情况。因此在原始的测序数据中，会有部分错误率较高的序列甚至是并非来源于样本的杂质序列。这些高错误率或者非样本来源的杂质序列将占用大量测序仪比对计算的时间，即为了判断测序仪所输出读段的质量，需要将读段与参考序列进行比对，这会耗费大量的计算时间，降低测序仪的工作效率。另外，不可避免的是，在很多情况下，是不存在参考序列的。即使存在参考序列，很多情况下仍然很难区分低质量的读段是由于测序过程的错误造成的数据波动还是生物样本特征造成的数据波动。例如某些读段中的一部分序列有可能随机比对到参考基因组上，在样本情况复杂、涉及不同物种来源的测序中，这些随机比对上的错误序列可能对不同物种测序序列的拆分以及后续的生信分析造成影响。由此，针对三代测序技术例如Genocare测序平台，急需新的测序读段质量评价方法，尤其是不需要通过比对处理就能够对测序读段质量进行评价的方法。发明人发现通过采用上述确定读段质量的方法，能够有效地对单分子双色测序的读段进行质量评价。

在本发明的第五方面，本发明提出了一种计算设备，根据本发明的实施例，该计算设备包括：处理器和存储器；所述存储器，用于存储计算机程序；所述处理器，用于执行所述计算机程序以实现前面第一方面至第三方面所述的方法。

在本发明的第六方面，本发明提出了一种计算机可读存储介质，所述存储介质包括计算机指令，当所述指令被计算机执行时，使得所述计算机实现如前面第一方面至第三方面所述的方法。

需要说明的是，在本发明中，在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，“计算机可读介质”可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(RAM)，只读存储器(ROM)，可擦除可编辑只读存储器(EPROM或闪速存储器)，光纤装置，以及便携式光盘只读存储器(CDROM)。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。本发明描述的各种计算机可读存储介质可代表用于存储信息的一个或多个设备和/或其它机器可读存储介质。术语“机器可读存储介质”可包括但不限于无线信道和能够存储、包含和/或承载指令和/或数据的各种其它介质。

应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

在本发明的第七方面，本发明提出了一种测序装置，其包括：前面第五方面所述的计算设备；或者前面第六方面所述的计算机可读存储介质。

实施例

在本实施例中，选择读长(rl)、平均轮空长度(emp_cyc)、平均连续读取轮数(rep_RG)、平均亮度(int_mean)、重复比例(rep_perc)、平均信噪比(snr_mean)、模板距离(temp_dist)作为特征(或者特征参数、特征或参数)，采用人类gDNA测序数据作为训练集和测试集。这些数据完成了与人类基因组hg19的比对，发明人从比对得到的sam文件中获取每一条序列是否成功比对上的信息。这些数据取自3台不同的测序仪的7次不同的上机实验，数据来自11条不同的测序通道，共取了其中63万条序列作为数据集。其中随机拆分70％作为训练集，30％作为测试集。所有训练和测试数据均来源于16样本测序模式下获得的测序数据(即每条通道一种样本，所有样本不带标签index，测序过程不测标签)。

在获得训练集和测试集后，分别按照本文中所描述的方法确定相关特征的取值，并采用Voigt2_Delta函数拟合得到特征的取值-比对读段或非比对读段的概率分布图。

其中，Voigt2_Delta函数的数学表达式为：

voigt_Delta2＝b

+b

在Voigt2_Delta函数中，x表示各特征取值作为自变量，比对读段或非比对读段的概率作为因变量。

图4显示了采用Voigt2_Delta函数拟合得到的特征的取值-比对读段或非比对读段的概率分布图示例(采用平均亮度)。

通过一系列特征的取值-比对读段或非比对读段的概率分布图，可以得到先验概率、条件概率的初始值，接下来采用朴素贝叶斯模型对这些初始值进行监督式训练，从而得到了经过优化的模型，该模型可以用于预测测序读段被分类为比对序列的概率。接下来，可以通过P(y

发明人在实时测序和分析软件TGGS上实现了质量分数值输出的功能。在原始测序结果fastq文件中以及比对结果sam文件中，每条序列均会标注出按照模型计算出的质量分数。

为了证明本发明的读段质量判断方法在区分比对读段与非比对读段时的效果。在本实施例中，所有展示数据均来源于16样本测序模式下获得的测序数据(即每条通道一种样本，所有样本不带标签，测序过程不测标签)。分别展示了人类gDNA混样测试集数据的结果、大肠杆菌样本的结果、不同来源的人类基因组样本的结果。其中，图5显示了人类gDNA混样测试集数据中实际比对读段与实际非比对读段的质量分数值分布，图6显示了人类gDNA混样测试集数据在不同的质量分数过滤阈值之下，对应的比对率及序列过滤和损失比例，图7显示了人类gDNA混样测试集数据中质量分数与对应的唯一比对读段的平均错误率的关联，图8显示了肠杆菌测序数据中实际比对读段与实际非比对读段的质量分数值分布；图9显示了大肠杆菌测序数据在不同的质量分数过滤阈值之下，对应的比对率及序列过滤和损失比例；图10显示了大肠杆菌测序数据中质量分数与对应的唯一比对读段的平均错误率的关联；图11显示了不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后非比对读段过滤比例的统计结果；图12不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后比对读段损失比例的统计结果；图13不同来源的人基因组样本测序数据中，使用不同质量分数过滤阈值进行过滤后比对率的统计结果。通过这些数据图可以看出，测序读段的质量判断方法所得到的的质量分数能够较好地区分开比对和非比对读段，并可以在控制比对读段损失的情况下有效过滤非比对读段，并提升比对率。另外，根据图7和图10可以看出本发明的判断方法与对应序列的错误率呈负相关关联。

由此，采用本发明质量分数评价体系对于测序数据的分析的优势至少在于：

a)能够在进行比对之前筛除掉很大比例的无法比对上的序列以及高错误率序列，大幅减少比对分析和后续生信分析的耗时。

b)降低高错误率和非比对读段的输出比例，帮助减少错误序列对测序数据分析结果的影响。例如当测序样本有多种物种来源时(例如以检测病原微生物感染为目的而采集的人类体液样本等)，如果存在大量高错误率、实际应该是非比对读段的测序序列，将很难用比对结果将其区分开。这是由于样本中存在多物种时，非比对读段随机比对上不同序列的概率会增加。因此在比对之前根据质量分数进行一步过滤能够去除掉很大一部分这样的高错误率非比对读段，从而减少这些序列在后续生信分析过程中对分析结果造成的背景噪音影响。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。