CN112365930B - 一种为基因数据库确定最佳序列比对阈值的方法 - Google Patents

Abstract

一种为基因数据库确定最佳序列比对阈值的方法，包括：1)获取蛋白质序列；2)从蛋白质序列中移除被包括在基因数据库中的序列，创建假基因数据集；3)将基因数据库中的蛋白质序列划分子类，作为真基因数据集；4)合并假基因数据集与真基因数据集，针对任意一条蛋白质序列，模拟高通量测序所产生的特定长度的DNA序列，得到模拟数据集；5)进行序列比对，对比对阈值进行取值；6)判定序列比对结果，计算真阳性、错配、假阳性、假阴性、真阴性的数量；7)计算灵敏度、准确度和马修斯相关系数；8)以相似度为X轴，E值为Y轴，灵敏度、准确度或马修斯相关系数为Z轴，绘制三维曲面图；9)在三维曲面图确定基因数据库的最佳序列比对阈值。

Description

一种为基因数据库确定最佳序列比对阈值的方法

技术领域

本发明属于生物技术领域，涉及一种结合了模拟数据集构建、序列比对效果评估等手段，为基因数据库确定最佳序列比对阈值的方法。

背景技术

近年来，高通量测序技术迅速发展。因其具有通量大、精确度高、信息量丰富等优点，高通量测序技术被广泛应用于微生物生态学研究，探究复杂微生物群体的多样性、群落结构及生态作用。尤其地，高通量测序技术解决了大部分微生物无法被分离培养的问题，为这部分微生物的代谢潜能及生态功能研究提供了有力的技术手段。

利用高通量测序技术获取微生物群体的宏基因组数据，将其在基因数据库中进行序列比对，是从人工系统或自然环境样品中检测特定基因的常用手段，能够克服大部分微生物无法被分离培养的局限性，弥补了纯培养技术的不足。利用高通量测序技术及序列比对方法，研究者能够探究人工系统(如污水生物脱氮反应器等)或自然生态系统(如农田土壤、深海热泉等)中特定基因的多样性及时空分布规律。例如，将宏基因组数据在抗性基因数据库ARDB及CARD(the Comprehensive Antibiotic Research Database)中进行序列比对，以检测抗生素抗性基因；在肽酶数据库MEROPS中进行序列比对，以检测蛋白质水解酶编码基因；在碳水化合物活性酶数据库CAZy中进行序列比对，以检测碳水化合物活性酶编码基因等。

通过序列比对的方法从海量的宏基因组数据中检测特定基因，其比对准确性及比对速度高度依赖于所选定的比对阈值，即相似度和E值。在不适宜的比对阈值下检测特定基因，会造成比对结果不可靠或比对速度很慢，如过低的相似度或过高的E值会使得某些非目标基因被错误地提取出来，假阳性率随之上升。

然而，目前为基因数据库确定序列比对阈值，往往基于经验主义而缺乏科学的数据计算及分析。具体来说，序列比对阈值的确定主要以同类基因或其他基因的已有的研究结果为参考，为基因数据库选定普遍认为较高的相似度及较低的E值，以控制假阳性率，或者直接采用序列比对工具(如CEMERA)的默认比对阈值。在这种情况下，确定的相似度和E值具有普适性而缺乏针对性，难以使基因数据库表现出最佳的比对性能。因此，确立一种为基因数据库快速、科学地确定最佳序列比对阈值(相似度及E值)的方法，对人工系统或自然环境样品中特定基因的序列鉴定及丰度、多样性分析至关重要。

发明内容

针对目前普遍基于经验主义为基因数据库选定最佳序列比对阈值(相似度及E值)的问题，本发明首次提出，利用蛋白质序列数据库(例如，Swiss-Prot蛋白质序列数据库)的全部蛋白质序列构建模拟数据集，将模拟序列在基因数据库中作序列比对，并以灵敏度、准确度以及马修斯相关系数(Matthews correlation coefficient，MCC)评估序列比对效果，从而为基因数据库快速、科学地确定最佳序列比对阈值(即，相似度和E值)的关键技术。

本发明提供的一种为基因数据库确定最佳序列比对阈值(相似度及E值)的方法，其包括：

步骤1)，从第一数据库中获取蛋白质序列并下载到本地；

步骤2)，从步骤1)所获取的蛋白质序列中移除已被包括在所述基因数据库中的序列，将余下的蛋白质序列作为假基因数据集；

步骤3)，根据功能将所述基因数据库中的蛋白质序列划分成子类，作为真基因数据集；

步骤4)，合并所述假基因数据集与所述真基因数据集，针对其中任意一条蛋白质序列，在第1～N个氨基酸中，随机选取任意一个位置往后切分形成若干条长度为N个氨基酸的序列，以模拟高通量测序所产生的预定长度的DNA序列，得到包含若干条模拟序列的模拟数据集；

步骤5)，用所述模拟数据集中的模拟序列在所述基因数据库中作序列比对，以等差数列对比对阈值进行取值，得到不同比对阈值下每条模拟序列的最佳匹配序列，所述比对阈值包括相似度和E值；

步骤6)，对每一次序列比对的结果进行判定，计算真阳性(True positive，TP)、错配(Cross-talking，CT)、假阳性(False positive，FP)、假阴性(False negative，FN)、真阴性(True negative，TN)的数量；

步骤7)，根据所述真阳性、所述错配、所述假阳性、所述假阴性以及所述真阴性的数量，分别计算每一次序列比对的灵敏度、准确度以及马修斯相关系数(Matthewscorrelation coefficient，MCC)；

步骤8)，以所述相似度为X轴，所述E值为Y轴，以所述灵敏度、所述准确度或所述马修斯相关系数为Z轴，分别绘制三维曲面图；

步骤9)，在所述三维曲面图上选点，当所述相似度及所述E值下所述灵敏度、所述准确度以及所述马修斯相关系数均具有超出预定值的值时，将此时的所述相似度和所述E值确定为所述基因数据库的最佳序列比对阈值。

优选地，所述第一数据库为UniProt蛋白质数据库的Swiss-Prot蛋白质序列数据库(https://www.uniprot.org/uniprot/？query＝reviewed:yes)。

优选地，将步骤3)中所得到的所述真基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为真序列，这些真序列能够模拟高通量测序所产生的DNA序列。将步骤2)中所得到的所述假基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为假序列，这些假序列能够模拟高通量测序所产生的DNA序列。在步骤6)，判定所述真阳性、所述错配、所述假阳性、所述假阴性以及所述真阴性的依据如下：假设一条序列是真序列且属于所述真基因数据集中的第一子类：当将其与所述基因数据库比对后判定为真序列，且归类到所述第一子类中时，则判定结果为真阳性，当将其与所述基因数据库比对后判定为真序列，但归类到所述第一子类以外的子类中时，则判定结果为错配，当将其与所述基因数据库比对后判定为假序列，即未比对上时，则判定结果为假阴性；假设一条序列是假序列：当将其与所述基因数据库比对后判定为真序列并归类到任意子类中，则判定结果为假阳性；当将其与所述基因数据库比对后判定为假序列，即未比对上时，则判定结果为真阴性。

优选地，所述灵敏度、所述准确度以及所述马修斯相关系数的计算公式分别如下：

其中，TP是所述真阳性的数量，FN是所述假阴性的数量，FP是所述假阳性的数量，TN是所述真阴性的数量。

优选地，在步骤8)，使用Matlab软件来绘制所述三维曲面图。

本发明利用蛋白质数据库尤其是Swiss-Prot蛋白质序列数据库的全部蛋白质序列构建模拟数据集，并将模拟序列在基因数据库中作序列比对，得到不同比对阈值下每条模拟序列的最佳匹配序列，这使得最终选定的比对阈值更符合实际情况，适用于高通量测序所产生的DNA序列。

本发明在模拟数据集的构建过程中，采用了特定的随机切分方法。具体来说，合并假基因数据集与真基因数据集，对其中任意一条蛋白质序列，在第1～N个氨基酸中，随机选取一个位置，以N个氨基酸的长度往后切分形成若干条长度为N个氨基酸的序列，并去掉首尾长度不足N个氨基酸的短序列，从而得到模拟数据集。采用这种随机切分方法，既使得模拟序列具有随机性，又充分利用了Swiss-Prot蛋白质序列数据库的蛋白质序列，还保证了模拟数据集不过于庞大，从而提高了比对阈值选定过程的效率。

本发明对每一次序列比对的结果进行判定，得到TP、CT、FP、FN、TN的数量，从而分别计算灵敏度、准确度及MCC，以这三个指标评估序列比对效果(越趋近于1则越优越)，以获取最佳序列比对阈值(相似度及E值)，此时序列比对具有较高的灵敏度、准确度及MCC，即基因数据库具有优越的比对性能，比对结果准确、可靠。

与现有技术相比，本发明的有益效果是：

利用本发明提供的方法，能够为基因数据库快速、科学地确定最佳序列比对阈值(即，相似度和E值)，使基因数据库表现出优越的比对性能，较之经验主义方法更符合实际情况，适用于高通量测序所产生的DNA序列，为人工系统与自然生态系统中特定基因的序列鉴定及数据分析，提供必要的分析手段。

附图说明

图1是本发明的为基因数据库确定最佳序列比对阈值的方法的流程图。

图2是实施例中以相似度为X轴，以E值为Y轴，以灵敏度为Z轴绘制的三维曲面图。

图3是实施例中以相似度为X轴，以E值为Y轴，以准确度为Z轴绘制的三维曲面图。

图4是实施例中以相似度为X轴，以E值为Y轴，以MCC为Z轴绘制的三维曲面图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和具体实施方式对本发明进一步详细说明。

图1示出了本发明所提供的为基因数据库确定最佳序列比对阈值的方法的流程图。具体而言，该方法主要包括以下步骤。

步骤1)，从第一数据库中获取蛋白质序列并下载到本地。该第一数据库优选为UniProt蛋白质数据库的Swiss-Prot蛋白质序列数据库(https://www.uniprot.org/uniprot/？query＝reviewed:yes)。

步骤2)，从步骤1)所获取的蛋白质序列中移除已被包括在基因数据库中的序列，将余下的蛋白质序列作为假基因数据集。

步骤3)，根据功能将基因数据库中的蛋白质序列划分成子类，作为真基因数据集。

步骤4)，合并上述得到的假基因数据集与真基因数据集，针对其中任意一条蛋白质序列，在第1～N个氨基酸中，随机选取一个位置，以N个氨基酸的长度往后切分形成若干条长度为N个氨基酸的序列，并去掉首尾长度不足N个氨基酸的短序列，通过上述随机切分方法，产生若干条长度为N的氨基酸序列，以模拟高通量测序所产生的特定长度的DNA序列，得到包含若干条模拟序列的模拟数据集。

步骤5)，用步骤4)得到的模拟数据集中的模拟序列在基因数据库中作序列比对，以等差数列对比对阈值进行取值，得到不同比对阈值下每条模拟序列的最佳匹配序列，比对阈值包括相似度和E值。

步骤6)，对每一次序列比对的结果进行判定，计算真阳性(True positive，TP)、错配(Cross-talking，CT)、假阳性(False positive，FP)、假阴性(False negative，FN)、真阴性(True negative，TN)的数量。将步骤3)中得到的真基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为真序列，这些真序列能够模拟高通量测序所产生的DNA序列。将步骤2)中得到的假基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为假序列，这些假序列能够模拟高通量测序所产生的DNA序列。如表1所示，在步骤6)，判定真阳性、错配、假阳性、假阴性以及真阴性的依据如下：假设一条序列是真序列且属于真基因数据集中的第一子类(例如，子类A)：当将其与基因数据库比对后判定为真序列，且归类到子类A中时，则判定结果为真阳性，当将其与基因数据库比对后判定为真序列，但归类到第一子类以外的子类(例如，子类B)中时，则判定结果为错配，当将其与基因数据库比对后判定为假序列，即未比对上时，则判定结果为假阴性；假设一条序列是假序列：当将其与基因数据库比对后判定为真序列并归类到任意子类中，则判定结果为假阳性；当将其与基因数据库比对后判定为假序列，即未比对上时，则判定结果为真阴性。应当理解的是，第一子类可以选自真基因数据集中的任意一个子类。

表1

步骤7)，根据真阳性、错配、假阳性、假阴性以及真阴性的数量，分别计算每一次序列比对的灵敏度、准确度以及马修斯相关系数(Matthews correlation coefficient，MCC)，灵敏度、准确度以及马修斯相关系数的计算公式分别如下：

其中，TP是真阳性的数量，FN是假阴性的数量，FP是假阳性的数量，TN是真阴性的数量。

步骤8)，以相似度为X轴，以E值为Y轴，以灵敏度、准确度或马修斯相关系数为Z轴，分别绘制三维曲面图。优选使用Matlab软件来绘制三维曲面图。

步骤9，在步骤8)所绘制的三维曲面图上选点，当在某一相似度和E值下灵敏度、准确度以及马修斯相关系数均超出预定值时，将此时的相似度和E值确定为基因数据库的最佳序列比对阈值。

实施例

以下以群体感应基因数据库为例，详细说明确定最佳序列比对阈值(即，相似度及E值)的过程。

步骤1)，在UniProt蛋白质数据库的Swiss-Prot蛋白质序列数据库(https://www.uniprot.org/uniprot/？query＝reviewed:yes)中，将全部蛋白质序列下载到本地，共557134条蛋白质序列。

步骤2)，从步骤1)得到的557134条蛋白质序列中，移除已被包括在群体感应基因数据库中的蛋白质序列(共245条)，将余下的556889条蛋白质序列作为假群体感应基因数据集。针对假群体感应基因数据集中的蛋白质序列，在序列名称后标注“F”，以序列F4HRV8为例，标注如下：

步骤3)，将群体感应基因数据库中的蛋白质序列，根据功能类别划分成11个子类，作为真群体感应基因数据集。这11个子类包括：酰化高丝氨酸内酯类(Acyl-homoserinelactones，AHL)的合成酶编码基因LuxI、LuxM、HdtS及降解酶编码基因Acylase、Lactonase，霍乱弧菌I类自诱导物(Cholerae autoinducer-1，CAI-1)的合成酶编码基因CqsA，可扩散因子(Diffusible signal factor，DSF)的合成酶编码基因RpfF，环二鸟苷酸(c-di-GMP)的合成酶编码基因DGC及降解酶编码基因PDE，以及2型自诱导物(Autoinducer-2，AI-2)的合成酶编码基因LuxS及降解酶编码基因LsrG。对真群体感应基因数据集中的蛋白质序列，在序列名称后标注由大写字母“T”和基因名称组成的标签，如下表2所示。

表2

步骤4)，合并上述得到的假群体感应基因数据集与真群体感应基因数据集，对其中任意一条蛋白质序列，在第1～50个氨基酸中，随机选取一个位置，以50个氨基酸的长度往后切分，形成若干条长度为50个氨基酸的序列，并去掉首尾长度不足50个氨基酸的短序列。通过上述随机切分方法，产生3485393条长度为50个氨基酸的序列，以模拟高通量测序所产生的长度为150bp的DNA序列，包括3451824条标注为“F”的序列及33569条标注为“T”的序列，即为模拟数据集。

步骤5)，用模拟数据集中的模拟序列在群体感应基因数据库中作序列比对，共完成2500次比对，得到不同比对阈值下每条模拟序列的最佳匹配序列。以等差数列对相似度和E值进行取值，相似度在0～100％内进行取值且公差为2％，Lg(E值)在0～(-10)内进行取值且公差为-0.2。

步骤6)，对2500次序列比对的结果进行判定，分别计算TP、CT、FP、FN、TN的数量。

步骤7)，针对2500次序列比对，基于TP、CT、FP、FN、TN的数量，分别计算灵敏度、准确度及MCC。

步骤8)，使用Matlab软件，以相似度为X轴，以E值为Y轴，以灵敏度、准确度或MCC为Z轴，绘制三维曲面图，分别如图2～4所示。

步骤9)，根据图2～4，当相似度为50％、E值为1e-5时，序列比对的灵敏度、准确度及MCC均超出了0.950的预定值，其中灵敏度＝0.986，准确度＝0.957，MCC＝0.971，表明群体感应基因数据库表现出优越的比对性能，此时的相似度及E值作为最终确定的最佳序列比对阈值。

最后，所属领域的普通技术人员应当理解：以上列举的仅是本发明的具体实施例子，并非旨在说明本公开的范围(包括权利要求)被限于为群体感应基因数据库选定最佳序列比对阈值(相似度及E值)。因此，凡在本发明的精神和原则之内，所做的任何省略、修改、等同替换、改进等为特定基因数据库选定最佳序列比对阈值(相似度及E值)的方法，均应包含在本发明的保护范围之内。

Claims (5)

1.一种为基因数据库确定最佳序列比对阈值的方法，其包括：

步骤1)，从第一数据库中获取蛋白质序列并下载到本地；

步骤2)，从步骤1)所获取的蛋白质序列中移除已被包括在所述基因数据库中的序列，将余下的蛋白质序列作为假基因数据集；

步骤3)，根据功能将所述基因数据库中的蛋白质序列划分成子类，作为真基因数据集；

步骤4)，合并所述假基因数据集与所述真基因数据集，针对其中任意一条蛋白质序列，在第1～N个氨基酸中，随机选取任意一个位置往后切分形成若干条长度为N个氨基酸的序列，以模拟高通量测序所产生的具有预定长度的DNA序列，得到包含若干条模拟序列的模拟数据集；

步骤5)，用所述模拟数据集中的模拟序列在所述基因数据库中作序列比对，以等差数列对比对阈值进行取值，得到不同比对阈值下每条模拟序列的最佳匹配序列，所述比对阈值包括相似度和E值；

步骤6)，对每一次序列比对的结果进行判定，计算真阳性、错配、假阳性、假阴性、真阴性的数量；

步骤7)，根据所述真阳性、所述错配、所述假阳性、所述假阴性以及所述真阴性的数量，分别计算每一次序列比对的灵敏度、准确度以及马修斯相关系数；

步骤8)，以相似度为X轴，以E值为Y轴，以灵敏度、准确度或马修斯相关系数为Z轴，分别绘制三维曲面图；

步骤9)，在所述三维曲面图上选点，当在某一相似度和E值下灵敏度、准确度以及马修斯相关系数均超出预定值时，将此时的相似度和E值确定为所述基因数据库的最佳序列比对阈值。

2.根据权利要求1所述的为基因数据库确定最佳序列比对阈值的方法，其特征在于：所述第一数据库为UniProt蛋白质数据库的Swiss-Prot蛋白质序列数据库。

3.根据权利要求1所述的为基因数据库确定最佳序列比对阈值的方法，其特征在于：

将步骤3)中得到的所述真基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为真序列，

将步骤2)中得到的所述假基因数据集中的蛋白质序列通过步骤4)处理后得到的模拟序列作为假序列，

在步骤6)，判定所述真阳性、所述错配、所述假阳性、所述假阴性以及所述真阴性的依据如下：

假设一条序列是真序列且属于所述真基因数据集的第一子类：当将其与所述基因数据库比对后判定为真序列，且归类到所述第一子类中时，则判定结果为真阳性，当将其与所述基因数据库比对后判定为真序列，但归类到所述第一子类以外的子类中时，则判定结果为错配，当将其与所述基因数据库比对后判定为假序列，即未比对上时，则判定结果为假阴性；

假设一条序列是假序列：当将其与所述基因数据库比对后判定为真序列并归类到任意子类中，则判定结果为假阳性；当将其与所述基因数据库比对后判定为假序列，即未比对上时，则判定结果为真阴性。

4.根据权利要求1所述的为基因数据库确定最佳序列比对阈值的方法，其特征在于：

在步骤7)，所述灵敏度、所述准确度以及所述马修斯相关系数的计算公式分别如下：

其中，MCC是马修斯相关系数，TP是所述真阳性的数量，FN是所述假阴性的数量，FP是所述假阳性的数量，TN是所述真阴性的数量。

5.根据权利要求1或2所述的为基因数据库确定最佳序列比对阈值的方法，其特征在于：在步骤8)，使用Matlab软件来绘制所述三维曲面图。

Images