## Significance A and B for protein ratios

Significance A/B最先是发表于2008年Nature Biotechnology期刊上，MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification，这篇文章主要是介绍Maxquant这款用于蛋白组定量分析软件的，非常有名，而其附录中作者提到了如何通过protein ratio来计算显著性（P值），附录PDF

The assumption that the ratio distribution is normal in a case where no differential regulation takes place is reasonable, since the protein ratios are obtained as averages of many SILAC peptide ratios. In the limit of a large number of SILAC peptides per protein a normal distribution can be assumed due to the central limit theorem

#### Significance A

• The base of the logarithm is irrelevant for our calculation，对每个ratio作log转换，至于log2还是log10都可以，不影响结果
• To make a robust and asymmetrical estimate of the standard deviation of the main distribution we calculate the 15.87, 50, and 84.13 percentiles r-1, r0, and r1，为了估计一个健壮以及不对称的标准差，作者先计算了15.87，50，84.34三个百分位的值，然后根据下述公式做标准化获得一个Z值，其中`r1-r0``r0-r-1`则是分布左右两边的标准差。因此当ratio值大于ro时，标准差则为`r1-r0`；而当ratio值小于ro时，标准差则为`r0-r-1`，具体公式如下图所示：

• Under the null hypothesis of normal tails of the log ratio distribution the probability of obtaining a value this large or larger is，将ratio做标准化转化后，这时的Z值则符合标准正态分布了，那么其对应的概率计算则就变得很简单了，类似于标准的假设检验方法了，按照作者意思，计算在正态分布中大于Z值的概率（标准正态分布曲线下大于Z值的那部分面积），可以理解为假设检验中的单尾检验？具体公式如下，其实就是个标准正态分布的概率累计函数

#### Significance B

Significance B的计算过程类似于Significance A，或者说是改进版？

• 作者发现那些高丰度蛋白相比低丰度蛋白，其中未发生变化的蛋白的ratio分布更加集中，如下图所示，那么高丰度蛋白的ratio也就相对会更加显著

• 考虑到这个丰度对结果影响，作者在Significance A算法的基础上，将protein根据intensity分成N个bin，每个bin至少包含300个蛋白，最后对每个bin进行一遍Significance A计算，然后将结果整合在一起

#### 代码实现

``````get_significance <- function(ratio){
ratio <- log2(as.numeric(ratio))
order_ratio <- ratio[order(ratio)]
quantiletmp <- quantile(order_ratio, c(0.1587,0.5,0.8413))
rl <- as.numeric(quantiletmp[1])      #对应公式中的r-1
rm <- as.numeric(quantiletmp[2])      #对应公式中的r0
rh <- as.numeric(quantiletmp[3])      #对应公式中的r1
p <- unlist(lapply(ratio, function(x){
if (x > rm){
z <- (x-rm)/(rh-rm)
pnorm(z,lower.tail = F)
}else{
z <- (rm-x)/(rm-rl)
pnorm(z,lower.tail = F)
}
}))
}
``````

``````p <- get_significance(data)
``````