CN106529204A - 一种基于半监督学习的交联质谱多谱排序方法 - Google Patents

Abstract

本发明提供一种基于半监督学习的交联质谱多谱排序方法，包括：1)分别对每个谱图进行单谱匹配和排序，得到对应的最优的交联二肽单谱匹配结果；提取当前每一个肽谱匹配结果的多谱匹配特征向量，其中包括SVM分数、母离子误差比例特征和修饰比例特征等动态特征；2)在所得到的交联二肽匹配结果中，取FDR在预设的FDR阈值以内的属于正样本的结果构建正样本库，取所有负样本的结果构建反样本库；基于新的训练样本更新各个多谱匹配特征向量；3)训练SVM分类器；4)用本轮训练后的SVM分类器对所有交联二肽结果进行重打分；5)根据预设的迭代条件判断是否继续进行迭代，结束迭代时基于当前SVM分数输出多谱排序结果。本发明的多谱排序方法灵敏度高且性能稳定。

Description

一种基于半监督学习的交联质谱多谱排序方法

技术领域

本发明涉及生物信息学技术，具体地说，本发明涉及交联质谱的鉴定技术。

背景技术

目前，化学交联质谱技术已成为研究蛋白质结构与蛋白间相互作用的主流技术。现有技术中，由于交联二肽的搜索空间巨大，在计算方法方面的研究相对滞后，在大规模数据库下鉴定交联蛋白质存在精度低的问题。由于单谱内的打分函数的归一化效果较弱，不适合多谱之间肽谱匹配情况的比较。因此，通常需要使用多谱排序算法来改善交联肽段鉴定的灵敏度，目前存在两类多谱排序算法，一类是计算E-value(E-value是用于度量肽谱随机匹配的可能性的参量)，然后基于E-value进行排序，另一种是采用机器学习方法进行排序。下面简要介绍这两种方法，然后分析这两种方法的不足。

第一种方法的代表是pLink 1.0软件。如前文所述，E-value一般用来度量肽谱随机匹配的可能性，它反映的是匹配统计显著性。在pLink 1.0软件中，E-value为给定谱图对应的打分不低于t的错误肽谱匹配的数目的期望，其中，t为当前肽谱匹配的细打分分数。这种方法虽然能够实现多谱排序，但对交联肽段鉴定的灵敏度的改善效果有限。

第二种方法的代表软件有Protein Prospector、MXDB、Kojak、pLink 2.0等。其中，Protein Prospector和MXDB两款软件均采用监督学习中的支持向量机(SVM)模型进行交联肽段鉴定的多谱排序工作。其中，Protein Prospector使用了两个特征构建特征向量，其中一个是为了刻画匹配较差的那条肽段的匹配情况的特征，另一个是与匹配谱峰的覆盖率相关的特征。MXDB软件使用了15个特征，其中有7对特征是分别刻画交联肽段中两条肽段各自匹配情况的特征，剩下一个特征是考虑两条肽段匹配情况的综合分数。Kojak软件将搜索阶段中打分第一名的结果作为第三方软件Percolator的输入，Percolator软件利用半监督学习中的SVM模型来对所有的鉴定结果进行重打分。pLink 2.0软件也是利用半监督学习中的SVM进行多谱排序。相对于第一种方法，上述第二种多谱排序方法的灵敏度更高，速度更快，但是它存在不稳定的缺点，即在不同的数据集上灵敏度波动较大。

因此，当前迫切需要一种灵敏度高且性能稳定的多谱排序解决方案。

发明内容

本发明的任务是提供一种灵敏度高且性能稳定的多普排序解决方案。

本发明提供了一种基于半监督学习的交联质谱多谱排序方法，包括下列步骤：

1)对于参与多谱排序的多个谱图，分别对每个谱图进行单谱匹配和排序，得到对应的最优的交联二肽单谱匹配结果；提取当前每一个肽谱匹配结果的多谱匹配特征向量，构成所述多谱匹配特征向量的特征包括：动态特征和静态特征；所述动态特征包括SVM分数、母离子误差比例特征和修饰比例特征；

2)根据预设的迭代条件判断是否进行迭代，如果是，进入步骤3)，如果否，则输出基于当前SVM分数的多谱排序结果；

3)选择训练样本构建本次迭代的训练样本库，其中，在所得到的交联二肽匹配结果中，取FDR在预设的FDR阈值以内的属于正样本的结果构建正样本库，取所有负样本的结果构建反样本库；然后，基于所选择的训练样本更新更新各个交联二肽单谱匹配结果的多谱匹配特征向量；

4)利用本次迭代的训练样本库训练SVM分类器；

5)用本轮训练后的SVM分类器对所有交联二肽结果进行重打分，然后返回步骤2)，不断迭代直至输出基于当前SVM分数的多谱排序结果。

其中，所述步骤1)中，所述最优的交联二肽单谱匹配结果是单谱内排序所使用的细打分分数第一名的匹配结果。

其中，所述步骤1)中，所述交联二肽单谱匹配结果包括包括谱图标识、肽序列、交联位点和对应的分数。

其中，所述步骤1)中，所述静态特征包括：单谱内排序所使用的细打分分数、α肽匹配谱峰强度比例、β肽匹配谱峰强度比例、α最长序列标签长度比例、β最长序列标签长度比例、较短子肽段长度、分数差和一级谱的母离子可信度。

其中，所述步骤2)中，所述预设的迭代条件是迭代次数是否超过阈值，或者是当前的特征差是否小于预设的阈值。

其中，所述步骤3)包括下列子步骤：

31)计算各个交联二肽匹配结果的FDR；

32)取FDR在预设的FDR阈值以内的属于正样本的结果构建正样本库，取所有负样本的结果构建反样本库，进而获得本次迭代的训练样本库；

33)基于本次迭代的训练样本库，更新各个多谱排序特征向量中的母离子误差比例特征和修饰比例特征。

其中，所述步骤31)包括下列子步骤：

311)将各张谱图对应的交联二肽匹配结果按得分排序；

312)按交联二肽的假发现率控制公式计算假发现率FDR＝(U-F)/T；其中U表示一条肽段来自反样本库，另一条来自正样本库的肽谱匹配结果的个数，F表示两条肽段均来自反样本库的肽谱匹配结果的个数，T表示两条肽段均来自正样本库的肽谱匹配结果的个数；

313)对交联二肽匹配结果序列的FDR值进行平滑处理，使得FDR值保持单调递增。

其中，所述步骤33)包括下列子步骤：

331)计算每个正样本母离子误差；

332)将母离子误差划分为若干个小区间中，并计算每个小区间中正样本的个数；

333)采用向下累计的方式计算各区间的累计频率；

334)统计各种可变修饰在正样本中的比例；

335)将母离子误差对应区间的累积频率作为母离子误差比例特征的值，将所含可变修饰对应的修饰比例的乘积作为修饰比例特征。

与现有技术相比，本发明具有下列技术效果：

1、本发明提供了一种灵敏度高且性能稳定的用于交联蛋白鉴定的多普排序方法。

2、本发明能够稳定提升了大规模数据库下交联蛋白鉴定的灵敏度。

3、本发明特别适合用于基于质谱技术的规模化交联谱图的鉴定。

4、实验表明，在相同的数据集下，本发明的鉴定灵敏度高于pLink 1.0、pLink 2.0和Kojak软件。

附图说明

以下，结合附图来详细说明本发明的实施例，其中：

图1示出了根据本发明一个实施例所提供的基于半监督学习的多谱排序方法的流程图。

具体实施方式

图1示出了根据本发明一个实施例所提供的基于半监督学习的多谱排序方法的流程图，参考图1，该方法包括下列步骤：

步骤1：分别对每个谱图进行单谱匹配和排序，得到最优(通常是得分第一名)的交联二肽段单谱匹配结果。本实施例中，对于样品同一批次打谱所得到的所有谱图，分别获取每个谱图匹配到的分数第一的候选肽段。一张谱图和候选交联二肽进行匹配，构成一个肽谱匹配。第一名结果指的是打分第一名的肽谱匹配结果，这个肽谱匹配结果包括谱图标识、肽序列、交联位点、分数等信息。上述匹配打分可以使用常用的单谱内排序的细打分分数，本实施例中，采用KSDP分数。

这样，在匹配完成后，每个谱图都有一个第一名的肽谱匹配结果，对于这个肽谱匹配结果，提取11维特征向量作为其多谱匹配特征向量。11维特征向量中的每个特征均进行归一化，以便于后续的数据处理。表1示出了上述交联二肽的肽谱匹配结果的11个特征。

表1

在具体实现上，可以先生成交联二肽的理论碎片离子,然后将理论碎片离子与实验谱图中的实验谱峰进行匹配，对于每一个交联二肽的搜索结果(即匹配结果)，计算出上述11个特征。其中SVM分数在第一次迭代时没有上一轮的预测分，因此可以直接使用预设的初始值替代。

在一个实施例中，上述11维特征中，SVM-Score，PrecursorError，ModRatio是动态特征，在迭代过程中不断变化，另外8维特征是静态特征，始终保持不变。特征归一化操作在迭代前进行，只对静态特征(8维)进行归一化。归一化是将同一种特征的(不同结果的)不同值进行归一化。

根据一个实施例，归一化方法如下：在所有第一名的交联二肽结果中，获取每个特征的最大值和最小值；利用最小最大值归一化方法将特征值转化到0到1之间。公式如(1)所示。

其中X为原始特征值，X′为归一化之后的特征值，X_min为当前特征值对应的最小特征值，X_max为当前特征值对应的最大特征值。

步骤2：根据预设的迭代条件判断是否进行迭代，如果是，进入步骤3，如果否，则输出基于当前SVM分数的多谱排序结果。预设的迭代条件可以是迭代次数是否超过阈值，也可以是当前的特征差是否小于预设的阈值，还可以是特征差连续数据迭代均小于预设的阈值。

步骤3：选择训练样本，基于所选择的训练样本更新各个多谱排序特征向量中的母离子误差比例特征和修饰比例特征。

在一个实施例中，本步骤3包括下列子步骤：

步骤31，将各张谱图对应的肽谱匹配结果按得分排序。这里的得分指的是SVM分数，但存在特例，在第一次迭代时尚没有SVM分数，因此使用KSDP分数进行排序。

步骤32，按交联二肽的假发现率控制公式计算假发现率FDR。FDR＝(U-F)/T。U表示一条肽段来自反库，另一条来自正库的肽谱匹配结果的个数，F表示两条肽段均来自反库的肽谱匹配结果的个数，T表示两条肽段均来自正库的肽谱匹配结果的个数。

步骤33，基于步骤31的排序，从后往前扫描所有结果，并记录最低的FDR值。若当前FDR值小于最低FDR值，则更新最低FDR值，将当前FDR值作为最低值，当前结果的q-value为更新后的最低FDR值。若当前FDR值大于最低FDR值，则最低FDR值保持不变，当前结果的q-value为最低FDR值。q-value是用于对FDR值进行平滑处理的中间变量。在上述子步骤33中，用q-value替代FDR值，并保持q-value值单调递增。这样能够抑制FDR计算公式直接获得的FDR值的抖动，进而减小误差。

步骤34，取FDR在1％内的正库结果作为正样本，所有反库结果作为负样本。此处1％为预设的FDR阈值。

步骤35，计算每个正样本母离子误差。

步骤36，将母离子误差划分为若干个小区间中，并计算每个小区间中正样本的个数。

步骤37，采用向下累计的方式计算各区间的累计频率。对数据进行分组后，分布在各组内的数据个数称为频数，各组频数与全部频数之和的比值称为频率。累积频率指的是某一数值以下或某一数值以上的频率之和。向下累积表示从变量值大的一方向变量值小的一方累加。频率的最终累积值为100％。

步骤38，统计各种可变修饰在正样本中的比例。

上述步骤中，取FDR在1％内的正库结果作为正样本，所有反库结果作为负样本。样本的选取不受累积频率的影响，但是反过来，累计频率受正样本分布的影响。其中，步骤38是为了后面计算ModRatio特征。一般认为正样本中比例高的修饰更加可信，因此可通过迭代使得排序结果更加可信。

步骤39，将母离子误差对应区间的累积频率作为母离子误差比例特征的值。对于每一个交联二肽匹配结果，将所含可变修饰对应的修饰比例的乘积作为修饰比例(ModRatio)特征。

步骤4：在特征向量更新完毕后，利用所选取的训练样本库(包括正库和反库)，训练线性SVM分类器。本实施例中SVM分类器可参考文献：Kall,L.,J.D.Canterbury,J.Weston,W.S.Noble and M.J.MacCoss,Semi-supervised learning for peptideidentification from shotgun proteomics datasets.Nat Methods,2007.4(11):p.923-5.，对SVM分类器的具体细节和训练方法不再赘述。

步骤5：用本轮训练后的SVM分类器对所有交联二肽结果进行重打分。这里的所有交联二肽结果指的是步骤1中得出的所有谱图各自对应的第一名的交联二肽结果。返回步骤2，判断是否进入下一轮迭代。

上述步骤2至5不断循环，即可不断修正多谱排序结果。

一般用户设定的FDR值是固定的。本实施例能够保证在FDR相同的情况下，获得更多的肽谱匹配结果(多谱排序结果就是最终给用户的肽谱匹配结果)。在多谱排序结果中，排在前面的结果表示高可信，它们被报告给用户。常用FDR来估计错误率。比如方法1使得前100个结果中有5个是错误的(错误率为5％)，前200个结果中有20个错误的，方法2使得前100个结果中只有2个是错误的，前200个结果中只有10个是错误的。那么就认为方法2比方法1好，因为在相同的FDR控制下(5％)，能够鉴定到更多的结果(100VS 200)。本实施例的鉴定结果基于SVM分数由高到低排列，越前面的错误率越少。

在实际测试中，相对于现有的用E-value做后处理的方案，在相同的FDR的前提下，本实施例的方法在人库上获得了95.3％的鉴定率，而现有的E-value方案的鉴定率只有74.4％。可以看出，在相同的FDR的前提下，本发明能够显著提升谱图的鉴定率。

上述实施例中，所采用的各个特征的定义可以参考文献：Trnka,M.J.,P.R.Baker,P.J.Robinson,A.L.Burlingame and R.J.Chalkley,Matching cross-linked peptidespectra:only as good as the worse identification.Mol Cell Proteomics,2014.13(2):p.420-34.；以及Wang,J.,V.G.Anania,J.Knott,J.Rush,J.R.Lill,P.E.Bourne andN.Bandeira,Combinatorial approach for large-scale identification of linkedpeptides from tandem mass spectrometry spectra.Mol Cell Proteomics,2014.13(4):p.1128-36.。

进一步地，在一个优选实施例中，还给出了11维特征的具体计算方法。下面分别列出。

(1)SVM分数(SVM-Score)

第一轮迭代时，不存在上一轮的预测分，所有样本的SVM-Score特征值均为零。从第二轮迭代开始，SVM-Score为上一轮的预测分，也是当前得分。每一轮迭代结束后，会产生新的SVM-Score，并将新的SVM-Score作为当前得分。SVM-Score综合了其它十个特征，对正负样本的区分能力是所有特征中最强的。

(2)KSDP分数(KSDP)

KSDP是搜索阶段采用的细打分分数，主要用来进行单谱内排序。KSDP综合考虑了理论离子匹配个数，匹配离子的连续性，共存性，匹配谱峰强度等因素，对正负样本的区分能力较强。KSDP特征值在迭代过程中保持不变。

(3)α肽匹配谱峰强度比例和β肽匹配谱峰强度比例(AlphaInt and BetaInt)

AlphaInt表示α肽匹配谱峰强度与谱图中总强度之比，而BetaInt表示β肽匹配谱峰强度与谱图中总强度之比。其中谱峰强度均指平滑后的谱峰强度，其强度等于原始谱峰强度开根号。交联二肽的碎片离子种类多，按碎裂次数可以将其划分为一次碎裂和二次碎裂。一次碎裂形成的离子只有一个断裂位点，如果断裂发生在α肽上，就认为这是α肽的理论碎片离子，而如果断裂发生在β肽上，就认为这是β肽的理论碎片离子。二次碎裂形成的离子有两个断裂位点，如果两次断裂都发生在α肽上，就认为这是α肽的理论碎片离子，而如果两次断裂都发生在β肽上，就认为这是β肽的理论碎片离子。如果一次断裂发生在α肽上，另一次断裂发生在β肽上，就认为这是双肽的理论碎片离子。因此，理论碎片离子按碎裂发生的肽段划分，可以分为三类，依次为α肽离子，β肽离子，双肽离子。

如果实验谱峰与α肽离子匹配上了，就将这根实验谱峰的强度加在α肽上。如果实验谱峰与β肽离子匹配上了，就将这根实验谱峰的强度加在β肽上。如果实验谱峰与双肽离子匹配上了，就将这根实验谱峰的强度的一半加在α肽上，另一半加在β肽上。如果一根实验谱峰同时被多种离子匹配上了，那么只有碎片离子误差绝对值最小的匹配被保留下来。

(4)α最长序列标签长度与α肽段长度之比和β最长序列标签长度与β肽段长度之比(AlphaTag and BetaTag)

序列标签长度指的是连续匹配上离子的断裂位点个数。在计算序列标签时，仅考虑b，y离子。比如α肽的长度为10，α肽上连续匹配上的b离子的断裂位点个数最多为3个，连续匹配上的y离子的断裂位点个数最多为5个，那么α最长序列标签长度为5，AlphaTag为0.5。特别地，若α肽序列与β肽序列完全一样，则将AlphaTag和BetaTag分别乘以一个衰减系数μ，μ的值在0和1之间。

(5)较短子肽段长度(LogLen)

LogLen表示将较短的子肽段的长度取自然对数。比如α肽段长度为10，β肽段长度为5，那么LogLen特征的值为In(5)≈1.61。

(6)分数差(DiffScore)

DiffScore特征量化了第一名结果比第二名结果的匹配情况好的程度。DiffScore特征的计算公式为：

其中Score1表示第一名结果的细打分分数，Score2表示第二名结果的细打分分数。DiffScore特征对正负样本的区分能力较强。

(7)pParse编号(pParseNum)

pParse是pFind团队开发的一款软件，其功能是母离子质量与电荷的校准、混合谱的导出。pParse导出的二级谱中每张谱的标题格式为“谱图名称.scan号.scan号.电荷数目.pParse编号.dta”。比如，一张二级谱的标题为“UTP_B.2569.2569.3.2.dta”，表示谱图的名称是UTP_B，scan号为2569，电荷数目为3，pParse编号为2。scan号与保留时间有关，pParse编号与母离子质量、电荷有关。一级谱的母离子质量和电荷越可信，pParse编号值越小。pParse编号的取值范围为0～6，0表示可信度最高，6表示可信度最低。scan号相同但pParse编号不同的谱图的碎片离子信息完全相同，不同之处在于母离子信息，即母离子质量和电荷。pParseNum从二级谱的标题中提取，利用了一级谱的信息。

(8)母离子误差比例(PrecursorError)

PrecursorError这个特征的计算分为四步。第一步，计算每个正样本母离子误差。第二步，将母离子误差划分为若干个小区间中，并计算每个小区间中正样本的个数。第三步，采用向下累计的方式计算各区间的累计频率。对数据进行分组后，分布在各组内的数据个数称为频数，各组频数与全部频数之和的比值称为频率。累积频率指的是某一数值以下或某一数值以上的频率之和。向下累积表示从变量值大的一方向变量值小的一方累加。频率的最终累积值为100％。第四步，将母离子误差对应区间的累积频率作为PrecursorError的值。

(9)修饰比例(ModRatio)

ModRatio表示肽段所含的各种修饰在正样本中的比例的乘积。比如α肽上含2个氧化修饰Oxidation[M]，β肽上含1个磷酸化修饰Phospho[S]，Oxidation[M]在正样本中的比例为0.1，Phospho[S]在正样本中的比例为0.05，那么其ModRatio的值为0.1*0.1*0.05，即0.0005。

最后应说明的是，以上实施例仅用以描述本发明的技术方案而不是对本技术方法进行限制，本发明在应用上可以延伸为其它的修改、变化、应用和实施例，并且因此认为所有这样的修改、变化、应用、实施例都在本发明的精神和教导范围内。

Claims (8)

1.一种基于半监督学习的交联质谱多谱排序方法，其特征在于，包括下列步骤：

1)对于参与多谱排序的多个谱图，分别对每个谱图进行单谱匹配和排序，得到对应的最优的交联二肽单谱匹配结果；提取当前每一个肽谱匹配结果的多谱匹配特征向量，构成所述多谱匹配特征向量的特征包括：动态特征和静态特征；所述动态特征包括SVM分数、母离子误差比例特征和修饰比例特征；

2)选择训练样本构建本次迭代的训练样本库，其中，在所得到的交联二肽匹配结果中，取FDR在预设的FDR阈值以内的属于正样本的结果构建正样本库，取所有负样本的结果构建反样本库；然后，基于所选择的训练样本更新更新各个交联二肽单谱匹配结果的多谱匹配特征向量；

3)利用本次迭代的训练样本库训练SVM分类器；

4)用本轮训练后的SVM分类器对所有交联二肽结果进行重打分；

5)根据预设的迭代条件判断是否继续进行迭代，如果是，则返回步骤2)进行新的迭代，如果否，则输出基于当前SVM分数的多谱排序结果。

2.根据权利要求1所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤1)中，所述最优的交联二肽单谱匹配结果是单谱内排序所使用的细打分分数第一名的匹配结果。

3.根据权利要求2所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤1)中，所述交联二肽单谱匹配结果包括谱图标识、肽序列、交联位点和对应的分数。

4.根据权利要求3所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤1)中，所述静态特征包括：单谱内排序所使用的细打分分数、α肽匹配谱峰强度比例、β肽匹配谱峰强度比例、α最长序列标签长度比例、β最长序列标签长度比例、较短子肽段长度、分数差和一级谱的母离子可信度。

5.根据权利要求1所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤5)中，所述预设的迭代条件是迭代次数是否超过阈值，或者是当前的特征差是否小于预设的阈值。

6.根据权利要求1所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤2)包括下列子步骤：

21)计算各个交联二肽匹配结果的FDR；

22)取FDR在预设的FDR阈值以内的属于正样本的结果构建正样本库，取所有负样本的结果构建反样本库，进而获得本次迭代的训练样本库；

23)基于本次迭代的训练样本库，更新各个多谱排序特征向量中的母离子误差比例特征和修饰比例特征。

7.根据权利要求6所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤21)包括下列子步骤：

211)将各张谱图对应的交联二肽匹配结果按得分排序；

212)按交联二肽的假发现率控制公式计算假发现率FDR＝(U-F)/T；其中U表示一条肽段来自反样本库，另一条来自正样本库的肽谱匹配结果的个数，F表示两条肽段均来自反样本库的肽谱匹配结果的个数，T表示两条肽段均来自正样本库的肽谱匹配结果的个数；

213)对交联二肽匹配结果序列的FDR值进行平滑处理，使得FDR值保持单调递增。

8.根据权利要求6所述的基于半监督学习的交联质谱多谱排序方法，其特征在于，所述步骤23)包括下列子步骤：

231)计算每个正样本母离子误差；

232)将母离子误差划分为若干个小区间中，并计算每个小区间中正样本的个数；

233)采用向下累计的方式计算各区间的累计频率；

234)统计各种可变修饰在正样本中的比例；

235)将母离子误差对应区间的累积频率作为母离子误差比例特征的值，将所含可变修饰对应的修饰比例的乘积作为修饰比例特征。