在转录组分析中，差异分析是必不可少的一步。那什么是差异分析呢？差异分析的结果又该怎么解读？以《GEO数据库转录组芯片数据处理与R分析：以GSE65682为例》一文中的数据集（GSE65682）为例，今天就让我们一起来深入了解一下，同时认识与它紧密相关的几个关键指标与图形，以及怎么做差异分析。
1. 转录组差异分析基础知识
1.1 什么是转录组？
简单来说，转录组是指一个活细胞在特定时间和环境下，所能转录出来的所有RNA的总和。因此，转录组分析就是为了了解基因的表达调控机制和功能。
1.2 差异分析是什么？为什么要做差异分析？
转录组差异分析是一种生物信息学分析方法，旨在比较不同样本或条件下转录组数据的差异，以揭示基因表达的变化。从所有基因中找出不同样本中表达具有差异的基因，了解这些差异基因如何影响生物体的生理功能和疾病发生发展。
1.3 差异分析结果解读
行名为基因名称，logFC表示基因在两个不同条件（例如处理组与对照组）之间的表达倍数的对数值。正值表示该基因为上调基因（处理组相对于对照组），负值表示该基因为下调基因。AveExpr表示该基因在所有样本中表达水平的平均值。PValue是统计学检验变量，代表差异显著性，一般pvalue小于0.05代表具有显著性差异，但可根据具体情况适当调整。adj.P.Val为校正后的P值，通过FDR校正调整原始P值，以控制多重比较中的假阳性率。校正后的P值较小的基因更有可能是真正的差异表达基因。singnificant代表差异显著性标志，down表示下调基因，up表示上调基因，no表示无显著性差异基因。
1.4 如何判断一个基因是否是差异基因呢？
通常会综合考虑logFC、pvalue或adj.P.Val这三个指标。一般来说，一个基因要被判定为差异基因，需要同时满足logFC符合设定的阈值（如logFC绝对值大于2）、pvalue小于0.05或adj.P.Val小于0.05。当然这些阈值要根据实际情况做相应的调整，常见的logFC会选择1，1.2，1.5等，pvalue常见的有0.05,0.01等。
注：pvalue和adj.P.Val选择一个就行了，从严格意义上讲选择adj.P.Val会更好，但是有时候选adj.P.Val的话，一个差异基因也筛选不出来，要根据实际情况做相应的变动
1.5 用火山图呈现差异结果
做完差异分析后通常要做火山图，来呈现差异分析的结果。
火山图是一种常用的可视化工具，它以logFC值为横坐标，以 -log10(padj)为纵坐标。在火山图中，每个点代表一个基因，横坐标表示基因表达的变化倍数，纵坐标表示差异的显著性。左上角蓝色区域为下调基因，右上角红色区域为上调基因，中间灰色区域为无显著性差异基因。通过火山图，我们可以快速直观地识别出哪些基因是显著差异基因，以及它们是上调还是下调。
2. 差异分析（基于limma）
2.1. 数据准备
我们将使用《GEO数据库转录组芯片数据处理与R分析：以GSE65682为例》一文中数据处理后的表达矩阵进行差异分析演示。
# 加载必要的包library(limma) library(dplyr)导入处理好的表达矩阵： expr_data <- read.csv("GSE65682_exp_unique.csv",row.names = 1)
表达矩阵如下图所示，每一行是一个基因，每一列是一个样本，中间数字的是基因表达量。
导入处理好的分组信息表：
group <- read.csv("GSE65682_group.csv") sample_group <- group %>%mutate(group = ifelse(group == "Alive", "Sepsis",ifelse(group == "Dead", "Sepsis","Healthy"))) group_info <- sample_group
样本分组数据包含样本名称和分组信息。第一列样品名称（sample），第二列为分组名称（group）。
提取表达矩阵和分组因子
# 从 group_info 中提取唯一的样本名称/标识符 sample_ids <- unique(group_info$sample)# 提取 expr_data 中对应的列 expr_data_subset <- expr_data[, colnames(expr_data) %in% sample_ids]# 从 group_info 中提取分组因子 group_factor <- factor(group_info$group[match(colnames(expr_data_subset), group_info$sample)])
根据匹配的表达矩阵和分组因子，进行 limma 分析
# 创建设计矩阵 design <- model.matrix(~0 + group_factor) colnames(design) <- levels(group_factor)# 将表达矩阵转换为 voom 对象 # v <- voom(expr_data_subset, design)# 拟合线性模型，这次使用正确的数值型矩阵 fit <- lmFit(expr_data_subset, design)# 设置对比（假设有两个组：group1 和 group2） cont.matrix <- makeContrasts(Healthy - Sepsis, levels = design)# 进行对比测试 fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2)# 获取结果 results <- topTable(fit2,n = Inf, adjust = "fdr")## 导出所有的差异结果 nrDEG = na.omit(results) ## 去掉数据中有NA的行或列 diffsig <- nrDEG # write.csv(diffsig, "all.limmaOut.csv") ## 我们使用|logFC| > 0.5，padj < 0.05（矫正后P值） foldChange = 0.5 padj = 0.05 ## 筛选出所有差异基因的结果 All_diffSig <- diffsig[(diffsig$adj.P.Val < padj & (diffsig$logFC>foldChange | diffsig$logFC < (-foldChange))),] write.csv(All_diffSig, "GSE65682_all.filterdegs.csv") ##输出差异基因数据集
在得到差异分析的结果后，我们接下来使用火山图进行可视化
# 加载必要的可视化包 library(ggplot2) library(ggrepel) library(ggthemes)# 准备火山图数据 volcano_data <- nrDEG volcano_data$gene <- rownames(volcano_data)# 添加差异表达状态列 volcano_data$singnificant <- "NO" volcano_data$singnificant[volcano_data$logFC > foldChange & volcano_data$adj.P.Val < padj] <- "UP" volcano_data$singnificant[volcano_data$logFC < -foldChange & volcano_data$adj.P.Val < padj] <- "DOWN"# 选择要标记的基因（前10个上调和下调基因） top_up <- volcano_data %>% filter(singnificant == "UP") %>% arrange(adj.P.Val) %>% head(5)top_down <- volcano_data %>% filter(singnificant == "DOWN") %>% arrange(adj.P.Val) %>% head(5)top_genes <- rbind(top_up, top_down)volcano_plot <- ggplot(volcano_data, aes(x = logFC, y = -log10(adj.P.Val), color = singnificant)) +geom_point(alpha = 0.6, size = 1.5) +scale_color_manual(values = c("DOWN" = "blue", "NO" = "grey", "UP" = "red"),labels = c("Down", "Not significant", "Up")) +geom_vline(xintercept = c(-foldChange, foldChange), linetype = "dashed", alpha = 0.5) +geom_hline(yintercept = -log10(padj), linetype = "dashed", alpha = 0.5) +geom_text_repel(data = top_genes, aes(label = gene), size = 3, max.overlaps = 20, show.legend = FALSE) +labs(x = "log2 Fold Change", y = "-log10(padj)") +theme_minimal() +theme(legend.position = "right")print(volcano_plot)# 保存火山图 ggsave("volcano_plot.png", volcano_plot, width = 10, height = 8)
火山图结果：
