CN111599413B - 一种测序数据的分类单元组分计算方法 - Google Patents

Abstract

本发明涉及一种测序数据的分类单元组分计算方法。本发明基于“测序读出序列的次分类单元频率”指标及其计算框架，用于衡量序列比对结果中分类单元误比对的情况，能够有效去除分类单元组分计算中的假阳性结果，提高组分计算的特异性和准确性。同时，本发明还通过剔除异常分类单元后重新统计的策略实现了误比对序列向真实组分结果的回归，有效校正了分类单元丰度的定量结果。

Description

一种测序数据的分类单元组分计算方法

技术领域

本发明涉及生信分析领域，具体涉及一种测序数据的分类单元组分计算方法。

技术背景

感染类疾病是由病原微生物引起的一类病症，感染源种类繁多且病患众多，给全球各国的公共卫生带来重大影响。根据世界卫生组织的数据，以2016年为例，仅下呼吸道感染即在全球导致约300万人死亡。与此同时，感染类疾病的盲目治疗所带来的抗生素滥用问题也有日益严重的趋势。而对于感染病原体的精准检测是解决上述问题的最重要一环。

感染类疾病病原体检测的传统手段是微生物培养，但培养具有检测时间长、灵敏度较低的缺点。聚合酶链式反应(以下简称PCR)方法检测时间较短、灵敏度较高，但每次只能检测一种病原体。基于测序技术的病原体检测直接对样本的全部DNA进行测序分析，具有检测范围广、灵敏度高的特点。

纳米孔测序技术是最近几年兴起的新一代测序技术。纳米孔测序技术弥补了第二代测序平台的劣势，不仅测序片段的读长较二代测序高出一到两个数量级，而且建库和测序时间短。此外，测序设备小巧便携，数据能够实时获取并进行后续分析，较好地解决了测序场地的限制，以及报告反馈的延误。因此这一技术非常适宜于感染微生物病原体检测的应用。本技术领域常规的纳米孔测序的物种组分计算流程如下：

1.测序运行过程中使用ONT MinKNOW软件实时收集原始测序数据；

2.使用ONT Albacore或ONT Guppy软件转换原始电信号数据，生成碱基序列；

3.使用Minimap2软件，进行基于hg38人类参考基因组的宿主序列去除；

4.使用What′s In My Pot？(WIMP)软件计算物种组分，最后进行物种丰度过滤。

其中，使用WIMP软件进行物种组分计算的过程为：

1.使用Centrifuge软件进行序列比对；

2.根据每一条测序读出序列的比对情况判断其所属物种；

3.统计每个物种的支持该物种的测序读出序列数，计算出物种的绝对丰度和相对丰度；

4.对物种结果进行用户定义的丰度过滤(例如使用1％的相对丰度阈值)。

但上述测序数据的常规分析方法存在物种结果假阳性较高(特异性较低)的缺陷，对于病原体结果准确性有很大影响。如何合理地去除序列比对过程中引入的物种假阳性，是现有技术中亟需解决的技术问题。

有鉴于此，提出本发明。

发明内容

本发明需要解决的核心问题是，如何通过数据分析方法尽量去除序列比对过程所引入的假阳性分类单元结果。考虑到测序数据的序列比对过程中，由于亲缘关系邻近的分类单元的基因组之间存在一定比例的相似序列，来源于某一分类单元的测序读出序列可能误比对到其他邻近分类单元的基因组，进而引起分类单元组分计算的错误。面对误比对现象，如果仅通过依次评估每条测序读出序列的比对情况来确定某分类单元是否存在，假阳性结果总会得到部分保留，本发明开创性地采取一种包含了比对状况整体统计(summary-based analysis)的计算框架来确定分类单元组分结果的真实与否。

现有的分类单元组分计算方法只是使用丰度筛选(例如去除相对丰度1％以下的分类单元)的策略来进行消极的假阳性去除，而并没有构建一种通过评估比对结果的整体性分布规律来判断误比对引入的假阳性分类单元的积极策略。

本发明考虑到误比对问题源于亲缘关系邻近的分类单元的基因组之间存在的相似序列。因此，真阳性比对结果大多并不来源于分类单元之间的相似序列，而假阳性比对结果则主要来源于分类单元之间的相似序列。那么一定可以通过一些指标或指标的组合来体现出真阳性比对与假阳性比对在整体统计上所具有的差异。

基于上述原理，首先本发明发现，如果以优先比对到的分类单元为单位，将测序读出序列分为不同小组，真阳性分类单元组的测序读出序列的专一比对比例相对较高，而假阳性分类单元组的测序读出序列的专一比对比例相对较低。之后通过数据探究，本发明发现与“专一比对比例”指标原理相同的“测序读出序列的次分类单元频率”指标具有更好的区分度(见图2)。于是本发明构建了一套包含“测序读出序列的分类单元频率”计算方法和基于“测序读出序列的次分类单元频率”指标来衡量序列比对的分析结果中假阳性分类单元结果的量化方法。这一分类单元层面的筛选方法能够有效去除宏基因组分类单元组分计算中的假阳性分类单元结果，提高特异性和准确性。最后本发明还通过“剔除异常分类单元后重新统计”的策略实现误比对的测序读出序列向真实组分结果的回归，在提高分类单元结果特异性的同时，有效校正了分类单元丰度的定量结果。

因此，本发明的第一目的是提供一种测序数据的分类单元组分计算方法及其系统。

本发明的第二目的是提供一种降低测序数据生信分析中假阳性分类单元结果的方法及其系统。

基于上述目的，本发明提供如下技术方案：

本发明提供一种测序数据生信分析方法，其特征在于，包括如下步骤：

步骤1)测序数据比对步骤；

步骤2)按分类单元分组步骤；

步骤3)测序读出序列的次分类单元频率统计步骤；

在一些实施方式中，所述步骤1)，采用保留非专一比对结果的比对软件对测序读出序列进行序列比对，优选的，所述软件为BLASTN软件。

在一些实施方式中，所述步骤2)，基于步骤1)比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组，即每组测序读出序列所优先支持的分类单元相同。

在一些实施方式中，所述步骤3)，针对步骤2)的测序读出序列分组，统计每组测序读出序列的次分类单元频率。在一些具体实施方式中，所述步骤2)的具体步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元(例如所比对上的物种集合)，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，即为该组测序读出序列的次分类单元频率。

本发明提供一种降低测序数据生信分析中假阳性分类单元结果的方法，其特征在于，包括如下步骤：

步骤1)测序数据比对步骤；

步骤2)按分类单元分组步骤；

步骤3)测序读出序列的次分类单元频率统计步骤；

步骤4)假阳性分类单元排除步骤；

在一些实施方式中，所述步骤1)，采用保留非专一比对结果的比对软件对测序读出序列进行序列比对，优选的，所述软件为BLASTN软件。

在一些实施方式中，所述步骤2)，基于步骤1)比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组，即每组测序读出序列所优先支持的分类单元相同。

在一些实施方式中，所述步骤3)，针对步骤2)的测序读出序列分组，统计每组测序读出序列的次分类单元频率。在一些具体实施方式中，所述步骤2)的具体步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元(例如所比对上的物种集合)，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，即为该组测序读出序列的次分类单元频率。

在一些实施方式中，所述步骤4)，将每组的次分类单元频率数值与次分类单元频率阈值进行比较，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除，所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果(alignment)进行剔除。

在一些实施方式中，所述测序数据来源于泌尿感染宏基因组测序数据。

在一些实施方式中，所述次分类单元频率阈值为15-30％，优选的，为20％。

在一些实施方式中，所述次分类单元频率阈值也可以通过如下方法计算：使用一定量样本作为训练集，通过对比传统培养和/或PCR的鉴定结果，确认常规生信分析中的真阳性与假阳性结果；重新进行生信分析，将每个样本的测序读出序列按照比对结果优先支持的分类单元分为不同组，即每组测序读出序列优先支持同一分类单元；统计每组读出序列所比对上的分类单元的分类单元频率，并得出该组测序读出序列的次分类单元频率；统计优先支持的分类单元在传统培养和/或qPCR结果中为真阳性的每组测序读出序列的次分类单元频率，再统计优先支持的分类单元在培养结果或qPCR结果中为假阳性的每组测序读出序列的次分类单元频率，即获得可区分两者的阈值。

本发明还提供一种降低测序数据生信分析中假阳性分类单元结果的系统，其特征在于，所述系统包括如下模块：

模块1)测序数据比对模块；

模块2)按分类单元分组模块；

模块3)测序读出序列的次分类单元频率统计模块；

模块4)假阳性分类单元排除模块；

在一些实施方式中，所述模块1)，采用保留非专一比对结果的比对软件对测序读出序列进行序列比对，优选的，所述软件为BLASTN软件。

在一些实施方式中，所述模块2)，基于模块1)获得的比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组，即每组测序读出序列所优先支持的分类单元相同。

在一些实施方式中，所述模块3)，针对模块2)的测序读出序列分组，统计每组测序读出序列的次分类单元频率。

在一些实施方式中，模块3)行使的具体步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，即为该组测序读出序列的次分类单元频率。

在一些实施方式中，所述模块4)，将每组的次分类单元频率数值与次分类单元频率阈值进行比较，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除，所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果(alignment)进行剔除。

在一些实施方式中，所述测序数据来源于泌尿感染宏基因组测序数据。

在一些实施方式中，所述次分类单元频率阈值为15-30％，优选的，为20％。

在一些实施方式中，上述降低测序数据生信分析中假阳性分类单元结果的方法或模块针对的测序数据来自于第二代测序平台或第三代测序平台；优选的，来自于Illumina、ION TORRENT、PacBio、Roche、Helicos、ABI或纳米孔测序平台；更优选的，来自于纳米孔测序平台。

在一些实施方式中，上述降低测序数据生信分析中假阳性分类单元结果的方法或模块测序数据为基因组测序数据；优选的为宏基因组测序数据；更优选的，为泌尿感染宏基因组测序数据。

本发明还提供一种测序数据的分类单元组分计算方法，所述方法包括如下步骤：

步骤1)测序数据比对步骤；

步骤2)按分类单元分组步骤；

步骤3)测序读出序列的次分类单元频率统计步骤；

步骤4)假阳性分类单元排除步骤；

步骤5)分类单元丰度统计步骤。

在一些实施方式中，所述步骤1)，采用保留非专一比对结果的比对软件对测序读出序列进行序列比对，优选的，所述软件为BLASTN软件。

在一些实施方式中，所述步骤2)，基于步骤1)比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组，即每组测序读出序列所优先支持的分类单元相同。

在一些实施方式中，所述步骤3)，针对步骤2)的测序读出序列分组，统计每组测序读出序列的次分类单元频率。在一些具体实施方式中，所述步骤2)的具体步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，即为该组测序读出序列的次分类单元频率。

在一些实施方式中，所述步骤4)，将每组的次分类单元频率数值与次分类单元频率阈值进行比较，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除，所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果(alignment)进行剔除。

在一些实施方式中，所述步骤5)，对剔除异常分类单元之后的比对结果，重新按照比对结果优先支持的分类单元对测序读出序列进行分组，并统计每个分组的序列数(即该分类单元的绝对丰度)和所占总读出序列数的比例(即该分类单元的相对丰度)。

在一些实施方式中，所述测序数据来源于泌尿感染宏基因组测序数据。

在一些实施方式中，所述步骤4)中的次分类单元频率阈值，可以是本领域对于特定样本类型所已知的经验值，通常为15-30％，优选的，为20％；

在一些实施方式中，所述次分类单元频率阈值也可以通过如下方法计算：使用一定量样本作为训练集，通过对比传统培养和/或PCR的鉴定结果，确认常规生信分析中的真阳性与假阳性结果；重新进行生信分析，将每个样本的测序读出序列按照比对结果优先支持的分类单元分为不同组，即每组测序读出序列优先支持同一分类单元；统计每组读出序列所比对上的分类单元的分类单元频率，并得出该组测序读出序列的次分类单元频率；统计优先支持的分类单元在传统培养和/或qPCR结果中为真阳性的每组测序读出序列的次分类单元频率，再统计优先支持的分类单元在培养结果或qPCR结果中为假阳性的每组测序读出序列的次分类单元频率，即获得可区分两者的阈值。

本发明还提供一种测序数据的分类单元组分计算系统，其特征在于，所述系统包括如下模块：

模块1)测序数据比对模块；

模块2)按分类单元分组模块；

模块3)测序读出序列的次分类单元频率统计模块；

模块4)假阳性分类单元排除模块；

模块5)分类单元丰度统计模块

在一些实施方式中，所述模块1)，采用保留非专一比对结果的比对软件对测序读出序列进行序列比对，优选的，所述软件为BLASTN软件。

在一些实施方式中，所述模块2)，基于模块1)获得的比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组，即每组测序读出序列所优先支持的分类单元相同。

在一些实施方式中，所述模块3)，针对模块2)的测序读出序列分组，统计每组测序读出序列的次分类单元频率。

在一些实施方式中，所述模块3)行使的具体步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，即为该组测序读出序列的次分类单元频率。

在一些实施方式中，所述步骤4)，将每组的次分类单元频率数值与次分类单元频率阈值进行比较，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除；优选的，所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果(alignment)进行剔除。

在一些实施方式中，所述测序数据来源于泌尿感染宏基因组测序数据。

在一些实施方式中，所述步骤4)中的次分类单元频率阈值，可以是本领域对于特定样本类型所已知的经验值，通常为15-30％，优选的，为20％。

在一些实施方式中，所述次分类单元频率阈值也可以通过如下方法统计：使用一定量样本作为训练集，通过对比传统培养和/或PCR的鉴定结果，确认常规生信分析中的真阳性与假阳性结果；重新进行生信分析，将每个样本的测序读出序列按照比对结果优先支持的分类单元分为不同组，即每组测序读出序列优先支持同一分类单元；统计每组读出序列所比对上的分类单元的分类单元频率，并得出该组测序读出序列的次分类单元频率；统计优先支持的分类单元在传统培养和/或qPCR结果中为真阳性的每组测序读出序列的次分类单元频率，再统计优先支持的分类单元在培养结果或qPCR结果中为假阳性的每组测序读出序列的次分类单元频率，即获得可区分两者的阈值。

在一些实施方式中，所述模块5)，对剔除异常分类单元之后的比对结果，重新按照比对结果优先支持的分类单元对测序读出序列进行分组，并统计每个分组的序列数(即该分类单元的绝对丰度)和所占总读出序列数的比例(即该分类单元的相对丰度)。

在一些实施方式中，上述分类单元组分计算方法或模块中的测序数据来自于第二代测序平台或第三代测序平台；优选的，来自于Illumina、ION TORRENT、PacBio、Roche、Helicos、ABI或纳米孔测序平台；更优选的，来自于纳米孔测序平台。

在一些实施方式中，上述分类单元组分计算方法或模块测序数据为基因组测序数据；优选的为宏基因组测序数据；更优先的为泌尿感染宏基因组测序数据。

本发明有益技术效果：

1.本发明是对常规的物种组分计算方法的改进，开创性地提出基于序列比对结果的整体统计量进行物种筛选的方法，提供一种新的生信分析方法。

2.本发明首次通过引入“测序读出序列的次分类单元频率”的计算，解决了原有的常规物种组分计算方法难以解决的误比对引入的假阳性分类单元结果的去除问题，有效提升了病原体检测的准确性和特异性。

3.本发明的计算框架独立于具体测序平台的选择，可适用于第二代测序技术和第三代测序技术等多种平台的测序数据，并可应用于不同来源或不同物种的检测样本。

附图说明

图1：通过序列比对进行分类单元组分计算的理想情况和现实情况；

图2：使用36例泌尿测试样本探究“专一比对比例”与“次分类单元频率”对于真阳性和假阳性分类单元结果的区分度情况；

图3：常规方法的物种检出结果与培养及qPCR验证的一致性情况(绝对丰度阈值为100条序列)；

图4：本发明方法的物种检出结果与培养及qPCR验证的一致性情况(绝对丰度阈值为100条序列)；

图5：常规方法的物种检出结果与培养及qPCR验证的一致性情况(绝对丰度阈值为200条序列)；

图6：本发明方法的物种检出结果与培养及qPCR验证的一致性情况(绝对丰度阈值为200条序列)。

具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。

部分术语定义

除非在下文中另有定义，本发明具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本发明。

如本发明中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由...组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。

本发明中的术语“大约”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的土10％，优选土5％。

此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本发明描述的实施方案能以不同于本发明描述或举例说明的其它顺序实施。

以下术语或定义仅仅是为了帮助理解本发明而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。

本发明中的术语“测序读出序列”：英文为“read”或“reads”，指测序平台读出的一条或一组核酸序列。

本发明中的术语“比对结果”：英文为“alignment”，指一条测序读出序列与一条参考序列之间的对应结果，一条测序读出序列可以同时具有多个比对结果。

本发明中的术语“分类单元”：英文为“taxon”，指共享一些相同特征的一群生物体，如原生动物门(Protozoa)、灵长目(Primates)、金黄色葡萄球菌(Staphylococcusaureus)、肠道沙门氏菌肠道亚种(Salmonella enterica subsp.enterica)等都分别是一个分类单元。不同的分类单元可能具有不同的分类级别(如原生动物对应“门”这个分类级别，灵长目对应“目”这个分类级别，金黄色葡萄球菌对应“种”这个分类级别，肠道沙门氏菌肠道亚种对应“亚种”这个分类级别)，也可能具有相同的分类级别(如金黄色葡萄球菌、表皮葡萄球菌、溶血葡萄球菌、人葡萄球菌都对应“种”这个分类级别)。物种、菌株是两类特殊的分类单元，也是本发明的分类单元组分计算方法优选关注的类别。

本发明中的术语“物种”：英文为“species”，物种是一类特殊的分类单元，指可以交配并繁衍后代的一群生物体。

本发明中的术语“(一组分类单元)互斥”，指在该组分类单元中任意选择两个分类单元A和B，满足分类单元A既不包含分类单元B，也不被分类单元B所包含。例如，“大肠杆菌、肠道沙门氏菌、肺炎克雷伯菌”这一组三个分类单元互斥；而“克雷伯菌、肺炎克雷伯菌”这两个分类单元不互斥。

本发明中的术语“一条测序读出序列比对到某分类单元”、“一条测序读出序列比对上某分类单元”、“一条测序读出序列支持某分类单元”：指该条测序读出序列的比对结果中包含了来自该分类单元的参考序列。

本发明中的术语“一条测序读出序列优先比对到某分类单元”、“一条测序读出序列优先比对上某分类单元”、“一条测序读出序列优先支持某分类单元”：指该分类单元是该条测序读出序列的比对结果所对应的所有互斥的分类单元中比对分值最高的。比对分值最高的分类单元的判定方法是将该条测序读出序列的比对结果按照各自对应的分类单元进行分组并对比对分值求和，将比对分值之和最高的分类单元判定为比对分值最高的分类单元。

本发明中的术语“专一比对”：指一条测序读出序列的所有比对结果所对应的参考序列都来自同一分类单元。

本发明中的术语“非专一比对”：指一条测序读出序列的比对结果中同时包含了来自两个或两个以上互斥的分类单元的参考序列。

本发明中的术语“误比对”：指实际来自某一分类单元的测序读出序列的比对结果包含了来自于既不包含前述分类单元、也不被前述分类单元所包含的另一个分类单元的参考序列。需要说明的是，这里的比对结果在计算上没有错误，碱基一致率可以很高，但比对上的分类单元与样本的来源分类单元事实不符。这种误比对一般发生在基因组亲缘关系比较邻近的分类单元之间。

本发明中的术语“某分类单元在一组测序读出序列中的分类单元频率(taxonfrequency)”：指在给定的一组测序读出序列中，支持该分类单元的测序读出序列占序列总数的比例。

本发明中的术语“一组测序读出序列的次分类单元频率(minor taxonfrequency)”：指每条测序读出序列优先支持的分类单元都相同的一组测序读出序列中，数值第二高的分类单元频率。即如果一组测序读出序列中每条序列优先支持的分类单元都相同，则寻找被该组中至少一条测序读出序列所比对上的所有分类单元，对于找出的这组分类单元，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，其中数值最高的必然是该组序列优先支持的分类单元所对应的分类单元频率(100％)，而全部百分比中数值第二高的，即为该组测序读出序列的次分类单元频率(minor taxon frequency)。

本发明中的术语“假阳性分类单元”：指分类单元组分计算结果为阳性，但事实上不存在于样本中的分类单元。“假阳性分类单元”的一种特殊情况是“假阳性物种”。

本发明中的术语“异常分类单元”：指在本发明方法中通过“测序读出序列的次分类单元频率”指标被判为异常的分类单元。“异常分类单元”的一种特殊情况是“异常物种”。

本发明中的“测序数据的分类单元组分计算方法”：通常是指对测序下机的数据通过生信分析获得测序数据中各种分类单元是否出现及所占比例情况的计算方法，本发明优选的是基于“测序读出序列的次分类单元频率”指标来衡量序列比对结果中误比对情况，进而获得测序数据的分类单元组分情况，该方法能够有效去除分类单元组分计算中的假阳性结果。可以理解的是，本发明通过引入“测序读出序列的次分类单元频率”的计算，解决了原有的常规分类单元组分计算方法难以解决的误比对引入的假阳性分类单元的去除问题，有效提升了病原体检测的特异性和准确性；同时，本发明的计算框架独立于具体测序平台的选择，并不受限于测序平台，可适用于第二代测序技术和第三代测序技术等多种平台的测序数据，本发明仅是优选于纳米孔测序平台；而且该计算框架也独立于具体测序数据来源，本领域可以理解，本发明的计算框架是针对任何同源序列的误比对，因此，序列来源并不会限制本发明的应用，除了本发明优选的宏基因组数据来源外，其他基因组或基因数据来源同样适于本发明。

本发明整体的技术思路，示例性的解释如下，但不作限制：

1)序列比对使用能保留非专一比对结果的比对软件获得测序读出序列的原始比对结果：

为了达到保留非专一比对的所有结果的需求，本发明使用诸如BLASTN软件的Megablast方法进行序列比对。

2)测序读出序列的分类单元频率的计算：

序列比对后，对于一组测序读出序列，统计该组测序读出序列所比对上的所有互斥分类单元中，比对上每个分类单元的测序读出序列数占测序读出序列总数的比例。

举例：如果共有4条测序读出序列，原始比对结果分别为：第1条和第2条都专一比对上大肠杆菌，第3条专一比对上肺炎克雷伯菌，第4条同时比对上大肠杆菌、肺炎克雷伯菌和产气克雷伯菌。那么这4条测序读出序列中大肠杆菌、肺炎克雷伯菌、产气克雷伯菌的分类单元频率分别为75％、50％、和25％。

3)测序读出序列的次分类单元频率的计算：

当给定一组所优先支持的分类单元皆相同的测序读出序列时，该组测序读出序列所比对上的所有互斥分类单元的分类单元频率中最大的，必然是所优先支持的分类单元的分类单元频率——100％，而这些互斥分类单元所对应的分类单元频率中排行第二高的频率，则为该组测序读出序列的次分类单元频率。根据实际经验，如果该组测序读出序列优先支持的分类单元为真阳性，其专一比对比例一般较高，而次分类单元频率一般较低。

4)基于“测序读出序列的次分类单元频率”指标来判定假阳性分类单元：

将样本的测序读出序列按照优先支持的分类单元情况分为不同小组，即每组测序读出序列的优先支持分类单元相同。然后统计每组测序读出序列的次分类单元频率。如果该组序列的次分类单元频率数值大于阈值(阈值计算方法见下文)，则判定为异常分类单元。

5)次分类单元频率阈值的计算：

次分类单元频率阈值可以是本领域知晓的经验值，比如15-30％，优选的，20％左右。本发明意识到对于不同类型的感染类疾病，次分类单元频率阈值可能因病原体类别的差异而有一定差异，可针对不同类型的疾病预先使用上述方法确定次分类单元频率阈值。示例性的，可以通过如下方法进行计算：使用一定数量的样本作为训练集，通过培养方法获得样本的分类单元鉴定结果。使用常规组分计算方法获得样本的分类单元组分结果。对于培养结果与常规生信结果中不一致的分类单元结果，通过qPCR方法进行鉴定，找出生信结果中的真阳性与假阳性结果。重新进行数据分析，将每个样本的测序读出序列按照比对结果优先支持的分类单元分为不同组，即每组测序读出序列优先支持同一分类单元。然后统计每组读出序列所比对上的分类单元的分类单元频率，并得出该组测序读出序列的次分类单元频率。统计优先支持的分类单元在培养结果或qPCR结果中为真阳性的每组测序读出序列的次分类单元频率，再统计优先支持的分类单元在培养结果或qPCR结果中为假阳性的每组测序读出序列的次分类单元频率，然后得出一个可以最大程度区分两者的阈值。

6)剔除异常分类单元后重新统计，实现误比对序列向真实组分结果的回归，有效校正分类单元丰度的定量结果的方法：

分析获得样本的异常分类单元列表后，将原始比对结果文件中比对上异常分类单元的所有比对结果进行剔除。注意是以比对结果(alignment)为单位进行剔除，而不是以测序读出序列为单位进行剔除。然后重新进行所有分类单元的丰度统计，得出分类单元组分计算结果。

下面结合具体实施例来阐述本发明。

实施例1发明设计

本发明不考虑a)采样、建库、测序过程引入样本的污染，和b)条码拆分错误引入样本的污染。因为前者作为实验操作引入的污染，可以通过在运营中建立阴性对照等实验方式进行污染排查，这一手段不在本发明的讨论范围；后者一方面可以通过选用区分效果更好的条码体系来解决(不在本发明的讨论范围中)，另一方面通过一些不在本发明讨论范围的定量阳性对照实验，可以得到错误引入比例的经验值，并用于丰度筛选以解决此假阳性。

1、序列比对时，一条测序读出序列可能同时比对上多个互斥的分类单元，参见图1，在理想情况下，如果测序得到的每条read都单一地比对上某一分类单元，本发明即可直接计算出每个分类单元的支持reads数，得出分类单元组分结果。

但在现实中，本发明发现部分reads可能不只比对上某一分类单元，而是同时比对上多个互斥的分类单元。

本发明知道，对于测序得到的单条read来说，它事实上只来源于某一个分类单元。这条read之所以同时比对上多个互斥的分类单元，不是因为比对软件的比对出错，而是因为亲属关系邻近的分类单元的基因组之间存在大量的序列相似性。

2、一条实际来源于某分类单元的read可能只比对到另外的分类单元，而未比对到原分类单元。

一条测序read进行序列比对后的可能情况有如下：

a)一条实际来源于分类单元A的read可能专一地比对上原分类单元A；

b)一条实际来源于分类单元A的read可能同时比对上原分类单元A和另外的一种或几种分类单元，且比对分值最高的分类单元仍然是原分类单元A；

c)一条实际来源于分类单元A的read可能同时比对上原分类单元A和另外的一种或几种分类单元，且比对分值最高的物种不是原分类单元A；

d)一条实际来源于分类单元A的read可能比对上一种或几种分类单元，但没有比对上原分类单元A。

e)一条实际来源于分类单元A的read可能没有比对上任何分类单元。

表1为发酵乳杆菌的纯菌Spike-in实验的序列比对结果，其展示了一条来源于分类单元A但并未比对上分类单元A参考序列的read。发酵乳杆菌纯菌的测序结果中，有31684条reads优先比对上真实来源的分类单元发酵乳杆菌，有48条reads优先比对上样本中并不存在的分类单元口乳杆菌，且48条中有29条reads专一比对到了这一并不存在的分类单元口乳杆菌。

表1：发酵乳杆菌的纯菌Spike-in实验的序列比对结果

	物种名	优先比对reads数	相对丰度	专一比对reads比例
					1	Lactobacillus_fermentus	31684	99.295％	85.179％
2	Lactobacillus_oris	48	0.150％	60.417％

3、常规的组分计算无法排除部分假阳性分类单元结果

计算分类单元组分的过程一般包括序列比对、结果筛选和统计，最终输出组分计算结果。常规方法中的筛选过程一般包括以下几类：

首先，alignment可能需要进行筛选：

-基于比对分值(BLASTN软件使用bitscore作为比对分值)的绝对大小

-基于比对分值的相对大小

-基于比对分值的排序名次(例如每条read只保留分值最高的1个alignment)

然后，read可能需要进行筛选：

-基于比对分值的大小

-基于是否专一比对到某一分类单元

最后，分类单元结果可能需要进行筛选：

-基于支持该分类单元的reads的绝对丰度和相对丰度

-基于已知的病原体列表等。

从前一部分的实验数据可知，实际来源于某一分类单元A的reads不一定优先比对到原分类单元A，甚至可能不比对上原分类单元A。这说明无法找出一种方法来准确判断每条read所实际对应的分类单元。即使丢弃所有非专一比对的reads，只保留专一比对上某一分类单元的reads，根据上面提到的发酵乳杆菌(Lactobacillus fermentus)纯菌Spike-in实验结果，结果中仍然会输出假阳性的分类单元：口乳杆菌(Lactobacillus oris)。

有鉴于此，本发明提出基于大量reads比对结果的整体统计量进行物种筛选的方法，即不是只根据每条read自身的比对数据来筛选该read的比对信息，而是把根据计算规则划为同组的所有reads的比对情况放在一起综合考虑作为背景，并用于保留和去除比对信息。

具体实验方法如下：

1)序列生成：ONT GridION测序平台产生的数据，通过ONT Guppy软件将电信号转换为碱基信号，得到fastq格式的序列信息。

2)样本拆分：使用ONT Guppy软件，根据文库的barcode序列将序列拆分成属于不同样本的序列集合，同时去除接头序列。

3)序列质控：统计序列的长度和质量值(quality score)，对于Nanopore长读长数据，过滤掉长度小于500bp或平均测序质量值小于8的序列。

4)宿主去除：使用minimap2软件，将过滤后的序列比对到人类基因组(hg38版本)，保留未比对上人类基因组的序列，进行下一步分析。

5)序列比对：使用BLASTN软件进行序列比对(参数设置如下：-task megablast-word_size 28-gapopen 0-gapextend 2-penalty-4-reward 1-max_hsps 1-max_target_seqs 50)。

6)比对结果的筛选：对于每条read，找出比对结果中的最高比对分值，去除该条read的比对分值低于最高比对分值的95％的比对结果。

7)序列按优先支持分类单元情况进行分组：以分类单元为单位，将reads分为优先支持不同分类单元的不同分组。

8)统计异常分类单元：计算每组reads的次分类单元频率。如果该组reads的次分类单元频率数值大于阈值(根据泌尿项目得出的经验值为20％)，则判定为异常分类单元。

9)剔除异常分类单元的比对结果：分析获得样本的异常分类单元列表后，将原始比对结果文件中比对上异常分类单元的所有比对结果进行剔除。

10)计算分类单元丰度：对剔除异常分类单元之后的比对结果，重新对优先支持不同的分类单元的reads进行分组，并统计每组的reads条数(该分类单元的绝对丰度)及占所有有比对结果的reads的百分比(该分类单元的相对丰度)。

11)通过用户定义的绝对丰度阈值和相对丰度阈值进行分类单元结果筛选。

实施例2临床实验验证

本发明收集泌尿感染患者的114例尿液样本，对每个样本进行了微生物培养和PCR检测，以微生物培养和PCR检测的综合结果判断样本中是否存在某分类单元。其中有36例样本用于计算次分类单元频率阈值；余下的78例样本用于计算常规生信分析方法与本发明的新方法的分类单元结果的表现，以说明新方法相比原有常规方法的结果提升效果，具体如下：

1、本发明取36例泌尿感染患者样本作为训练集，通过培养方法获得样本的分类单元鉴定结果。使用常规组分计算方法获得样本的分类单元组分结果。对于培养结果与常规生信结果中不一致的分类单元结果，通过qPCR方法进行鉴定，找出生信结果中的真阳性与假阳性结果。重新进行数据分析，将每个样本的测序读出序列按照比对结果优先支持的分类单元分为不同组，即每组测序读出序列优先支持同一分类单元。然后统计每组读出序列所比对上的分类单元的分类单元频率，并得出该组测序读出序列的次分类单元频率。统计优先支持的分类单元在培养结果或qPCR结果中为真阳性的每组测序读出序列的次分类单元频率，再统计优先支持的分类单元在培养结果或qPCR结果中为假阳性的每组测序读出序列的次分类单元频率，然后得出一个可以最大程度区分两者的阈值。

结果如图2所示，最终确定泌尿感染样本的次分类单元频率阈值为20％。

2、使用ONT GridION测序平台进行感染样本的测序，每一批下机数据在完成碱基序列生成、样本划分、质量控制、和宿主序列去除之后，分别使用常规生信分析方法与本发明方法来进行物种层级的分类单元组分计算。

常规生信分析的组分计算流程和本发明的组分计算流程所使用的参考序列数据库(nt)和丰度筛选阈值都相同，所不同的是常规分析方法使用Centrifuge软件进行序列比对，而本发明的分析方法使用BLASTN软件进行序列比对，然后通过“测序读出序列的次分类单元频率”指标排除异常物种并重新统计结果。具体实验步骤参见实施例1。

图3至图6展示本发明和常规生信分析两种方法的结果，可以看到，通过本发明的异常分类单元剔除方法，组分计算的特异性和一致率都有了明显改进。对于绝对丰度阈值为100条序列的结果，常规方法检出结果的特异性为20.6％，一致率为87.5％；本发明的新方法的检出结果的特异性为29.2％，一致率为90.8％。对于绝对丰度阈值为200条序列的结果，常规方法检出结果的特异性为25.8％，一致率为88.0％；本发明的新方法的检出结果的特异性为36.4％，一致率为92.0％。因此，新方法的特异性和准确性具有明显提升。

以上对本申请具体实施方式的描述并不限制本申请，本领域技术人员可以根据本申请作出各种改变或变形，只要不脱离本申请的精神，均应属于本申请所附权利要求的范围。

Claims (15)

1.一种降低测序数据生信分析中假阳性分类单元结果的方法，其特征在于，所述方法包括如下步骤：

步骤1），测序数据比对步骤；

步骤2），按分类单元分组步骤；

所述步骤2），基于步骤1）比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组；

步骤3），测序读出序列的次分类单元频率统计步骤；

所述步骤3），针对步骤2）的测序读出序列分组，统计每组测序读出序列的次分类单元频率；

步骤4）假阳性分类单元排除步骤；

所述步骤4），将每组的次分类单元频率数值与次分类单元频率阈值进行比较，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除；所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果进行剔除。

2.根据权利要求1所述的降低测序数据生信分析中假阳性分类单元结果的方法，其特征在于，所述步骤3）中次分类单元频率的统计步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，为该组测序读出序列的次分类单元频率。

3.根据权利要求2所述的降低测序数据生信分析中假阳性分类单元结果的方法，其特征在于，所述步骤1）中采用保留非专一比对结果的比对软件对测序读出序列进行序列比对。

4.根据权利要求3所述的降低测序数据生信分析中假阳性分类单元结果的方法，其特征在于，所述比对软件为BLASTN。

5.一种测序数据的分类单元组分计算方法，其特征在于，所述方法包括权利要求1-4任一所述步骤，并进一步包括：

步骤5），分类单元丰度统计步骤；

所述步骤5），对步骤4）剔除异常分类单元之后的比对结果，重新按照比对结果优先支持的分类单元对测序读出序列进行分组，并统计每个分组的序列数和所占总读出序列数的比例。

6.根据权利要求1-5任一所述的方法，其特征在于，所述测序数据来自于第二代测序平台或第三代测序平台。

7.根据权利要求6所述的方法，其特征在于，所述测序数据来自于Illumina、IONTORRENT、PacBio、Roche、Helicos、ABI或纳米孔测序平台。

8.根据权利要求1-5任一所述的方法，其特征在于，所述测序数据为基因组测序数据。

9.根据权利要求8所述的方法，其特征在于，所述测序数据为宏基因组测序数据。

10.根据权利要求9所述的方法，其特征在于，所述测序数据为泌尿感染宏基因组测序数据。

11.根据权利要求10所述的方法，其特征在于次分类单元频率阈值为15-30%。

12.一种降低测序数据生信分析中假阳性分类单元结果的系统，其特征在于，所述系统包括如下模块：

模块1），测序数据比对模块；

模块2），按分类单元分组模块；

模块3），测序读出序列的次分类单元频率统计模块；

模块4），假阳性分类单元排除模块；

所述模块1），采用保留非专一比对结果的比对软件对测序读出序列进行序列比对；

所述模块2），基于模块1）获得的比对结果，按照比对结果优先支持的分类单元进行测序读出序列的分组；

所述模块3），针对模块2）的测序读出序列分组，统计每组测序读出序列的次分类单元频率；

所述模块4），比较每组的次分类单元频率数值与次分类单元频率阈值，若大于次分类单元频率阈值，则确定该组测序读出序列所优先支持的分类单元为异常分类单元并剔除；所述剔除为，获得样本的全部异常分类单元后，将原始比对结果文件中比对上异常分类单元的所有比对结果进行剔除。

13.根据权利要求12所述的降低测序数据生信分析中假阳性分类单元结果的系统，其特征在于，所述模块3）中次分类单元频率的统计步骤如下：对于每组测序读出序列，找出至少被一条该组测序读出序列所比对上的所有互斥分类单元，对于找出的分类单元集合，分别计算每个分类单元在该组测序读出序列中的支持序列数占该组测序读出序列总数的百分比，这些百分比中数值第二大的，为该组测序读出序列的次分类单元频率。

14.一种测序数据的分类单元组分计算系统，其特征在于，所述系统包括权利要求12或13所述系统，并进一步包括：

模块5），分类单元丰度统计模块；

所述模块5），对剔除异常分类单元之后的比对结果，重新按照比对结果优先支持的分类单元对测序读出序列进行分组，并统计每个分组的序列数和所占总读出序列数的比例。

15.根据权利要求12-14任一所述系统，其特征在于，所述测序数据来自于泌尿感染宏基因组测序数据。

Images