CN111537659A

CN111537659A - 一种筛选生物标志的方法

Info

Publication number: CN111537659A
Application number: CN202010492720.XA
Authority: CN
Inventors: 刘超; 郭一洁; 宫鹏云; 李春燕
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2020-06-03
Filing date: 2020-06-03
Publication date: 2020-08-14

Abstract

本发明公开了一种筛选生物标志的方法，其包括：以非标记定量分析方法，分别获得两组生物样品中不同生物标志类物质的名称和对应的相对强度值，并将不同生物标志类物质的相对强度值进行归一化处理，对两组生物样品中每一种生物标志类物质对应的归一化后的平均强度值分别进行t检验，计算出p‑value，并采用BH法对p‑value进行校正，以计算得到的FDR值作为q‑value，其小于或等于显著性水平，则该q‑value对应的生物标志类物质为差异表达。通过BH方法校正p‑value，在检验出尽可能多的候选差异表达蛋白质的同时将错误发现率控制在一个可以接受的范围内，从而提高了检验准确性。

Description

一种筛选生物标志的方法

技术领域

本发明涉及生物技术与蛋白质组学技术领域，具体而言，涉及一种筛选生物标志的方法。

背景技术

生物体内含有多种物质，其中有些物质发生了差异性表达，需要通过一定的技术手段来获得差异性表达的信息，例如，定量蛋白质组学通过对质谱仪产生的大量质谱数据进行分析，比较不同样本中蛋白质相应肽段的信号强度，从而对肽段对应的蛋白质进行相对定量。定量蛋白质组学中，差异蛋白质的检测是一项重要的研究目标。

专利号为ZL201310397694.2的专利中公开了一种检测差异表达蛋白质的方法，包括：肽谱匹配、可信度评价、肽段信号提取、肽段比值计算、蛋白质比值计算、统计学分析。其中，肽谱匹配对实际二级谱图与理论二级谱图进行匹配打分，取可信度高的匹配结果，从实际得到的一级谱图中提取每个肽段在各种样品中的信号，计算相同肽段在不同样品中的信号比值，并将肽段比值整合为相应蛋白质比值，并给出蛋白质比值的置信区间，最终通过统计学分析确定差异蛋白。专利号为ZL200910045221.X也公开了一种利用蛋白质芯片技术鉴定差异表达蛋白的新方法，包括：将实验样本进行蛋白质芯片的实验操作，每个样本得到一系列蛋白质离子峰的数据；根据实验设计的情况对样本进行分组，得到不同样本间的差异蛋白峰；最终对差异蛋白峰进行蛋白质预测。

但是，以上方法和其他已有的利用蛋白质质谱数据进行定性定量分析的方法，均存在的主要问题是差异蛋白质的检测结果的准确性无法保证。

鉴于此，特提出本发明。

发明内容

本发明的目的在于提供一种筛选生物标志的方法，以改善上述问题。

本发明是这样实现的：

本发明的实施例提供了一种筛选生物标志的方法，其包括：以非标记定量分析方法，分别获得两组生物样品中不同生物标志类物质的名称和对应的相对强度值，并将不同生物标志类物质的相对强度值进行归一化处理；其中，两组生物样品对应的生物标志类物质的种类相同。

对两组生物样品中每一种生物标志类物质对应的归一化后的平均强度值分别进行t检验，以通过以下公式获得t-value。

其中，m和n分别为两组生物样品中的样品个数，s(i)为标准差，s₀为预设的0～2范围内的最小差异倍数值，a为自由度，

和

分别为两组生物样品中的第i种相同生物标志类物质的归一化后的平均强度值。

通过t分布的累积分布函数来检索观察t-value的绝对值的累积概率，从而计算出p-value，采用Benjamini-Hochberg法对p-value进行校正，计算得到FDR值。如果某一个p-value所对应的FDR值大于排序的前一位p值所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值；反之则保留计算的FDR值。最终得到的值记为q-value。q-value小于或等于显著性水平，则该q-value对应的生物标志类物质为差异表达的生物标志类物质。

需要说明的是，一些实施方式中，上述筛选生物标志的方法不以疾病的诊断和治疗为目的。

本发明具有以下有益效果：通过将不同生物标志类物质的相对强度值进行归一化处理，以提高后续t检验的准确度，进一步计算出p-value并通过BH方法校正p-value，在检验出尽可能多的候选差异表达蛋白质的同时将错误发现率控制在一个可以接受的范围内，从而提高了检验准确性。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明一种实施方式的流程图；

图2为实施例1的酵母溶菌液和纯化蛋白质浓度；

图3为实施例1的样品制备信息；

图4为实施例1质谱图；

图5为实施例1质谱图重构的色谱曲线；

图6为显示样品蛋白质强度分布的直方图；

图7为显示样品之间相关性的散点图；

图8为显示差异蛋白结果的火山图；

图9为sp|P44015|VAC2_YEAST在火山图中位置显示的火山图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将对本发明实施例中的技术方案进行清楚、完整地描述。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。

下面对本发明提供的一种筛选生物标志的方法进行具体说明。

本发明的一些实施方式提供了一种筛选生物标志的方法，其包括：

S1、以非标记定量分析方法，分别获得两组生物样品中不同生物标志类物质的名称和对应的相对强度值，并将不同生物标志类物质的相对强度值进行归一化处理；其中，两组生物样品对应的生物标志类物质的种类相同。

具体地，一些实施方式中，生物标志类物质包括但不限于基因转录或翻译产物、大分子蛋白质和小分子代谢物中的任意一种，例如，生物标志类物质为蛋白质，生物标志为差异表达的蛋白质。

非标记定量(label-free)蛋白质组学技术是通过液质联用技术对蛋白质酶解肽段进行质谱分析，该技术不需要使用昂贵的稳定同位素标签做内部标准，只需分析大规模鉴定蛋白质时所产生的质谱数据，比较不同样品中相应肽段的信号强度，从而对肽段对应的蛋白质进行相对定量。

经典的蛋白质组学流程首先用胰酶将蛋白质酶切成肽段，然后经过液相色谱分离，根据肽段的疏水性、离子强度、等电点的差异进行分离洗脱，接着利用电喷雾离子源产生气相离子.质谱仪根据它们的质荷比将其分离并记录，这一阶段称为母离子扫描.然后按照不同质谱仪设定的规则，通常选择丰度最高的5～9个母离子依次进行二级碎裂，得到二级碎片离子谱图，再返回至母离子扫描，进入下一循环。这样的一种质谱数据获取方式被称为数据依赖获取(data dependent acquisition,DDA)模式，也是最为常见的一种模式。在一级质谱图中，每一个母离子都包括三维信息，即液相色谱保留时间、质荷比以及离子强度。由于电喷雾离子化得到的母离子信号强度与离子浓度存在相关性，每个母离子都是离子化的肽段，因此在一级谱图中提取出鉴定肽段对应的离子峰强度即可以反映出肽段的丰度，常用来表示离子峰强度的参数有峰高度、峰面积等。

具体地，一些实施方式中，生物标志类物质为蛋白质时，非标记定量分析方法是使用蛋白质定量软件进行非标记定量分析，导入数据库文件和含有蛋白质质谱数据的raw文件，数据依赖性采集方法将给定质荷比范围内的离子碎裂，提取母离子谱图，并在保留时间内积分得到峰面积，通过比较样品的峰面积完成定量。在通过质谱仪扫描记录的肽段强度来确定蛋白质的相对强度时，影响定量准确度的主要因素是具有相似质荷比信息的共洗脱离子的信号干扰。

一些实施方式中，蛋白质的质谱数据是将液相色谱法分离得到的肽段通过通过质谱仪检测所有离子的质荷比信息，并记录其相应保留时间上的信号强度，最终得到结果序列数据，并记录在结果文件raw文件中。

一些实施方式中，在用液相色谱法分离之前，采用胰蛋白酶将蛋白质酶解成肽段。即需要在进行质谱分析之前，需要制备样品以去除复杂混合物中不感兴趣的蛋白质，降低样品的复杂性。由于肽段比整个蛋白质更能有效离子化，通常使用胰蛋白酶将蛋白质酶解成肽段，并用液相色谱法分离。

一些实施方式中，还需要确定生物样品中蛋白质种类，以对应匹配不同蛋白对应的相对强度值，具体过程为利用质谱法进行蛋白质鉴定和序列分析，对蛋白质的鉴定通常使用搜库方法实现。将实际得到的二级谱图与数据库中理论二级谱图进行肽谱匹配，确定样品中含有的肽段。通过对匹配得到的候选肽段进行匹配打分，取得分最高的匹配结果，最终得到对蛋白质的鉴定结果。

为了后续方便对蛋白质的数据进行处理，读取定量结果文件，将其蛋白名和强度值信息存入二维矩阵，以使得后续对蛋白质数据的处理能够通过对矩阵操作实现。

进一步地，为了提高筛选结果的准确性，生物标志类物质为蛋白质时，在将蛋白质的相对强度值归一化处理之前或之后，可对不同蛋白质的蛋白名和对应的相对强度值的数据进行过滤，即将蛋白名和相对强度值存入二维矩阵，然后过滤来自反库的蛋白质和有潜在污染的蛋白质，其中，来自反库的蛋白质蛋白名前缀为“REV_”，有潜在污染的蛋白质蛋白名前缀为“CON_”。

可选地，采用循环倒序删除的方法，倒序遍历蛋白名，判断是否含有目标前缀，若符合条件则删除该数据行强度值以及蛋白名信息。该步骤完成后得到含有蛋白名和强度值的二维矩阵。首次删除操作完成后，无论之后的元素怎么移动，之前的元素对应的索引不发生变化。倒序遍历元素方法避免空间回收导致索引越界、定位错误的情况发生。

进一步地，若在归一化之前进行过滤，为了提高筛选结果的准确性还包括对过滤后的数据进行缺失值处理。缺失值产生的原因包括蛋白质仅在个别样品中存在，或者浓度低于质谱仪检测下限，或者蛋白质质谱图无法得到正确匹配。具体的操作方式可采取将含有0值的蛋白质数据过滤的方法，例如，倒序遍历相对强度值，判断是否含有0值，若符合条件则删除该行。

需要说明的是，还可以通过其他方法对缺失值进行处理，即判断技术重复实验对同一个蛋白质的定量值的个数。用户自行设置一个值，如果定量值中缺失值的个数大于等于该值，则删除该蛋白质数据；如果小于该值，则根据正态分布在缺失处进行插补。

在对蛋白质数据进行过滤和缺失值处理后，能够提高蛋白质数据的有效性，再通过归一化处理，其目的是为了数据方便处理所做的线性变换，不改变原始数据的数值排序而将其映射到想要的范围内。而在具体的蛋白质非标定量分析中，由于质谱仪输出的强度值高达11个数量级，因此会导致后期数据处理出现较大的误差，而无法获得有效的结果。因此，一些实施方式中，归一化处理为将不同蛋白质的相对强度值换算为2为底数的对数的强度值，便于后续处理。该步骤完成后得到含有蛋白名和取完以2为底的强度值的二维矩阵。

最小二乘法做归一化，归一化本质上是一种线性变化，同时不改变原始数据的排序。当数据比较集中时，样本的方差小，归一化后数据分布会分散；而当数据分布较分散时，样本方差较大，归一化后数据会被集中到更小的范围内。绘制蛋白质强度相关性分析图，可以通过数据分布情况判断是否需要进行归一化处理。假设两个技术重复实验的数据之间的线性关系为f(x)＝ax+b，最小二乘法的思想计算总误差的平方和为：

∈＝∑(f(x_i)-y_i)²＝∑(ax_i+b-y_i)²

不同的a、b会导致不同的总误差平方和，根据多元微积分的知识，当

时取最小值，解方程组即可得到a、b的值。

进一步地，本发明的实施方式是为了比较两组生物样品进而筛选中差异性表达的生物标志例如差异表达的蛋白质。因此，为了提高结果的准确性，两组生物样品均为同一类物质，例如，均为酵母溶菌液，其含有相同种类的蛋白质，大部分蛋白质的比例均为1：1，少量蛋白质的含量差异较大，属于差异性表达。因此，在进行筛选过程中，每组生物样品的样品数量至少为3个，优选3个。即对每组生物样品均进行至少2次或3次的技术重复性实验，每次进行的均为一个样品。例如，对其中一组的生物样品做3次技术重复性实验，记为1A、1B、1C；对另外一组生物样品也做3次技术重复性实验，记为2A、2B、2C。

进一步，一些实施方式中，还对各样品进行相关性分析；计算得到pearson系数并绘制相关性分析散点图，横、纵坐标分别为相比较的两次技术重复实验的归一化的强度值，图中点表示一个蛋白质。pearson系数用以估计两样本的相关性，定义为两样本的协方差和两者标准差乘积的比值。pearson系数取值范围在-1和1之间，pearson系数的绝对值越接近1则说明两样品的相关性越好。

总体相关系数ρ定义为两个变量X、Y之间的协方差和两者标准差乘积的比值，如下：

估算样本点协方差和标准差，可以得到样本相关系数(即样本皮尔森相关系数)，常用r表示：

R还可以由(Xi，Yi)样本点的标准分数均值估计得到与上式等价的表达式：

其中:

σ_X

为Xi样本的标准分数、样本均值和样本标准差，n为样本数量。

S2、对两组生物样品中每一种生物标志类物质对应的归一化后的平均强度值分别进行t检验，以通过以下公式获得t-value；

其中，m和n分别为两组生物样品中的样品个数，s(i)为标准差，s0为预设的0～2范围内的最小差异倍数值，a为自由度，

和

以生物标志类物质是蛋白质为例，上述t检验通过比较t-value的绝对值与临界值的大小判断该蛋白质是否为差异表达蛋白质。当t-value的绝对值大于临界值例如0.05时，有充分的理由判断蛋白质为差异表达。

S3、通过t分布的累积分布函数来检索观察t-value的绝对值的累积概率，从而计算出p-value，采用Benjamini-Hochberg法对p-value进行校正，计算得到FDR值。如果某一个p-value所对应的FDR值大于排序的前一位p值所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值；反之则保留计算的FDR值。最终得到的值记为q-value。q-value小于或等于显著性水平，则该q-value对应的生物标志类物质为差异表达的生物标志类物质。

具体地，使用t分布的累积分布函数(cumulative distribution function,CDF)来检索观察t统计量的绝对值的累积概率，从而计算出p-value。累积分布函数输入t-value，输出点小于输入值的概率。查表得到对应的p-value。该步骤完成后得到含有蛋白名、取完以2为底的强度值和p-value的二维矩阵。

假设检验也叫显著性检验，是以小概率反证法的逻辑判断假设是否成立的统计学方法。本发明实施方式原假设H0为两组数据来自同一分布，即两份生物样品之间没有没有显著差异。相应备择假设H1为两组数据来自不同分布，也就是说两份生物样品之间存在显著差异，存在差异蛋白。

假设检验存在两种错误。第一类错误即假阳性错误，就是在假设检验作推断结论时，拒绝了实际上是正确的原假设。犯第一类错误的概率用α表示。第一类错误是针对原假设而言的，α就是事先规定的允许犯第一类错误的概率值，如规定α＝0.05，意味着在某特定总体抽样，100次拒绝H0的假设检验中，最多有5次允许发生第一类错误。与此相应，推断正确的可能性为1-α，1-α又称为可信度。第二类错误即假阴性错误，即接受实际上是不成立的H0。就是无效假设原本是不正确的，但所算得的统计量不足以拒绝它，错误地得出了无差别的结论。第二类错误是针对备择假设而言的，其概率值用β表示。β值的大小一般未知。

p-value是拒绝原假设H0的最小显著性水平。这里的p-value是指由H0成立时的检验统计量出现在由样本计算出来的检验统计量的末端或更末端处的概率值。显著性水平是指当原假设为正确时错误拒绝的概率，用α表示。它是公认的小概率事件的概率值，必须在每一次统计检验之前确定，一般取0.01或0.05，表示正确接受原假设的概率为99％或95％。p-value是拒绝原假设H0的最小显著性水平。显著性水平是指当原假设为正确时错误拒绝的概率，用α表示，一般取0.05，表示正确接受原假设的概率为95％。当p-value<α时，得到结论：拒绝原假设H0，接受备择假设H1。这样下结论的理由是：在原假设成立的条件下，备择假设发生是小概率事件，即现有样本信息不支持原假设因而拒绝它。当p-value>α时，接受原假设。若p-value＝α，认为没有充分的理由拒绝或接受原假设，通常增加样本容量继续做假设检验，直至得到严格大于或小于结果。

进一步地，采用多重假设检验筛选出差异蛋白。每个蛋白质都需要做一次假设检验来判断是否为差异表达的蛋白质。当假设检验的次数大于等于2时，每个单独的假设检验都具有I型错误，即错误拒绝原假设，犯I型错误的概率会随着假设检验的个数迅速增加。蛋白质组学希望尽量多地识别出差异表达地蛋白质，并且能够容忍和允许总的拒绝原假设事件中发生少量地错误识别。即需要在错误发现率和总的拒绝次数之间寻找一种平衡，在检验出尽可能多的候选差异表达蛋白质的同时将错误发现率控制在一个可以接受的范围。通常需要做多重假设检验对每个p-value进行校正，控制错误发现率(False DiscoveryRate,FDR)小于等于显著性水平。蛋白质组学常用的显著性水平的取值为0.05和0.01，表示显著性结果的5％或1％为假阳性。

因此，一些实施方式中，采用Benjamini-Hochberg对p-value进行校正包括：将p-value按照升序排列，使用公式FDR＝(n/i)p-value计算每个p-value对应的FDR值，其中，n为p-value个数，i为p-value排序后的序号，计算出来FDR值。如果某一个p-value所对应的FDR值大于排序的前一位p值所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值；反之则保留计算的FDR值。最终得到的值记为q-value。筛选小于或等于显著性水平的q-value，筛选出的q-value对应的蛋白质为差异表达的蛋白；优选地，显著性水平为0.05。

进一步地，为了更为直观地展示差异表达的生物标志类物质(蛋白质)，一些实施方式还包括：

S4、对筛选结果绘制火山图，火山图横坐标为生物标志类物质对应的归一化后的平均强度值，纵坐标为q-value取以10为底的对数再取负数，每个点表示一个生物标志类物质。

横坐标越接近0表示该蛋白质在两样品之间的差异越小，纵坐标越大表示p-value越小，即是差异性表达的蛋白质的可能性越大。火山图中可将筛选的差异表达的蛋白质标记为其他颜色进行区别。

进一步地，一些实施方式，还包括输出差异生物标志类物质的信息列表。列表信息包括生物标志类物质的名称、火山图横坐标、火山图纵坐标、q-value和是否差异表达。

本发明的一些实施方法还提供了一种筛选差异蛋白的方法，如图1所示，其具体包括以下步骤：

步骤1、搜集样品。在进行质谱分析之前，需要制备样品以去除复杂混合物中不感兴趣的蛋白质，降低样品的复杂性。由于肽段比整个蛋白质更能有效离子化，通常使用胰蛋白酶将蛋白质酶解成肽段，并用液相色谱法分离。

步骤2、采集样品的质谱数据。液相色谱法分离得到的肽段通过质谱仪检测所有离子的质荷比信息，并记录其相应保留时间上的信号强度，最终得到结果序列数据用于确定样品的原始蛋白质成分，记录在结果文件raw文件中。

步骤3、鉴定样品中蛋白质种类。利用质谱法进行蛋白质鉴定和序列分析，对蛋白质的鉴定通常使用搜库方法实现。将实际得到的二级谱图与数据库中理论二级谱图进行肽谱匹配，确定样品中含有的肽段。通过对匹配得到的候选肽段进行匹配打分，取得分最高的匹配结果。最终得到对蛋白质的鉴定结果。

步骤4、对蛋白质进行定量分析。使用蛋白质定量软件进行无标记定量分析，导入raw文件和数据库文件。数据依赖性采集方法将给定质荷比范围内的离子碎裂，提取母离子谱图，并在保留时间内积分得到峰面积或者提取色谱峰最高点强度，通过比较样品的面积或强度完成定量。

步骤5、读取定量结果文件，将其蛋白名和强度值信息存入二维矩阵，后续对蛋白质数据的处理通过对矩阵操作实现。

步骤6、根据蛋白名前缀对步骤1中矩阵进行过滤。非标记定量蛋白质组学通常过滤来自反库的蛋白质、有潜在污染的蛋白质，其中，来自反库的蛋白质蛋白名前缀为“REV_”，有潜在污染的蛋白质蛋白名前缀为“CON_”。这里使用循环倒序删除的方法，倒序遍历蛋白名，判断是否含有目标前缀，若符合条件则删除该数据行强度值以及蛋白名信息。该步骤完成后得到含有蛋白名和强度值的二维矩阵。

步骤7、对缺失值进行处理，并对矩阵中的强度值进行归一化处理。采取将含有0值的蛋白质数据过滤的方法。倒序遍历强度值，判断是否含有0值，若符合条件则删除该行。由于质谱仪输出的强度值高达11个数量级，取其以2为底的对数，便于后续处理。该步骤完成后得到含有蛋白名和取完以2为底的强度值的二维矩阵。

步骤8、绘制样品各技术重复实验蛋白质强度分布直方图。横坐标为取过以2为底的对数的强度值，划分为等间隔的若干个区间，区间个数可由用户设定；纵坐标为某强度区间内蛋白质的个数。强度分布直方图显示数据大致服从某一分布，作为后续假设检验方法选择的依据。

步骤9、对样品进行相关性分析，计算得到pearson系数并绘制相关性分析散点图。横、纵坐标分别为相比较的两次技术重复实验的取过以2为底的对数的强度值，图中点表示蛋白质。

步骤10、对步骤7所得矩阵强度值进行分组。双样本t检验之前，将两份生物样品数据存放在一个变量中，同时对同一份生物样品的技术重复实验设置相同的名字。

步骤11、对步骤10分组后样品进行双样本t检验，得到p-value并对应存入矩阵。t检验通过检查来自两个样品的平均值之间的标准误差来确定它们是否具有显著的差异。原假设为两个样本具有相同的均值。相应备择假设为两个样品具有不同的均值。使用t分布的累积分布函数来检索观察t统计量的绝对值的累积概率，从而计算出p-value。累积分布函数输入t-value，输出点小于输入值的概率。查表得到对应的p-value。该步骤完成后得到含有蛋白名、取完以2为底的强度值和p-value的二维矩阵。

步骤12、进行多重假设检验筛选出差异蛋白。使用的p-value校正方法为BH(Benjamini-Hochberg)，将p-value按照升序排列，使用公式FDR＝(n/i)p-value计算每个p-value对应的FDR值，其中，n为p-value个数，i为p-value对应序号数。将计算出来的FDR值作为新p-value，如果某一个p-value所对应的FDR小于前一位p-value所对应的FDR值，选择使用与它前一位相同的值，反之则保留计算的FDR值。校正后的p-value即q-value，q-value小于等于显著性水平，则该q-value对应的蛋白质为差异表达的蛋白质。

步骤13、对筛选结果绘制火山图，直观展示差异表达的蛋白质。火山图横坐标为log₂(sample1技术重复实验的均值/sample2技术重复实验的均值)，纵坐标为调整后的p-value取以10为底的对数再取负数，每个点表示一个蛋白质。横坐标越接近0表示该蛋白质在两样品之间的差异越小，纵坐标越大表示p-value越小，即是差异性表达的蛋白质的可能性越大。火山图中将步骤13筛选的差异表达的蛋白质标记为红色或其他颜色。

步骤14、输出差异蛋白的信息列表。列表信息包括蛋白质的蛋白名、火山图横坐标、火山图纵坐标、q-value和是否差异表达。

以下结合实施例对本发明的特征和性能作进一步的详细描述。

实施例1

本实施例以2015年在圣路易斯市的会展中心举行的美国质谱年会上，ABRF(TheAssociation of Biomolecular Resource Facilities)报告了由其组织的蛋白质组学数据分析评测实际操作过程为例。

iPRG原始数据来源于2份相同的经过酶切的酵母溶菌液分别人工混入少量不同浓度的纯化的蛋白质，这样处理过的2份样品中，绝大多数蛋白质的比例都是1:1，只有少量蛋白质的比例不同。如图2所示，其中，A、B、C、D、E、F分别代表人工混入的不同浓度纯化的蛋白质，Sample1、Sample2表示两份生物样品，数字表示差异蛋白在样品中的浓度，单位为fmol，颜色深浅表示浓度高低，颜色越深即浓度越大。蛋白质来源信息如表1所示。接着分别做三次重复实验，一共得到6个RAW文件。一级谱图、二级谱图、母离子及电荷状态等信息都存储于RAW文件中。2份样品分别记为1、2号，三次重复实验分别记为A、B、C，6个RAW文件分别记为1A、1B、1C、2A、2B、2C。实验之间相互独立，6个RAW文件之间两两可比。RAW文件由iPRG官方统一提供。具体样本制备信息如图3所示。

表1纯化蛋白质名称及来源

本实施例数据分析目标是找到人工加入的浓度有差异的蛋白质。

一、定量分析的参数设置

下面介绍本实施例在非标记定量时的参数配置。

每组均使用Label-free quantification(LFQ)定量方法。Intensity作为蛋白质原始强度值，LFQ Intensity则是将原始强度值在三次重复实验样本之间进行矫正，以消除处理、上样、预分、仪器等造成的样本间误差。选择使用LFQ Intensity进行后续数据处理。

每组添加Variable modifications和Fixed modifications如表2所示。

表2修饰参数设置

表2设置的固定修饰为Carbamidomethyl(C)，指的是修饰位点位于肽段含有羧基的C端，修饰基团为Carbamidomethyl，即半胱氨酸上发生烷基化修饰。可变修饰Acetyl(Protein N-term)、Deamidation(N)分别指的是位于蛋白质含有氨基的N端发生乙酰化、脱酰胺修饰，Oxidation(M)为甲硫氨酸上发生氧化修饰。这些都是人体内常见的氨基酸。

添加数据库文件。数据库文件中包含核苷酸序列或氨基酸序列的信息，氨基酸常用大写字母表示，每一个大写字母代表一个氨基酸或多个氨基酸。数据库文件用来生成肽段搜索空间。

使用Match between runs方法。表示只匹配相同或相邻碎片的特征。理想情况下，每个RAW文件都可以得到足够的信息，可以从测量得到的二级质谱中检索肽段序列。然而，有时可用信息不足，测得的二级质谱要么不足以识别序列，要么甚至没有被测得。为了得到一个已识别的特征，可以根据质荷比信息，通过在质量和保留时间窗口内的匹配，从另一个RAW文件中获得二级质谱和肽段序列信息。为此，首先要将保留时间对齐。除了已经被二级质谱数据库搜索引擎测序和鉴定的肽外，这增加了可用于量化的肽的数量。

二、定量算法

1.重构色谱曲线

RAW文件中质谱本身不能反映离子信号从无到有再消失的过程，如图4所示，通过提取母离子/碎片离子的在各保留时间处的峰值，绘制色谱曲线，反映信号在质谱仪中的强度变化过程。将每条肽段的同位素峰在色谱曲线中的强度加和，作为原始的肽段信号强度，以对色谱曲线进行重构，如图5所示。

2.强度归一化

理论上，归一化可以消除肽段信号在实验时受到的不可避免的干扰。由于构成蛋白质的肽段是在相邻几个扫描中分别传播的，首先需要确定每条肽段的归一化系数，对一个蛋白质的每条肽段赋予归一化系数后再将强度进行加和，作为蛋白质归一化后的强度值。

当样品个数增多时，H(N)计算起来非常复杂，因此可以应用Fast label-free算法，将每个样品当成一个节点，两个样品中同时鉴定到的肽段作为连接这两点的一条边，两点之间只要有一条边相连即可进行计算。参数中设置每个节点至少有三个、平均有六个相邻节点。使用Fast LFQ算法，通过减少样本两两比较的数量来确保合理的计算时间。

设置“LFQ min.ratio count”为2，表示当一个蛋白质在两个样品中被鉴定到的相同肽段数目大于或等于2时，才将这些肽段用于计算比值，而当一个蛋白质中的所有肽段都不符合条件时，将该蛋白质的强度记为0。这个参数的值越大，得到的结果越精确。

三、数据处理流程

将ProteinName和LFQ intensity存为二维矩阵。如表3所示。

表3二维矩阵

对来自反库的蛋白质、有潜在污染的蛋白质在矩阵中进行过滤。具体操作为：

根据ProteinName进行数据过滤。过滤掉来自反库的蛋白质(前缀为“REV_”)数据行，过滤掉有潜在污染的蛋白质(前缀为“CON_”)数据行。

对缺失值进行处理：

对于一个蛋白质的定量结果，根据每份生物样品的若干次技术重复实验中含有的缺失值个数进行过滤。这里采取过滤掉蛋白质强度含有0的所有数据行。也就是说，对于一个蛋白质的定量结果，6个raw文件中的强度值中有一个0值就在二维矩阵中将该蛋白质数据行删除。

对LFQ intensity数据取以2为底的对数，蛋白质强度比值的计算便由乘除法变成了加减法运算。

绘制蛋白质强度分布直方图，如图6所示，表示各样品中蛋白质相对丰度的分布情况。

观察样品间的相关性，绘制所有样品的相关性散点图。如图7所示，横(从左至右)、纵(从上至下)坐标依次为sample1_A、sample1_B、sample1_C、sample2_A、sample2_B、sample2_C，一个散点图每个点表示一个蛋白质，数字为pearson系数，系数值越接近1表示两样品相关性越强。

对所得矩阵强度值进行分组。双样本t检验之前，将两份生物样品数据存放在一个变量中，同时对同一份生物样品的技术重复实验设置相同的名字。

对两组生物样品中每一种生物标志类物质对应的归一化后的平均强度值分别进行t检验，以通过以下公式获得t-value；

其中，m和n分别为两组生物样品中的样品个数，s(i)为标准差，s₀为0，a为自由度，

和

通过t分布的累积分布函数来检索观察t-value的绝对值的累积概率，从而计算出p-value，采用Benjamini-Hochberg对p-value进行校正，即将p-value按照升序排列，使用公式FDR＝(n/i)p-value计算每个p-value对应的FDR值，其中，n为p-value个数，i为p-value排序后的序号，计算出来FDR值。如果某一个p-value所对应的FDR值大于排序的前一位p值所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值；反之则保留计算的FDR值。最终得到的值记为q-value。筛选小于或等于显著性水平的q-value，筛选出的q-value对应的蛋白质为差异表达的蛋白；其中，选择显著性水平为0.05。

控制被错误拒绝的比例FDR(false discovery rate)＝5％。FDR指错误拒绝原假设的个数占所有被拒绝的原假设个数的比例的期望值，是假设检验错误率的控制指标，也是筛选出的差异蛋白的评价指标。FDR取0.05时，最多有5％的差异蛋白是错误的。

绘制火山图(Volcano plot)。火山图将所有检测到的蛋白质的差异显著性进行可视化展示。本文使用Difference(sample1-sample2)作为横坐标，即log₂(ratio)，其中ratio为两份生物样品3次技术重复实验的均值的比值，-log(p-value)作为纵坐标绘制火山图如图8所示。p-value越小，-log(p-value)越大，则拒绝原假设的显著性越大，同时蛋白质强度比值越大，log₂(ratio)越远离0，则蛋白质表达差异越显著。火山图中具有文字标记的点即为筛选出的差异蛋白。

输出差异蛋白的信息列表。如表4所示，列表信息包括蛋白质的蛋白名、火山图横坐标、火山图纵坐标、q-value和是否差异表达。第二列Difference表示该蛋白质在火山图中的横坐标，第三列Pvalue(-log10)表示该蛋白质在火山图中的纵坐标，第四列qvalue为校正后的p-value，第五列Significant_Difference分别用TRUE或FALSE标明是否为差异表达蛋白质。

表4差异蛋白质列表

四、结果显示

1.定量结果显示

鉴定得到3170个蛋白质，过滤掉来自反库或有污染的47个蛋白质，用于定量分析的蛋白质有3123个。其中，来自反库的蛋白质体现为Protein IDs的“REV_”前缀。有污染的蛋白质体现为Protein IDs的“CON_”前缀。

2.差异蛋白结果显示

差异蛋白ID如表5所示，图8火山图中文字标记处的空心方块处显示差异蛋白位置。

表5差异蛋白ID

3.手工验证

以蛋白质sp|P44374|SFG2_YEAST为例，查看该蛋白质的任意两条肽段的单同位素峰在6次实验中对应保留时间处的原始强度，如表6所示。

表6查看RAW文件中肽段强度

观察上表可以发现，这两条肽段在1号、2号样品的三次重复实验中的强度变化不大，而在两组生物样品中存在一个数量级的差异。因此，用实施例方法找差异蛋白是可行的。

4.结果正确性分析

如表7所示，本次数据分析所得结果存在漏报，而不存在错报。

表7差异蛋白正确ID

5.漏掉的差异蛋白分析

在火山图中查找并定位未被筛选出来的蛋白质sp|P44015|VAC2_YEAST和sp|P44683|PGA4_YEAST，如图9所示，sp|P44015|VAC2_YEAST可以在火山图中定位到。由于其p-value比较大，且其Difference(sample1-sample2)较小，即在两组样本中强度相差不大，无法通过火山图被筛选出来。而sp|P44683|PGA4_YEAST在火山图中找不到，但存在于定量结果中。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种筛选生物标志的方法，其特征在于，其包括：

以非标记定量分析方法，分别获得两组生物样品中不同生物标志类物质的名称和对应的相对强度值，并将不同生物标志类物质的相对强度值进行归一化处理；其中，两组生物样品对应的生物标志类物质的种类相同；

对两组生物样品中每一种生物标志类物质对应的归一化后的平均强度值分别进行t检验，并通过以下公式获得t-value；

和

分别为两组生物样品中的第i种相同生物标志类物质的归一化后的平均强度值；

通过t分布的累积分布函数来检索观察t-value的绝对值的累积概率，从而计算出p-value，采用Benjamini-Hochberg法对p-value进行校正，以计算得到的FDR值作为校正后的q-value，q-value小于或等于显著性水平，则该q-value对应的生物标志类物质为差异表达的生物标志类物质。

2.根据权利要求1所述的筛选生物标志的方法，其特征在于，所述生物标志类物质选自基因转录或翻译产物、大分子蛋白质和小分子代谢物中的任意一种，优选地，所述生物标志类物质为蛋白质，所述生物标志为差异表达的蛋白质。

3.根据权利要求1所述的筛选生物标志的方法，其特征在于，所述生物标志类物质为蛋白质时，在将蛋白质的相对强度值归一化处理之前或之后，对不同蛋白质的蛋白名和对应的相对强度值的数据进行过滤，所述过滤为将蛋白名和相对强度值存入二维矩阵，然后过滤来自反库的蛋白质和有潜在污染的蛋白质，其中，来自反库的蛋白质蛋白名前缀为“REV_”，有潜在污染的蛋白质蛋白名前缀为“CON_”；

优选地，对所有数据行采用循环倒序删除的方法，倒序遍历蛋白名，判断是否含有目标前缀，若符合条件则删除该数据行强度值以及蛋白名信息。

4.根据权利要求3所述的筛选生物标志的方法，其特征在于，还包括对过滤后的数据进行缺失值处理，所述缺失值处理包括将相对强度值为0值的蛋白质数据进行过滤，倒序遍历相对强度值，判断是否含有0值，若符合条件则删除该行；

优选地，归一化处理为将不同蛋白质的相对强度值换算为2为底数的对数的强度值。

5.根据权利要求2所述的筛选生物标志的方法，其特征在于，每组生物样品的样品数量至少为3个，优选3个，优选地，绘制两组生物样品中各样品的技术重复性实验的蛋白强度分布直方图，横坐标为归一化后的强度值，划分等间隔的多个区间，纵坐标为对应强度值的区间内蛋白质的个数，通过蛋白强度分布直方图的相似性来判断两组生物样品的蛋白质的含量的相似性；

优选地，对各样品进行相关性分析；计算得到pearson系数并绘制相关性分析散点图，横、纵坐标分别为相比较的两次技术重复实验的归一化的强度值，图中点表示蛋白质。

6.根据权利要求1所述的筛选生物标志的方法，其特征在于，采用Benjamini-Hochberg法对p-value进行校正包括：将p-value按照升序排列，使用公式FDR＝(n/i)p-value计算每个p-value对应的FDR值，其中，n为p-value个数，i为p-value排序后的序号，如果某一个p-value所对应的FDR值大于排序的前一位p值所对应的FDR值，则放弃公式计算出来的FDR值，选用与它前一位相同的值；反之则保留计算的FDR值，最终得到的值记为q-value，筛选小于或等于显著性水平的q-value，筛选出的q-value对应的蛋白质为差异表达的蛋白；优选地，显著性水平为0.05。

7.根据权利要求6所述的筛选生物标志的方法，其特征在于，对筛选结果绘制火山图，火山图横坐标为生物标志类物质对应的归一化后的平均强度值，纵坐标为q-value取以10为底的对数再取负数，每个点表示一个生物标志类物质。

8.根据权利要求7所述的筛选生物标志的方法，其特征在于，还包括输出差异生物标志类物质的信息列表，列表信息包括生物标志类物质的名称、火山图横坐标、火山图纵坐标、q-value和是否差异表达。

9.根据权利要求1～8任一项所述的筛选生物标志的方法，其特征在于，所述生物标志类物质为蛋白质时，非标记定量分析方法是使用蛋白质定量软件进行非标记定量分析，导入数据库文件和含有蛋白质质谱数据的raw文件，数据依赖性采集方法将给定质荷比范围内的离子碎裂，提取母离子谱图，并在保留时间内积分得到峰面积，通过比较样品的峰面积完成定量。

10.根据权利要求9所述的筛选生物标志的方法，其特征在于，蛋白质的质谱数据是将液相色谱法分离得到的肽段通过通过质谱仪检测所有离子的质荷比信息，并记录其相应保留时间上的信号强度，最终得到结果序列数据，并记录在结果文件raw文件中，优选地，在用液相色谱法分离之前，采用胰蛋白酶将蛋白质酶解成肽段。