差异分析、显著性标记及统计作图的自动实现R代码示例

差异分析、显著性标记及统计作图的自动实现R代码示例

前篇对差异分析中，标注显著性水平“*”或“abc”的方法概念作了简要的描述。同时留了一个思考问题，自写一个R代码，将多重比较+统计作图这两步串联起来，实现差异分析、显著性水平的标注以及统计图的自动完成。白鱼小编不才，今天给大家分享两种实现过程。

两种方法一种对应了参数检验，单因素ANOVA+ Tukey HSD多重比较；另一种对应了非参数的方法，Kruskal-Wallis检验+非参数的多重比较Behrens-Fisher检验。作图数据和R脚本的网盘链接（提取码pmbr）：https://pan.baidu.com/s/12ejK5gu1m3-C04_WQO1qsA示例数据“test_data.txt”：sample，样本名称：group，各分组变量；value，各因变量。接下来，我们期望基于各组比较各因变量在组间的差异水平，并作图直观地展示出来。

单因素ANOVA + Tukey HSD

如果数据满足方差分析的前提假设，正态性、方差齐性等，那么方差分析肯定就是首选了。先执行单因素 ANOVA 比较整体差异，再执行多重比较查看两两差异，这个方法是公认的。

以下自写R代码整合了单因素ANOVA（整体的多组间比较）、Tukey HSD多重比较（事后两两比较）以及ggplot2作图，便于自动完成一系列的需要。#构建函数，统计各组均值、标准差，执行 ANOVA 以及 Tukey HSD 获得显著性
library(reshape2)
library(multcomp)
library(ggplot2)

aov_tukey <- function(data, groups, values, p = 0.05) {
stat_anova <- NULL
abc_list <- NULL

#统计检验
for (i in groups) {
for (j in values) {

#单因素 ANOVA，整体差异
dat <- data[c(i, j)]
names(dat) <- c('class', 'var')
dat$class <- factor(dat$class)
fit <- aov(var~class, dat)
p_value <- summary(fit)[[1]][1,5]

#单因素 ANOVA 结果整理
if (p_value < 0.001) sig <- '***'
else if (p_value >= 0.001 & p_value < 0.01) sig <- '**'
else if (p_value >= 0.01 & p_value < 0.05) sig <- '*'
else sig <- ''
stat_anova <- rbind(stat_anova, c(paste(i, j, sep = '/'), p_value, sig))

#Tukey HSD 检验（multcomp 包），多重比较
tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(class = 'Tukey')), sig = p, decreasing = TRUE)

#cld() 自动得到了显著性“abc”水平，提取显著性标记“abc”
sig <- data.frame(tuk$mcletters$Letters, stringsAsFactors = FALSE)
names(sig) <- 'sig'
sig$class <- rownames(sig)
sig$var <- j

#均值标准差统计
abc_ij <- cbind(aggregate(dat$var, by = list(dat$class), FUN = mean), aggregate(dat$var, by = list(dat$class), FUN = sd)[2])
names(abc_ij) <- c('class', 'mean', 'sd')
abc_ij$class <- as.character(abc_ij$class)

#合并结果
abc_ij$group <- i
abc_ij <- merge(abc_ij, sig, by = 'class')
abc_ij <- abc_ij[c(4, 6, 1, 2, 3, 5)]
abc_list <- rbind(abc_list, abc_ij)
}
}

#ggplot2 作图，柱状图
plot_bar <- ggplot(data = abc_list, aes(x = class, y = mean)) +
geom_col(aes(fill = class), color = NA, show.legend = FALSE) +
geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
facet_grid(var~group, scale = 'free' , space = 'free_x') +
geom_text(aes(label = sig, y = mean + sd + 0.3*(mean+sd))) +
labs(x = '', y = '')

#ggplot2 作图，箱线图
data <- melt(data[c(groups, values)], id = groups)
data <- melt(data, id = c('variable', 'value'))
data <- data[c(3, 1, 4, 2)]
names(data) <- c('group', 'var', 'class', 'value')

value_max <- aggregate(data$value, by = list(data$group, data$var), FUN = max)
names(value_max) <- c('group', 'var', 'value')
value_max <- merge(value_max, abc_list[c(-4, -5)], by = c('group', 'var'))

plot_box <- ggplot(data = data, aes(x = class, y = value)) +
geom_boxplot(aes(fill = class), show.legend = FALSE) +
facet_grid(var~group, scale = 'free' , space = 'free_x') +
geom_text(data = value_max, aes(label = sig, y = value + 0.3*value)) +
labs(x = '', y = '')

#return
stat_anova <- data.frame(stat_anova, stringsAsFactors = FALSE)
names(stat_anova) <- c('group', 'anova_pvalue', 'sig')
stat_anova$anova_pvalue <- as.numeric(stat_anova$anova_pvalue)

list(stat = list(anova = stat_anova, tukey = abc_list), plot = list(plot_bar = plot_bar, plot_box = plot_box))
}好了，接下来我们调用以上命令执行。这里由于只是演示，所以暂且忽略示例数据本身是否满足方差分析的前提假设，直接进入方差分析及多重比较过程。实际的数据分析中，切记一定要考虑这点的。##调用函数运行
dat <- read.table('test_data.txt', header = TRUE)

#groups 指定分组列，可指定多列；values 指定变量列，可指定多列；p 指定显著性水平
#这里作为测试，就先假设这些数据满足方差分析的前提条件了
stat_result <- aov_tukey(data = dat, groups = c('group1', 'group2'), values = c('value1', 'value2', 'value3', 'value4'), p = 0.05)

#差异分析统计结果
stat_result$stat$anova
stat_result$stat$tukey

#显著性标记“abc”作图结果
stat_result$plot$plot_bar
stat_result$plot$plot_box多组间方差分析统计，以及事后两两的多重比较结果，已经将显著性水平“*”和“abc”标注其中。和作图结果一起，如下所示。

Kruskal-Wallis + Behrens-Fisher

当原始数据不满足方差分析的条件时，可以考虑转化数据（如log转换），看转化后的数据是否满足。

或者更换为非参数的方法，这里展示一个针对于非参数检验的方法示例，先执行Kruskal-Wallis检验比较整体差异，再执行Behrens-Fisher的非参数多重比较查看两两差异。在下文开始前，这里允许我吐槽一件事。本人的很多经验学自《R语言实战 第二版》，它的154页有这一段话。因此，对于非参数的检验方法，我也是一直在使用wilcoxon检验+p值校正的方法，代替多重比较。然后，还是大神师妹提醒的（就在今天，你们信吗？真的太及时了，非常感谢），《R语言实战 第一版》，它的154页居然是这样的…..好吧，原来非参数的多重比较有现成的方法（本篇中使用的“Behrens-Fisher”就是npmc()里的一种方法，还有另一种“Steel”，详见函数帮助），看到这一幕的我差点吐血……那至于为什么第二版把这个函数删除了呢？可能是npmc包在3.0及以后版本的R中已经不更新的原因吧，所以见不着了。既然如此，考虑到之前对非参数检验中，使用wilcoxon检验+p值校正来代替多重比较的方法已经做了描述，因此本篇就不再讨论它（事实上，我已经写好了，然后经师妹一提醒，就赶紧修改了）。接下来，我们就使用npmc()函数，完成非参数的多重比较过程。以下自写R代码整合了Kruskal-Wallis检验（整体的多组间比较）、Behrens-Fisher非参数多重比较（事后两两比较）以及ggplot2作图，便于自动完成一系列的需要。这里Behrens-Fisher检验得到的是两组间的p值，对于显著性“abc”的标注，尝试通过p值再结合中位数水平手动判断，大体方法如下（参考前述）。（1）首先根据中位数大小，将各组由高往低排排序，中位数最高的组标注为“a”；（2）将中位数最高的组与第二高的组相比，若差异显著，则第二组标注为“b”；若不显著，继续比较其与中位数第三高的组的差异；（3）若中位数最高的组与第二高的组不显著，中位数第二高的组与第三高的组显著，则第二高的组就同样标注为“a”，第三高的组标注为“b”；若中位数最高的组与第二高的组不显著、中位数第二高的组与第三高的组不显著，但中位数最高的组与第三高的组显著，则第二高的组就标注为“ab”，第三高的组标注为“b”；（4）然后以标注为“b”的组的中位数为标准，以此类推，继续循环往后比较，直到最小中位数的组被标记，且比较完毕为止。不过总感觉这种方法有点怪怪的（不确定中位数在这里是否有代表性），所以，下面的代码虽可供参考，但实际应用时慎用！！#构建函数，统计各组中位数，执行 Kruskal-Wallis 以及 Behrens-Fisher 获得显著性，并手动基于中位数高低和 p 值判断“abc”水平
library(reshape2)
library(npmc)
library(ggplot2)

krus_BF <- function(data, groups, values, p = 0.05) {
stat_kruskal <- NULL
stat_BF <- NULL
abc_list <- NULL

#统计检验
for (i in groups) {
for (j in values) {

#Kruskal-Wallis 检验，整体差异
dat <- data[c(i, j)]
names(dat) <- c('class', 'var')
dat$class <- factor(dat$class)
fit <- kruskal.test(var~class, dat)
p_value <- fit$p.value

#获取 Kruskal-Wallis 检验 p 值
if (p_value < 0.001) sig <- '***'
else if (p_value >= 0.001 & p_value < 0.01) sig <- '**'
else if (p_value >= 0.01 & p_value < 0.05) sig <- '*'
else sig <- ''
stat_kruskal <- rbind(stat_kruskal, c(paste(i, j, sep = '/'), p_value, sig))

#Behrens-Fisher 检验（npmc 包），非参数多重比较
BF_test <- npmc(dat, df = 2, alpha = p)
BF <- BF_test$test$BF
info <- BF_test$info

BF$group <- paste(i, j, sep = '/')
for (l in 1:nrow(BF)) {
x <- strsplit(as.vector(BF[l,'cmp']), '-')[[1]]
BF[l,'class1'] <- rownames(info)[which(info$group.index == x[1])]
BF[l,'class2'] <- rownames(info)[which(info$group.index == x[2])]
}
stat_BF <- rbind(stat_BF, BF[c(12:14, 2:11)])

#Behrens-Fisher 检验 p 值矩阵获取
stat_p <- BF[c(13, 14, 10)]
stat_p1 <- stat_p[c(2, 1, 3)]
names(stat_p1)[1:2] <- c('class1', 'class2')
stat_p <- rbind(stat_p, stat_p1)
stat_p <- dcast(stat_p, class1~class2, value.var = 'p.value.2s')
rownames(stat_p) <- stat_p$class1

#计算各组中位数
abc_ij <- aggregate(dat$var, by = list(dat$class), FUN = mean)
abc_ij$group <- i
abc_ij$var <- j
abc_ij <- abc_ij[c(3, 4, 1, 2)]
names(abc_ij) <- c('group', 'var', 'class', 'median')

#根据中位数和 p 值，手动添加显著性标记“abc”
abc_ij <- abc_ij[order(abc_ij$median, decreasing = TRUE), ]
class <- as.vector(abc_ij$class)

m = 1; n = 2; l = 1
abc_ij[m,'sig'] <- letters[l]
while (n < length(class)) {
for (n in (m+1):length(class)) {
if (stat_p[class[m],class[n]] < p) {
abc_ij[n,'sig'] <- letters[l+1]
if (n - m != 1) {
for (x in (m+1):(n-1)) {
if (stat_p[class[x],class[n]] < p) abc_ij[x,'sig'] <- letters[l]
else abc_ij[x,'sig'] <- paste(letters[l], letters[l+1], sep = '')
}
}
l = l + 1
m = n
break
}
}
}
abc_ij[(m:n),'sig'] <- letters[l]

abc_list <- rbind(abc_list, abc_ij)
}
}

#ggplot2 作图，箱线图
data <- melt(data[c(groups, values)], id = groups)
data <- melt(data, id = c('variable', 'value'))
data <- data[c(3, 1, 4, 2)]
names(data) <- c('group', 'var', 'class', 'value')

value_max <- aggregate(data$value, by = list(data$group, data$var), FUN = max)
names(value_max) <- c('group', 'var', 'value')
value_max <- merge(value_max, abc_list[-4], by = c('group', 'var'))

plot_box <- ggplot(data = data, aes(x = class, y = value)) +
geom_boxplot(aes(fill = class), show.legend = FALSE) +
facet_grid(var~group, scale = 'free' , space = 'free_x') +
geom_text(data = value_max, aes(label = sig, y = value + 0.3*value)) +
labs(x = '', y = '')

#return
stat_kruskal <- data.frame(stat_kruskal, stringsAsFactors = FALSE)
names(stat_kruskal) <- c('group', 'kruskal_pvalue', 'sig')
stat_kruskal$kruskal_pvalue <- as.numeric(stat_kruskal$kruskal_pvalue)

stat_BF[which(stat_BF$p.value.2s < 0.001),'sig.2s'] <- '***'
stat_BF[which(stat_BF$p.value.2s >= 0.001 & stat_BF$p.value.2s < 0.01),'sig.2s'] <- '**'
stat_BF[which(stat_BF$p.value.2s >= 0.01 & stat_BF$p.value.2s < 0.05),'sig.2s'] <- '*'
stat_BF[which(stat_BF$p.value.2s >= 0.05),'sig.2s'] <- ''

list(data = value_max, stat = list(kruskal = stat_kruskal, BF_p = stat_BF, BF_abc = abc_list), plot = plot_box)
}接下来我们调用以上命令执行。和上述示例数据一致。
##调用函数运行
dat <- read.table('test_data.txt', header = TRUE)

#groups 指定分组列，可指定多列；values 指定变量列，可指定多列；p 指定显著性水平
stat_result <- krus_BF(data = dat, groups = c('group1', 'group2'), values = c('value1', 'value2', 'value3', 'value4'), p = 0.05)

#差异分析统计结果
stat_result$stat$kruskal
stat_result$stat$BF_p
stat_result$stat$BF_abc

#显著性标记“abc”作图结果
stat_result$plot整体的非参数检验（Kruskal-Wallis检验），已经将显著性水平“*”标注其中。Behrens-Fisher非参数多重比较结果，已经将显著性水平“*”标注其中。对于各列数值的含义，可通过?npmc查看帮助文档。这个是根据中位数的高低和Behrens-Fisher的p值，手动标注的“abc”水平。以及箱线图。友情链接