CN103776891A - 一种检测差异表达蛋白质的方法 - Google Patents

Abstract

本发明涉及一种检测差异表达蛋白质的方法，面向定量蛋白质组学中的基于一级谱图信息的标记和非标记的相对定量数据分析，包括肽谱匹配、可信度评价、肽段信号提取、肽段比值计算、蛋白质比值计算、统计学分析，根据某蛋白质在两种或多种样品中对应的质谱信号强度比值判断其是否是差异表达蛋白质。对于近百GB的规模的质谱实验采集的数据，快速地自动化分析，对不同蛋白质在质谱仪中的信号尽可能精准地提取蛋白质信号；从统计学意义上确定蛋白质差异表达，并对结果的准确性进行评价。

Description

一种检测差异表达蛋白质的方法

技术领域

本发明涉及定量蛋白质组学领域，特别涉及一种检测差异表达蛋白质的方法。

背景技术

定量蛋白质组学主要通过生物质谱技术研究复杂生物样品中蛋白质的表达情况，其中一项重要的研究目标是在不同生物样品中检测差异表达蛋白质。

对于某些重要的蛋白质，如与癌症相关的核心岩藻糖蛋白质，在癌症患者体内的表达量要远远高于正常人。该差异表达的核心岩藻糖蛋白质可以作为肝癌早期诊断标志物。在其他诸如生物信号传导、细胞衰老以及翻译后修饰等生命过程中，差异表达蛋白质都起到了主导作用。据报道，2011年，美国药品研究与制造商协会（PhRMA）的成员公司总共投入495亿美元研究经费，其中绝大部分经费用于寻找疾病标志物，即检测和确认患者体内与正常人体内差异表达的蛋白质。这些差异表达蛋白质与疾病有直接关系，可以测量相关人群这些蛋白质的表达量，用于诊断是否患病。

传统的生物化学方法，如蛋白质印迹法，一次实验需要耗费数天的时间，却只能检测一个或数个高丰度蛋白质在两个样品中表达量的差异。生物体内含有数以千计的蛋白质，使用这类方法一一进行检测需要消耗大量的人力和物力。而生物质谱技术具有灵敏度好、动态线性范围大、通量高等优点，可以一次性分析数千蛋白质。德国和瑞士的科学家报道说使用生物质谱技术数个小时内即可检测到大肠杆菌体内的5000个蛋白质，占大肠杆菌全部蛋白质的90%以上；可以检测到人体内1万个以上蛋白质，并且数目还在增加。

生物质谱技术的核心功能是将蛋白质在生物体中的表达量的信息“数字化”，即对某复杂生物样品采集质谱数据，而后使用计算技术分析这些数据，还原蛋白质的表达量信息。2009年，当今国际上蛋白质组学3大主要引领者之一，瑞士联邦理工学院分子系统生物学研究所的Ruedi Aebersold教授在期刊《Nature Methods》上发表了一篇文章，认为使用计算技术分析质谱数据是整个定量蛋白质组学中最关键、最有挑战性的研究问题。这里面的难点包括：（1）一次质谱实验采集的数据可以达到近百GB的规模，需要有方法快速地自动化分析；（2）不同蛋白质在质谱仪中的信号可能会相互重叠、相互干扰，需要专门的算法去干扰，尽可能精准地提取蛋白质信号；（3）需要从统计学意义上确定哪些蛋白质是差异表达的，并对结果的准确性进行评价。

目前已有的面向蛋白质质谱数据的分析方法，最主要的问题是对信号重叠、信号干扰的处理能力较弱，得出的定性、定量结果中有相当比例都是不准确的，并且几乎没有算法对蛋白质的定量结果的准确性进行评价。这导致使用该方法检测的差异表达蛋白质无法完全信任，需要后续使用多种传统的生物化学手段再进行一一验证，还要耗费了大量的人力和物力，且增长了研究周期。该问题一直是限制标志物检测研究取得突破性进展的瓶颈。

目前蛋白质定量的主流技术是，面向定量蛋白质组学中的基于一级谱图信息的标记和非标记的相对定量数据分析，广泛应用于生物标志物发现、临床诊断、生物信号传导过程以及翻译后修饰研究等领域。

发明内容

为了解决上述问题，本发明的目的在于提出一种检测差异表达蛋白质的方法，面向定量蛋白质组学中的基于一级谱图信息的标记和非标记的相对定量数据分析，根据某一蛋白质在两种或多种生物样品中对应的质谱信号强度比值判断其是否是差异表达蛋白质。

本发明公开了一种检测差异表达蛋白质的方法，包括：

步骤1，对质谱数据进行预处理，用于将该蛋白质的原始二进制质谱数据转换为文本格式，并建立索引；

步骤2，对二级谱图进行肽谱匹配，确定样品中含有的肽段，用于对二级谱图与蛋白质数据库中的记录的肽段进行匹配打分，取高可信的匹配结果；或者直接从二级谱图推测肽段序列，取高可信的结果；

步骤3，提取每个肽段在两种生物样品和多种生物样品中的信号，以多个同位素曲线的形式进行表示；

步骤4，将相同肽段在不同样品中的信号之间建立对应关系，计算其表达量差异的肽段比值和置信区间；

步骤5，将肽段比值归并为蛋白质比值，并给出蛋白质比值的置信区间；

步骤6，确定差异表达的蛋白质。

所述的检测差异表达蛋白质的方法，肽段信号以多个同位素色谱曲线的形式进行表示。

所述的检测差异表达蛋白质的方法，基于干扰最小的同位素色谱曲线计算肽段比值，具体方法为局部最小一乘法，并计算比值的置信区间。

所述的检测差异表达蛋白质的方法，将肽段比值归并为蛋白质比值采用核密度估计方法。

所述的检测差异表达蛋白质的方法，步骤2还包括：

步骤11，对每张二级谱图进行处理，只保留强度最大的前200个谱峰；

步骤12，对输入的每张二级谱图，在蛋白质数据库中寻找最相似的肽段。

所述的检测差异表达蛋白质的方法，步骤2还包括：

步骤21，对肽谱匹配按照打分从高到低进行排序；

步骤22，控制肽谱匹配的假发现率。

所述的检测差异表达蛋白质的方法，步骤3还包括：

步骤31，读取肽段排序列表，并对一级谱图进行预处理，同时设定必要的数据分析参数，包括要分析的质谱数据、数据的类型、所属物种的蛋白质数据库；

步骤32，提取肽段信号，以多个同位素色谱曲线的形式表示。

步骤321，对每个肽谱匹配，计算每个肽段的理论同位素分布；

步骤322，根据理论同位素分布，在鉴定到该肽段的二级谱图前后2分钟保留时间范围内的一级谱图上确定实际同位素峰；

对于这个范围内的某张一级谱图，如果某理论同位素峰的质荷比正负10ppm范围内有谱峰，则记录下来；如果10ppm范围内有多根谱峰，那么每次选其中一个谱峰，和理论同位素分布计算余弦夹角，取值最高的那个组合为实际的同位素峰；

步骤323，对于每一个实际的同位素峰，沿保留时间把它们连成一条曲线，表示该肽段从没有信号到有信号再到信号消失的过程。

先根据鉴定到肽段的二级谱图的扫描号，找到距离该二级谱图最近的一级谱图。而后以该一级谱图所在扫描号为基准，在这条曲线上向前寻找起始点，向后寻找终止点，当曲线的上某点的强度低于最高曲线强度的10%时，停止；该曲线在原来的基础上就小了一些；再以此曲线为基础，寻找起始点和终止点附近的极值点，将该极值点与起始点或终止点之间的部分删除；

步骤324，如果是标记定量实验，那么肽段有轻、重标两种形式。

所述的检测差异表达蛋白质的方法，步骤4还包括：计算肽段在两个样品中对应信号的强度差异，以比值的形式表示；取干扰最小的曲线计算比值。

所述的检测差异表达蛋白质的方法，步骤4还包括：

步骤41，对于一个肽段，用向量

和

表示计算的肽段同位素标记的轻、重标形式的理论同位素丰度；n和m分别表示肽段轻、重标形式的同位素峰数目；l₁＝(l_1,1,...,l_1,k)表示肽段轻标形式单同位素色谱曲线，k是色谱曲线跨越的一级谱图的数目；l₂＝(l_2,1,...,l_2,k)表示肽段轻标形式第一同位素色谱曲线，依次类推，l_n＝(l_n,1,...,l_n,k)表示肽段轻标形式第n同位素色谱曲线；类似的，h₁＝(h_1,1,...,h_1,k)表示肽段重标形式单同位素色谱曲线，h_m＝(l_m,1,...,h_m,k)肽段重标形式第m同位素色谱曲线；

步骤42，对上一步的各同位素曲线进行归一化，以确保肽段轻、重标记间各取任一同位素曲线都可用于计算肽段比值；

归一化方法如下： $L_1=l_1\·\frac{{\mathrm{\Σt}}_i^l}{t_1^l}=(l_{\mathrm{1,1}},...,l_{1,k}),$ $L_2=l_2\·\frac{{\mathrm{\Σt}}_i^l}{t_2^l}=(l_{\mathrm{2,1}},...,l_{2,k}),$ ……， $L_n=l_n\·\frac{{\mathrm{\Σt}}_i^l}{t_n^l}=(l_{n,1},...,l_{n,k});$ $H_1=h_1\·\frac{{\mathrm{\Σt}}_j^h}{t_1^h}=(h_{\mathrm{1,1}},...,h_{1,k}),$ $H_2=h_2\·\frac{{\mathrm{\Σt}}_j^h}{t_2^h}=(h_{2,1},...,h_{2,k}),$ ……，

L₁，……L_i，……，L_n表示归一化后肽段轻标形式的各同位素曲线，H₁，……，H_j，……，H_m表示归一化后肽段重标形式的各同位素曲线；

步骤43，计算m×n个比值，分别是使用H_j和L_i计算比值，其中i＝1,...,n，j＝1,...,m；

所述的检测差异表达蛋白质的方法，步骤43中计算比值包括如下步骤：

给定X=(x₁,…x_k)，Y=(y₁,…,y_k)表示分别从肽段轻标形式和重标形式中各取了一条同位素曲线。

步骤431）以两个曲线的中心为基础，取三个点为局部曲线：中点左边的点，中点，中点右边的点；即取

计算

和

的余弦夹角t；

步骤432）

和

分别向两端延伸，即取 $\tilde{X}=(x_{\frac{k}{2}-2},x_{\frac{k}{2}-1},x_{\frac{k}{2}},x_{\frac{k}{2}+1},x_{\frac{k}{2}+2}),$

$\tilde{Y}=(y_{\frac{k}{2}-2},y_{\frac{k}{2}-1},y_{\frac{k}{2}},y_{\frac{k}{2}+1},y_{\frac{k}{2}+2}),$ 再计算

和的余弦夹角t；

步骤433）如果t大于一个阈值（根据经验，设定为0.8），就继续（2）中的步骤，否则停止；

步骤434）延伸停止后，使用过原点的最小一乘法

计算

和

的比值a。

步骤435）计算比值的置信区间，计算

当样本大小n→∞时，有

那么，可得a的1-α的区间估计是

其中u_α/2是标准正态分布N(0,1)的“上α/2分位点”；令 $\hat{f}\left(0\right)=\frac{1}{\mathrm{nw}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}K\left(\frac{{\hat{e}}_i}{w}\right),$ 选择核函数 $K\left(x\right)=\frac{3}{4}(1-x^2)I_{\left|x\right|\≤1},$ 窗宽 $w=1.66n^{-1/5}\widehat{\mathrm{\σ}},$ 其中 $\widehat{\mathrm{\σ}}=\sqrt{s_n^2}.$

步骤44，计算的m×n个比值分别对应m×n个置信区间，取置信区间最小的那个作为肽段的比值，该置信区间作为肽段比值的置信区间。

所述的检测差异表达蛋白质的方法，步骤5还包括：

步骤51，使用核密度估计的方法推断蛋白质的比值；

步骤52，使用了高斯核，并假设每个肽段的比值服从一个高斯分布：

$f_{\mathrm{pep}}\left(x\right)=\frac{1}{N}\·\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_i}{e}^{\frac{{(x-r_i)}^2}{2{{\mathrm{\σ}}_i}^2}}$

其中，r_i是第i个肽段比值，σ_i是其标准差，N是其所对应的蛋白质所鉴定或定量的肽段的总数；

步骤53，定义F(x)＝∑f_pep(x)，根据F(x)拟合一个高斯曲线，其形式如下：

$f_{\mathrm{pro}}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{pro}}}{e}^{\frac{{(x-R)}^2}{2{\mathrm{\σ}}_{\mathrm{pro}}^2}}$

其中，R就是这个蛋白质的比值，σ_pro是其标准差,根据这个标准差，就可以计算R的置信区间。

所述的检测差异表达蛋白质的方法，步骤6还具体包括如下步骤：

步骤61，对所有比值取log2变换，并通过加或减一个常数，使其中值为0；

步骤62，拟合所有比值的分布；

步骤63，根据拟合的分布计算每个比值的p-value；

步骤64，计算假发现率，报告假发现率小于1%的那些蛋白质比值作为可信的显著差异表达蛋白质，并以报表的形式展示给用户。

本发明的主要内容和特点包括：

（1）支持多种实验技术，包括定量蛋白质组学中的基于一级谱图信息的标记，如细胞培养中稳定同位素标记氨基酸（Stable Isotope Labeling byAmino acids in Cell culture，SILAC）、¹⁵N体内代谢、¹⁸O化学标记等技术以及非标记技术。

（2）数据分析全流程自动化，定性、定量紧密结合，可以一次性分析近万个蛋白质，加快了研究进度，缩短了研究周期。

（3）设计了一整套方法处理蛋白质和肽段信号重叠问题，尽可能排除干扰，定量结果准确性好。

（4）德国马普学会生物化学研究所Matthias Mann实验室开发了一种类似的方法MaxQuant，该方法发表在领域内权威期刊《Nature Biotechnology》上。本发明的方法支持的实验技术比MaxQuant多至少2种，同硬件条件下分析速度提高接近10倍，相同数据集下定量结果准确性（以标准差为准）提高约26%，报告的差异表达蛋白质的错误率减小了3倍。

(5)本发明在北京生命科学研究所等国内多家单位进行了试用，从统计学意义上确定蛋白质差异表达，并对结果的准确性进行评价，得到的结果明显优于各单位之前使用的方法，解决了实际问题，取得了预期的效果。

附图说明

图1显示了本发明的实施过程；

图2显示了本发明的六大核心模块；

图3展示了一个受干扰的肽段信号。

具体实施方式

本发明是通过分析质谱仪采集的蛋白质在不同样品中的数据，来检测哪些是差异表达的，具体来说，面向定量蛋白质组学中的基于一级谱图信息的标记和非标记的相对定量数据分析，根据某蛋白质在两种或多种生物样品中对应的质谱信号强度比值判断其是否是差异表达蛋白质。

如图1所示为本发明的完整实施过程，也是标记定量的实验流程，主要包括：

（1）取要进行差异蛋白质检测的两份不同生物样品，样品1为某生物正常组织，并进行了稳定同位素标记处理；样品2是该生物病变组织。

对其中一份样品使用细胞培养中稳定同位素标记氨基酸（Stable IsotopeLabeling by Amino acids in Cell culture，SILAC）或¹⁵N体内代谢等技术进行处理，而后，将两份样品混合，混合样品中含有近万个蛋白质。也可以使用非标记定量的技术，在后续分别采集两份样品的质谱数据，直接分析。目前的主流技术是将蛋白质酶切为肽段，所以，后续质谱仪采集的是肽段的信号，在数据分析阶段将进行肽段到蛋白质的归并。

（2）使用质谱仪采集上述混合样品的数据。一般来说，质谱仪采集的数据由两部分组成：一部分是一级谱图，记录了肽段的强度信息，用于肽段的定量分析；另一部分是二级谱图，记录了肽段的序列信息，用于肽段的定性分析。

对质谱数据进行预处理，用于是将该蛋白质的原始二进制质谱数据转换为文本格式，并建立索引，方便后续处理；

确定样品中含有的肽段，用于对二级谱图与蛋白质数据库中的记录的肽段进行匹配打分，取高可信的匹配结果；或者直接从二级谱图推测肽段序列，取高可信的结果；

（3）分析相关数据，确定差异表达蛋白质。这一步骤是本发明有所创新之处，也是目前生物质谱技术的发展瓶颈。

分别提取肽段在样品1和样品2中的信号，以多个同位素曲线的形式进行表示；

将相同肽段在不同样品中的信号之间建立对应关系，以信号强度的比值表示表达量的差异，计算该比值和对应的置信区间；

将肽段比值归并为蛋白质比值，并给出蛋白质比值的置信区间；大部分蛋白质的比值接近1:1，只有与病变有直接关系的蛋白质比值差异较大，即本发明要检测的差异表达蛋白质；

确定差异表达的蛋白质。

图2显示了本发明的六大核心模块（对应以下六个主要步骤）：1肽谱匹配、2可信度评价、3肽段信号提取、4肽段比值计算、5蛋白质比值计算、6统计学分析。这六个模块协同工作，检测差异表达蛋白质。

具体实施步骤有：

步骤1：对二级谱图进行肽谱匹配，确定样品中含有的肽段。输入n张二级谱图，就可以得到n个肽段，所以称为“肽谱匹配”。

步骤11：对每张二级谱图进行处理，只保留强度最大的前200个谱峰。

步骤12：对输入的每张二级谱图，在蛋白质数据库中寻找最相似的肽段。可以使用蛋白质数据库搜索软件，如pFind（Wang,L.H.,et al.,Rapid CommunMass Spectrom,2007；Fu,Y.,et al.,Bioinformatics,2004）完成该步骤。也可以使用从头测序软件，如pNovo（Chi,H.,et al.,J Journal ofProteome Research,2010），直接从二级谱图推测肽段序列。

步骤2：步骤1中得到的n个肽谱匹配可能会有错误匹配，本步骤需要从中挑取相对可信的那部分。

步骤21：对步骤1中得到n个肽谱匹配按照打分从高到低进行排序。

步骤22：使用哈佛大学生物医学院Gygi教授开发的诱饵库方法，控制肽谱匹配的假发现率（False Discovery Rate，FDR），一般取假发现率低于1%的部分，即挑选打分排在前m位的肽谱匹配（0<m<n），这部分的错误的数目不超过1%。这部分肽段输给步骤3用于提取定量信息。

步骤3：对每个肽段，分别提取其在两个样品中的信号。

步骤31：读入步骤2中的肽段列表，并对一级谱图进行预处理，同时设定必要的数据分析参数，包括要分析的质谱数据、数据的类型、所属物种的蛋白质数据库。

步骤32：提取肽段信号，以多个同位素色谱曲线的形式表示，如图3所示。

步骤321：对每个肽谱匹配，使用emass算法（Rockwood,A.L.,et al.,J Am Soc Mass Spectrom,2006）计算每个肽段的理论同位素分布。

步骤322：根据理论同位素分布，在鉴定到该肽段的二级谱图前后2分钟保留时间范围（该范围可由用户设定）内的一级谱图上确定实际同位素峰。对于这个范围内的某张一级谱图，如果某理论同位素峰的质荷比正负10ppm（该偏差可由用户设定）范围内有谱峰，则记录下来。如果10ppm范围内有多根谱峰，那么每次选其中一个谱峰，和理论同位素分布计算余弦夹角，取值最高的那个组合为实际的同位素峰。

步骤323：对于每一个实际的同位素峰，沿保留时间（时间从小到大）把它们连成一条曲线，表示该肽段从没有信号到有信号再到信号消失的过程。该曲线包含了真正的肽段的信号，但是，在它的起始和结束阶段，可能会有其他信号。所以，还需要在这个曲线上确定真正的肽段信号“起始点”和“终止点”，所用方法如下：先根据鉴定到肽段的二级谱图的扫描号，找到距离该二级谱图最近的一级谱图。而后以该一级谱图所在扫描号为基准，在这条曲线上向前（扫描号变小）寻找“起始点”，向后（扫描号变大）寻找“终止点”，当曲线的上某点的强度低于最高曲线强度的10%时，停止。该曲线在原来的基础上就小了一些。再以此曲线为基础，寻找起始点和终止点附近的极值点，将该极值点与起始点或终止点之间的部分删除。

步骤324：如果是标记定量实验，那么肽段有轻、重标两种形式。这两种形式的肽段信号在一个一级谱图文件中。为每个肽段的轻、重标两种形式的信号建立一一映射，称为一个“信号对”（signal pair）。

如果是非标记定量实验，那么肽段在两次实验中各有信号。为每个肽段在两次实验中的信号建立一一映射，称为一个“信号对”。

步骤4：计算肽段在两个样品中对应信号的强度差异，以比值的形式表示。在本发明中，肽段的信号是以多个同位素曲线的形式进行表示的。发明中包含了一个打分算法，可以判断各曲线受到的干扰情况，取干扰最小的曲线计算比值。

步骤41：对于一个肽段，用向量和

表示步骤3中计算的肽段同位素标记的轻、重标形式的理论同位素丰度。n和m分别表示肽段轻、重标形式的同位素峰数目。l₁＝(l_1,1,...,l_1,k)表示肽段轻标形式单同位素色谱曲线，k是色谱曲线跨越的一级谱图的数目；l₂＝(l_2,1,...,l_2,k)表示肽段轻标形式第一同位素色谱曲线，依次类推，l_n＝(l_n,1,...,l_n,k)表示肽段轻标形式第n同位素色谱曲线。类似的，h₁＝(h_1,1,...,h_1,k)表示肽段重标形式单同位素色谱曲线，h_m＝(l_m,1,...,h_m,k)肽段重标形式第m同位素色谱曲线。

步骤42：对上一步的各同位素曲线进行归一化，以确保肽段轻、重标记间各取任一同位素曲线都可用于计算肽段比值。归一化方法如下： $L_1=l_1\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_1^l}=(l_{\mathrm{1,1}},...,l_{1,k}),$ $L_2=l_2\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_2^l}=(l_{\mathrm{2,1}},...,l_{2,k}),$ ……， $L_n=l_n\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_n^l}=(l_{n,1},...,l_{n,k});$ $H_1=h_1\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_j^h}{t_1^h}=(h_{\mathrm{1,1}},...,h_{1,k}),$ $H_2=h_2\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_j^h}{t_2^h}=(h_{2,1},...,h_{2,k}),$ ……，

L₁，……L_i，……，L_n表示归一化后肽段轻标形式的各同位素曲线，H₁，……，H_j，……，H_m表示归一化后肽段重标形式的各同位素曲线。

步骤43：我们计算m×n个比值，分别是使用H_j和L_i计算比值，其中i＝1,...,n，j＝1,...,m。每个比值计算使用的方法如下：

给定X=(x₁,…x_k)，Y=(y₁,…,y_k)表示分别从肽段轻标形式和重标形式中各取了一条同位素曲线。

（1）以两个曲线的中心为基础，取三个点为局部曲线：中点左边的点，中点，中点右边的点；即取

计算

和的余弦夹角t；

（2）

和

分别向两端延伸，即取 $\tilde{X}=(x_{\frac{k}{2}-2},x_{\frac{k}{2}-1},x_{\frac{k}{2}},x_{\frac{k}{2}+1},x_{\frac{k}{2}+2}),$

$\tilde{Y}=(y_{\frac{k}{2}-2},y_{\frac{k}{2}-1},y_{\frac{k}{2}},y_{\frac{k}{2}+1},y_{\frac{k}{2}+2}),$ 再计算和

的余弦夹角t；

（3）如果t大于一个阈值（根据经验，设定为0.8），就继续（2）中的步骤，否则停止；

（4）延伸停止后，使用过原点的最小一乘法

计算和

的比值a。

（5）计算比值的置信区间。计算

当样本大小n→∞时，有

那么，可得a的1-α的区间估计是

其中u_α/2是标准正态分布N(0,1)的“上α/2分位点”。令 $\hat{f}\left(0\right)=\frac{1}{\mathrm{nw}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}K\left(\frac{{\hat{e}}_i}{w}\right),$ 选择核函数 $K\left(x\right)=\frac{3}{4}(1-x^2)I_{\left|x\right|\&le;1},$ 窗宽 $w=1.66n^{-1/5}\widehat{\mathrm{\&sigma;}},$ 其中

一般情况下，取α为0.0025，那么计算的是比值的“99.75%置信区间”，它表示每个比值最可能的取值范围。

步骤44：上步计算的m×n个比值分别对应m×n个置信区间，我们取置信区间最小的那个作为肽段的比值，该置信区间作为肽段比值的置信区间。一般情况下可以认为该比值是在肽段轻、重标形式中各取了干扰最小的同位素色谱曲线得到的。如图3所示，肽段的重标形式的第一、二、三同位素色谱曲线受到了干扰，但单同位素色谱曲线没有受到干扰（箭头所示）。该曲线和肽段轻标形式求得的比值置信区间最小。

步骤5：计算蛋白质的比值和置信区间。

使用核密度估计的方法推断蛋白质的比值。使用了高斯核，并假设每个肽段的比值服从一个高斯分布：

$f_{\mathrm{pep}}\left(x\right)=\frac{1}{N}\&CenterDot;\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_i}{e}^{\frac{{(x-r_i)}^2}{2{{\mathrm{\&sigma;}}_i}^2}}$

其中，r_i是第i个肽段比值，σ_i是其标准差，N是其所对应的蛋白质所鉴定或定量的肽段的总数。定义F(x)＝∑f_pep(x)，它是一个不规则的曲线。可以根据F(x)拟合一个高斯曲线，其形式如下：

$f_{\mathrm{pro}}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{pro}}}{e}^{\frac{{(x-R)}^2}{2{\mathrm{\&sigma;}}_{\mathrm{pro}}^2}}$

其中，R就是这个蛋白质的比值，σ_pro是其标准差。根据这个标准差，就可以计算R的置信区间。也就是说，可以认为R是服从正态分布的，R～N(R^*,σ_pro)，其中R^*是R的真实值，R的99.75%置信区间为为[R-3·σ_pro,R+3·σ_pro]。

步骤6，将一次实验得到的所有蛋白质的比值进行以下处理，最终得出显著差异的蛋白质。

步骤61，对所有比值取log2变换，并通过加或减一个常数，使其中值为0。

步骤62，拟合所有比值的分布，一般服从正态分布。

步骤63，根据拟合的分布计算每个比值的p-value。

步骤64，使用的是目前领域内比较公认的Benjamini-Hochberg(BH)方法计算假发现率（False Discovery Rate，FDR），报告假发现率小于1%的那些蛋白质比值作为可信的显著差异表达蛋白质，并以报表的形式展示给用户。

综上，本发明对于近百GB的规模的质谱实验采集的数据，快速地自动化分析，对不同蛋白质在质谱仪中的信号尽可能精准地提取蛋白质信号；从统计学意义上确定蛋白质差异表达，并对结果的准确性进行评价。

Claims (15)

1.一种检测差异表达蛋白质的方法，其特征在于，面向定量蛋白质组学中的基于一级谱图信息的标记和非标记的相对定量数据分析，根据某一蛋白质在两种或多种生物样品中对应的质谱信号强度比值判断其是否是差异表达蛋白质。

2.如权利要求1所述的检测差异表达蛋白质的方法，其特征在于，包括：

步骤1，对质谱数据进行预处理，用于是将该蛋白质的原始二进制质谱数据转换为文本格式，并建立索引；

步骤2，对二级谱图进行肽谱匹配，确定样品中含有的肽段，用于对二级谱图与蛋白质数据库中的记录的肽段进行匹配打分，取高可信的匹配结果；或者直接从二级谱图推测肽段序列，取高可信的结果；

步骤3，提取每个肽段在两种生物样品和多种生物样品中的信号，以多个同位素曲线的形式进行表示；

步骤4，将相同肽段在不同样品中的信号之间建立对应关系，计算其表达量差异的肽段比值和置信区间；

步骤5，将肽段比值归并为蛋白质比值，并给出蛋白质比值的置信区间；

步骤6，确定差异表达的蛋白质。

3.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，肽段信号以多个同位素色谱曲线的形式进行表示。

4.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，基于干扰最小的同位素色谱曲线计算肽段比值，具体方法为局部最小一乘法，并计算比值的置信区间。

5.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，将肽段比值归并为蛋白质比值采用核密度估计方法。

6.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，步骤2还包括：

步骤11，对每张二级谱图进行处理，只保留强度最大的前200个谱峰；

步骤12，对输入的每张二级谱图，在蛋白质数据库中寻找最相似的肽段。

7.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，步骤2还包括：

步骤21，对肽谱匹配按照打分从高到低进行排序；

步骤22，控制肽谱匹配的假发现率。

8.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，步骤3还包括：

步骤31，读取肽段排序列表，并对一级谱图进行预处理，同时设定必要的数据分析参数，包括要分析的质谱数据、数据的类型、所属物种的蛋白质数据库；

步骤32，提取肽段信号，以多个同位素色谱曲线的形式表示。

9.如权利要求8所述的检测差异表达蛋白质的方法，其特征在于，步骤32还包括：

步骤321，对每个肽谱匹配，计算每个肽段的理论同位素分布；

步骤322，根据理论同位素分布，在鉴定到该肽段的二级谱图前后2分钟保留时间范围内的一级谱图上确定实际同位素峰；

对于这个范围内的某张一级谱图，如果某理论同位素峰的质荷比正负10ppm范围内有谱峰，则记录下来；如果10ppm范围内有多根谱峰，那么每次选其中一个谱峰，和理论同位素分布计算余弦夹角，取值最高的那个组合为实际的同位素峰；

步骤323，对于每一个实际的同位素峰，沿保留时间把它们连成一条曲线，表示该肽段从没有信号到有信号再到信号消失的过程。

10.如权利要求9所述的检测差异表达蛋白质的方法，其特征在于，步骤323还包括：

先根据鉴定到肽段的二级谱图的扫描号，找到距离该二级谱图最近的一级谱图。而后以该一级谱图所在扫描号为基准，在这条曲线上向前寻找起始点，向后寻找终止点，当曲线的上某点的强度低于最高曲线强度的10%时，停止；该曲线在原来的基础上就小了一些；再以此曲线为基础，寻找起始点和终止点附近的极值点，将该极值点与起始点或终止点之间的部分删除；

步骤324，如果是标记定量实验，那么肽段有轻、重标两种形式。

11.如权利要求8所述的检测差异表达蛋白质的方法，其特征在于，步骤4还包括：计算肽段在两个样品中对应信号的强度差异，以比值的形式表示；取干扰最小的曲线计算比值。

12.如权利要求11所述的检测差异表达蛋白质的方法，其特征在于，步骤4还包括：

步骤41，对于一个肽段，用向量

和

表示计算的肽段同位素标记的轻、重标形式的理论同位素丰度；n和m分别表示肽段轻、重标形式的同位素峰数目；l₁＝(l_1,1,...,l_1,k)表示肽段轻标形式单同位素色谱曲线，k是色谱曲线跨越的一级谱图的数目；l₂＝(l_2,1,...,l_2,k)表示肽段轻标形式第一同位素色谱曲线，依次类推，l_n＝(l_n,1,...,l_n,k)表示肽段轻标形式第n同位素色谱曲线；类似的，h₁＝(h_1,1,...,h_1,k)表示肽段重标形式单同位素色谱曲线，h_m＝(l_m,1,...,h_m,k)肽段重标形式第m同位素色谱曲线；

步骤42，对上一步的各同位素曲线进行归一化，以确保肽段轻、重标记间各取任一同位素曲线都可用于计算肽段比值；

归一化方法如下： $L_1=l_1\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_1^l}=(l_{\mathrm{1,1}},...,l_{1,k}),$ $L_2=l_2\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_2^l}=(l_{\mathrm{2,1}},...,l_{2,k}),$ ……， $L_n=l_n\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_i^l}{t_n^l}=(l_{n,1},...,l_{n,k});$ $H_1=h_1\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_j^h}{t_1^h}=(h_{\mathrm{1,1}},...,h_{1,k}),$ $H_2=h_2\&CenterDot;\frac{{\mathrm{\&Sigma;t}}_j^h}{t_2^h}=(h_{2,1},...,h_{2,k}),$ ……，L₁，……L_i，……，L_n表示归一化后肽段轻标形式的各同位素曲线，H₁，……，H_j，……，H_m表示归一化后肽段重标形式的各同位素曲线；

步骤43，计算m×n个比值，分别是使用H_j和L_i计算比值，其中i＝1,...,n，j＝1,...,m；

步骤44，计算的m×n个比值分别对应m×n个置信区间，取置信区间最小的那个作为肽段的比值，该置信区间作为肽段比值的置信区间。

13.如权利要求12所述的检测差异表达蛋白质的方法，其特征在于，步骤43中计算比值包括如下步骤：

给定X=(x₁,…x_k)，Y=(y₁,…,y_k)表示分别从肽段轻标形式和重标形式中各取了一条同位素曲线。

步骤431，以两个曲线的中心为基础，取三个点为局部曲线：中点左边的点，中点，中点右边的点；即取计算和

的余弦夹角t；

步骤432，

和

分别向两端延伸，即取 $\tilde{X}=(x_{\frac{k}{2}-2},x_{\frac{k}{2}-1},x_{\frac{k}{2}},x_{\frac{k}{2}+1},x_{\frac{k}{2}+2}),$

$\tilde{Y}=(y_{\frac{k}{2}-2},y_{\frac{k}{2}-1},y_{\frac{k}{2}},y_{\frac{k}{2}+1},y_{\frac{k}{2}+2}),$ 再计算

和

的余弦夹角t；

步骤433）如果t大于一个阈值（根据经验，设定为0.8），就继续（2）中的步骤，否则停止；

步骤434）延伸停止后，使用过原点的最小一乘法

计算

和的比值a。

步骤435，计算比值的置信区间，计算

当样本大小n→∞时，有

那么，可得a的1-α的区间估计是

其中u_α/2是标准正态分布N(0,1)的“上α/2分位点”；令 $\hat{f}\left(0\right)=\frac{1}{\mathrm{nw}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}K\left(\frac{{\hat{e}}_i}{w}\right),$ 选择核函数 $K\left(x\right)=\frac{3}{4}(1-x^2)I_{\left|x\right|\&le;1},$ 窗宽 $w=1.66n^{-1/5}\widehat{\mathrm{\&sigma;}},$ 其中 $\widehat{\mathrm{\&sigma;}}=\sqrt{s_n^2}.$

14.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，步骤5还包括：

步骤51，使用核密度估计的方法推断蛋白质的比值；

步骤52，使用了高斯核，并假设每个肽段的比值服从一个高斯分布：

$f_{\mathrm{pep}}\left(x\right)=\frac{1}{N}\&CenterDot;\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_i}{e}^{\frac{{(x-r_i)}^2}{2{{\mathrm{\&sigma;}}_i}^2}}$

其中，r_i是第i个肽段比值，σ_i是其标准差，N是其所对应的蛋白质所鉴定或定量的肽段的总数；

步骤53，定义F(x)＝∑f_pep(x)，根据F(x)拟合一个高斯曲线，其形式如下：

$f_{\mathrm{pro}}\left(x\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{pro}}}{e}^{\frac{{(x-R)}^2}{2{\mathrm{\&sigma;}}_{\mathrm{pro}}^2}}$

其中，R就是这个蛋白质的比值，σ_pro是其标准差,根据这个标准差，就可以计算R的置信区间。

15.如权利要求2所述的检测差异表达蛋白质的方法，其特征在于，步骤6还具体包括如下步骤：

步骤61，对所有比值取log2变换，并通过加或减一个常数，使其中值为0；

步骤62，拟合所有比值的分布；

步骤63，根据拟合的分布计算每个比值的p-value；

步骤64，计算假发现率，报告假发现率小于1%的那些蛋白质比值作为可信的显著差异表达蛋白质，并以报表的形式展示给用户。

Images