代谢组学中差异代谢物的识别仍然是一个巨大的挑战，并在代谢组学数据分析中发挥着突出的作用。由于分析、实验和生物的模糊性，代谢组学数据集经常包含异常值，但目前可用的差异代谢物识别技术对异常值很敏感。作者这里提出了一种基于权重的具有稳健性火山图方法，助于从含有离群值的代谢组数据中更加准确鉴定差异代谢物。
作者将提出的方法与其它9种代谢组中常用的方法(包含t-test)进行比较，发现有较低的misclassification error rates(MERs)和较高的AUC，证明了该方法的优越性。
安装加载数据集
Github：https://github.com/nishithkumarpaul/Rvolcano
library(devtools)
install_github("nishithkumarpaul/Rvolcano")
library(Rvolcano)
data("DummyFullData")
dim(DummyFullData)
[1] 40 40
数据集由40行代谢物和40个样本组成，其中（癌症和健康组各20生物学重复）
计算foldchange
fcval <- foldChngCalc(DummyFullData,nSampG1 = 20,nSampG2 = 20)
head(fcval)
## 查看该函数
edit(foldChngCalc)
function (data, nSampG1, nSampG2)
{
g1fold <- apply((data[, 1:nSampG1]), 1, weightedMean)
g2fold <- apply((data[, (nSampG1 + 1):(nSampG1 + nSampG2)]),
1, weightedMean)
foldChange <- g2fold - g1fold
return(foldChange)
}
我们查看下foldChngCalc该函数,可以看到主要是通过apply函数计算了weightedMean值，即加权平均数, 而不是我们常用的mean(算术平均数)。
算术平均数又称均值，是统计学中基本的平均指标，就是简单的把所有数加起来然后除以个数。
加权平均数即加权平均值，是将各数值乘以相应的权数，然后加总求和得到总体值，再除以总的单位数。
通过一个例子计算下，观察weightedMean函数和mean有啥区别，,我们创建一个正态分布的向量集x和一个含有较高离群值的向量集，分布计算一下结果。
> set.seed(123)
> x<-c(rnorm(200,3,1))
> x1<-c(rnorm(180,3,1),rnorm(20,60,3)) #contain 20 outliers#
> weightedMean(x)
[1] 2.957601
> mean(x)
[1] 2.99143
> weightedMean(x1)
[1] 3.062994
> mean(x1)
[1] 8.712792
## 查看该函数的计算方法
function (a, b = 0.2)
{
w <- NULL
pw <- NULL
for (i in 1:length(a)) {
w[i] <- exp(-(b/2) * (mad(a)^-2) * (a[i] - median(a))^2)
pw[i] <- w[i] * a[i]
}
tpw <- sum(pw)
tw <- sum(w)
wm <- tpw/tw
dd <- c(tpw, tw, wm)
return(wm)
}
结果发现，不含离群值时，加权与算术平均数都一样，而含有离群值时候，结果容易形成误差，算术平均数易受极端数据的影响，极端值的出现，会使平均数的真实性受到干扰，加权差异有助于在进行计算时考虑更多数据，以便获得最准确的结果。
计算Pvalue
使用p.valcalc函数进行差异分析，使用robust vesion的t-test获取p值，公式可以查看函数源码，用到了weightedVar加权方差。
## 添加Pvalue
pval <- NULL
for (i in 1:dim(DummyFullData)[1]){
pval[i] <- p.valcalc(DummyFullData[i,1:20],DummyFullData[i,21:40])
}
##查看该源码
edit(p.valcalc)
function (x, y)
{
t.stat <- (weightedMean(x) - weightedMean(y))/sqrt(((length(x) -
1) * weightedVar(x) + (length(y) - 1) * weightedVar(y)) *
(1/(length(x) + length(y) - 2)) * ((1/length(x)) + 1/length(y)))
pval <- 2 * pt(abs(t.stat), (length(x) + length(y) - 2),
lower.tail = FALSE)
return(pval)
}
绘制火山图
log2foldchange和Pvalue都有了，直接使用RobVolPlot函数绘图
##Robust volcano plot
RobVolPlot(fcval,pval,cexcutoff = 0.7,cexlab = 0.8,
main = 'Volcano Plot',
xlab = 'log2fc',
ylab = '-log(Pvalue)',
plimit = 0.05,fclimit = 2)
edit(RobVolPlot)
ggplot2绘图
已经得到了log2foldchange和Pvalue数据了，我们就可以自己拿着数据绘制图形了
library(ggrepel)
## 获取数据
result <- data.frame(log2foldchange = fcval,p_value = pval)
result$threshold = factor(ifelse(result$p_value < 0.05 & abs(result$log2foldchange) >= 1,
ifelse(result$log2foldchange>= 1 ,'Up','Down'),'NoSignificance'),
levels=c('Up','Down','NoSignificance'))
##差异代谢物
re <- subset(result,p_value < 0.05)
ggplot(result,aes(log2foldchange, -log10(p_value)))+
geom_point(aes(color = threshold), size = 5) +
scale_colour_manual(values = c('#fb81f5', '#6868ed', 'gray40'), labels = c('Enriched metabolites', 'Depleted metabolites', 'No diff')) +
# 添加x轴标签
xlab(bquote(Log[2] ~ '(fold change)')) +
# 添加y轴标签
ylab(bquote(-Log[10] ~ italic('Pvalue')))+
# 添加标签：
geom_text_repel(data = re, aes(color = threshold,label = rownames(re)),max.overlaps = 20,
size = 4.8 ,fontface = "bold",arrow = arrow(length=unit(0.01, "npc")),force = 1, max.iter = 3e3,
box.padding = unit(0.6, 'lines'), segment.color = 'black', show.legend = FALSE)+
theme_bw()+
# 调整主题和图例位置：
theme(panel.grid = element_blank(),
legend.position = c(0.01,0.25),
legend.justification = c(0,1)
)+
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey",size = 1)+
geom_vline(xintercept = c(-1.2,1.2), linetype = "dashed", color = "grey",size = 1)
参考资料
