作者将提出的方法与其它9种代谢组中常用的方法(包含t-test)进行比较，发现有较低的misclassification error rates(MERs)和较高的AUC，证明了该方法的优越性。
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)
参考资料
