CN102495127B - 一种基于概率统计模型的蛋白质二级质谱鉴定方法 - Google Patents

Abstract

本发明公开了一种基于概率统计模型的蛋白质二级质谱鉴定方法，该方法首先虚拟酶解蛋白质数据库序列，并根据肽段的质量数对酶解后的肽段建立肽段数据库和肽段数据库索引；然后根据待分析实验图谱中母离子的核质比在肽段数据库中找出符合要求的候选肽段，并对找到的所有候选肽段产生符合要求的理论图谱；然后对待分析实验图谱进行去同位素和去噪处理；对处理后的待分析实验图谱和每张候选肽段的理论图谱进行匹配打分，选择分值最高的候选肽段作为此实验图谱的鉴定结果；最后针对所有实验图谱鉴定结果进行整体假阳性控制。该方法鉴定有效质谱的数量和蛋白质肽段数量均高于目前现有算法，且可动态选峰，运行速度快。

Description

一种基于概率统计模型的蛋白质二级质谱鉴定方法

技术领域

本发明涉及蛋白质二级质谱鉴定领域，特别涉及一种基于概率统计模型的蛋白质二级质谱鉴定方法。

背景技术

随着基质辅助激光解吸(matrix-assisted laser desorption ionization，MALDI)和电喷雾(Electrospray Ionization，ESI)两种软电离技术的出现，使生物质谱能够较少地引入杂质和保持肽段分子的完整性，从而使生物质谱大规模应用于蛋白质分析。目前，生物质谱已成为蛋白质组研究的支撑技术之一，其主要是利用串联质谱(LC MS/MS)来分析蛋白质样品。在蛋白质组的生物信息学研究中，质谱数据处理是十分重要的研究内容，其任务是从带有复杂噪声或者部分信息缺失的数据中推断样品的蛋白质组成。数据库搜索是质谱数据处理的主要方法，其基本过程如图1所示，即将实验图谱和数据库中的理论酶切图谱进行比对、打分，选择分值最高的匹配作为搜索结果后选肽段。

蛋白质二级质谱鉴定涉及到诸多方面的内容，其主要涉及到母离子价态的确定、效质谱峰的选取和匹配打分模型。目前针对鉴定结果整体质量控制的方法主要是应用随机数据库方法对整体鉴定结果进行阳性率控制，其基本思想是：先针对真实蛋白质数据库和实验数据集构建一个随机数据库，然后同时或者分别搜索真实蛋白质数据库和新构建的随机数据库，通过随机数据库肽段匹配来模拟正常数据库中的随机匹配，从而估计正常数据库中随机匹配的特征分布，确定不同过滤标准，Kall’s于2008年在Proteome上公开了一种方法，具体是采用如下公式来得到整体数据集的假阳性率(False Discovery Rate，FDR)。

$\mathrm{FDR}=\frac{N_R}{N_N}$

目前蛋白质二级质谱鉴定算法根据匹配打分模型大致可以分为两类：解释型模型和概率统计模型。其中著名的商业软件SEQUEST的算法是解释型模型，而另一个商业软件Mascot的算法是概率统计模型。另外还有一些免费的鉴定算法，例如比较有影响力的基于统计模型的算法有X！Tandem和OMSSA。其中X！Tandem用的是超几何模型，OMSSA用的是泊松分布模型。这些基于统计模型的算法中考虑的是实验质谱峰匹配与不匹配，并没有考虑质谱峰的连续匹配情况，更较少考虑到质谱峰峰强的概率模型。在基于解释模型的算法中，其中Sequest考虑了离子连续匹配和峰强。但它统一把峰强分别定义为三个值：50(b和y离子)、25(b，y离子脱水和脱氨离子)和10(a离子)，没有充分体现实验质谱离子的特征。

因此，研究一种能大大提高蛋白质有效质谱和蛋白质肽段数量的二级质谱鉴定方法具有很高的理论和实际应用价值。

发明内容

本发明的主要目的在于克服现有技术的缺点与不足，提供一种基于概率统计模型的蛋白质二级质谱鉴定方法，该方法鉴定有效质谱和蛋白质肽段的数量均高于现有算法。

本发明的目的通过以下的技术方案实现：一种基于概率统计模型的蛋白质二级质谱鉴定方法，具体包括以下步骤：

(1)虚拟酶解蛋白质数据库序列，并根据肽段的质量数对酶解后的肽段建立肽段数据库和肽段数据库索引；

(2)根据待分析实验图谱中母离子的核质比在步骤(1)所述的肽段数据库中找出符合要求的候选肽段，并对找到的所有候选肽段产生符合要求的理论图谱；

(3)对待分析实验图谱进行去同位素和去噪处理；

(4)将步骤(3)得到的待分析实验图谱和步骤(2)中得到的每张候选肽段的理论图谱进行匹配打分，选择分值最高的候选肽段作为此实验图谱的鉴定结果；

(5)针对所有实验图谱鉴定结果进行整体假阳性控制。

所述步骤(1)具体包括以下步骤：

(1-1)读取待分析二级质谱样本中物种蛋白质序列库文件的一条蛋白质序列；

(1-2)根据用户设定的蛋白酶，找到蛋白质序列中的酶切位点，在符合规则的酶切位点产生断裂，从而产生没有漏切位点的肽段或存在漏切位点的断裂肽段；

(1-3)计算步骤(1-2)所得到的各个虚拟酶切后肽段的质量数，根据每个氨基酸的分子量计算每个肽段的质量数；

(1-4)将肽段信息写入肽段数据库中以该肽段取整后质量数命名的文件中；

(1-5)读取下一条蛋白质序列，重复步骤(1-2)-(1-4)，直到所有的蛋白序列被酶解和存入肽段数据库；

(1-6)按文件名的数字从小到大读出文件中的肽段信息，每读一个文件，按照文件中肽段的质量数从小到大进行排序，然后存入database.ind文件中，同时，以1da为单位对所有肽段建立查找索引database.index，其查找索引包括以下信息：其质量数，这些肽段在database.ind文件中的开始位置，该区间内的肽段的个数。

所述步骤(1-3)中在计算质量数之前首先对每个氨基酸的质量建立索引，其对20个氨基酸的索引和翻译后修饰的索引方法如下：

(1-3-1)启用一个与ASCII码等同大小的数组，该数组的下标与氨基酸单字母简写的ASCII码数值一致，其数组中保存氨基酸的质量数；

(1-3-2)把单字母表示氨基酸的肽段序列中每个字母依次转换成其对应ASCII码的数值，然后根据氨基酸索引表的数值计算每条虚拟酶解后的肽段的质量数。

所述步骤(2)在肽段数据库中找出符合要求的候选肽段的具体步骤是：

(2-1-1)加载步骤(1-6)中的database.index文件信息到内存数组index，读取待分析二级质谱的母离子核质比值和电荷信息，并计算其母离子去电荷后的质量数；

(2-1-2)根据容许的质量误差和步骤(2-1)所述的质量数在index数组中查找相应肽段在文件database.ind中的开始位置和行数，然后加载此区间内的所有肽段信息；

(2-1-3)根据用户所采用质谱仪的精确度，对步骤(2-1-2)加载到内存的肽段进行进一步的筛选，作为此待分析二级质谱的候选肽段。

所述步骤(2)中产生符合要求的理论图谱具体包括以下步骤：

(2-2-1)在肽段的上部和下部分别表示出b和y离子的类型，其中b离子是从左边往右标记，y离子是从右往左进行标记，模拟经过离子碎裂过程候选肽段可能产生的理论碎片b、y离子；

(2-2-2)在步骤(2-2-1)产生b、y离子基础上，根据下面两条的规则产生b、y离子丢水和丢氨情况下的碎片离子：

(2-2-2-1)如果在b离子或y离子中含有S，T，E和D四种氨基酸中的一种，那么产生该离子对应的丢水碎片离子b-H₂O和y-H₂O；

(2-2-2-2)如果在b离子或y离子中含有R，K，Q和N四种氨基酸中的一种，那么产生该离子对应的丢氨碎片离子b-NH₃和y-NH₃；

(2-2-3)产生二价的碎片离子，其规则为：待分析二级质谱母离子价态是一价的则考虑一价的碎片离子峰，当母离子价态＞＝2时，并且对应的碎片离子中包含R、K和H三种氨基酸其中一种时，则考虑二价碎片离子峰，在理论图谱中产生其对应的二价的碎片离子峰。

所述步骤(3)对待分析实验图谱进行去同位素的具体步骤是：

(3-1-1)初始化，将三个比较峰的mz值和其强度，全部设为0，设三个峰mz值是：mz_1＝0，mz_2＝0，mz_3＝0，其峰强对应是mz_1_in＝0，mz_2_in＝0和mz_3_in＝0，并设置用于存储非同位峰的保留峰容器，并已知测量质量误差m；

(3-1-2)读取一个峰的信息，将其作为第三个峰的值，即mz_3和mz_3_in；然后判断第三个峰是否是前两个峰的同位素峰，

(3-1-2-1)如果以下三个条件任意一个条件成立，则认为是同位素峰，

a.|mz_3-mz_2-1|＜＝m并且mz_2_in＞mz_3_in；

b.|mz_3-mz_1-1|＜＝m并且mz_1_in＞mz_3_in；

c.|mz_2-mz_1|＜＝m并且mz_2_in＞mz_3_in，此为相同的峰信息，记录误差，

然后进入步骤(3-1-2-3)；

(3-1-2-2)如果步骤(3-1-2-1)中各条件均不成立，则认为目前进入第三位置的峰不是同位素峰，将其作为保留峰存入保留峰容器中，然后进入步骤(3-1-2-3)；

(3-1-2-3)三个峰向前平移一位，空出第三个峰的位置，即将第一个峰的信息去掉，第三个和第二个峰的信息分别作为新的第二个和第一个峰的信息；

(3-1-3)逐个读取下一个峰的信息，重复步骤(3-1-2)，直到处理完一张二级质谱图所有峰信息，最后得到的保留峰容器中的峰即为去同位素峰后的非同位素峰。

所述步骤(3)，对待分析实验图谱进行去噪处理，即除去信号峰中的噪声峰，具体步骤是：

(3-2-1)首先选取局部最强峰，包括以下步骤：

(3-2-1-1)根据步骤(3-1-3)得到的去同位素后的离子峰，找到全局最强峰，然后以此峰为中心，分别向左右各平移50Da，形成一个搜索窗口，在这100Da范围内挑选离子峰强度排名前n位的峰，然后记录这n个峰的信息；

(3-2-1-2)以已搜索区域为中心，再分别向左右各平移50Da，在左右各形成1个搜索区域，在这100Da范围内挑选离子峰强度排名前n位的峰，然后记录这n个峰的信息；

(3-2-1-3)重复进行(3-2-1-1)和(3-2-1-2)两步，直到该质谱文件所有的质荷比信息被提取完成；

(3-2-2)根据步骤(3-2-1-1)得到的全局最强峰，搜索峰值大于等于全局最强峰峰值*0.33的峰，作为全局相对高峰，判断这些峰是否已记录在步骤(3-2-1)中，是则不做处理，否则记录峰的信息；

(3-2-3)将选取的局部最强峰和全局相对高峰进行合并，得到最终选取的用于鉴定的峰。

所述步骤(4)待分析实验图谱和理论图谱进行打分的具体步骤如下；

(4-1)待分析实验图谱和理论图谱进行匹配打分的具体步骤是：

(4-1-1)逐个读取峰信息判断理论图谱和选峰后的实验图谱是否匹配，如果理论图谱与实验图谱对应峰的核质比之差小于等于质谱仪的测量误差，则认为这两个峰匹配，之后记录其匹配的信息；

(4-1-2)设E为产生的理论碎片的个数；K为理论图谱和选峰后的实验图谱匹配个数，Q代表随机匹配概率事件，i为随机匹配概率，i＝0.01*n，P为在E个理论峰中有K个峰匹配的概率，则P由下面二项式分布概率密度函数计算：

$\left\{\begin{array}{c}Q=i+\mathrm{factor}\\ P=\left(\begin{array}{c}K\\ E\end{array}\right)Q^K{(1-Q)}^{E-K}\end{array}\right.;$

其中：

(4-2)待分析实验图谱和理论图谱进行连续匹配打分的具体步骤是：

设E1为理论图谱产生的理论连续匹配个数；K1为实验图谱实际连续匹配的个数；B_factor为背景值，B_factor＝统计大量实验图谱连续匹配的平均值/统计大量对应理论图谱连续匹配的平均值，Q1反映了某一图谱在步骤(4-1)匹配情况下连续匹配的概率，P1是在E1个理论连续匹配个数中实际存有K1个连续匹配的概率，由下面二项式分布概率密度函数计算：

$\left\{\begin{array}{c}Q1=B\_\mathrm{factor}*K/E\\ P1=\left(\begin{array}{c}K1\\ E1\end{array}\right){Q1}^{K1}{(1-Q1)}^{E1-K1}\end{array}\right.;$

所述待分析实验图谱和理论图谱连续匹配个数具体是指图谱中两两连续匹配的对数；

(4-3)对匹配峰强度信息进行分析，求得强度因子，具体步骤是：

设M_I为统计所有实验图谱中某两个氨基酸产生的峰大于等于最强峰的33％的个数，设M_E为期望总的离子的个数，则两个氨基酸中间的断裂概率Yi通过下式得到：

Yi＝M_I/M_E；

进而得到强度因子Infactor为：

Infactor＝(1+Ym+Bm))/(1+0.155*m)；

其中Ym＝∑Y_i；Bm＝∑B_i；Ym、Bm分别为实验图谱强度大于全局最强峰的33％的匹配峰Yi和Bi分值之和；m_p为一张实验图谱中强度大于最强峰的33％的匹配的个数；其中0.155是理论平均匹配值；

(4-4)有机组合上述步骤(4-1)，(4-2)和(4-3)的打分方法，采用下面公式得到肽段的得分：

PEP-S＝Infactor*(-10)*log₁₀ ^(P*P1)；

(4-5)对计算的PEP-S分数去除背景值，首先设在真实库和随机库统计概率相等时的背景值为其在某种情况下的背景值B_B，背景值B_B是经过贝叶斯网络学习得到的；

然后得到去背景值PEP-S_M：

PEP-S_M＝PEP-S-B_B；

(4-6)取出下一个肽段，重复执行步骤(4-1)-(4-5)，直到符合此图谱母离子误差的所有肽段均被打分处理；

(4-7)对此图谱所有候选肽段的得分PEP-S_M进行排序，值最大的作为当前图谱的鉴定结果。

所述步骤(5)针对所有实验图谱鉴定结果进行整体假阳性控制，具体包括以下步骤：

(5-1)统计待分析实验图谱所有二级图谱中的鉴定结果，取出肽段得分最小值和最大值；

(5-2)按分数统计待分析图谱的鉴定结果，统计最小值和最大值之间大于各得分中的正库肽段和随机库肽段的个数，设大于等于某一个分值时，其鉴定结果在真实库中的个数为N_N，在随机数据库中的个数为N_R，按照下述公式计算每个分值为阀值时FDR的值，

$\mathrm{FDR}=\frac{N_R}{N_N}$

(5-3)将产生的FDR＜＝0.01的值最小的分值作为待分析实验图谱的整体阀值；

(5-4)以步骤(5-3)得到的整体阀值过滤待分析实验图谱的鉴定结果，如果鉴定分数小于此阀值则被过滤掉，其结果作为待分析实验图谱最终鉴定结果。

所述步骤(3-2-1)和步骤(4-1)中的n取值范围为3～6。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明主要对生物质谱产生的二级质谱数据进行解释和鉴定，其鉴定有效质谱的数量和蛋白质肽段数量均高于目前的常用的国外商业软件的算法。目前现有技术中鉴定的有效质谱的数量和蛋白质肽段数量按从小到大顺序为：sequest，omissa(NCBI肽段开发)，X！tandom，ProteinPilot，Mascot。其中Mascot鉴定最多，本鉴定方法结果优于Mascot。

2、本发明方法在选取有效质谱峰与以前方法有了很大不同。此算法采取在全局最高峰的左右各50da动态的选取离子峰，并保留大于全局最高峰33％的峰，此选峰方法即保留了局部峰的信息又保留全局峰的信息，并且还加入了一个动态选峰机制。

3、本发明方法的打分模型主体是基于二项式分布统计模型，但加入了一些别的统计元素的全新打分模型。其方法与前人的方法不同，前人的统计方法只考虑了峰的匹配和不匹配因素。本方法不仅考虑了离子匹配和不匹配，还考虑了离子的连续匹配情况、离子峰的强度信息。在考虑强度信息时，本方法是基于实验质谱统计分析来给离子峰的强度一个统计值，完全有区别以前的人为设定值的方法。

4、本发明方法在实现过程中建立了多重索引技术方法，极大的提高了软件运行的速度。其中氨基酸和翻译后修饰索引的方法可以同时设置230钟翻译后修饰，与目前鉴定软件仅能设置10种翻译后修饰不同。

附图说明

图1是现有技术中二级质谱鉴定的基本流程图；

图2是本发明方法的主流程图；

图3是本发明方法中蛋白质虚拟酶解示意图；

图4是原始4个峰的去同位素执行过程中三个峰和保留峰的状态改变过程；

图5是本发明实施例中选取局部最强峰时的示意图；

图6是本发明实施例中选取全局相对高峰时的示意图；

图7是经过图5和图6所示选峰后得到的结果图；

图8是候选肽段可能产生的理论碎片b、y离子原理示意图；

图9是待分析实验图谱和理论图谱的匹配标注图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

如图1所示，一种基于概率统计模型的蛋白质二级质谱鉴定方法，包括以下步骤：

(1)虚拟酶解蛋白质数据库序列，并根据肽段的质量数对酶解后的肽段建立肽段数据库和肽段数据库索引；

(2)根据待分析实验图谱中母离子的核质比在步骤(1)所述的肽段数据库中找出符合要求的候选肽段，并对找到的所有候选肽段产生符合要求的理论图谱；

(3)对待分析实验图谱进行去同位素和去噪处理；

(4)将步骤(3)得到的待分析实验图谱和步骤(2)中得到的每张候选肽段的理论图谱进行匹配打分，选择分值最高的候选肽段作为此实验图谱的鉴定结果；

(5)针对所有实验图谱鉴定结果进行整体假阳性控制。

本发明一种蛋白质数据库序列进行虚拟酶解并对酶解后肽段建立肽段数据库索引方法，包括：

步骤1)读取质谱分析样本(待分析二级质谱的样本)的物种蛋白质序列库文件的一条蛋白质序列。

步骤2)按照表1根据用户设定蛋白酶和容许的漏切位点个数对此蛋白质序列进行虚拟理论酶切。目前大部分Trypsin酶进行蛋白质酶解实验，从表第一行可知Trypsin酶是对蛋白质C-Term敏感的，也就是说蛋白质序列C端可能会被切掉一个氨基酸；其酶切位点KR，也就是说其酶在序列的K和R氨基酸上发生酶切作用；其限制酶切位点是P，也就是说序列K和R氨基酸上发生酶切时，如果其后面一个氨基酸是P则不能发生酶切作用。

表1蛋白酶酶切位点表

蛋白质酶	敏感端	酶切位点	限制酶切位点
				Trypsin	C-Term	KR	P
Arg-C	C-Term	R	P
				Asp-N	N-Term	D
Asp-N_ambic	N-Term	DE
				Chymotrypsin	C-Term	FLWY	P
CNBr	C-Term	M

其步骤2)详细过程是：

A根据表1找到蛋白质序列中包含符合上面规则的酶切位点；

B在符合规则的酶切位点产生断裂，产生没有漏切位点的肽段；

C产生存在漏切位点的断裂肽段。

图2给出了蛋白质虚拟Trypsin酶解示意图。将一个蛋白质序列按照上述规则虚拟酶解，最后得到了7个肽段。

步骤3)计算每个虚拟酶切后肽段的质量，根据每个氨基酸的分子量计算每个肽段的质量数；由于计算肽段质量数计算频率高，在计算质量数之前首先对每个氨基酸的质量建立索引，如表2所示，其对20个氨基酸的索引和翻译后修饰的索引方法如下：

A启用一个数组与ASCII码等同大小的数组(大小为250)；

B这个数组的下标与氨基酸单字母简写的ASCII码数值一致，其数组中保存氨基酸的质量数。除了20氨基酸的位置放置没有修饰的氨基酸，其它位置(大概有230)个可以处理翻译后修饰，此方法可以同时处理230种修饰。

表220个氨基酸索引表

数组	氨基酸简写	数组值	化学组成
				AA(1)		14.00307	N
AA(2)		15.99491	O
				AA(3)		1.007825	H
AA(4)		12	C
				AA(65)	A	71.037114	H(5)C(3)NO
AA(66)	B	115.02694	H(5)C(4)NO(3)
				AA(67)	C	103.0092	H(5)C(3)NOS
AA(68)	D	115.026943	H(5)C(4)NO(3)
				AA(69)	E	129.04259	H(7)C(5)NO(3)
AA(70)	F	147.06841	H(9)C(9)NO
				AA(71)	G	57.02146	H(3)C(2)NO
AA(72)	H	137.05891	H(7)C(6)N(3)O
				AA(73)	I	113.08406	H(11)C(6)NO
AA(75)	K	128.09496	H(12)C(6)N(2)O
				AA(76)	L	113.084064	H(11)C(6)NO
AA(77)	M	131.040485	H(9)C(5)NOS
				AA(78)	N	114.042927	H(6)C(4)N(2)O(2)
AA(80)	P	97.052764	H(7)C(5)NO
				AA(81)	Q	128.058578	H(8)C(5)N(2)O(2)
AA(82)	R	156.101111	H(12)C(6)N(4)O
				AA(83)	S	87.032028	H(5)C(3)NO(2)
AA(84)	T	101.047679	H(7)C(4)NO(2)
				AA(86)	V	99.068414	H(9)C(5)NO
AA(87)	W	186.079313	H(10)C(11)N(2)O
				AA(89)	Y	163.063329	H(9)C(9)NO(2)

之后，把单字母表示氨基酸的肽段序列中每个字母依次转换成其对应ASCII码的数值，根据氨基酸索引表的数值计算肽段的质量，例如：假设有一个肽段为ACD，那么肽段ACD的ASCII码数值是65、67、68，那么其肽段的质量数为数组AA下标为65、67、68的值的和并加上水的分子量，因为肽段有C(H)和N端(OH)，即：

2*AA(3)+AA(2)+AA(65)+AA(67)+AA(68)

＝2*1.007825+15.99491+71.037114+103.0092+115.026943＝307.0838

根据氨基酸索引表计算出每条虚拟酶解后的肽段的质量数。

步骤4)把计算得到质量数的肽段放入肽段数据库，以每1da为单位将所有酶解后肽段分别存入相应的文件中。把肽段的质量数取整，例如307.0838取整为307，之后将肽段的信息在质量数取整的文件中末尾追加保存，即在文件名为307的文件末尾追加一行存入肽段的信息。按照上面方法把每条肽段放入肽段数据库。

步骤5)读取下一条蛋白质序列，重复步骤2)，3)，4)，直到所有的蛋白序列被酶解和存入肽段数据库。

步骤6)合并每1da为单位文件的肽段信息并对其建立索引文件：按文件名的数字从小到大读出文件中的肽段信息，每读一个文件，按照文件中肽段的质量数从小到大进行排序，之后从小到大顺序存入database.ind文件中，并删除每个读取肽段信息文件。例如文件名为1000的文件保存着质量数为1000da-1001da的所有肽段的信息，读取其文件的肽段信息，并排序，之后将排序后肽段信息存入database.ind文件中，并删除文件名为1000的文件。将database.ind每行存入一个肽段，其文件格式如下表3所示。与此同时，按照1da对酶解所有肽段建立查找索引database.index，其查找索引记录下信息：1.第一列保存其质量数，例如1000，表示质量数位为1000da-1001da肽段，第二列是这些肽段在database.ind文件开始位置，第三列是酶解肽段在1000da-1001da的个数，即1000da-1001da肽段在database.ind文件中的行数。根据database.index可以知道1000da-1001da在文件database.ind中的位置。其结构如表4。

表3database.ind索引表

表4database.index索引表

肽段质量数索引编号	文件开始位置	肽段数量
			1005	0	2
1064	56	2
			1089	224	2
1106	282	2
			1117	340	4

本发明提供了根据待分析二级质谱母离子查找候选肽段的方法，实现方法步骤：

步骤1)加载database.index文件信息到内存数组index，读取待分析二级质谱的母离子mz值和电荷信息，并计算其母离子去电荷后的质量数，例如有一个mz＝2100.2，charge＝2的母离子信息，其去电荷的质量数为mz*2-2＝4198.2。

步骤2)根据容许的质量误差查找index数组记录并读取相应肽段信息，假设质量误差为0.01.4198.2-0.01＝4198.1和4198.2+0.01＝4198.3，4198.1和4198.3取整都为4198da，查找index数组找到其在文件database.ind中的开始位置和行数，在此位置开始顺序读取相应的行数加入内存中，即加载了4198～4199Da的所用肽段信息。

步骤3)对内存加载肽段进行步的筛选，由于内存加载的是1da为单位的所有肽段，然而在实际中，根据质谱仪型号的不同，其精确度也不同，例如Orbitrap质量误差为0.1，因此要对以上加载到内存肽段按照用户实验质谱类型进行进一步的筛选，即按照质谱质量误差为0.1，筛选出质量数范围在4198.1～4198.3Da的肽段，作为此待分析二级质谱的候选肽段。

本发明提供了对待分析实验图谱快速去同位素和去噪处理(选峰)的方法，如图4所示，其方法包括以下步骤：

步骤1)待分析二级实验图谱进行去同位素峰处理。理论上同位素峰之间质荷比mz相差1且同位素峰之间的峰强受自然界同位素丰度控制，例如自然界C12丰度高于C13的丰度，其质谱峰的高度也高于C13。自然界中稳定同位素中，低分子量的基本上丰度都占其丰度的最高位。在质谱中，一个同位素峰群中，第一峰基本上应该是最高峰。实际质谱仪的测量中，由于质谱仪都存在测量误差。根据质谱仪类型不同，其测量的精确度也不同，例如LTQ质谱仪的测量误差是0.5Da。由于一张质谱的系统误差一样，也就是说同位素峰总是向右或向左偏离理论值，因此认为两个峰mz1和mz2符合|mz1-mz2-1|＜0.25da既为同位素峰。基于以上人类的认知，再结合质谱实际情况，这里提出以下快速去同位素峰的方法：

A.初始化，将三个比较峰的mz值和其强度，全部设为0(假设三个峰mz值是：mz_1＝0，mz_2＝0，mz_3＝0，其峰强对应是mz_1_in＝0，mz_2_in＝0和mz_3_in＝0，并设置保留峰容器(用于存储非同位峰)；

B.读取一个峰的信息，假设mz_curr＝245，in_curr＝80，测量质量误差m＝0.25.

a.把目前的峰放入第三个峰的位置，即mz_3＝mz_curr，mz_3_in＝mz_curr；

b.第三个峰与第一个峰和第二个峰比较，判断是否是前两个峰的同位素峰。即：如果以下三个条件任意一个条件成立，认为是同位素峰，

a.|mz_3-mz_2-1|＜＝m并且mz_2_in＞mz_3_in

b.|mz_3-mz_1-1|＜＝m并且mz_1_in＞mz_3_in

c.|mz_2-mz_1|＜＝m并且mz_2_in＞mz_3_in(此为相同的峰信息，记录误差)

执行三个峰向前平移一位，空出第三个峰的位置即：

mz_1＝mz_2，mz_1_in＝mz_2_in

mz_2＝mz_3，mz_2_in＝mz_3_in

否则，认为目前进入第三位置的峰不是同位素峰，将其作为保留峰存入保留峰容器中，并三个峰向前平移一位，空出第三个峰的位置。即：

mz_1＝mz_2，mz_1_in＝mz_2_in

mz_2＝mz_3，mz_2_in＝mz_3_in

C.逐个读取下一个峰的信息，重复2)直到处理完一张二级质谱图所有峰信息，最后得到的保留峰容器中的峰即为去同位素峰后的非同位素峰。

例如：图4给出了原始4个峰的去同位素执行过程中三个峰和保留峰的状态改变过程。

步骤2)对去同位素峰后的非同位素峰(保留峰)进行去噪处理，去噪处理就是除去信号峰的噪声峰，以免影响鉴定结果，在质谱数据处理中也叫做选峰机制，既选择非噪声信号的峰。此发明中的选峰的规则是：

(1)选取局部最强峰：对去同位素后的离子峰，每100Da选取六个局部最高强度峰，这样选取离子峰是因为氨基酸分子的平均质量是100Da，而且本方法考虑六种离子类型的峰(b，y，b-H₂O，y-H₂O，b-NH₃，y-NH₃)的匹配情况，这样选择的目的是能够选取到局部最强峰。此选峰方法区别与其他选峰方法，具体方法如下：

(a)找到全局最强峰，在图5中为红色的实线。

(b)在最强峰度值对应的质荷比为952.449Da分别向左右各平移50Da视为窗口1(见图5实线所包括区域)，也就是对窗口1中的质荷比值在902.449Da到1002.449Da范围内的质荷比信息按其峰的强度从高到底进行排序，挑选离子峰强度排名前六位的峰的信息(峰强度最强的六个峰)。

(c)在窗口1的左边界向左平移50Da，窗口1的右边界向右平移50Da(图5中的虚线所包括窗口范围内)定义为窗口2，在窗口2范围内的所有峰的信息按峰强度从高到低进行排序，挑选出排名前六位的质荷比信息。

(d)重复进行(c)两步，直到该质谱文件所有的质荷比信息被提取完成。

(2)选取全局最强峰：选取大于等于全局最高峰强度33％的峰，同时是上一步局部选峰信息没有选到的峰进行保留，具体步骤如下：

(a)找出全局最强峰，如图6中所示，得到峰的强度值为maxintense＝6602.7。

(b)计算全局最强峰的阈值，threshold＝maxintense*0.33，在这里阈值等于2178.9。

(c)在图6窗口中大于阈值(图6中水平点划线)以上的峰强信息未被在步骤(1)中选取为强峰的，在这里要把这些峰加载到全局最强峰中。

图7为对一张二级质谱执行(1)和(2)两步后最终选取峰的信息图。

本发明提出了一种产生候选肽段的理论图谱方法规则，其方法包含以下内容：

步骤1)产生候选肽段可能产生的理论碎片b、y离子：由于肽段中肽键位置键能最低，最容易断裂，因此b、y离子一般在二级质谱峰中大量存在。例如肽段SIEFYEDYYGVK经过离子碎裂(CID)过程产生基本的碎片离子(b离子、y离子)如图8所示，在肽段的上部和下部分别表示出了b和y离子的类型，其中b离子是从左边往右标记，比如b1为S、b2为SI、b3为SIE以此类推；而y离子是从右往左进行标记，y1为K、y2为VK、y3为GVK以此类推。

步骤2)在步骤1)产生b、y离子基础上，产生b、y离子丢水(b-H₂O，y-H₂O)和丢氨(b-NH₃，y-NH₃)情况下的碎片离子。然而不是每一个b、y离子都会产生丢水(b-H₂O，y-H₂O)和丢氨的碎片离子。根据下面两条的规则产生丢水和丢氨碎片离子：

(1)如果在b离子或y离子中含有S、T、E和D四种氨基酸中的一种，那么产生该离子对应的丢水碎片离子(b-H₂O和y-H₂O)。

(2)如果在b离子或y离子中含有R、K、Q和N四种氨基酸中的一种，那么该离子要考虑其对应的丢氨碎片离子(b-NH₃和y-NH₃)。

例如肽段SIEFYEDYYGV的前四个b离子及丢水丢氨的离子信息见表5，表5中给出了上面所述四种离子类型，在这里只产生了四种离子类型，其表中值为0的是没有产生相应丢失情况。表5中的第一列数据代表b1、b2、b3、b4的核质比，第二列核质比值都为零，其原因是从肽段前四个氨基酸SIEF中包含在R、K、Q和N这四种氨基酸的任意一种，按上规则(2)，不产生这些离子丢水情况。所，第三列的b离子丢水信息，从肽段中可以发现它们中都包含着S氨基酸，所以产生丢水离子。

表5一价理论质谱信息

步骤3)产生二价的碎片离子，其规则为：待分析二级质谱母离子价态是一价的则考虑一价的碎片离子峰，当母离子价态＞＝2时，并且对应的碎片离子中包含R、K和H三种氨基酸其中一种时，则考虑二价碎片离子峰。

假设待分析二级质谱母离子价态为3，则要考虑表5中肽段SIEF的二价离子(见表6)，因为在这几个类型离子中没有包含R、K和H氨基酸中的任何一种，因此不产生二价的碎片离子(其值都为0)。

表6肽段SIEF的二价碎片离子

本发明提出了一种待分析选峰后的实验二级质谱与候选肽段的理论图谱进行匹配打分的方法，其打分方法包含以下内容：

步骤1)候选肽段理论图谱与选峰后的实验二级质谱进行匹配打分，其打分采用以下公式进行匹配打分

$\left\{\begin{array}{c}Q=i+\mathrm{factor}\\ P=\left(\begin{array}{c}K\\ E\end{array}\right)Q^K{(1-Q)}^{E-k}\end{array}\right.$

公式中i＝0.06，由于每100da选取6个峰，其随机匹配概率为0.06；其中factor等于实验图谱选取大于33％峰个数/实验质谱的峰范围，其原理与i相同；E为产生的理论碎片的个数；K为实验图谱匹配峰的个数，Q代表随机匹配概率事件，P由二项式分布概率密度函数计算，含义是在E个理论峰中有K个峰匹配的概率。

匹配步骤1)之前包括以下过程：

(1)对某张特定二级质谱进行去同位素峰和选峰。按照上述规则对其进行去同位素峰和选峰，并获得factor因子的值。

(2)取出查找肽段数据库的候选肽段集中取出一个候选肽段，按上述产生理论质谱的步骤2、步骤3产生理论图谱。

二级质谱匹配打分包括以下过程：

A.逐个读取峰信息查找理论图谱和选峰后的实验图的匹配是否，如果理论图谱与实验图谱对应峰的核质比之差小于等于质谱仪的测量误差，则认为这两个峰匹配，之后记录其匹配的信息；本实施例中质谱仪的测量误差为0.5。

B.统计步骤3)中理论图谱和选峰后的实验图匹配个数K，并且统计理论图谱峰的个数E，根据上面公式计算P值。

例如100峰选峰后的信息如图9所示，其中误差容限(m＝0.5)，对于一个候选肽段SIEFYEDYYGVK理论质谱理论图谱如表5、6，每个理论质谱的mz值通过查找实验图谱看是否在实验图谱中有匹配。查找后其匹配情况如表7、8所示，带括号的mz值表示理论图谱与实验图谱有匹配，其理论图谱在实验图谱的匹配表示图，其中有标识峰为实验图谱的匹配峰。

同时作出匹配标注图，如图9所示，通过匹配标注图可以看出那些离子匹配成功，得到匹配后的表格。

表7一价离子匹配

b	b-NH3	b-H2O	y	y-NH3	y-H2O
						88.03931	0	70.02871	147.1128	130.0863	0
201.1234	0	183.1128	(246.1812)	229.1547	0
						(330.166)	0	(312.1554)	(303.2027)	286.1762	0
(477.2344)	0	(459.2238)	(466.266)	(449.2395)	0
						(640.2977)	0	(622.2871)	(629.3293)	(612.3028)	0
769.3403	0	751.3297	(744.3563)	(727.3298)	726.3457
						(884.3672)	0	(866.3566)	(873.3989)	856.3724	(855.3883)
(1047.431)	0	(1029.42)	(1036.462)	(1019.436)	(1018.452)
						(1210.494)	0	(1192.483)	(1183.531)	1166.504	(1165.52)
(1267.515)	0	(1249.505)	(1312.573)	(1295.547)	(1294.563)
						(1366.584)	0	(1348.573)	1425.657	1408.631	1407.647

表8二价离子匹配表

b++	b++-NH3	b++-H2O	y++	y++-NH3	y++-H2O
						0	0	0	74.06004	65.54679	0
0	0	0	123.5942	115.081	0
						0	0	0	152.105	143.5917	0
0	0	0	233.6366	(225.1234)	0
						0	0	0	315.1683	306.6551	0
0	0	0	372.6818	364.1685	363.6765
						0	0	0	437.2031	428.6898	428.1978
0	0	0	(518.7347)	510.2215	509.7294
						0	0	0	592.2689	583.7557	583.2636
0	0	0	(656.7902)	648.277	647.7849
						0	0	0	713.3323	704.819	704.327

计算其值为：统计以上匹配信息：其为K＝65，E＝26，factor＝0，i＝0.06按照公式2计算匹配概率P1

步骤2)对步骤1)查找到的匹配情况进行连续匹配分析并进行统计打分。在考虑连续匹配的情况中，为了简化模型和提高速度，本方法只考虑两个离子的连续的情况，如果有三个离子连续，例如，b1、b2、b3匹配，这里可把三个连续匹配转换成两个连续的匹配也就是b1和b2、b2和b3，多个连续匹配都可以依此理转换成两个连续匹配。统计理论的理论连续匹配的情况和实验图谱实际连续匹配的情况，采用以下公式计算连续匹配得分：

$\left\{\begin{array}{c}Q1=B\_\mathrm{factor}*K/E\\ P1=\left(\begin{array}{c}K1\\ E1\end{array}\right){Q1}^{K1}{(1-Q1)}^{E1-K1}\end{array}\right.;$

设E1为理论图谱产生的理论连续匹配个数；K1为实验图谱实际连续匹配的个数；B_factor为背景值，通过统计大量的鉴定结果获取，B_factor＝统计大量实验图谱连续匹配的平均值/统计大量对应理论图谱连续匹配平均值，B_factor反映了实际连续匹配的概率，其统计后值为0.00843，其中K、E上步骤1)中的K和E，Q1反映了某一图谱在步骤1)匹配情况下连续匹配的概率，P1由二项式分布概率密度函数计算，含义是在E1个理论连续中实际存有K1个连续匹配的概率。

根据连续匹配方法，图9中匹配标注图可以看(y2，y3，y4，y5，y6，y7，y8，y9，y10)，(b3，b4，b5)，(b7，b8，b9，b10，b11)都属于三个以上连续的情况，比如(b3，b4，b5)也可以等价于两个连续b3和b4连续，b4和b5连续，同理，(y2，y3，y4，y5，y6，y7，y8，y9，y10)可以等价于8个连续，依次类推，表7和表8有57个理论连续即E1＝57，实验图谱连续匹配为17个既K1＝17；按公式3计算连续匹配的概率P2。

步骤3)对步骤1)查找到的匹配情况进行峰强度因素考虑，即对匹配峰强度信息进行考虑。此发明中的峰强度考虑是通过大量统计分析实际鉴定图谱鉴定结果来获得峰强度概率表，这里认为碎片离子峰强度产生与b、y离子的物理化学属性有关，也就是说某一碎片离子可能总是出现比较高的峰。从化学角度看，20种氨基酸与b、y离子强度与两边断裂离子的氨基酸化学键有关，表7中的峰强度概率由以下公式获取：

Yi＝M_I/M_E

其中M_I为统计所有实验图谱中某两个氨基酸产生的峰大于等于最强峰的33％的个数；M_E为期望总的离子的个数，Yi为两个氨基酸中间断裂概率，计算强度因子可以通过查表4得出某张二级质谱超过33％峰匹配峰的因子，并按照下面公式累加每个匹配峰(＞＝33％的峰)进行累加得到Ym和Bm。

Ym＝∑Y_i

Bm＝∑B_i

进而得到强度因子Infactor为：

Infactor＝(1+Ym+Bm))/(1+0.155*m_p)；

其中Ym、Bm分别为实验图谱强度大于全局最强峰的33％的匹配峰Yi和Bi分值之和；m_p为一张实验图谱中强度大于最强峰的33％的匹配的个数；其中0.155是理论平均匹配值。

表9峰强概率统计表

例如上面匹配中有多个实验图谱峰全局最高峰的33％，并且有多个匹配上理论图谱，通过肽段序列可以知道每个匹配峰在那两个氨基酸位置断裂，通过表9可以查出相应的概率值，按上面公式计算得到infactor＝1.3813。

步骤4)有机组合上述步骤1)，2)和3)的打分方法，采用下面公式得到肽段的得分：

PEP-S＝Infactor*(-10)*log₁₀ ^(P*P1)；

例如肽段SIEFYEDYYGV上述结果计算得到PEP-S＝1.3813*257.8791。

步骤5)：对步骤4)计算的PEP-S分子去除背景值，PEP-S分数很可能与肽段长度、母离子价态等因素有相关性，也就是说，有可能PEP-S一价母离子产生的分数都大于二价母离子产生的分数，所以导致二价母离子的图谱被很少的鉴定。采用以下公式去背景值：

PEP-S_M＝PEP-S-B_B；

设在真实库和随机库统计概率相等时的背景值为其在某种情况下的背景值B_B，经过贝叶斯网络学习，各种情况背景值如表10所示。

表10背景值表

例如：肽段SIEFYEDYYGV没有修饰位点其背景值是8，没有漏切位点其背景值为8，实验图谱母离子为2价其背景值为15，肽段质量数的16.4102，去除背景值：

PEP-S_M＝257.8791*1.3813-16.4102-8-8-15＝308.7982。

步骤6)取出下一个肽段，重复执行前面五步，直到符合此图谱母离子误差的所有肽段均被打分处理。对上面实验图谱候选肽段重复前面五步得到肽段得分，如表11所示。

表11二级图谱(编号100)候选肽段打分结果

步骤7)对此图谱所有候选肽段的得分PEP-S_M从大到小进行排序，排名第一位的即为此图谱的鉴定结果。表格中VNDEIEIVGIK即为上述一张二级图谱的鉴定结果。

本专利提供了一种对整体鉴定结果假阳性率控制的方法，上述过程是对一张二级质谱进行分析，得到一张二级质谱鉴定结果，对于待分析的每张二级质谱重复上序过程，就得到每张二级质谱的鉴定结果。此结果用以下方法进行整体的质量控制，其步骤包括：

步骤1)统计待分析实验图谱所有二级图谱中的鉴定结果，取出肽段得分最小值和最大值；

步骤2)按分数统计待分析图谱的鉴定结果，统计最小值和最大值之间大于各得分中的正库肽段和随机库肽段的个数，例如大于等于某一个分值其鉴定结果在真实库中的个数N_N和在随机数据库中的个数N_R，按照下述公式计算每个分值为阀值时FDR的值，

$\mathrm{FDR}=\frac{N_R}{N_N}$

步骤3)按得分值从小到大寻找每个分值，直到找到FDR＜＝0.01时，该分值为其待分析图谱的整体阀值；

步骤4)根据第三步找到整体阀值，以此阀值过滤待分析图谱的鉴定结果，也就是说小于此阀值结果被过滤掉，其结果作为待分析实验图谱最终鉴定结果。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (8)

1.一种基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，包括以下步骤：

（1）虚拟酶解蛋白质数据库序列，并根据肽段的质量数对酶解后的肽段建立肽段数据库和肽段数据库索引；

（2）根据待分析实验图谱中母离子的核质比在步骤（1）所述的肽段数据库中找出符合要求的候选肽段，并对找到的所有候选肽段产生符合要求的理论图谱；

（3）对待分析实验图谱进行去同位素和去噪处理；去噪处理，即除去信号峰中的噪声峰，具体步骤是：

（3-2-1）首先选取局部最强峰，包括以下步骤：

（3-2-1-1）根据得到的去同位素后的离子峰，找到全局最强峰，然后以此峰为中心，分别向左右各平移50Da，形成一个搜索窗口，在这100Da范围内挑选离子峰强度排名前n位的峰，然后记录这n个峰的信息；

（3-2-1-2）以已搜索区域为中心，再分别向左右各平移50Da，在左右各形成1个搜索区域，在这100Da范围内挑选离子峰强度排名前n位的峰，然后记录这n个峰的信息；

（3-2-1-3）重复进行（3-2-1-1）和（3-2-1-2）两步，直到该质谱文件所有的质荷比信息被提取完成；

（3-2-2）根据步骤（3-2-1-1）得到的全局最强峰，搜索峰值大于等于全局最强峰峰值*0.33的峰，作为全局相对高峰，判断这些峰是否已记录在步骤（3-2-1）中，是则不做处理，否则记录峰的信息；

（3-2-3）将选取的局部最强峰和全局相对高峰进行合并，得到最终选取的用于鉴定的峰；

（4）将步骤（3）得到的待分析实验图谱和步骤（2）中得到的每张候选肽段的理论图谱进行匹配打分，选择分值最高的候选肽段作为此实验图谱的鉴定结果；具体步骤如下；

（4-1）将待分析实验图谱和理论图谱进行匹配打分的具体步骤是：

（4-1-1）逐个读取峰信息判断理论图谱和选峰后的实验图谱是否匹配，如果理论图谱与实验图谱对应峰的核质比之差小于等于质谱仪的测量误差，则认为这两个峰匹配，之后记录其匹配的信息；

（4-1-2）设E为产生的理论碎片的个数；K为理论图谱和选峰后的实验图谱匹配个数，Q代表随机匹配概率事件，i为随机匹配概率，i=0.01*n，P为在E个理论峰中有K个峰匹配的概率，则P由下面二项式分布概率密度函数计算：

$\left\{\begin{array}{c}Q=i+\mathrm{factor}\\ P=\left(\begin{array}{c}K\\ E\end{array}\right)Q^K{(1-Q)}^{E-K}\end{array}\right.;$

（4-2）将待分析实验图谱和理论图谱进行连续匹配打分的具体步骤是：

设E1为理论图谱产生的理论连续匹配个数；K1为实验图谱实际连续匹配的个数；B_factor为背景值，B_factor=统计大量实验图谱连续匹配的平均值/统计大量对应理论图谱连续匹配的平均值，Q1反映了某一图谱在步骤（4-1）匹配情况下连续匹配的概率，P1是在E1个理论连续匹配个数中实际存有K1个连续匹配的概率，由下面二项式分布概率密度函数计算：

$\left\{\begin{array}{c}Q1=B\_\mathrm{factor}*K/E\\ P1=\left(\begin{array}{c}K1\\ E1\end{array}\right){Q1}^{K1}{(1-Q1)}^{E1-K1}\end{array}\right.;$

所述待分析实验图谱和理论图谱连续匹配个数具体是指图谱中两两连续匹配的对数；

（4-3）对匹配峰强度信息进行分析，求得强度因子，具体步骤是：

设M_I为统计所有实验图谱中某两个氨基酸产生的峰大于等于最强峰的33%的个数，设M_E为期望总的离子的个数，则两个氨基酸中间的断裂概率Yi通过下式得到：

Yi=M_I/M_E；

进而得到强度因子Infactor为：

Infactor=(1+Ym+Bm）)/(1+0.155*m_p)；

其中Ym=∑Y_i；Bm=∑B_i；Ym、Bm分别为实验图谱强度大于全局最强峰的33%的匹配峰Yi和Bi分值之和；m_p为一张实验图谱中强度大于最强峰的33%的匹配的个数；其中0.155是理论平均匹配值；

（4-4）有机组合上述步骤（4-1），（4-2）和（4-3）的打分方法，采用下面公式得到肽段的得分：

PEP-S=Infactor*（-10）*log₁₀ ^（P*P1）；

（4-5）对计算的PEP-S分数去除背景值，首先设在真实库和随机库统计概率相等时的背景值为其在某种情况下的背景值B_B，背景值B_B是经过贝叶斯网络学习得到的；

然后得到去背景值PEP-S_M：

PEP-S_M=PEP-S—B_B；

（4-6）取出下一个肽段，重复执行步骤（4-1）-（4-5），直到符合此图谱母离子误差的所有肽段均被打分处理；

（4-7）对此图谱所有候选肽段的得分PEP-S_M进行排序，值最大的作为当前图谱的鉴定结果；

（5）针对所有实验图谱鉴定结果进行整体假阳性控制。

2.根据权利要求1所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（1）具体包括以下步骤：

（1-1）读取待分析二级质谱样本中物种蛋白质序列库文件的一条蛋白质序列；

（1-2）根据用户设定的蛋白酶，找到蛋白质序列中的酶切位点，在符合规则的酶切位点产生断裂，从而产生没有漏切位点的肽段或存在漏切位点的断裂肽段；

（1-3）计算步骤（1-2）所得到的各个虚拟酶切后肽段的质量数，根据每个氨基酸的分子量计算每个肽段的质量数；

（1-4）将肽段信息写入肽段数据库中以该肽段取整后质量数命名的文件中；

（1-5）读取下一条蛋白质序列，重复步骤（1-2）-（1-4），直到所有的蛋白序列被酶解和存入肽段数据库；

（1-6）按文件名的数字从小到大读出文件中的肽段信息，每读一个文件，按照文件中肽段的质量数从小到大进行排序，然后存入一个数据库文件database.ind中，同时，以1da为单位对所有肽段建立查找索引database.index，其查找索引包括以下信息：其质量数，这些肽段在database.ind文件中的开始位置，该区间内的肽段的个数。

3.根据权利要求2所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（1-3）中在计算质量数之前首先对每个氨基酸的质量建立索引，其对20个氨基酸的索引和翻译后修饰的索引方法如下：

（1-3-1）启用一个与ASCII码等同大小的数组，该数组的下标与氨基酸单字母简写的ASCII码数值一致，其数组中保存氨基酸的质量数；

（1-3-2）把单字母表示氨基酸的肽段序列中每个字母依次转换成其对应ASCII码的数值，然后根据氨基酸索引表的数值计算每条虚拟酶解后的肽段的质量数。

4.根据权利要求2所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（2）在肽段数据库中找出符合要求的候选肽段的具体步骤是：

（2-1-1）加载步骤（1-6）中的database.index文件信息到内存数组index，读取待分析二级质谱的母离子核质比值和电荷信息，并计算其母离子去电荷后的质量数；

（2-1-2）根据容许的质量误差和步骤（2-1-1）所述的质量数在index数组中查找相应肽段在文件database.ind中的开始位置和行数，然后加载此区间内的所有肽段信息；

（2-1-3）根据用户所采用质谱仪的精确度，对步骤（2-1-2）加载到内存的肽段进行进一步的筛选，作为此待分析二级质谱的候选肽段。

5.根据权利要求1所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（2）中产生符合要求的理论图谱具体包括以下步骤：

（2-2-1）在肽段的上部和下部分别表示出b和y离子的类型，其中b离子是从左边往右标记，y离子是从右往左进行标记，模拟经过离子碎裂过程候选肽段可能产生的理论碎片b、y离子；

（2-2-2）在步骤（2-2-1）产生b、y离子基础上，根据下面两条的规则产生b、y离子丢水和丢氨情况下的碎片离子：

（2-2-2-1）如果在b离子或y离子中含有S，T，E和D四种氨基酸中的一种，那么产生该离子对应的丢水碎片离子b-H₂O和y-H₂O；

（2-2-2-2）如果在b离子或y离子中含有R，K，Q和N四种氨基酸中的一种，那么产生该离子对应的丢氨碎片离子b-NH₃和y-NH₃；

（2-2-3）产生二价的碎片离子，其规则为：待分析二级质谱母离子价态是一价的则考虑一价的碎片离子峰，当母离子价态>=2时，并且对应的碎片离子中包含R、K和H三种氨基酸其中一种时，则考虑二价碎片离子峰，在理论图谱中产生其对应的二价的碎片离子峰。

6.根据权利要求1所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（3）对待分析实验图谱进行去同位素的具体步骤是：

（3-1-1）初始化，将三个比较峰的mz值和其强度，全部设为0，设三个峰mz值是：mz_1=0,mz_2=0,mz_3=0，其峰强对应是mz_1_in=0,mz_2_in=0和mz_3_in=0，并设置用于存储非同位峰的保留峰容器，并已知测量质量误差m；

（3-1-2）读取一个峰的信息，将其作为第三个峰的值，即mz_3和mz_3_in；然后判断第三个峰是否是前两个峰的同位素峰，

（3-1-2-1）如果以下三个条件任意一个条件成立，则认为是同位素峰，

a.|mz_3-mz_2-1|<=m并且mz_2_in>mz_3_in；

b.|mz_3-mz_1-1|<=m并且mz_1_in>mz_3_in；

c.|mz_2-mz_1|<=m并且mz_2_in>mz_3_in，此为相同的峰信息，记录误差，然后进入步骤（3-1-2-3）；

（3-1-2-2）如果步骤（3-1-2-1）中各条件均不成立，则认为目前进入第三位置的峰不是同位素峰，将其作为保留峰存入保留峰容器中，然后进入步骤（3-1-2-3）；

（3-1-2-3）三个峰向前平移一位，空出第三个峰的位置，即将第一个峰的信息去掉，第三个和第二个峰的信息分别作为新的第二个和第一个峰的信息；

（3-1-3）逐个读取下一个峰的信息，重复步骤（3-1-2），直到处理完一张二级质谱图所有峰信息，最后得到的保留峰容器中的峰即为去同位素峰后的非同位素峰。

7.根据权利要求1所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（5）针对所有实验图谱鉴定结果进行整体假阳性控制，具体包括以下步骤：

（5-1）统计待分析实验图谱所有二级图谱中的鉴定结果，取出肽段得分最小值和最大值；

（5-2）按分数统计待分析图谱的鉴定结果，统计最小值和最大值之间大于各得分中的正库肽段和随机库肽段的个数，设大于等于某一个分值时，其鉴定结果在真实库中的个数为N_N，在随机数据库中的个数为N_R，按照下述公式计算每个分值为阀值时FDR的值，

$\mathrm{FDR}=\frac{N_R}{N_N}$

（5-3）将产生的FDR<=0.01的值最小的分值作为待分析实验图谱的整体阀值；

（5-4）以步骤（5-3）得到的整体阀值过滤待分析实验图谱的鉴定结果，如果鉴定分数小于此阀值则被过滤掉，其结果作为待分析实验图谱最终鉴定结果。

8.根据权利要求1所述的基于概率统计模型的蛋白质二级质谱鉴定方法，其特征在于，所述步骤（3-2-1）和步骤（4-1）中的n取值范围为3～6。

Images