嘿，针对你的微阵列基因表达数据问题，我来给你详细拆解一下：

一、t-test与Ebayes方法对微阵列数据的适用性 #

首先明确：t-test和Ebayes（尤其是limma包实现的Ebayes）都适用于微阵列数据，但两者的适用场景和统计效能有明显差异：

• 传统独立样本t-test：可以用，但你的每组样本量只有4个，这种小样本场景下，t-test的方差估计会非常不稳定，容易出现假阳性或假阴性结果，所以不是最优选择。
• Ebayes（limma包实现）：这是微阵列数据分析的黄金标准方法！它会利用所有基因的表达信息来收缩每个基因的方差估计，完美解决小样本下的方差不稳定问题，能显著提升统计检验的效能，更适配你的数据情况。

二、具体操作步骤（R语言） #

假设你已经通过read.table()或read.csv()将数据集读入R环境，命名为expr_data。

方法1：独立样本t-test（小样本慎用） #

1. 分组提取表达值并计算倍数变化

# 分离铜暴露组与对照组的表达列
cu_group <- expr_data[, 2:5]
control_group <- expr_data[, 6:9]

# 计算每个基因的log2倍数变化（Cu组均值 - Control组均值）
expr_data$log2FC <- rowMeans(cu_group) - rowMeans(control_group)

2. 逐基因执行t-test获取P值

# 自定义t-test函数（采用Welch t-test，不假设两组方差相等）
calc_t_pvalue <- function(cu_vals, ctrl_vals) {
  t_result <- t.test(cu_vals, ctrl_vals, alternative = "two.sided", var.equal = FALSE)
  return(t_result$p.value)
}

# 应用到每行基因
expr_data$p_value <- apply(cbind(cu_group, control_group), 1, function(row) {
  cu <- row[1:4]
  ctrl <- row[5:8]
  calc_t_pvalue(cu, ctrl)
})

方法2：Ebayes（limma包，推荐使用） #

limma是微阵列数据分析的经典工具包，Ebayes是其核心方差收缩方法。

1. 安装并加载limma包

if (!require("limma")) {
  install.packages("limma")
}
library(limma)

2. 构建实验设计矩阵与对比矩阵

# 定义分组因子（Control作为基准组）
group <- factor(c(rep("Cu", 4), rep("Control", 4)), levels = c("Control", "Cu"))

# 构建设计矩阵
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)

# 构建对比矩阵，指定比较Cu组 vs Control组
contrast_mat <- makeContrasts(Cu - Control, levels = design)

3. 拟合模型并执行Ebayes方差收缩

# 拟合线性模型
fit <- lmFit(expr_data[, 2:9], design)
# 应用对比矩阵
fit_contrast <- contrasts.fit(fit, contrast_mat)
# 执行Ebayes方差收缩
fit_ebayes <- eBayes(fit_contrast)

4. 提取结果到原数据集

# 提取所有基因的结果（不排序）
results <- topTable(fit_ebayes, number = nrow(expr_data), adjust = "fdr", sort.by = "none")

# 将log2倍数变化、P值、校正后P值合并到原数据
expr_data$log2FC <- results$logFC
expr_data$p_value <- results$P.Value
expr_data$adj_p_value <- results$adj.P.Val # FDR校正后的P值，火山图更常用

三、绘制火山图 #

用ggplot2包绘制标准火山图：

if (!require("ggplot2")) {
  install.packages("ggplot2")
}
library(ggplot2)

# 定义显著差异基因的阈值（可根据需求调整）
expr_data$significance <- ifelse(abs(expr_data$log2FC) > 0.5 & expr_data$adj_p_value < 0.05,
                                 "Significant", "Not Significant")

# 绘制火山图
ggplot(expr_data, aes(x = log2FC, y = -log10(adj_p_value), color = significance)) +
  geom_point(alpha = 0.7, size = 1.5) +
  scale_color_manual(values = c("Significant" = "red", "Not Significant" = "gray")) +
  # 添加阈值参考线
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "black") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
  labs(x = "log2(Fold Change)", y = "-log10(Adjusted P-value)",
       title = "Volcano Plot: Copper Exposure vs Control", color = "Significance") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

内容的提问来源于stack exchange，提问作者A.Carter