setwd("E:\\20220927宏基因组教学\\02后期分析\\05willcox")
library(doBy)
gene <- read.table('table.l5.relative-SE.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
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
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, 'table.l5.relative-SE.wilcox.txt', sep = '\t', row.names = FALSE, quote = FALSE)
#############################################################################################
这段代码的作用是进行Wilcoxon秩和检验，以比较两组样本中的基因表达差异，并输出显著差异的基因列表。
- 首先使用setwd()函数将当前工作目录设置为指定路径。
- 使用read.table()函数读取文件table.l5.relative-SE.txt和group.txt，其中table.l5.relative-SE.txt包含基因表达数据，group.txt包含每个样本所属的组别信息。
- 定义一个空数据框result，用于存储显著差异的基因。
- 使用for循环遍历每行基因表达数据。
- 将当前行的基因表达数据转换为数据框gene_n，并将第一列（基因名称）的变量名设置为‘gene’。
- 将样本名作为一列变量加入gene_n，并将gene_n与group数据框根据样本名进行合并，得到每个样本所属的组别信息。
- 将组别信息转换为因子变量（factor）。
- 使用wilcox.test()函数进行Wilcoxon秩和检验，计算两组样本之间的p值。
- 如果p值<0.05，则使用summaryBy()函数计算两组样本的基因表达平均值和中位数，并将结果存储到stat数据框中。
- 将基因名称和stat中的数据组合为一行，加入到结果数据框result中。
- 循环结束后，将result转换为数据框，并设置列名。
- 使用p.adjust()函数对p值进行多重检验校正，得到p_adjust列。
- 使用write.table()函数将结果数据框写入文件table.l5.relative-SE.wilcox.txt中。