CN103439441B - 一种基于子集错误率估计的肽鉴定方法 - Google Patents

Abstract

本发明涉及一种基于子集错误率估计的肽鉴定方法，其步骤包括：1）用质谱仪分析需要鉴定的肽样品以生成串联质谱；2）把所述串联质谱搜索包含目标肽序列的目标-诱饵蛋白质数据库，将得到的肽鉴定结果按照打分从高到低排序；3）给定打分阈值x，采用迁移FDR方法估计打分高于x的第k类肽鉴定子集的错误率FDR_k(x)；4）通过调整打分阈值x找到的x最小值，使得估计出的FDR_k(x)小于给定的错误率控制水平α，所得打分高于x的第k类肽鉴定结果即为可以接受的可靠鉴定结果。本发明通过迁移FDR方法估计子集错误率，进而通过该子集错误率得到肽的可靠鉴定结果，具有更高的鉴定准确率。

Description

一种基于子集错误率估计的肽鉴定方法

技术领域

本发明属于蛋白质分析技术领域，具体涉及一种基于子集错误率估计的肽鉴定方法。

背景技术

众所周知，绝大多数生物的遗传信息保存在DNA中。DNA通过转录过程生成信使RNA，而信使RNA又通过翻译过程生成蛋白质，从而实现了遗传信息由DNA到RNA再到蛋白质的传递，这一过程也被称为生命的中心法则。在从RNA翻译生成蛋白质的过程中，20种氨基酸以肽键顺序相连所形成的链状分子被称为肽，而其中分子量达到一定级别的肽则被称为蛋白质。大多数蛋白质在翻译形成后，会在蛋白质中的某些氨基酸上增加某种功能团(如在蛋白质的N末端加入乙酰)，或增加了其他的蛋白质或肽，或改变了氨基酸的化学性质或结构，这一过程被称为发生了化学修饰，由于该过程发生在前述的翻译过程后，因此在蛋白质氨基酸上所发生的变化也被称为蛋白质翻译后修饰。

液相色谱与质谱仪联用，并结合数据库搜索计算是目前蛋白质组学中鉴定蛋白质及其翻译后修饰的常用方法。在这种方法中，通过液相色谱与质谱仪联用可以得到蛋白质样品的实验串联质谱。实验串联质谱的获取过程包括：蛋白质样品首先被选定的蛋白酶水解，形成肽混合物；肽混合物通过液相色谱进行分离，不同物理化学性质的肽先后从色谱柱中流出；从色谱柱中流出的肽不断进入质谱仪；肽在质谱仪中被离子化，具有特定质量电荷比的肽离子在能量作用下碎裂形成碎片离子，碎片离子被分离和检测形成肽碎片离子谱；通过以上过程便得到蛋白质的实验串联质谱。在得到实验串联质谱后就可以从实验串联质谱中鉴定肽的氨基酸序列，进而鉴定蛋白质。

从实验串联质谱中鉴定肽的氨基酸序列时通常采用数据库搜索计算的方法。在计算过程中，数据库中所保存的蛋白质序列被模拟水解成候选肽，然后再将候选肽理论碎裂，生成理论串联质谱；将模拟计算得到的多个理论串联质谱依次与前述液相色谱与质谱仪联用所得到的实验串联质谱相比较，根据相似度进行打分，得分最高的理论质谱对应的肽就是实验质谱的鉴定结果。如果生成实验串联质谱的肽序列存在于数据库中的话，就可能将其鉴定出来。为了鉴定发生翻译后修饰的蛋白质，一种常见的基于串联质谱的鉴定方法是在数据库搜索时指定一些可变修饰类型，然后在生成候选肽时同时考虑发生和不发生指定修饰的情况，当候选肽中有多个可能的修饰位点时考虑所有可能的组合。

在基于质谱的蛋白质组学研究中，一次蛋白质组质谱实验通常能够产生数千至百万规模的串联质谱。通过数据库搜索鉴定这些质谱图，就产生了数目巨大的有待确认的肽鉴定结果。然而，由于谱图信号差、存在未知修饰、以及打分算法的缺陷等原因，这些结果的一部分(往往是大部分)是不正确的。所以，需要根据鉴定分值对鉴定结果进行过滤以及FDR(FalseDiscovery Rate，中文可翻译为假发现率或者错误发现率，参见参考文献：Benjamini,Y.and Y.Hochberg,Controlling the false discovery rate:a practical and powerful approach to multipletesting.Journal of the Royal Statistical Society,Series B(Methodological),1995.57(1):p.289-300.)的估计和控制。目前最常用和有效的肽鉴定FDR估计方法是目标-诱饵库搜索方法(参见文献：Elias,J.E.and S.P.Gygi,Target-decoy search strategy for increased confidence inlarge-scale protein identifications by mass spectrometry.Nat Methods,2007.4(3):p.207-14.)。在这种方法中，通过搜索诱饵蛋白质序列(如目标库的反序列构成的数据库)来获得错误的鉴定，而FDR就用分值阈值之上的诱饵肽鉴定数量除以目标肽鉴定数量来估计。当鉴定数量较大时，这种目标-诱饵库搜索方法可以有效的估计肽鉴定FDR。但是如果鉴定数量较少的话，这种FDR估计方法就不准确了(参见文献：Huttlin,E.L.,et al.,Prediction of error associated withfalse-positive rate determination for peptide identification in large-scale proteomics experimentsusing a combined reverse and forward peptide sequence database strategy.J Proteome Res,2007.6(1):p.392-8.)。

在当前的肽和蛋白质鉴定研究中，各种形式的肽(例如带有各种修饰的肽和非修饰的肽、特异性酶切肽和非特异性酶切肽等等)通常都是混合在一起鉴定分析的，从最初的蛋白质样品处理到最后的FDR估计都是如此。但最后经常是只强调或报道鉴定结果的某个子集(例如带有某种修饰的肽，或者特异性酶切肽，或者含有某种特定氨基酸的肽等等)，造成估计的FDR可能严重高于或低于真实的子集FDR。这种混合鉴定分析有多方面的原因。比如，质谱仪和鉴定软件通常都无法有效识别不同种类的肽。并且，对所有鉴定结果进行整体的FDR估计也是容易计算的，并且这也是现有蛋白质组学软件所提供的功能。各种肽的混合鉴定，降低了数据分析的复杂度，但却使不同种类的肽的特性遭到掩盖。这种“不分青红皂白”的混合可能造成非常严重的后果，甚至错误的结论。但这一点尚未引起领域的重视，也没有清晰的认识。一种解决办法是将感兴趣的肽鉴定子集分离出来，单独用目标-诱饵库方法估计FDR。但是如果子集较小的话，这种分离估计FDR的做法就不准确了。子集越小，误差越大。极端情况下，如果子集只包含一个鉴定，那么根本无法用目标-诱饵库方法来估计该子集的FDR，从而无法实现肽的可靠鉴定。

发明内容

本发明的目的是提供一种基于子集错误率估计的肽鉴定方法，通过优化的子集错误率估计得到肽的可靠鉴定结果。

本发明采用的技术方案如下：

一种基于子集错误率估计的肽鉴定方法，其步骤包括：

1)用质谱仪分析需要鉴定的肽样品以生成串联质谱；

2)把所述串联质谱搜索包含目标肽序列的目标-诱饵蛋白质数据库，将得到的肽鉴定结果按照打分从高到低排序；

3)给定打分阈值x，采用迁移FDR方法估计打分高于x的第k类肽鉴定子集的错误率FDR_k(x)；

4)通过调整打分阈值x找到x的最小值，使得估计出的FDR_k(x)小于给定的错误率控制水平α，所得打分高于x的第k类肽鉴定结果即为可靠鉴定结果。

进一步地，采用所述迁移FDR方法得到的FDR_k(x)的计算公式为：

${\mathrm{FDR}}_k\left(x\right)=\frac{\mathrm{FDR}\left(x\right)}{\mathrm{FDR}\left(x\right)+\frac{{\mathrm{\λ}}_k\left(x\right)}{{\mathrm{\γ}}_k\left(x\right)}(1-\mathrm{FDR}\left(x\right))},$

其中，FDR(x)为打分在阈值x之上的所有目标鉴定的整体FDR，λ_k(x)为在鉴定正确并且打分大于x的条件下，一张谱图被鉴定为第k类肽的概率；γ_k(x)为在鉴定错误并且打分大于x的条件下，一张谱图被鉴定为第k类肽的概率。

进一步地，所述λ_k(x)的估计方法为：

其中，N(x)是打分大于x的目标鉴定的数量，N_k(x)是打分大于x且为第k类肽目标鉴定的数量；进而得到FDR_k(x)估计的一个中间形式：

${\mathrm{FDR}}_k\left(x\right)=\frac{N\left(x\right)}{N_k\left(x\right)}{\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right).$

进一步地，采用数据驱动的方法，将目标-诱饵库搜索结果中所有诱饵序列匹配作为训练数据，用于估计所述γ_k(x)。将所述γ_k(x)用x的线性函数来近似，即

给定一次数据库搜索的结果和任意的x值，计算在所有打分大于x的诱饵鉴定中第k类肽鉴定所占的比例，x值和该比例值构成线性函数的一个观测样本；通过改变x的值生成一个训练样本数据集，利用该训练样本数据集，使用最小二乘回归的方法估计线性函数中的系数a和b，并分别用和表示，则最终FDR_k(x)的估计为

进一步地，采用下式估计所述FDR(x)：

其中D(x)是打分在x之上的所有匹配诱饵序列的肽鉴定的数目；进而得到FDR_k(x)的最终估计结果：

以上过程可以对不同的k进行计算，从而得到全面的鉴定结果。

本发明的基于子集错误率估计的肽鉴定方法，提出了一种新的、优化的方式来估计子集错误率，进而通过该子集错误率可得到肽的可靠鉴定结果，比现有技术中的肽鉴定方法具有更高的准确率。当子集很小时，传统方法极不准确，甚至无法应用。本发明创造性的采用从整体FDR推算子集FDR的方法，避免了直接完全从数据估计子集FDR的传统方法，对于规模较小的子集，具有极大的优势。

附图说明

图1是本实施例的基于子集错误率估计的肽鉴定方法的步骤流程图；

图2为Mascot软件产生的诱饵修饰肽比例γ_k(x)与x之间的近似线性关系的示意图；

图3为SEQUEST软件产生的诱饵修饰肽比例γ_k(x)与x之间的近似线性关系的示意图；

图4为pFind软件产生的诱饵修饰肽比例γ_k(x)与x之间的近似线性关系的示意图。

具体实施方式

下面通过具体实施例和附图，对本发明做详细的说明。

图1是本实施例的基于子集错误率估计的肽鉴定方法的步骤流程图。如该图所示，首先用质谱仪分析需要鉴定的肽样品以生成串联质谱；然后把所述串联质谱搜索包含目标肽序列的目标-诱饵蛋白质数据库，并将得到的肽鉴定结果按照打分从高到低排序；然后给定打分阈值x，采用迁移FDR方法估计打分高于x的第k类肽鉴定子集的错误率FDR_k(x)；再通过调整打分阈值x找到的x最小值，使得估计出的FDR_k(x)小于给定的错误率控制水平α，所得打分高于x的第k类肽鉴定结果即为可以接受的可靠鉴定结果。下面更详细地说明该方法的各个步骤。

假设我们将一组串联质谱图搜索一个目标-诱饵蛋白质序列库。数据库搜索完成之后，每一个谱图匹配上了一个(得分最高的)肽段，并有相应的打分值(不妨假设匹配越好分值越高)，该匹配称作一个肽鉴定。我们只对某类(用“第k类”表示，比如带有某种特定修饰、或者含有某种特定的氨基酸等等)肽鉴定兴趣，并希望根据打分阈值估计和控制这些修饰鉴定的FDR。用于问题建模的符号表示如表1所示。

表1.符号表示

注意，表1中定义的FDR叫做贝叶斯FDR，当数据规模巨大并且独立同分布时，贝叶斯FDR与原始定义的FDR是等价的(参见文献：Benjamini,Y.and Y.Hochberg,Controlling thefalse discovery rate:a practical and powerful approach to multiple testing.Journal of the RoyalStatistical Society,Series B(Methodological),1995.57(1):p.289-300.)。

使用贝叶斯规则，可以得到FDR_k(x)与FDR(x)之间的一个关系式：

$\begin{array}{c}{\mathrm{FDR}}_k\left(x\right)\\ =P\left(F\right|I_k,X>x)\\ =\frac{P(F,I_k|X>x)}{P\left(I_k\right|X>x)}\\ =\frac{P\left(I_k\right|F,X>x)P\left(F\right|X>x)}{P\left(I_k\right|F,X>x)P\left(F\right|X>x)+P\left(I_k\right|T,X>x)P\left(T\right|X>x).}\\ =\frac{P\left(F\right|X>x)}{P\left(F\right|X>x)+\frac{P\left(I_k\right|T,X>x)}{P\left(I_k\right|F,X>x)}(1-P\left(F\right|X>x))}\\ =\frac{\mathrm{FDR}\left(x\right)}{\mathrm{FDR}\left(x\right)+\frac{{\mathrm{\λ}}_k\left(x\right)}{{\mathrm{\γ}}_k\left(x\right)}(1-\mathrm{FDR}\left(x\right))}\end{array}---\left(1\right)$

在实际中，所有鉴定结果的全集通常足够大，所以可以用目标-诱饵库方法准确地估计FDR(x)。然而，由第k类肽构成的鉴定子集有可能非常小，以至于无法准确估计其FDR。幸运的是，公式(1)意味着，如果能够计算λ_k(x)和γ_k(x)，那么就可以从FDR(x)间接地推算出FDR_k(x)，而不必使用不充足的数据直接估计FDR_k(x)。这里，从数据直接估计FDR_k(x)的一种计算方法为N_k(x)/D_k(x)，当D_k(x)很小时，这个估计显然会有很大的偏差。

可以有多种方法估计λ_k(x)和γ_k(x)。

λ_k(x)是在某鉴定为正确且打分大于x的前提下，该鉴定为第k类肽鉴定的概率。λ_k(x)的值跟谱图数据中由第k类肽产生的谱图所占的比例直接相关。这个比例可以根据经验设定，比如，如果认为某种修饰含量非常低，可以设定一个很小的值(例如0.0001)给λ_k(x)。

γ_k(x)是在鉴定错误且打分大于x的前提下，一个谱图被鉴定为第k类肽的概率。这个概率跟候选肽搜索空间中第k类肽所占的比例高度相关，并且从某种意义上说，就应该等于这个比例。所以，一种估计γ_k(x)的方法就是：将数据库中的蛋白质模拟水解，生成所有候选肽，计算其中第k类肽所占的比例，作为γ_k(x)的值。

但是以上方法都有一些缺点，实验表明下面给出的估计λ_k(x)和γ_k(x)的方法更为有效。

1.估计λ_k(x)

估计方法为：

其中，N(x)是打分大于x的目标鉴定的数量，N_k(x)是打分大于x且为第k类肽目标鉴定的数量。用公式(2)替换公式(1)中的λ_k(x)，得到

$\begin{array}{c}{\mathrm{FDR}}_k\left(x\right)\\ =\frac{\mathrm{FDR}\left(x\right)}{\mathrm{FDR}\left(x\right)+\frac{\frac{N_k\left(x\right)(1-{\mathrm{FDR}}_k\left(x\right))}{N\left(x\right)(1-\mathrm{FDR}\left(x\right))}}{{\mathrm{\γ}}_k\left(x\right)}(1-\mathrm{FDR}\left(x\right))}\\ =\frac{{\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right)}{{\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right)+\frac{N_k\left(x\right)}{N\left(x\right)}(1-\mathrm{FD}{R}_k\left(x\right))}\end{array}---\left(3\right)$

进一步有,

$\begin{array}{c}{\mathrm{FDR}}_k\left(x\right)({\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right)+\frac{N_k\left(x\right)}{N\left(x\right)}(1-{\mathrm{FDR}}_k\left(x\right)))\\ ={\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right)\end{array}---\left(4\right)$

或者

$\begin{array}{c}\frac{N_k\left(x\right)}{N\left(x\right)}{\mathrm{FDR}}_k\left(x\right)(1-{\mathrm{FDR}}_k\left(x\right))\\ ={\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right)(1-{\mathrm{FDR}}_k\left(x\right))\end{array}---\left(5\right)$

从公式(5)便得到FDR_k(x)估计的一个中间形式:

${\mathrm{FDR}}_k\left(x\right)=\frac{N\left(x\right)}{N_k\left(x\right)}{\mathrm{\γ}}_k\left(x\right)\mathrm{FDR}\left(x\right).---\left(6\right)$

这样，公式右边就剩下γ_k(x)是未知的。

2.估计γ_k(x)

如前所述，γ_k(x)跟候选肽搜索空间中第k类肽所占的比例高度相关。但是候选肽搜索空间是不可见的，并且不同的搜索引擎可能有不同枚举候选肽的方式。并且候选肽比例是常数，但γ_k(x)可能与x的值相关。这里用数据驱动的方法来估计γ_k(x)。在目标-诱饵库搜索方法中，所有诱饵序列匹配都是错误的鉴定。这些鉴定构成了天然的训练数据，可用于估计γ_k(x)。根据对实际数据的观察，γ_k(x)可以用x的线性函数来近似，即

给定一次数据库搜索的结果和任意的x值，很容易计算在所有打分大于x的诱饵鉴定中第k类肽鉴定所占的比例。x的值和计算所得的比例值就构成了线性函数的一个观测样本。通过改变x的值，就可以生成一个训练样本数据集：{(x_i,y_i|i＝1,2,…,n,y_i＝D_k(x_i)/D(x_i))}。利用这个训练集，使用最小二乘回归的方法，就可以估计线性函数中的两个系数(用和表示)了。图2至图4分别为Mascot软件、SEQUEST软件和pFind软件产生的诱饵修饰肽比例γ_k(x)与x之间的近似线性关系的示意图，说明了公式(7)的合理性。

最终，FDR_k(x)的估计为

该估计称作第k类肽鉴定的迁移FDR(Transferred FDR)，其含义是指该FDR是从整体FDR转换过来的，而不是直接从数据中估计出来的。

可以用任何有效的方法估计FDR(x)，最简单的估计就是：

所以，FDR_k(x)的一个简单估计就是：

如果打分阈值x很大的话，第k类肽的诱饵鉴定必然很少，计算出的比例值会不稳定，发生波动。这种波动就是造成直接FDR估计不准确的原因。这就是为什么不能直接用打分在x之上的诱饵鉴定来估计FDR_k(x)的原因。实际上，最好不要用大的x值来拟合公式(7)中的函数。对于大的x值，γ_k(x)的值应该用拟合的函数外推出来。此外，a和b的值不仅取决于质谱数据集，也取决于搜索的数据库和搜索参数，所以，每次搜索都应该重新估计a和b的值。

为了直观的理解迁移FDR的效果和计算过程，此处给出一个例子。假设在一次实验中，目的是鉴定磷酸化肽段，所以指定磷酸化为可变修饰。首先用质谱仪分析需要鉴定的肽样品以生成串联质谱，然后把串联质谱搜索包含目标肽序列的目标-诱饵蛋白质数据库，并将得到的肽鉴定结果按照打分从高到低排序。通过搜索目标-诱饵库搜索后，打分在阈值x＝37之上的目标鉴定数量为N(x)＝3249(其中磷酸化鉴定数量N_k(x)＝44个)和3个诱饵鉴定(全部为非修饰鉴定)。那么，整体FDR可估计为如果分开估计磷酸化鉴定的话，其FDR为0/44＝0。通过用所有诱饵鉴定拟合公式(7)，得到因此，根据公式(8)，迁移FDR为通过检查3249个目标鉴定，发现其中5个是错误鉴定，包括1个磷酸化鉴定。因此，磷酸化鉴定的实际错误率为1/44≈0.0227，跟迁移FDR估计值(0.0222)非常接近，但跟估计的整体FDR值(0.00092)相差甚远。

在可以计算FDR_k(x)之后，就可以根据指定的FDR控制水平，过滤出符合条件的第k类肽鉴定鉴定结果。假设指定的FDR控制水平为α(常用α值为0.01或0.05)，并且打分越高鉴定越可靠。那么，先设定打分阈值x为一个足够低的值，将有所有打分高于x的鉴定结果提取出来，用上述方法计算FDR_k(x)，如果计算出的FDR_k(x)值高于α，则提高阈值x的值(比如加一个足够小的步长σ，即x＝x+σ)，然后重新计算FDR_k(x)。重复上述过程，直到FDR_k(x)<α。这时，所有打分高于x的第k类肽鉴定鉴定结果即是满足错误率(FDR)小于α的肽鉴定结果。

需注意的是，以上整个过程可以对不同的k进行计算，再将结果合并，从而得到所有想要的鉴定结果。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求所述为准。

Claims (6)

1.一种基于子集错误率估计的肽鉴定方法，其步骤包括：

1)用质谱仪分析需要鉴定的肽样品以生成串联质谱；

2)把所述串联质谱搜索包含目标肽序列的目标-诱饵蛋白质数据库，将得到的肽鉴定结果按照打分从高到低排序；

3)给定打分阈值x，采用迁移FDR方法估计打分高于x的第k类肽鉴定子集的错误率FDR_k(x)；采用所述迁移FDR方法得到的FDR_k(x)的计算公式为：

${\mathrm{FDR}}_k\left(x\right)=\frac{\mathrm{FDR}\left(x\right)}{\mathrm{FDR}\left(x\right)+\frac{{\mathrm{\&lambda;}}_k\left(x\right)}{{\mathrm{\&gamma;}}_k\left(x\right)}(1-\mathrm{FDR}\left(x\right))},$

其中，FDR(x)为打分在阈值x之上的所有目标鉴定的整体FDR，λ_k(x)为在鉴定正确并且打分大于x的条件下，一张谱图被鉴定为第k类肽的概率；γ_k(x)为在鉴定错误并且打分大于x的条件下，一张谱图被鉴定为第k类肽的概率；

4)通过调整打分阈值x找到x的最小值，使得估计出的FDR_k(x)小于给定的错误率控制水平α，所得打分高于x的第k类肽鉴定结果即为可靠鉴定结果。

2.如权利要求1所述的方法，其特征在于：所述λ_k(x)的估计方法为：

其中，N(x)是打分大于x的目标鉴定的数量，N_k(x)是打分大于x且为第k类肽目标鉴定的数量；进而得到FDR_k(x)估计的一个中间形式：

${\mathrm{FDR}}_k\left(x\right)=\frac{N\left(x\right)}{N_k\left(x\right)}{\mathrm{\&gamma;}}_k\left(x\right)\mathrm{FDR}\left(x\right).$

3.如权利要求2所述的方法，其特征在于：采用数据驱动的方法，将目标-诱饵库搜索结果中所有诱饵肽鉴定作为训练数据，用于估计所述γ_k(x)。

4.如权利要求3所述的方法，其特征在于：将所述γ_k(x)用x的线性函数来近似，即

给定一次数据库搜索的结果和任意的x值，计算在所有打分大于x的诱饵鉴定中第k类肽鉴定所占的比例，x值和该比例值构成线性函数的一个观测样本；通过改变x的值生成一个训练样本数据集，利用该训练样本数据集，使用最小二乘回归的方法估计线性函数中的系数a和b，并分别用和表示，则最终FDR_k(x)的估计为

5.如权利要求4所述的方法，其特征在于：采用下式估计所述FDR(x)：

其中D(x)是打分在x之上的所有匹配诱饵序列的肽鉴定的数目；进而得到FDR_k(x)的最终估计结果：

6.如权利要求1或2所述的方法，其特征在于，估计所述γ_k(x)的方法是：将数据库中的蛋白质模拟水解以生成所有候选肽，计算其中第k类肽所占的比例，作为γ_k(x)的值。