氨基酸可信度和修饰位点定位的评估方法

摘要

本发明提供了一种氨基酸可信度评估模型训练方法。该方法包括：根据包含待训练氨基酸的训练肽段产生所述待训练氨基酸的背景肽段集合；从所述训练肽段和所述待训练氨基酸提取多个特征；以所提取的多个特征作为输入向量，以所述待训练氨基酸是否正确作为输出，训练分类模型，得到氨基酸可信度评估模型。本发明获得氨基酸可信度评估模型可用于氨基酸可信度评估和修饰位点定位的评估，提高了氨基酸可信度评估的准确率并且改善了修饰位点定位的评估性能。

说明书

技术领域

本发明涉及生物技术领域，尤其涉及一种氨基酸和修饰位点定位的评估方法。

背景技术

质谱技术已经成为生物学家分析生物样品的常规手段，其中，肽段和蛋白质鉴定方法学已经成为关键的一环。目前，基于串联质谱数据的肽段鉴定方法主要分两类：数据库搜索方法和从头测序方法。数据库搜索方法严重依赖于数据库的质量，如果正确肽段不在数据库中，那么鉴定结果会出现错误。从头测序方法不依赖于数据库信息，直接从谱图中获取肽段序列，从而发现一些不在数据库中的新肽段，例如，包括突变以及意外修饰等。目前，有越来越多的从头测序算法，包括：SHERENGA，PEAKS，PepNovo，pNovo，pNovo+，UniNovo，Novor和支持意外修饰鉴定的Open-pNovo。

然而，由于从头测序没有使用数据库信息作为先验，难免会将非常相似的肽段序列报告出来，从而导致错误率非常高。根据文献报道，从头测序获得的高分结果中存在近40％的错误率，因此，从头测序领域中如何控制假发现率(False Discovery Rate，FDR)仍然是个急需解决的问题。

根据经验，从头测序结果中的肽段序列存在部分连续序列正确而剩余序列错误的现象，基于这一特征，可以评估肽段序列中每个氨基酸的可信度，将高可信的氨基酸组成的子序列提取出来作为序列标签，然后使用基于序列标签的方式去检索数据库，从而能够获得肽段序列的报告。然而，目前为止，在现有技术中并没有相关文献报道如何具体评估氨基酸的可信度，对于氨基酸可信度评估的准确性也没有进行深入的评测。

因此，需要对现有技术进行改进，以准确地评估氨基酸的可信度并进而降低从头测序中检测肽段序列的错误率。

发明内容

因此，本发明的目的在于克服上述现有技术的缺陷，提供一种氨基酸可信度的评估方法以及修饰位点定位的评估方法。

根据本发明的第一方面，提供了一种氨基酸可信度评估模型训练方法。该方法包括以下步骤：

步骤1：根据包含待训练氨基酸的训练肽段产生所述待训练氨基酸的背景肽段集合；

步骤2：从所述训练肽段和所述待训练氨基酸提取多个特征；

步骤3：以所提取的多个特征作为输入向量，以所述待训练氨基酸是否正确作为输出，训练分类模型，得到氨基酸可信度评估模型。

在本发明的氨基酸可信度评估模型训练方法中，步骤1包括：

对于所述待训练氨基酸，枚举预定长度的子序列，其中，该子序列包含所述待训练氨基酸以及所述训练肽段中的其它氨基酸；

从所述训练肽段中枚举质量等于所述子序列质量的氨基酸全排列；

将所述氨基酸全排列与所述训练肽段中的剩余序列进行拼接，以获得所述待训练氨基酸的背景肽段集合。

在本发明的氨基酸可信度评估模型训练方法中，步骤2包括：

计算所述训练肽段的肽谱匹配打分psm₁、谱峰强度匹配比例psm₂和谱峰匹配数目比例psm₃，分别作为第一特征、第二特征和第三特征；

计算所述待训练氨基酸的背景肽段集合中的最好背景肽段的肽谱匹配打分psm′₁、谱峰强度匹配比例psm′₂和谱峰匹配数目比例psm′₃，并计算所述训练肽段打分与最好的背景肽段打分的分差，表示为psm₁-psm′₁、psm₂-psm′₂和psm₃-psm′₃，分别作为第四特征、第五特征和第六特征，其中，所述最好背景肽段是所述待训练氨基酸的背景肽段集合中肽谱匹配打分最高的背景肽段；

计算所述待训练氨基酸的位置信息、类别信息以及所述训练肽段的长度信息，分别作为第七特征、第八特征和第九特征。

在本发明的氨基酸可信度评估模型训练方法中，在步骤3中，所述分类模型包括支持向量机、决策树，随机森林、贝叶斯网络中的任一项。

根据本发明的第二方面，提供了一种氨基酸可信度的评估方法。该评估方法包括：

步骤51：根据包含待评估氨基酸的原始肽段产生所述待评估氨基酸的背景肽段集合；

步骤52：从所述原始肽段和所述待评估氨基酸提取多个特征；

步骤53：将所提取的特征输入根据本发明的氨基酸可信度评估模型训练方法所得到的氨基酸可信度评估模型，以获得所述待评估氨基酸的可信度打分分布。

在本发明的氨基酸可信度的评估方法中，还包括：

将所述待评估氨基酸的可信度打分分布拟合为Gamma分布；

基于所述Gamma分布计算所述待评估氨基酸的假发现率：

其中，FAR表示待评估氨基酸的假发现率，p_w和p_r分别表示错误、正确氨基酸的先验概率，Γ(X|α_w,β_w)表示错误氨基酸中高于打分阈值X的分布面积，Γ(X|α_r,β_r)表示正确氨基酸中高于打分阈值X的分布面积，X表示待评估氨基酸的打分阈值，α_w,β_w表示错误氨基酸打分分布的Gamma参数，α_r,β_r表示正确氨基酸打分分布的Gamma参数。

根据本发明的第三方面，提供了一种修饰位点定位的评估方法。该评估方法包括：

对于给定的肽段序列枚举能够发生磷酸化修饰的候选修饰位点；

根据本发明的氨基酸可信度的评估方法获得每个候选位点发生磷酸化修饰的可信度打分。

在本发明的修饰位点定位的评估方法中，还包括采用以下公式计算每个候选修饰位点发生磷酸化修饰的概率：

其中，p_i表示第i个修饰位点的先验概率，s_i表示候选的磷酸化位点i的可信度打分，t_i表示候选的磷酸化位点i是否发生磷酸化，如果t_i等于1，则表示发生磷酸化；t_i等于0，表示没发生磷酸化。

与现有技术相比，本发明的优点在于：提出了利用机器学习方法来评估氨基酸的可信度，准确率较高；第一次提出了氨基酸水平的假发现率(False Amino-acid Rate，FAR)的概念，并将其用于从头测序领域的质量控制；将对氨基酸可信度的评估与对修饰位点定位的评估两个问题进行了统一，并且修饰位点定位的评估性能也得以改善。

附图说明

以下附图仅对本发明作示意性的说明和解释，并不用于限定本发明的范围，其中：

图1示出了根据本发明一个实施例的氨基酸可信度模型训练方法的流程图；

图2示出了根据本发明一个实施例的产生背景肽段的示意图；

图3示出了根据本发明的一个实施例的氨基酸可信度的评估方法的流程图；

图4示出了氨基酸可信度打分分布与Gamma分布的示意；

图5示出了根据本发明一个实施例的修饰位点定位的评估方法的流程图。

具体实施方式

为了使本发明的目的、技术方案、设计方法及优点更加清楚明了，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

简言之，本发明的氨基酸可信度的评估方法包括两个过程，第一个过程是利用机器学习的方法通过训练获得用于评估氨基酸可信度的模型，第二个过程是利用训练好的模型获得待评估氨基酸的可信度。图1示出了根据本发明一个实施例的氨基酸可信度评估模型训练方法的流程图。

步骤S110，选择训练样本

在此步骤中，将选择用于机器学习的训练样本，包括正样本和负样本。

在一个实例中，选择训练样本的过程包括：

步骤111：使用数据库搜索用于获取训练样本所使用的生物样品对应的数据集，将假发现率FDR≤1％的结果作为标注集。

数据集是指将生物样品(里面包含了非常多的肽段，需要解析这些肽段信息)放入质谱仪器得到的数据集合，该数据集合中通常包含几万张谱图，每个谱图对应一条肽段序列，需要计算机去解析每张谱图的肽段序列。

将质谱仪器产生的数据集使用数据库(即已知的基因库)进行搜索，是指对数据集在该数据库中进行匹配打分，以找到数据集中每个谱图对应的在数据库中存在且打分最好的肽段序列。

为了保证所获取的训练样本的准确性，选择假发现率FDR≤1％的搜索结果作为标注集。

步骤112，在同样的数据集上，使用从头测序软件进行测序，以在不存在数据库信息的情况下，直接通过谱图信息来找到肽段序列。

步骤113，对于从头测序获得的每个氨基酸，如果氨基酸与标注集上的氨基酸类型一致，则认为是正样本；否则认为是负样本。

通过此步骤S110，可以获得用于训练的正样本肽段序列以及负样本肽段序列，在下文中，将用于训练的肽段序列称为训练肽段或训练肽段序列，训练肽段序列包含的氨基酸称为待训练氨基酸。

步骤S120、产生待训练氨基酸的背景肽段。

对于训练肽段序列a₁a₂…a_l，假设参与训练的氨基酸a_i，i的取值为1至l，其中，l表示肽段的长度，一般取值在6到30之间。根据本发明的一个实施例，产生背景肽段的步骤包括：

步骤121：对于待训练的氨基酸a_i，枚举所有k长的子序列，其中，k表示枚举的子序列长度，一般取值2到5之间，如果k取值过长，则算法会较慢；

步骤122：假设k选择为3，则有三种子序列形式，即a_i-2a_i-1a_i，a_i-1a_ia_i+1和a_ia_i+1a_i+2；

步骤123：枚举质量等于这三种子序列质量的所有氨基酸排列，分别对应三个集合：S₁，S₂和S₃；

步骤124：将三个集合与训练肽段序列中的剩余序列进行拼接得到背景肽段：a₁…a_i-3S₁a_i+1…a_l，a₁…a_i-2S₂a_i+2…a_l和a₁…a_i-1S₃a_i+3…a_l，由此，对于S₁，S₂和S₃中的每个集合通过拼接得到了背景肽段，每个背景肽段均是一个集合，其包含非常多的背景肽段，例如，背景肽段a₁…a_i-3S₁a_i+1…a_l即是一个集合。

在此步骤中，产生背景肽段的目的是，为了后续通过对比训练肽段和背景肽段的谱特征来判断氨基酸正确与否。参见图2示意的产生背景肽段的过程，其中，在图2的谱图中，横坐标m/z表示质荷比，即质量除以电荷，纵坐标指示谱峰的强度(intensity)。由图2可知，假设正确肽段序列(即训练肽段序列)为AQPSK，需要判断第一个氨基酸A的正确性，枚举质量等于AQPS质量的所有氨基酸排列，例如有：QAPS，APSQ，APQ S，…，TQPG。将所有背景氨基酸排列与训练肽段的剩余序列进行拼接产生背景肽段，包括：QAPSK，APSQK，APQSK，…,TQPGK。根据图2的相应谱图可以看出，训练肽段序列AQPSK的谱峰匹配数目为7，即y4、y3、y2、y1、b2、b3和b4，用score为7表示，而背景肽段QAPSK的谱峰匹配数目是6(即y3、y2、y1、b2、b3和b4)，背景肽段TQPGK的谱峰匹配数目是2(即y1和b4)，这两个背景肽段的谱峰匹配数目均低于训练肽段序列AQPSK的谱峰匹配数目7，因此，训练肽段序列中的氨基酸A可信度会比较高。

应理解的是，图2仅示意性示出了产生背景肽段的过程和意义，本发明是根据下文将详细介绍的利用提取背景肽段的多个特征通过机器学习方法来评估待训练氨基酸的正确性。

步骤S130、提取并选择训练肽段和待训练氨基酸的特征

该步骤的目的在于从训练肽段和待训练氨基酸的背景肽段中选择能有效评估待训练氨基酸可信度的特征。

在一个实例中，提取的特征包括下列项中的至少两项：1)训练肽段的肽谱匹配打分；2)训练肽段的谱峰强度匹配比例；3)训练肽段的谱峰数目匹配比例；4)训练肽段与打分最好的背景肽段的肽谱匹配分差；5)训练肽段与打分最好的背景肽段的谱峰强度匹配比例差；6)训练肽段与打分最好的背景肽段的谱峰数目匹配比例差；7)氨基酸位置信息(例如，从1到肽段长度l)；8)氨基酸类别信息；9)肽段长度信息。

具体而言，提取上述九维特征的过程包括：

计算训练肽段a₁a₂…a_l的肽谱匹配打分psm₁，训练肽段的谱峰强度匹配比例psm₂和训练肽段的谱峰数目匹配比例psm₃，作为特征1、特征2和特征3；

计算所有背景肽段的肽谱匹配打分、谱峰强度匹配比例和谱峰匹配数目比例。找到肽谱匹配打分最高的背景肽段，其对应的三个打分分别表示为psm′₁、psm′₂和psm′₃，计算训练肽段打分与最好的背景肽段打分的分差，表示为psm₁-psm′₁、psm₂-psm′₂和psm₃-psm′₃，作为特征4、特征5和特征6；

计算待训练氨基酸a_i的位置信息、氨基酸类别信息(该类别信息用于指示氨基酸类型，总共有20种氨基酸，使用26个大写的英文字母表示，其中，去掉字母B、J、O、U、X、Z)和训练肽段的长度信息(即对于训练肽段a₁a₂…a_l，其肽段长度为l)，作为特征7、特征8和特征9。

步骤S140，训练分类模型以获得氨基酸可信度评估模型

利用获得的正、负样本通过机器学习的方法训练分类器，以获得训练好的分类模型，即氨基酸可信度评估模型，其中，正样本用于识别正确的氨基酸，负样本用于识别错误的氨基酸，训练时的输入是利用步骤S130的过程提取的正样本和负样本中氨基酸的九维特征向量，训练的输出是氨基酸为正确氨基酸的打分或为错误氨基酸的打分。

在此步骤中，分类器可采用支持向量机(SVM)或者诸如决策树，随机森林RF、贝叶斯网络等的其它类型。在一个实施例中，采用SVM训练分类时，使用SVM的径向基核函数或者用其它的核函数，甚至是非核函数的线性分类均可。

利用训练好的分类模型即可评估任意待测氨基酸的可信度，参见图3所示的评估氨基酸可信度的方法的流程图。该实施例以SVM分类模型为例进行介绍，具体包括：

步骤S310、对待评估原始肽段的谱图进行预处理。

此步骤的目的在于，在从头测序之前，从包含待评估氨基酸的原始肽段对应的谱图中去除大量同位素谱峰和噪音谱峰，以避免对从头测序算法造成干扰，例如，删除母离子附近的谱峰，以及删除中性失水，如失去水分子和氨分子的谱峰。

在一个实例中，对谱图进行预处理的过程包括：

在原始肽段对应的谱图中枚举电荷，根据两两谱峰质量差，查找所有同位素峰簇；根据同位素峰簇中两两峰的质量差，判断电荷；如果质量差约为则为+n电荷；根据电荷量，将单同位素峰转化为单电荷质量，去掉其它同位素峰；去掉谱图中的母离子峰和母离子失水、失氨峰。

步骤S320、产生待评估氨基酸的背景肽段。

根据与步骤S120类似的过程对于待测的原始肽段序列a₁a₂…a_l中的评估氨基酸a_i产生背景肽段。

步骤S330、提取并选择原始肽段和待评估氨基酸的特征

利用与步骤S130类似的过程提取并选择待评估氨基酸的特征，类似地，提取的特征包括下列项中的至少两项：1)原始肽段的肽谱匹配打分；2)原始肽段的谱峰强度匹配比例；3)原始肽段的谱峰数目匹配比例；4)原始肽段与打分最好的背景肽段的肽谱匹配分差；5)原始肽段与打分最好的背景肽段的谱峰强度匹配比例差；6)原始肽段与打分最好的背景肽段的谱峰数目匹配比例差；7)氨基酸位置信息(例如，从1到肽段长度l)；8)氨基酸类别信息；9)肽段长度信息。

步骤S340，利用训练好的氨基酸可信度评估模型获得待评估氨基酸的可信度打分分布。

在此步骤中，将提取的原始肽段和待评估氨基酸的特征输入至上述获得的氨基酸可信度评估模型，对于所有的待测氨基酸使用训练好的SVM模型进行打分并绘制打分分布，为便于后续描述，将对氨基酸可信度的打分，取名为SVM-Score。

通过此步骤的氨基酸的可信度打分分布，可以识别待评估的氨基酸是否为正确的氨基酸，例如，如果打分高于预定的阈值，则认为所评估的氨基酸的可信度较高。

步骤S350，使用Gamma拟合氨基酸的可信度打分分布。

此步骤的目的在于，为了更精确的获得氨基酸的可信度，可对获得的打分分布进一步处理。

由于SVM-Score分布与Gamma分布类似，在此实施例中，使用Gamma分布来拟合氨基酸的可信度打分分布。例如，可使用EM(期望最大化)方法结合Gamma分布来拟合。由于待评估的氨基酸的打分分布必然是两个分布(因为存在两个类别：正确氨基酸和错误氨基酸)，所以使用两个Gamma分布Γ(X|α_w,β_w)和Γ(X|α_r,β_r)分别拟合SVM-Score打分的正确，其中X表示SVM-Score，α_w,β_w表示错误结果分布的Gamma参数，α_r,β_r表示正确结果分布的Gamma参数。

图4示出了SVM-Score分布与Gamma分布的示意，横坐标表示分数值，纵坐标表示对应分数的氨基酸比例(Percentage)，其中“real data”表示包括所有正确、错误的氨基酸的实际打分分布，“real incorrect”和“real correct”分别表示错误氨基酸和正确氨基酸的实际打分分布，“estimated incorrect”和“estimated correct”分别表示使用EM算法估计出来的错误氨基酸和正确氨基酸Gamma分布。

步骤S360，计算待评估氨基酸的假发现率FAR：

其中，p_w表示错误氨基酸的先验概率，p_w×Γ(X|α_w,β_w)表示超过阈值的错误氨基酸数目，p_r表示正确氨基酸的先验概率，p_r×Γ(X|α_r,β_r)表示超过阈值的正确氨基酸数目，上述氨基酸的假发现率FAR(False>

根据本发明的另一方面，还可利用训练好的分类模型进行修饰位点定位的评估，参见图5所示，该实施例包括以下步骤：

步骤510，对于给定肽段序列枚举可以发生磷酸化修饰的候选修饰位点。

例如，对于给定的肽段序列WQSHTPPYAEK，该序列上发生了一个磷酸化修饰，假设磷酸化修饰可以发生在S、T、Y三个氨基酸上。对此三个候选修饰位点定位的具体过程是：

枚举WQSHTPPYAEK上所有可以发生磷酸化修饰的修饰位点：WQpSHTPPYAEK，WQSHpTPPYAEK和WQSHTPPpYAEK，其中“pS”表示S氨基酸上发生了磷酸化修饰，“pT”表示氨基酸T上发生了磷酸化修饰，“pY”表示氨基酸Y上发生的磷酸化修饰，发生磷酸化修饰的位点可统一采用“pX”来表示，即表示X氨基酸上发生了磷酸化修饰，在此实施例中X可以为S、T、Y三个氨基酸中的任意一个，可以将pX看作是一个新的氨基酸。

步骤S520，利用训练好的氨基酸可信度评估模型获得每个候选位点发生磷酸化修饰的可信度打分。

根据本发明提供的氨基酸可信度的评估方法计算pX这个新的氨基酸的可信度打分，在该实施例中，有三个候选修饰位点S、T、Y，则打分分别表示为s₁，s₂和s₃。

步骤S530，计算每个候选修饰位点发生磷酸化修饰的概率。

计算每个候选修饰位点发生磷酸化的概率。在此实施例中，使用贝叶斯公式来计算磷酸化发生在候选位点的概率，即：

其中，p_i表示第i个修饰位点的先验概率，s_i表示使用本发明的方法得到的候选的磷酸化位点i的可信度打分，t_i表示第i个位点是否发生磷酸化，如果t_i等于1，则表示发生磷酸化；t_i等于0，表示没发生磷酸化。

综上所述，本发明将氨基酸可信度评估与修饰位点定位两个问题进行了统一，将带修饰的氨基酸认为是新的氨基酸，因此，氨基酸可信度评估的方法也可应用于修饰点定位的评估。

本发明的方法可以实现为软件、硬件或软硬件结合的方式。为了进一步验证本发明的效果，发明人将本发明的方法实现为软件并将其将与目前仅有的两个支持氨基酸可信度评估的软件PEAKS和Novor进行了对比。结果表明，在三个真实数据集上，本发明的方法远好于现有的两个软件，例如，在FAR控制为5％情况下，本发明能够比性能较好的软件PEAKS多鉴定124.8％的氨基酸；在修饰位点定位方面，在三个磷酸化富集的数据集上，本发明的方法也优于非常流行的软件Ascore和phosphoRS，例如，在FAR控制为1％的情况下，本发明的方法能够比Ascore多鉴定67.5％的磷酸化位点，比phosphoRS多鉴定65.6％的磷酸化位点，同时覆盖了Ascore和phosphoRS的98％结果以及单独鉴定到自身的21％结果。

需要说明的是，虽然上文按照特定顺序描述了各个步骤，但是并不意味着必须按照上述特定顺序来执行各个步骤，实际上，这些步骤中的一些可以并发执行，甚至改变顺序，只要能够实现所需要的功能即可。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以包括但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。