引言
转录组测序（RNA-seq）已成为现代生物学研究中不可或缺的技术手段，而差异表达分析则是其中最关键的分析环节之一。
本文将基于R语言的DESeq2包，通过airway数据集实战演示「从原始count数据到结果解读」的完整分析流程，重点解决分析过程中常见的「基因注释、数据预处理、结果可视化」等关键问题。
技术亮点
本教程将深入讲解以下核心技术：
「核心分析包」：DESeq2、airway、ggplot2、pheatmap
「关键函数」：DESeqDataSetFromMatrix()、DESeq()、results()、vst()
「高级技巧」：基因ID转换策略、重复基因处理、独立过滤原理
「可视化方法」：火山图、热图、PCA图的专业绘制与解读
完整分析流程
1. 环境准备与数据加载
首先确保所有必要的R包已安装加载。这里采用「自动化检查安装」的方法，提高代码可重复性：
# 检查并安装必要的Bioconductor包
if (!require("BiocManager", quietly = TRUE)){
install.packages("BiocManager")
}
# 自动检查、安装、加载所需包
required_packages <- c("DESeq2", "airway", "ggplot2", "pheatmap", "ggrepel",
"dplyr", "tidyr", "reshape2", "AnnotationDbi", "org.Hs.eg.db", "biomaRt")
for (pkg in required_packages) {
if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
message("Installing package: ", pkg)
BiocManager::install(pkg, update = FALSE)
if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
stop("Failed to install package: ", pkg)
}
}
}
加载airway数据集并查看基本信息：
# 加载airway数据集
data("airway") # RangedSummarizedExperiment对象是Bioconductor生态系统中存储基因组数据的标准容器，包含三个核心组件：assay（计数矩阵）、rowData（基因信息）和colData（样本信息）。
raw_counts <- assay(airway) # 提取计数矩阵
col_data <- colData(airway) # 提取样本信息
# 查看数据基本信息
cat("计数矩阵维度:", dim(raw_counts), "\n")
cat("样本分组情况:\n")
print(table(col_data$dex))
airway数据集记录了地塞米松处理对人气道平滑肌细胞的基因表达影响，包含4个处理组和4个对照组。
2. 基因注释与ID转换
原始数据通常使用Ensembl ID作为基因标识，但我们需要更易读的Gene Symbol进行后续分析和结果解读。
离线转换方法（备选）
# 获取原始基因ID（Ensembl ID）
ensembl_ids <- rownames(raw_counts)
# 使用org.Hs.eg.db包进行ID转换
convert_ensembl_to_symbol <- function(ensembl_ids) {
# 去除Ensembl ID的版本号。Ensembl ID中的版本信息（如ENSG00000000003.15）会频繁更新，而去除版本号后使用稳定ID（ENSG00000000003）能提高注释成功率。
clean_ids <- sub("\\..*", "", ensembl_ids)
# 使用mapIds进行ID转换
symbol_map <- mapIds(org.Hs.eg.db,
keys = clean_ids,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first") # 多重映射时取第一个
return(symbol_map)
}
# 执行ID转换
offline_symbols <- convert_ensembl_to_symbol(ensembl_ids)
multiVals = "first"参数处理一个Ensembl ID对应多个Gene Symbol的情况，虽然可能不是最优选择，但在大多数情况下是可接受的折衷方案。
在线转换方法（推荐）
convert_ensembl_to_symbol_online <- function(ensembl_ids) {
# 连接到Ensembl数据库获取最新注释
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# 去除版本号
clean_ids <- sub("\\..*", "", ensembl_ids)
# 获取基因Symbol
results <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = clean_ids,
mart = mart)
# 创建映射
symbol_map <- setNames(results$hgnc_symbol, results$ensembl_gene_id)
# 将原始ID（带版本号）映射到Symbol
final_map <- symbol_map[sub("\\..*", "", ensembl_ids)]
names(final_map) <- ensembl_ids
return(final_map)
}
# 执行ID转换
online_symbols <- convert_ensembl_to_symbol_online(ensembl_ids)
离线方法使用本地注释数据库，速度更快且可重复性好。在线方法能获取最新注释但依赖网络。实际分析中建议先尝试在线方法。
3. 基因去重处理
一个常见的挑战是多个Ensembl ID可能对应同一个Gene Symbol，需要合理处理：
# 创建基因注释数据框
gene_annotation <- data.frame(
ensembl_id = ensembl_ids,
symbol = online_symbols,
stringsAsFactors = FALSE
)
# 去除没有对应Symbol的基因
annotated_genes <- gene_annotation[!is.na(gene_annotation$symbol), ]
# 处理重复Symbol：保留每个Symbol中表达量最高的Ensembl ID
resolve_duplicate_symbols <- function(count_matrix, gene_annotation) {
# 计算每个基因的平均表达量
gene_means <- rowMeans(count_matrix)
# 合并表达量信息并去重
unique_genes <- data.frame(
ensembl_id = rownames(count_matrix),
mean_expression = gene_means,
stringsAsFactors = FALSE
) %>%
inner_join(gene_annotation, by = "ensembl_id") %>%
filter(!is.na(symbol)) %>%
group_by(symbol) %>%
arrange(desc(mean_expression)) %>% # 按表达量降序排列
slice(1) %>% # 取每个组的第一个（表达量最高的）
ungroup()
return(unique_genes)
}
# 执行去重
unique_genes <- resolve_duplicate_symbols(raw_counts, gene_annotation)
dedup_raw_counts <- raw_counts[unique_genes$ensembl_id, ]
rownames(dedup_raw_counts) <- unique_genes$symbol
选择表达量最高的转录本通常能够捕获该基因最主要的活性形式，这在大多数生物学场景下是合理的选择。其他策略包括取平均值等，应根据研究目的选择。
4. 数据预处理与质控
创建DESeq2数据对象
# 从计数矩阵创建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(
countData = dedup_raw_counts,
colData = col_data,
design = ~ dex # 指定实验设计。`~ dex`表示我们要研究`dex`变量（处理条件）对基因表达的影响，这是后续统计检验的基础。如果存在批次效应，可以在design公式中加入批次变量，如~ batch + dex。
)
过滤低表达基因
# 过滤低表达基因
keep <- rowSums(counts(dds) >= 10) >= 3 # 我们的过滤标准（至少3个样本中计数≥10）在统计效能和基因保留之间取得了良好平衡
dds <- dds[keep,]
cat("过滤后保留", nrow(dds), "个基因\n")
「为什么要过滤低表达基因？」
「提高统计可靠性」：低计数基因方差估计不可靠，会影响离散度估计
「减少多重检验负担」：降低假阳性率
「减少技术噪音」：避免将技术误差误认为生物学信号
数据标准化与变换
# 方差稳定变换
vsd <- vst(dds, blind = FALSE) # vst（方差稳定变换）：计算速度更快，特别适用于大样本量。rlog（正则化对数变换）：计算速度较慢，但对于小样本量可能表现稍好。
「vst变换的作用」：
解决RNA-seq数据方差依赖均值的问题。vst（方差稳定变换）比简单的log2转换更适合RNA-seq数据。
使数据更接近正态分布，适合下游可视化
blind = FALSE考虑实验设计，能更好揭示组间差异
PCA分析检查数据质量
pca_data <- plotPCA(vsd, intgroup = "dex", returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))
ggplot(pca_data, aes(x = PC1, y = PC2, color = group)) +
geom_point() +
theme_bw(base_size = 8) +
xlab(paste0("PC1: ", percent_var[1], "% variance")) +
ylab(paste0("PC2: ", percent_var[2], "% variance")) +
ggtitle("PCA plot") +
scale_color_manual(values = c("untrt" = "blue", "trt" = "red"))
PCA图可以直观显示样本间的相似性和分组情况，是重要的质控步骤。
5. DESeq2差异表达分析
运行核心分析流程
# 运行DESeq2分析（包含标准化、离散度估计、模型拟合等所有步骤）
dds <- DESeq(dds)
# 提取分析结果
res <- results(dds,
contrast = c("dex", "trt", "untrt"), # 指定比较：处理组 vs 对照组
alpha = 0.05, # 显著性阈值，控制错误发现率(FDR)为5%
lfcThreshold = 0) # log2FC阈值
# 添加显著性标记
res_df <- as.data.frame(res)
res_df$significant <- ifelse(res_df$padj < 0.05 & abs(res_df$log2FoldChange) >= 1,
ifelse(res_df$log2FoldChange >= 1, "Up", "Down"),
"Not significant")
res_df <- na.omit(res_df)
结果统计与保存
# 筛选显著差异表达基因
sig_genes <- subset(res[order(res$padj), ], padj < 0.05 & abs(log2FoldChange) > 1)
cat("显著差异基因数量:", nrow(sig_genes), "\n")
cat("上调基因数量:", sum(sig_genes$log2FoldChange > 0, na.rm = TRUE), "\n")
cat("下调基因数量:", sum(sig_genes$log2FoldChange < 0, na.rm = TRUE), "\n")
# 保存结果
write.csv(res_df, "DESeq2_results.csv", row.names = TRUE)
6. 结果可视化与解读
火山图展示差异表达格局
volcano_plot <- ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = significant), alpha = 0.5) +
scale_color_manual(name = "",
values = c("Not significant" = "gray",
"Up" = "red", "Down" = "blue")) +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey") +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey") +
labs(title = "Volcano Plot of Differential Expression",
x = expression(Log[2]~"Fold Change"),
y = expression(-Log[10]~"P.adj")) +
theme_bw(base_size = 8)
# 添加Top基因标签
volcano_plot +
geom_label_repel(data = res_df[head(rownames(sig_genes), 10), ],
aes(label = rownames(res_df[head(rownames(sig_genes), 10), ])),
box.padding = 0.5,
max.overlaps = Inf)
热图展示基因表达模式
# 选择前50个显著差异基因
select_genes <- rownames(sig_genes)[1:min(50, nrow(sig_genes))]
mat <- assay(vsd)[select_genes, ]
# 绘制热图
pheatmap(mat,
scale = "row", # 对行进行Z-score标准化
show_rownames = TRUE, # 显示基因名
show_colnames = TRUE, # 显示样本名
cluster_rows = TRUE, # 对行聚类
cluster_cols = TRUE, # 对列聚类
annotation_col = as.data.frame(col_data["dex"]),
main = paste("Heatmap of Top", length(select_genes), "Significant DEGs"),
fontsize_row = 8,
color = colorRampPalette(c("blue", "white", "red"))(50)
)
总结
正确的实验设计公式是分析的前提
合理的基因注释和去重策略保证结果可靠性
严格的质量控制和适当的过滤提高分析灵敏度
多角度的可视化帮助理解和验证结果
「完整代码已分享，欢迎在实际研究中应用并根据具体需求调整参数。如果你在分析中遇到问题或有更好的实践经验，欢迎在评论区分享交流！」