526互联

R：Wilcoxon秩和检验

发布时间 2023-11-01 17:15:30作者: 王哲MGG_AI

Wilcoxon秩和检验，也被称为Mann-Whitney U检验，是一种非参数统计检验方法，用于比较两组独立样本的中位数是否显著不同。这个检验适用于以下情况：

  1. 数据分布不满足正态分布假设：与某些参数统计检验（如t检验）不同，Wilcoxon秩和检验不要求数据满足正态分布假设。因此，它适用于那些数据分布不明确或偏离正态分布的情况。

  2. 方差齐性假设不成立：与一些其他统计检验（如方差分析）不同，Wilcoxon秩和检验也不要求样本的方差相等。

  3. 两组数据是独立的：这个检验通常用于比较两组不相关的样本，例如，两组不同的治疗组的病人数据。

具体的Wilcoxon秩和检验步骤如下：

  1. 将两组数据合并并按大小顺序排列，同时记录每个数据的来源组。

  2. 为每个数据点分配一个秩次，这些秩次是数据在合并数据中的相对大小顺序。如果有秩次相同的数据，则取其平均秩。

  3. 对于每组数据，计算其秩和，即将该组数据中的所有秩次相加。

  4. 比较两个秩和（通常称为U1和U2），并使用统计方法来评估它们是否显著不同。这个统计检验的结果会产生一个p值，该p值表示两组数据中的中位数是否存在显著差异。

如果p值小于显著性水平（通常设置为0.05），则可以拒绝原假设，认为两组数据的中位数存在显著差异。如果p值大于显著性水平，则不能拒绝原假设，即无法确定两组数据的中位数是否显著不同。

Wilcoxon秩和检验是一种有力的工具，适用于分析非正态分布的数据，特别是在小样本情况下。它被广泛用于生物统计学、社会科学、医学研究和其他领域的数据分析中。

# 设置工作目录到指定路径
setwd("C:\\Users\\Administrator\\Desktop")

# 加载doBy库，它提供了许多数据分析功能
library(doBy)

# 从文件'otu_table2.txt'中读取基因数据
gene <- read.table('otu_table2.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

# 从文件'metadata3-1.txt'中读取样本分组数据
group <- read.table('metadata3-1.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

# 初始化一个空的结果数据框
result <- NULL

# 遍历gene数据框的每一行
for (n in 1:nrow(gene)) {
  
  # 从gene数据框提取第n行的数据并转置，以获得gene_n
  gene_n <- data.frame(t(gene[n,]))
  
  # 获取基因的名称（列名）
  gene_id <- names(gene_n)[1]
  names(gene_n)[1] <- 'gene'
  
  # 添加样本名称列
  gene_n$sample <- rownames(gene_n)
  
  # 将gene_n与group数据框合并，基于sample列
  gene_n <- merge(gene_n, group, by = 'sample', all.x = TRUE)
  
  # 将group列转换为因子类型
  gene_n$group <- factor(gene_n$group)
  
  # 使用Wilcoxon检验测试基因的表达差异
  p_value <- wilcox.test(gene~group, gene_n)$p.value
  
  # 如果p值不是NA且小于0.05，则继续进行统计
  if (!is.na(p_value) & p_value < 0.05) {
    
    # 使用summaryBy函数计算每个组的基因表达均值和中值
    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转换为数据框并命名其列
result <- data.frame(result)
names(result) <- c('gene_id', 'group1', 'mean1', 'median1', 'group2', 'mean2', 'median2', 'p_value')

# 使用BH方法进行p值校正
result$p_adjust <- p.adjust(result$p_value, method = 'BH')

# 将结果数据框保存到文件'table.l5.relative-SE.wilcox.txt'
write.table(result, 'table.l5.relative-SE.wilcox.txt', sep = '\t', row.names = FALSE, quote = FALSE)

上面的代码输出的文件只保留显著相关，如果不想只保留显著相关，可以试试下面的代码

setwd("C:\\Users\\Administrator\\Desktop")
library(doBy)

gene <- read.table('otu_table2.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('metadata3-1.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

# 初始化结果数据框
result <- data.frame(
  gene_id = rownames(gene),
  group1 = NA,
  mean1 = NA,
  median1 = NA,
  group2 = NA,
  mean2 = NA,
  median2 = NA,
  p_value = NA,
  stringsAsFactors = FALSE
)

for (n in 1:nrow(gene)) {
  gene_n <- data.frame(t(gene[n,]))
  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
  
  stat <- summaryBy(gene~group, gene_n, FUN = c(mean, median))
  
  result[n, ] <- 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$p_adjust <- p.adjust(result$p_value, method = 'BH')

write.table(result, 'table.l5.relative-SE.wilcox.txt', sep = '\t', row.names = FALSE, quote = FALSE)