R语言执行两组间差异分析Wilcox秩和检验
在进行两组数据间的差异分析时，我们通常会想到使用t检验。但若数据不满足执行t检验的参数假设（例如数据分布不符合正态性，变量在本质上就严重偏倚或呈现有序关系），无法使用t检验分析时，可以考虑使用非参数的方法来完成。
数据预处理及正态性假设检验
首先将上述两个数据表读入R中，并合并在一起，以及数据的正态分布检验。library(reshape2)
#读入文件，合并分组信息，数据重排
alpha <- read.table('alpha.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', 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(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组中，具有丰度差异的基因。
library(doBy) #使用其中的 summaryBy() 以方便按分组计算均值、中位数
#读取数据
gene <- read.table('gene.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', 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 = '\t', 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检验的结果会相对可靠一些。
友情链接