浅谈蛋白组的差异蛋白分析

蛋白质组学在现在众多组学中属于起步较晚的组学之一,一般跟代谢组并列提起,因为都是通过质谱仪进行检测(虽然代谢组不一定要用质谱)。我们通过质谱下机数据,然后经过软件进行图谱分析(查库),然后获得肽段/蛋白的丰度值,接下来可能就跟RNA-Seq等NGS数据一样,需要在蛋白表达谱(虽然由于仪器的限制,并未能检测到全部低丰度蛋白,鉴定到的蛋白数目还是处于N千数量级)中筛选出丰度发生显著变化的蛋白(这也就是常说的比较蛋白组学)

RNA-Seq相比蛋白组发展较早,从芯片到NGS测序,其差异分析流程经过众多学者的完善趋于成熟;而蛋白组,根据其检测方法的不同,外加一些谱图分析软件的差异,并未能向NGS那样有着统一的规范输出结果。从一些蛋白组应用相关的文献来看,还处于较为简单的分析,比如T检验;这里并不是说T检验不好,而是从现在普遍生物学重复只做3个的现象来看,T检验的效应还是不够的吧?

我之前在一篇博文Differential expression in proteomics中简单整理了下最近几年发表关于蛋白组差异分析方法的文献,我整理的还是有偏差的,我只对蛋白组Label-free技术相关文献做了整理,而iTRAQ技术的则未关注(这是因为我觉得iTRAQ实在有点贵,不如label-free应用范围广,特别是当临床上百样本的时候)。

一般来说,现在蛋白组原始数据的采集分为DDA( Data-dependent acquisition)和DIA( Data-independent acquisition),Label-free技术的蛋白定量分析属于DDA,差不多就是shutgun原理;质谱仪在数据采集时先采集MS1数据,再采集N次MS2数据;因此Label-free的定量分析时出现了两种策略:第一种是基于MS2数据进行定量,也就是文献中说的Spectra Count,其一般符合泊松分布,因其后续差异分析可以参照NGS的RNA-seq分析的一些方法,比如DESeq2,edgeR等;第二种是基于MS1数据进行定量,计算每个肽段在谱图上的积分,也就是先通过MS1数据识别肽段,然后图形积分计算丰度,其在一些文献中普遍认为符合正态分布(个人觉得蛮勉强的),因此后续分析一些文献建议类似于芯片分析的一些方法,如T检验,limma等。

除了上述几个思路外,在最近的一篇博文一个R包的故事中提到MSstats作者在对蛋白组定量数据处理时,采用了一种完全不同于上述方法的一种的思路,其整理了每个Protein Group下所有肽段的丰度值,当其对同一个Protein Group在不同样本处理下进行差异分析时,并不是以Protein Group的丰度作为数据(因为一些定量软件不仅会给出每个肽段的丰度值,并用其算法得出每个Protein Group的丰度),而是类似将同一个Protein Group下每个肽段的丰度值看做一次’取样’,比如如果某个Protein Group下,样本A有20个肽段,样本B有15个肽段,那么根据其logsum t-test算法,则T检验的自由度则为33。。是不是很独特的思路,其实有类似思路的文献,但见识的比较少。。

其实在做差异分析之前,还有一个步骤也是需要值得注意的,那就是需不需要进行空值填充,这是Label-free技术的缺陷所导致。一般在Label-free定量结果会看到不少蛋白在一些样本中的丰度值是缺失的,这很普遍。一些文献是采用过滤策略,比如大于50%样本有值的Protein Group保留(这在一些芯片分析过滤中比较常见);一些文献则会采用补空策略,比如用Multiple imputation算法等,但是由于label-free空值产生原因的不确定性,因此单一的补空算法无法准确对应到每个空值,因此有些研究者发现补空虽然能增加灵敏度,但会增加假阳性

下面我结合上述文献结论,以一个Label-free数据为例,分别使用T检验,Limma以及MSstats等方法进行蛋白组的差异分析

测试数据下载:labelfree.csv

首先做一步过滤处理,以50%原则过滤掉缺失值较多的Protein,这步主要针对后续T检验和limma的预处理,而MSstats则有其默认的预处理方式

library(ggplot2)
library(dplyr)

data <- read.csv(file = "labelfree.csv", sep = ",", header = T, row.names = 1)
rm <- apply(data, 1, function(x){
  sum(x == 0) > 3
})
df <- data[!rm,]
> nrow(df)
[1] 5325

接着则是做log转化,这步有没有必要呢?Label-free技术是用峰谱积分计算丰度的,其丰度的数量级在10的7次到10的11次不等,从下面每个样本的均值和标准差可看出,丰度值数据分布非常离散

Mean <- apply(df, 2, mean)
SD <- apply(df, 2, sd)
> Mean
    A1         A2         A3         B1         B2         B3 
2366243926 2351546722 2165457844 2175212680 2708233812 2290120209
> SD
     A1          A2          A3          B1          B2          B3 
13098498790 14392930144 14179115278 11151997873 12092333164 10344859562

那么做log转化有什么用呢,一般作用为消除异方差;减少自变量数量级不一致的情况;使非线性变量关系转化为线性关系;最主要的还是为了使数据的呈现方式接近我们所希望的假设,从而更好的进行统计推断(以上来自知乎)

下面左边是原始数据AB两组的散点图,右边是做了对数转化后的

DEP_analysis

那么后续差异分析是用log转化后的数据还是原始数据呢,对于limma和MSstats两个软件来说,是必须做对数转化(其实前者是默认已经做了log转化的数据作为输入,而后者则是根据参数来设定如何做log转化,总之两者的差异分析结果都是做了log转化后的数据所产生的)。但T检验则是用哪个数据作为输入呢,因为统计推断T检验是比较两个样本均值是否有差异,但对做了log转化后的数据做T检验则不是比较算数平均数了,而是几何平均数,至少还算有明确的含义?

引用一位小伙伴对对数转化的理解,我觉得蛮有道理的

T检验或者方差分析之类的有参检验需要数据符合方差齐性,一般随着实验的进行,由于仪器或者人工操作,误差有可能会越来越大,会导致数据波动越来越大。这个时候可以用对数变换把方差变大的趋势给调回来

T检验

在做Two-Sample T检验前,需要样本总体符合正态分布,并且方差齐性,以测试数据为例(未做对数转化),先做个F检验

fvalue <- apply(df, 1, function(x){
  a <- factor(c(rep("A",3), rep("B",3)))
  var.test(x~a)
})
>sum(unlist(lapply(fvalue, function(x){
  x$p.value > 0.05
})))
[1] 4605

可看出5325个蛋白有4605个是方差齐性的

接下来用t.test函数做成组T检验,双尾(默认),方差不齐性的用var.equal = FALSE

pvalue <- apply(df, 1, function(x){
  a <- factor(c(rep("A",3), rep("B",3)))
  fvalue <- var.test(x~a)
  if (fvalue$p.value > 0.05){
    t.test(x~a, var.equal = T)
  }else{
    t.test(x~a, var.equal = F)
  }
})

result_ttest <- data.frame(ID=row.names(df), 
                           Pvalue = as.numeric(unlist(lapply(pvalue, function(x) x$p.value))),
                           FC = log2(as.numeric(unlist(lapply(pvalue, function(x) x$estimate[1]/x$estimate[2]))))
                           )
> sum(result_ttest$Pvalue < 0.05 & abs(result_ttest$FC) > 1)
[1] 499

可看出差异蛋白数目为499个;如果是预处理时做了对数转化的话,代码同上,结果符合相同的阈值的差异蛋白有694个

Limma

Limma就比较熟悉了,如果就测试数据两组比较而言,limma需要先做log2转化,然后设定分组矩阵,最后用线性模型拟合再用Empirical Bayes test求出检验的P值

library(limma)

df1 <- log2(df + 1)
group_list <- factor(c(rep("A",3), rep("B",3)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(df1)

fit <- lmFit(df1, design)
fit <- eBayes(fit, trend=TRUE)
result_limma <- topTable(fit, coef=2,n=Inf)
>sum(result_limma$P.Value < 0.05 & abs(result_limma$logFC) > 1)
[1] 952

Limma的差异蛋白数目就相比以上就略多了,满足阈值的有952个

MSstats

我按照之前博文一个R包的故事中的代码(做了略微的修改),对测试数据的原始数据进行了分析(因为MSstats包可以对Maxquant结果的最原始数据开始处理),分析过程这里不列出了(主要代码有点长。。。),结果文件可看MSstats_DEP.csv,结果满足阈值的差异蛋白有513个

Summary

其实我为什么要挑选这几种方法呢,因为一篇发表于2017年Briefings in Bioinformatics期刊的一篇文章比较过上述3种方法外加Generalized linear model with a gamma distribution(GLM-Gamma)方法对label-free数据的差异蛋白分析的结果,最终结果是Limma的Empirical Bayes test在低样本重复的情况下相比其他优势非常明显,随着样本重复数的增加,优势才会有略微缩小(总之作者认为limma对现在的label-free数据结果来说,相比其他方法,最佳)

最后将上述对数转化的T检验的p值做个多重检验(FDR),跟Limma和MSstats的adj.pavlue(校正后的P值)做个VENN图看看交集如何

DEP_venn

从图中可看出,Limma和T检验的覆盖度较好,可能由于这两种分析方式过程以及思路比较近似吧(都是默认蛋白水平的丰度值都符合正态分布),但limma在一些地方处理的比T检验更加好点(可能是由于limma在了线性拟合的缘故?),比如当一个样本丰度均为0,而另一样本丰度各个重复的丰度值偏差较大,这种情况T检验的效用比较糟糕,FDR值能难达到阈值,但这种情况,我觉得可以认为这个Protein Group在两个样本中表达是有差异的!!!至于MSstats,我不确定是否是我没有完全了解其用法,使得这方法筛选出来的差异蛋白跟其他两种相差较大(理论上按照文献中所展示的,不应该是这样的,但是那文献也没给出其计算的具体过程。。。),根据其R包代码,其分析过程的思路确实跟其他的方法有点不同。。。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处

R语言之非参数检验篇
R语言中方差齐性检验丨数析学院
在统计学中为什么要对变量取对数