CN111462818B - 测序产量预测方法和建立测序产量预测模型的方法及装置 - Google Patents

Abstract

一种测序产量预测方法和建立测序产量预测模型的方法及装置，测序产量预测方法包括：获取测序过程中的测序数据，包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；根据测序数据和测序产量预测模型公式，预测测序完成时的测序产量；和输出预测出的测序完成时的测序产量。本发明能够有效预测测序平台的测序产量，通过实时监控测序过程中的测序仪参数变化，实现测序产量早预测、测序异常早发现，从而缩短生产交付周期，并及时向前端反馈测序进展。

Description

测序产量预测方法和建立测序产量预测模型的方法及装置

技术领域

本发明涉及测序技术领域，具体涉及一种测序产量预测方法和建立测序产量预测模型的方法及装置。

背景技术

联合探针锚定聚合技术（cPAS，Combinatorial Probe-Anchor Synthesis）是一些测序仪使用的测序技术。DNA文库环化后形成单链环状DNA分子，经过滚环复制形成DNA纳米球（DNA nano ball，DNB）。DNA纳米球通过DNA装载技术，固定在规则阵列（Pattern Array）的芯片上。经过先进的光刻和干法刻蚀技术，形成可以装载DNA纳米球的芯片。测序时，首先DNA分子锚定荧光探针在DNB上进行聚合酶促反应，每扩增一个碱基，高分辨率成像系统将采集光信号，然后经过数字化处理，将光信号转化为DNA碱基信息。这种测序仪具有更低的试剂消耗、更低的重复率、更高的数据产出和更高的测序准确度等优势。这种测序仪的测序产量主要通过下机后统计单个通道（Lane）的读长（Reads(M)）和碱基（Base(G)）来反映。

目前测序过程中，通常是等待下机后统计单个通道（Lane）的测序产量，并针对测序产量不足的情况，进行新一轮的加测。但是在测序过程中，针对最终测序产量的监控和预测方法少之又少。因此，开发一种行之有效的方法，在测序未完成时，预测测序仪单个通道（Lane）的最终测序产量，能够缩短生产交付周期，提高交付服务质量。

发明内容

本发明提供一种测序产量预测方法和建立测序产量预测模型的方法及装置，能够有效预测测序平台的测序产量，通过实时监控测序过程中的测序仪参数变化，实现测序产量早预测、测序异常早发现，从而缩短生产交付周期，并及时向前端反馈测序进展。

根据第一方面，一种实施例中提供一种测序产量预测方法，包括：

获取测序过程中的测序数据，该测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；

根据上述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：

（a）测序产量= a+ b·BIC+ c·ESR+d·ACCGRR± 置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数；

（b）测序产量= m+ n·BIC+ p·FIT+q·ESR+ r·光强± 置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数；和

输出预测出的测序完成时的测序产量。

在优选实施例中，上述a的取值为-510.983及其±50%范围内的数值；上述b的取值为-5.581及其±50%范围内的数值；上述c的取值为1474.396及其±50%范围内的数值；上述d的取值为498.429及其±50%范围内的数值。

在优选实施例中，上述m的取值为-972.058及其±50%范围内的数值；上述n的取值为7.471及其±50%范围内的数值；上述p的取值为142.704及其±50%范围内的数值；上述q的取值为982.381及其±50%范围内的数值；上述r的取值为-0.003及其±50%范围内的数值。

在优选实施例中，上述置信区间为标准差或其大于1的倍数。

在优选实施例中，上述测序过程中的测序数据是测序进行至30-120个循环时的测序数据，更优选测序进行至60个循环时的测序数据。

根据第二方面，一种实施例中提供一种测序产量预测装置，包括：

测序数据获取单元，用于获取测序过程中的测序数据，上述测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；

测序产量预测单元，用于根据上述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：

（a）测序产量= a+ b·BIC+ c·ESR+d·ACCGRR± 置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数；

（b）测序产量= m+ n·BIC+ p·FIT+q·ESR+ r·光强± 置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数；和

测序产量输出单元，用于输出预测出的测序完成时的测序产量。

根据第三方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的测序产量预测方法。

根据第四方面，一种实施例中提供一种建立测序产量预测模型的方法，包括：

获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，该测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；

对上述测序产量数据进行预处理以排除异常值；

根据上述测序产量数据和上述测序数据，基于线性回归模型建立测序产量预测模型；和

利用上述测序产量预测模型和上述测序产量数据，模拟计算测序产量预测值并统计残差。

在优选实施例中，上述方法还包括：

根据上述测序产量预测值的大小对上述测序产量预测值分类，并按类统计上述残差的标准差。

在优选实施例中，上述方法还包括：

根据上述测序产量预测模型和上述残差及输入的测序过程中的测序数据，预测测序完成时的测序产量；和

按时间周期统计新产生的测序产量数据的残差，并在监测到上述残差的标准差偏离原残差的标准差时，重新进行上述测序产量预测模型的建立。

根据第五方面，一种实施例中提供一种建立测序产量预测模型的装置，包括：

数据获取单元，用于获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，上述测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；

异常值排除单元，用于对上述测序产量数据进行预处理以排除异常值；

模型建立单元，用于根据上述测序产量数据和上述测序数据，基于线性回归模型建立测序产量预测模型；

预测值计算单元，用于利用上述测序产量预测模型和上述测序产量数据，模拟计算测序产量预测值并统计残差。

根据第六方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第四方面的建立测序产量预测模型的方法。

本发明的测序产量预测方法和建立测序产量预测模型的方法，能够有效预测测序平台的测序产量，通过实时监控测序过程中的测序仪参数变化，实现测序产量早预测、测序异常早发现，从而缩短生产交付周期，并及时向前端反馈测序进展。

附图说明

图1为本发明实施例中测序产量预测方法流程图；

图2为本发明实施例中测序产量预测装置结构框图；

图3为本发明实施例中一种建立测序产量预测模型的方法流程图；

图4为本发明实施例中另一种建立测序产量预测模型的方法流程图；

图5为本发明实施例中建立测序产量预测模型的装置结构框图；

图6为本发明实施例中测序产量与拆分率和ESR的相关性关系图；

图7为本发明实施例中测序产量与不同阶段参数值的相关性关系图；

图8至图11为本发明实施例中不同的数值参数之间的相关性关系图；

图12为本发明实施例中聚类的簇数量k（Number of cluster k）与模型解释力（Gap statistic）之间的关系曲线；

图13为本发明实施例中采用逐步线性回归建立模型的预测变量与模型性能的关系图；

图14为本发明实施例中采用逐步线性回归建立模型的模型显著性结果图；

图15为本发明实施例中得到的模型结论参数结果图；

图16为本发明实施例中测序产量预测模型的置信区间与预测准确率关系图；

图17至图19为本发明实施例中按照测序产量预测模型公式进行测序产量预测结果图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本发明能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。

另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。

如图1所示，本发明一种实施例中提供一种测序产量预测方法，包括如下步骤：

S101：获取测序过程中的测序数据，该测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强。

本发明实施例中，所获取的测序过程中的测序数据是指测序过程中的各种参数，并非测序读长（Reads）本身。这样的参数包括但不限于BIC、ESR和ACCGRR。其中BIC（basecall information content）是指可用于作碱基调取（basecalling）的DNA纳米球DNB比率；ESR（Effective Spot Rate）是指经过过滤最终测序所得的读长数（Total Reads），占碱基调取（Basecall）能识别出来的DNB数量（DNB Number）的比率，即ESR=TotalReads/DNBNumber×100%；ACCGRR（Average accumulative continued Good Reads Rate）是指连续平均累计好读长（Reads）占比，反应测序的连续测序质量。其他可能的参数还包括dnbnum、总ESR（totalESR）、FIT、光强（Intensity）、信噪比（SNR）和拆分率（Split Rate）等。其中，Dnbnum，或Dnbnum DNB数量是指以测序通道（Lane）为单位，可以正常产出数据的FOV数。totalESR，是指总ESR（total EffectiveSpot Rate，totalESR）。FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强（Intensity）分布越集中（一种碱基光强很高，其他三种碱基光强很低），其FIT值越高；光强（Intensity）分布越离散，FIT值越低。光强（Intensity）：数据来源rawInt，为经过图像处理后提取的光强值，对rawInt的每一列进行排序（descend），取最高1%~15%（百分位点）之间的光强求平均值，为原始光强估计值。信噪比（SNR）：以单个DNB SNR的计算为例，碱基A最大光强（max intensity）作为信号，碱基CGT则为背景，CGT光强的方差则为噪声，那么碱基A第一循环（cycle1）的SNR（CGT_dev）估算方法为：碱基A的信号取所有碱基A的平均值（A_mean），碱基CGT矩阵的方差（CGT_dev）则为噪声，那么A_SNR=A_mean/CGT_dev。拆分率（Split Rate）：拆分出来的标签序列（Barcode）能对应到所有（例如1-96号）标签序列的比率。标签拆分率为标签序列拆分率柱状图，纵坐标为拆分率，横坐标为标签序列编号，只显示拆分率大于1%的标签序列情况。这些参数与测序产量有一定相关性，并且不同测序产量与测序产量的相关性不同，即权重不同，因此，根据这些参数并依据一定的模型可以预测测序产量。

本发明的测序产量预测方法能够在测序全过程中任意时刻预测测序完成时的测序产量，例如对于测序循环（Cycle）总数是210个循环的测序，本发明的测序产量预测方法能够从第1至210个循环中任意一个循环预测测序完成时的测序产量。然而，在一个优选实施例中，在测序进行至30-120个循环时，获取测序数据进行测序产量预测；更优选，在测序进行至60个循环时，获取测序数据进行测序产量预测。一方面保证及时获知测序完成时的测序产量预测结果，特别是结果少于预期时可以及时加测数据；另一方面保证测序产量预测的准确性。过早进行测序产量预测虽然能够及时获知测序完成时的测序产量预测结果，但是准确性会不足；过晚进行测序产量预测虽然能够保证测序产量预测的准确性，但是不能及时加测数据失去了预测的意义。

S102：根据上述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：

（a）测序产量= a+ b·BIC+ c·ESR+d·ACCGRR± 置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数。

（b）测序产量= m+ n·BIC+ p·FIT+ q·ESR+ r·光强± 置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数。

上述测序产量是根据测序数据和测序产量通过模拟的形式建立的经验性公式，根据上述常数和/或系数不同，上述测序产量预测模型公式可以广泛适用于具有各种不同测序因素的测序过程，例如不同的测序平台、不同的测序仪器、不同的测序原理、不同的测序策略、不同测序地点和测序时间段等的测序过程中测序产量的预测。随着测序因素的变化，上述常数和/或系数会发生变化。

在本发明的一些实施例中，a的取值为-510.983及其±50%范围内的数值；b的取值为-5.581及其±50%范围内的数值；c的取值为1474.396及其±50%范围内的数值；d的取值为498.429及其±50%范围内的数值；置信区间为标准差或其大于1的倍数。其中，±50%表示比本数（例如-510.983）小或大50%的数值；置信区间为标准差或其大于1的倍，例如，置信区间可以是1倍的标准差，2倍的标准差，3倍的标准差。置信区间所取数值越大，预测准确性或可靠性越大。

在本发明的一些实施例中，m的取值为-972.058及其±50%范围内的数值；n的取值为7.471及其±50%范围内的数值；p的取值为142.704及其±50%范围内的数值；q的取值为982.381及其±50%范围内的数值；r的取值为-0.003及其±50%范围内的数值。

S103：输出预测出的测序完成时的测序产量。

通过步骤S102获得的测序完成时的测序产量被输出，能够有效预测测序平台的测序产量，通过实时监控测序过程中的测序仪参数变化，实现测序产量早预测、测序异常早发现，从而缩短生产交付周期，并及时向前端反馈测序进展。

对应于上述测序产量预测方法，本发明实施例还提供一种测序产量预测装置，如图2所示，该测序产量预测装置包括：测序数据获取单元201，用于获取测序过程中的测序数据，测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；测序产量预测单元202，用于根据上述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：（a）测序产量= a+ b·BIC+ c·ESR+d·ACCGRR± 置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数；（b）测序产量= m+ n·BIC+ p·FIT+ q·ESR+ r·光强± 置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数；和测序产量输出单元203，用于输出预测出的测序完成时的测序产量。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

因此，本发明实施例还提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如本发明的测序产量预测方法。

如图3所示，本发明一种实施例中提供一种建立测序产量预测模型的方法，包括：

S301：获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，该测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强。

本发明实施例中，测序产量数据来源于已完成测序，可以是任何测序平台和测序仪器在任何时间和地点产生的测序产量数据，这些测序产量数据中包含每一个测序（每个RUN）的测序产量信息，例如若干M的测序读长（Reads）和/或若干G的碱基（Base）数量。所采用的测序产量数据越多越好，可以多至成千上万个测序（即RUN）的测序产量数据。测序产量数据越多，越有利于得到精确的测序产量预测模型。在本发明一个实施例中，采用来自华大基因公司武汉信息化平台（http://172.28.25.87:8085/）的3000多个测序（即RUN）的测序产量数据。

这些测序产量数据都有对应的测序数据，这些测序数据是指测序过程中的各种参数，并非测序读长（Reads）本身。这样的参数包括但不限于BIC、ESR和ACCGRR。其中BIC（basecall information content）是指可用于作碱基调取（basecalling）的DNB比率；ESR（Effective Spot Rate）是指经过过滤最终测序所得的读长数（Total Reads），占碱基调取（Basecall）能识别出来的DNB数量（DNB Number）的比率，即ESR=Total Reads/DNB Number×100%；ACCGRR（Averageaccumulative continued Good Reads Rate）是指连续平均累计好读长（Reads）占比，反应测序的连续测序质量。其他可能的参数还包括dnbnum、总ESR（totalESR）、FIT、光强（Intensity）、信噪比（SNR）和拆分率（Split Rate）等。其中，Dnbnum，或Dnbnum DNB数量是指以测序通道（Lane）为单位，可以正常产出数据的FOV数。totalESR，是指总ESR（total Effective Spot Rate，totalESR）。FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强（Intensity）分布越集中（一种碱基光强很高，其他三种碱基光强很低），其FIT值越高；光强（Intensity）分布越离散，FIT值越低。光强（Intensity）：数据来源rawInt，为经过图像处理后提取的光强值，对rawInt的每一列进行排序（descend），取最高1%~15%（百分位点）之间的光强求平均值，为原始光强估计值。信噪比（SNR）：以单个DNB SNR的计算为例，碱基A最大光强（max intensity）作为信号，碱基CGT则为背景，CGT光强的方差则为噪声，那么碱基A第一循环（cycle1）的SNR（CGT_dev）估算方法为：碱基A的信号取所有碱基A的平均值（A_mean），碱基CGT矩阵的方差（CGT_dev）则为噪声，那么A_SNR=A_mean/CGT_dev。拆分率（SplitRate）：拆分出来的标签序列（Barcode）能对应到所有（例如1-96号）标签序列的比率。标签拆分率为标签序列拆分率柱状图，纵坐标为拆分率，横坐标为标签序列编号，只显示拆分率大于1%的标签序列情况。这些参数与测序产量有一定相关性，并且不同测序产量与测序产量的相关性不同，即权重不同，因此，根据这些参数并依据一定的模型可以预测测序产量。

S302：对上述测序产量数据进行预处理以排除异常值。

测序产量数据中的异常值是无效数据，这些无效数据的存在，会影响结果的准确性，因此需要将这些无效数据排除出去。这些异常值包括，例如：（a）数据异常情况：从芯片原理上判断为不可能达到的异常值；（b）极低数值：数值为非零的极低数值，可能由于数据拆分或数据写入导致的数据异常；（c）测序失败和/或测序进行中的测序数据。

S303：根据上述测序产量数据和上述测序数据，基于线性回归模型建立测序产量预测模型。

本发明实施例中，基于线性回归模型建立测序产量预测模型，这是针对现有测序参数来选择的，不排除选择一般线性模型、广义线性模型、混合模型、线性回归模型、曲线回归模型、非线性回归模型、逻辑回归模型、神经网络模型、贝叶斯网络模型等其他模型算法及机器学习建模方法来预测产量。

S304：利用上述测序产量预测模型和上述测序产量数据，模拟计算测序产量预测值并统计残差。

在本发明一个实施例中，预测测序产量=-510.983-5.581*BIC平均值+1474.396*ESR平均值+498.429*ACCGRR平均值，残差标准差是26.2。上述预测测序产量公式中，给出了测序数据（或参数）BIC、ESR和ACCGRR的系数，即各自的权重系数。

在本发明一个实施例中，预测测序产量=-972.058+7.471*BIC平均值+142.704*FIT平均值+982.381*ESR平均值-0.003*光强平均值±置信区间，其中置信区间可以根据需要选择，如果设置为53.34058（2倍标准差）则准确率为95%；设置为26.67029（标准差）则准确率为83%。上述预测测序产量公式中，给出了测序数据（或参数）BIC、FIT、ESR和光强的系数，即各自的权重系数。

如图4所示，本发明一种实施例中提供一种建立测序产量预测模型的方法，包括：

S401：获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，该测序数据至少包括BIC、ESR和ACCGRR。

S402：对上述测序产量数据进行预处理以排除异常值。

S403：根据上述测序产量数据和上述测序数据，基于线性回归模型建立测序产量预测模型。

S404：利用上述测序产量预测模型和上述测序产量数据，模拟计算测序产量预测值并统计残差。

S405：根据上述测序产量预测值的大小对上述测序产量预测值分类，并按类统计上述残差的标准差。

本发明一个实施例中，根据测序产量预测值的大小，将测序产量预测值分为三类：（a）测序产量预测值小于600M，（b）测序产量预测值介于600M与700M之间，（c）测序产量预测值大于700M。其中，测序产量预测值小于600M的残差标准差为50.99097；测序产量预测值介于600M与700M之间的残差标准差为50.99097；测序产量预测值大于700M的残差标准差为19.37197。

相应地，得到：

预测产量及其概率区间=-510.983-5.581*BIC+1474.396*ESR+498.429*ACCGRR±置信区间，其中置信区间可以自由选择，采用2倍标准差则产量预测准确率为96.88%，采用1倍标准差则为89.24%。

S406：根据上述测序产量预测模型和上述残差及输入的测序过程中的测序数据，预测测序完成时的测序产量。

S407：按时间周期统计新产生的测序产量数据的残差，并在监测到上述残差的标准差偏离原残差的标准差时，重新进行上述测序产量预测模型的建立。

由于残差会随着新产生的测序产量数据的增加而变化，因此，为了修正这种变化，需要随时监控残差的标准差，当残差的标准差偏离原残差的标准差时，需要重新进行上述测序产量预测模型的建立，即重新回到步骤S401，并重复步骤S401至步骤S405的过程。这种监控时间周期进行，其中时间周期例如可以是每一天、每一个星期、每一个月或每三个月等作为一个周期。所谓“残差的标准差偏离原残差的标准差”，是指相比原残差的标准差，残差标准差发生明显变化，可以人为确定偏差/变化的标准，例如可以认为相比原残差的标准差变化超过10%、20%或50%等，是残差的标准差偏离原残差的标准差的情况。

对应于上述建立测序产量预测模型的方法，本发明实施例还提供一种建立测序产量预测模型的装置，如图5所示，包括：数据获取单元501，用于获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，上述测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；异常值排除单元502，用于对上述测序产量数据进行预处理以排除异常值；模型建立单元503，用于根据上述测序产量数据和上述测序数据，基于线性回归模型建立测序产量预测模型；预测值计算单元504，用于利用上述测序产量预测模型和上述测序产量数据，模拟计算测序产量预测值并统计残差。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

因此，本发明实施例还提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如本发明的建立测序产量预测模型的方法。

以下详述本发明的建立测序产量预测模型的方法的思路，以便更好地理解本发明。

首先，本发明通过模拟探究了测序产量（Reads，M）与分类参数的关系，其中分类参数包括测序类型、建库类型、建库和测序试剂、测序仪器等，结果显示：不同测序类型，测序产量分布不同；而不同建库类型，测序产量差异不明显；不同批次建库和测序试剂，可以看出有些批次对应的测序产量波动比较大；不同测序仪器，可以看出部分测序仪器波动较大。这些分类参数对于测序产量结果没有呈现出较强的相关性。

其次，分析测序产量与拆分率（Split Rate）和ESR的相关性。结果如图6所示，结果显示；ESR与测序产量强相关，相关系数为0.96，测序产量与拆分率（Split Rate）弱相关，相关系数为0.2。

再次，分析测序产量与不同阶段参数值的相关性，取五个点的ESR值，即全部循环的平均值（average）、前25个循环（pre25cycle）、前100个循环（pre100cycle）、前125个循环（pre125cycle）、前150个循环（pre150cycle），结果如图7所示，结果显示：测序产量与各个阶段ESR值都显著相关，可以利用ESR值来预测下机产量。

最后，分析不同的数值参数之间的相关性，选取6个参数，即BIC、FIT、ESR、ACCGRR、Intensity和SNR，随机挑选4个测序（4个RUN）的数据，对以上参数进行相关性分析。结果如图8-11所示，结果显示：（a）ESR与intensityA强相关（>0.9）；（b）intensityA与ACCGRR、ESR与intensity强相关（0.7-0.8）；（c）FIT与SNRA、intensity与FIT、FIT与ESR相关，但四次结果波动较大（0.2-0.7）。

随后进行数据聚类——主成分分析并聚类，图12示出了聚类的簇数量k（Numberof cluster k）与模型解释力（Gap statistic）之间的关系曲线，示出了在聚类的簇数量k是3时，是一个较佳簇数量（Optimal numberof clusters）的转折点，是一个较优的选择。

之后，采用逐步线性回归建立模型，比较各个模型的数据，调整模型加入次序，比较获得的模型性能。结果如图13所示，结果显示，从模型1至6，随着预测变量的增加（从a至f），R方（R²）逐渐增加，表明模型的解释性能越来越好。并且，采用逐步线性回归建立模型，比较各个模型的数据，调整模型加入次序，比较获得的模型显著性（Sig.）。结果如图14所示，结果显示，从模型1至6，模型显著性均较好，体现为Sig.为0。

最终，得到模型结论参数，如图15所示。

得到可用的测序产量预测模型公式：

预测产量及其概率区间=-510.983-5.581*BIC+1474.396*ESR+498.429*ACCGRR±置信区间，其中置信区间可以自由选择，采用2倍标准差（52.4）则产量预测准确率为96.88%，采用1倍标准差（26.2）则产量预测准确率为89.24%，如图16所示。

测序到60个循环（cycle）时，按照上述测序产量预测模型公式进行测序产量预测，预测结果如图17至图19所示。

在另一个实施例中按照如下方法进行测序产量预测。

1、训练集&测试集：

使用random函数将用于分析的数据分为训练集（data_train.txt）和测试集（data_test.txt）；使用训练集（data_train.txt）数据进行建模，使用测试集（data_test.txt）的数据观察预测是否准确。

采用线性回归模型，进入的方法为向后删除，添加自变量：ESR_平均值（ESR_average），BIC_平均值（BIC_average），FIT_平均值（FIT_average），光强_平均值（intensity_average），SNR_平均值（SNR_average）。得到的模型参数如下表1所示：

表1

模型公式如下：

预测产量=-972.058+7.471*BIC平均值+142.704*FIT平均值+982.381*ESR平均值-0.003*光强平均值±置信区间，其中置信区间可以根据需要选择，如果设置为53.34058（2倍标准差）则准确率为95%；设置为26.67029（标准差）则准确率为83%。将上述公式代入测试集（data_test.txt）中，发现预测效果高于预期，表2所示。

表2

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (11)

1.一种测序产量预测方法，其特征在于，所述方法包括：

获取测序过程中的测序数据，所述测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；其中BIC是指能够用于作碱基调取的DNA纳米球DNB的比率，ESR是指经过过滤最终测序所得的读长数占碱基调取能识别出来的DNB数量的比率，ACCGRR是指连续平均累计好读长占比，FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强分布越集中，其FIT值越高，光强分布越离散，FIT值越低；

根据所述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：

（a）测序产量= a+ b·BIC+ c·ESR+ d·ACCGRR±置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数，置信区间为标准差或其大于1的倍数；

（b）测序产量= m+ n·BIC+ p·FIT+ q·ESR+ r·光强±置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数，置信区间为标准差或其大于1的倍数；和

输出预测出的测序完成时的测序产量。

2.根据权利要求1所述的方法，其特征在于，所述a的取值为-510.983在±50%范围内的数值；所述b的取值为-5.581在±50%范围内的数值；所述c的取值为1474.396在±50%范围内的数值；所述d的取值为498.429在±50%范围内的数值；

任选地，所述m的取值为-972.058在±50%范围内的数值；所述n的取值为7.471在±50%范围内的数值；所述p的取值为142.704在±50%范围内的数值；所述q的取值为982.381在±50%范围内的数值；所述r的取值为-0.003在±50%范围内的数值。

3.根据权利要求1所述的方法，其特征在于，所述测序过程中的测序数据是测序进行至30-120个循环时的测序数据。

4.根据权利要求3所述的方法，其特征在于，所述测序过程中的测序数据是测序进行至60个循环时的测序数据。

5.一种测序产量预测装置，其特征在于，所述装置包括：

测序数据获取单元，用于获取测序过程中的测序数据，所述测序数据包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；其中BIC是指能够用于作碱基调取的DNA纳米球DNB的比率，ESR是指经过过滤最终测序所得的读长数占碱基调取能识别出来的DNB数量的比率，ACCGRR是指连续平均累计好读长占比，FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强分布越集中，其FIT值越高，光强分布越离散，FIT值越低；

测序产量预测单元，用于根据所述测序数据和如下测序产量预测模型公式（a）或（b），预测测序完成时的测序产量：

（a）测序产量= a+ b·BIC+ c·ESR+ d·ACCGRR±置信区间，其中，a为常数，b、c和d分别为相应测序数据的系数，置信区间为标准差或其大于1的倍数；

（b）测序产量= m+ n·BIC+ p·FIT+ q·ESR+ r·光强±置信区间，其中，m为常数，n、p、q和r分别为相应测序数据的系数，置信区间为标准差或其大于1的倍数；和

测序产量输出单元，用于输出预测出的测序完成时的测序产量。

6.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1至4任一项所述的测序产量预测方法。

7.一种建立测序产量预测模型的方法，其特征在于，所述方法包括：

获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，所述测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；其中BIC是指能够用于作碱基调取的DNA纳米球DNB的比率，ESR是指经过过滤最终测序所得的读长数占碱基调取能识别出来的DNB数量的比率， ACCGRR是指连续平均累计好读长占比，FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强分布越集中，其FIT值越高，光强分布越离散，FIT值越低；

对所述测序产量数据进行预处理以排除异常值；

根据所述测序产量数据和所述测序数据，基于线性回归模型建立测序产量预测模型；和

利用所述测序产量预测模型和所述测序产量数据，模拟计算测序产量预测值并统计残差。

8.根据权利要求7所述的方法，其特征在于，所述方法还包括：

根据所述测序产量预测值的大小对所述测序产量预测值分类，并按类统计所述残差的标准差。

9.根据权利要求7所述的方法，其特征在于，所述方法还包括：

根据所述测序产量预测模型和所述残差及输入的测序过程中的测序数据，预测测序完成时的测序产量；和

按时间周期统计新产生的测序产量数据的残差，并在监测到所述残差的标准差偏离原残差的标准差时，重新进行所述测序产量预测模型的建立。

10.一种建立测序产量预测模型的装置，其特征在于，所述装置包括：

数据获取单元，用于获取至少一个已完成测序的测序产量数据，以及相应测序产量数据所对应的测序数据，所述测序数据至少包括BIC、ESR和ACCGRR，或BIC、FIT、ESR和光强；其中BIC是指能够用于作碱基调取的DNA纳米球DNB的比率，ESR是指经过过滤最终测序所得的读长数占碱基调取能识别出来的DNB数量的比率，ACCGRR是指连续平均累计好读长占比，FIT用来反映各碱基信号和噪声的分布情况，校正后各碱基的光强分布越集中，其FIT值越高，光强分布越离散，FIT值越低；

异常值排除单元，用于对所述测序产量数据进行预处理以排除异常值；

模型建立单元，用于根据所述测序产量数据和所述测序数据，基于线性回归模型建立测序产量预测模型；

预测值计算单元，用于利用所述测序产量预测模型和所述测序产量数据，模拟计算测序产量预测值并统计残差。

11.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求7至9任一项所述的建立测序产量预测模型的方法。

Images