统计 | R语言执行两组间差异分析Wilcox秩和检验

Posted 2021-04-18 soilhealth

两组间差异的非参数检验之Wilcox秩和检验在R中实现

在进行两组数据间的差异分析时，我们通常会想到使用。但若数据不满足执行t检验的参数假设（例如数据分布不符合正态性，变量在本质上就严重偏倚或呈现有序关系），无法使用t检验分析时，可以考虑使用非参数的方法来完成。

就两组数据的比较而言，wilcox秩和检验（或称Mann-Whitney U检验）是常见的非参数检验方法之一。 本文简介怎样在R中进行wilcox秩和检验，以实现两组间非参数差异分析。

本文使用的作图数据的网盘链接（提取码o8lr）：

https://pan.baidu.com/s/1b-1INL4HFrsIOvs_0UfByw

文件“alpha.txt”为某16S细菌群落测序所获得的部分alpha多样性指数数据，包含3列信息： sample，样本名称； observed_species和shannon分别为两种类型的alpha多样性指数。 文件“group.txt”为各样本分组信息，第一列（sample）为各样本名称； 第二列（group）为各样本的分组信息。

以上使用的示例数据与前文“”中的数据一致。 已知group3的shannon指数数据分布并不符合正态性，此时，若我们想比较group2和group3的shannon指数间是否存在显著差异，就不适合使用t检验（暂且不考虑对数据进行合理的转化后是否会满足t检验的参数假设），可采用非参数的方法（本文中介绍使用wilcox秩和检验）去实现。

数据预处理及正态性假设检验

首先将上述两个数据表读入R中，并合并在一起，以及数据的正态分布检验。

library(reshape2)
 
#读入文件，合并分组信息，数据重排
alpha <- read.table('alpha.txt', sep = '	', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group  <- read.table('group.txt', sep = '	', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
alpha <- melt(merge(alpha, group, by = 'sample'), id = c('sample', 'group'))
 
#选择要比较的分组（此处查看 group1 与 group2 在 shannon 指数上是否存在显著差异）
shannon_23 <- subset(alpha, variable == 'shannon' & group %in% c('2', '3'))
shannon_23$group <- factor(shannon_23$group)
head(shannon_23, 10)
 
#Shapiro-Wilk 检验数据是否符合正态分布（发现不符合正态分布）
tapply(shannon_23$value, shannon_23$group, shapiro.test)

选取的数据框“shannon_23”内容如下所示。第一列（sample），两组数据中所含样本名称；第二列（group），两组分组名称，且分组列已转化为因子类型；第三列（variable），alpha多样性指数shannon指数；第四列（value），shannon指数的数值。

通过Shapiro-Wilk检验得知数据分布不满足正态性。这里p值小于0.05表明数据违背了正态性分布的零假设。

Wilcoxon检验

不符合正态性前提的数据，无法应用t检验去比较差异。我们考虑使用非参数的方法作为替代，对于两组数据的比较，可使用wilcoxon检验。类似于t检验，根据样本间是否独立，wilcoxon检验分为wilcox秩和检验以及wilcox符号秩和检验。

wilcox秩和检验

  

假设样本间是相互独立的，直接使用wilcox秩和检验去处理。它是独立样本t检验的一种非参数替代方法。

此处使用的wilcox.test()与t检验t.test()的参数很相似。wilcox_test()中默认两组间相互独立（默认参数paired = FALSE），执行独立样本的wilcox秩和检验；同时默认的备择假设是双侧的（默认参数alternative = 'two.sided'），即执行双侧检验，可分别使用“alternative = 'less'”或“alternative = 'greater'”执行单侧wilcox检验。

##wilcox 秩和检验，我们执行了一个双侧检验
wilcox_test <- wilcox.test(value~group, shannon_23, paired = FALSE, alternative = 'two.sided')
wilcox_test
wilcox_test$p.value

由于p值（约为0.003）小于0.05，即拒绝了原假设（原假设两组间没有差异），group2和group3的shannon指数间存在显著不同。

wilcox符号秩和检验

  

假设样本间并非相互独立的，可考虑wilcox符号秩和检验，它是非独立样本t检验的一种非参数替代方法。例如，非独立组设计（dependent groups design）、重复测量设计（repeated measures design）等。尽管此时你选用独立样本的wilcox秩和检验方法也是可行的，这种分析方法本身并没问题（仅仅是在统计算法上存在一些不同，相较之下可能wilcox符号秩和检验会更合适一些）。

此时在wilcox.test()中设定参数“paired = TRUE”即可执行wilcox符号秩和检验。

##wilcox 符号秩和检验，我们执行了一个双侧检验
wilcox_test <- wilcox.test(value~group, shannon_23, paired = TRUE, alternative = 'two.sided')
wilcox_test
wilcox_test$p.value

根据p值（0.039，低于0.05）可知group2和group3的shannon指数间存在显著不同。

可视化展示

  

考虑作图将两组差异进行可视化展示。例如，一个简单的箱线图示例。

#boxplot() 箱线图
boxplot(value~group, data = shannon_23, col = c('blue', 'orange'), ylab = 'Shannon', xlab = 'Group', main = 'wilcox test: p-value = 0.00295')

Wilcox秩和检验的一个批处理示例

相较于参数分析的t检验，wilcox秩和检验不必事先验证数据分布的正态性，因此理论上来讲，只要是两组数据间的差异分析均可使用wilcox秩和检验去完成，因此其适用范围相较于t检验也更广泛。在数据量较大的情况下（可能会存在部分数据满足t检验分析的条件，而另一部分数据则不满足，无法做到全部使用t检验去实现），可以考虑使用循环逐一挑选分组后，直接执行wilcox秩和检验进行各两两分组间的差异分析。尽管这种方法比较“粗暴”，但也不失为一种备选方案。

如下将展示一个批处理示例。

网盘附件中提供了另一示例数据集“gene.txt”。表格中每一行为一种基因，每一列为一个样本，交叉区域为各基因在各样本中的相对丰度。接下来，我们期望通过wilcox秩和检验，找到在group1和group2组中，具有丰度差异的基因。

##wilcox 检验批处理示例
library(doBy) #使用其中的 summaryBy() 以方便按分组计算均值、中位数
 
#读取数据
gene <- read.table('gene.txt', sep = '	', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '	', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
result <- NULL
 
#使用循环，逐一对各基因进行两组间 wilcox 秩和检验
for (n in 1:nrow(gene)) {
gene_n <- data.frame(t(gene[n,subset(group, group %in% c(1, 2))$sample]))
gene_id <- names(gene_n)[1]
names(gene_n)[1] <- 'gene'

gene_n$sample <- rownames(gene_n)
gene_n <- merge(gene_n, group, by = 'sample', all.x = TRUE)

gene_n$group <- factor(gene_n$group)
p_value <- wilcox.test(gene~group, gene_n)$p.value
if (!is.na(p_value) & p_value < 0.05) {
stat <- summaryBy(gene~group, gene_n, FUN = c(mean, median))
result <- rbind(result, c(gene_id, as.character(stat[1,1]), stat[1,2], stat[1,3], as.character(stat[2,1]), stat[2,2], stat[2,3], p_value))
}
}
 
#输出统计结果
result <- data.frame(result)
names(result) <- c('gene_id', 'group1', 'mean1', 'median1', 'group2', 'mean2', 'median2', 'p_value')
result$p_adjust <- p.adjust(result$p_value, method = 'BH') #推荐加个 p 值校正的过程
write.table(result, 'gene.wilcox.txt', sep = '	', row.names = FALSE, quote = FALSE)

我们主要输出这些结果： gene_id，基因id； group1和group2，分别为所需比较的分组1和分组2的名称； mean1、median1、mean2、median2，分别为各基因在分组1、2中的平均丰度以及中位数数值； p_value，显著性p值，此处仅输出了p值低于0.05的结果（即只保留相对丰度在group1、2中具有显著差异的基因）； p_adjust，同时通过Benjamini方法校正p值（嗯嗯，这里的数据p值校正后，没有差异基因……）。

特别说明

  

既然参数检验的前提条件有些苛刻，自己的数据不一定都满足参数分析的条件，那么以后需要用到组间差异比较时，直接全部使用非参数的检验就不可以了？

虽然对全部数据直接使用非参数的检验方式理论上没啥问题，但还是有点粗暴了一些。两种方法（此处比较了t检验和wilcox秩和检验）毕竟还是有差别的，非参数方法普遍没有参数方法严格。对于符合参数检验条件的数据来讲，使用参数检验还有可能会鉴别出非参数检验鉴别不到的差异，此时需要特别关注。例如，某数据符合t检验的条件，t检验的p值显著，但wilcox检验p值不显著，那么这时t检验的结果会相对可靠一些。