RNA-seq（3）：用 DESeq2 做差异表达分析——以 airway 数据为例

以解牛之法析生信，观微雀之形览科研。

今天我们一起来学习使用DESeq2进行转录组的差异分析。

三、DEseq2差异分析及对比

使用数据：airway
工具：DESeq2
相关R包：DESeq2（分析）、ggplot2（绘图）

1. 工具介绍：DESeq2

DESeq2是 RNA-seq 差异表达分析 最主流的 R 包之一（发表于 2014 年，已被引用数万次），专门用于基于负二项分布模型的计数数据差异分析。它由 Bioconductor 项目维护，是转录组分析的标准工具。

输入数据：原始整数计数矩阵（不能是 TPM/FPKM），行=基因，列=样本

标准化方法： median-of-ratios 方法，通过计算每个基因的相对计数因子来校正文库大小差异，对异常值更稳健

差异检验： 基于 Wald 检验 或 似然比检验（LRT），自动估计基因离散度

处理小样本： 使用 经验贝叶斯方法 将基因间离散度信息“收缩”（shrinkage），使小样本（如 n=3）也能获得可靠结果

log2FC 收缩： 提供 apeglm 和 ashr 两种收缩方法，对低计数基因的 log2FC 进行适度压缩，提高结果稳定性

DESeq2说明文档

Bioconductor 手册：https://bioconductor.org/packages/devel/bioc/manuals/DESeq2/man/DESeq2.pdf

github仓库：https://github.com/thelovelab/DESeq2

使用说明：

# 如何使用DESeq(object,test=c("Wald","LRT"),fitType=c("parametric","local","mean","glmGamPoi"),sfType=c("ratio","poscounts","iterate"),betaPrior,full=design(object),reduced,quiet=FALSE,minReplicatesForReplace=7,modelMatrixType,useT=FALSE,minmu=if(fitType=="glmGamPoi")1e-06else0.5,parallel=FALSE,BPPARAM=bpparam())

2. 下载及加载包

# # 第一次安装可以选中ctrl+shift+c取消注释，运行安装# 安装biocManager# if (!require("BiocManager", quietly = TRUE))# install.packages("BiocManager")# # 差异分析核心DESeq2包# BiocManager::install(c("DESeq2", "apeglm"), update = FALSE)# # 数据处理# BiocManager::install(c("tidyverse", "tibble"), update = FALSE)# # 绘图包# BiocManager::install(c("ggplot2", "EnhancedVolcano", "ggrepel", "pheatmap", "ggpubr", "RColorBrewer"), update = FALSE)# 1.2.3 加载包，确定已经安装了包后再加载# 如果报错或者提示没有安装，返回上一步library(DESeq2)library(apeglm)library(tidyverse)library(ggplot2)library(EnhancedVolcano)library(ggrepel)library(pheatmap)library(ggpubr)

3. 示例数据加载：

airway数据：

if(!require("BiocManager",quietly=TRUE))install.packages("BiocManager")options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")# 切换中科大镜像下载BiocManager::install("airway")#下载airway数据，2014年发表library(airway)#加载包data(airway)#加载airway数据集airway#查看数据集是否加载成功

Bioconductor 页面：https://bioconductor.org/packages/release/data/experiment/html/airway.html

原始 GEO 项目：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778

airway 数据集来源于 GEO 项目 GSE52778，是该项目的子集（包含 4 个细胞系各 1 对处理/对照样本，共 8 个样本），数据处理方式与原始研究一致。

注：如果加载不成功，也可以手动下载，但是此步需要进一步处理才可使用。
地址：https://bioconductor.org/packages/release/data/experiment/html/airway.html
或者使用下载的GSE52778，但是也需要适当的处理才可使用。

4. 标准的 DESeq2 差异分析流程

可从airway中：获取原始整数技术矩阵并构建分组信息——过滤——构建 DESeq2 对象——运行 DESeq2，进行差异基因的比对。

# 标准的分析流程# 2.1 提取数据library(airway)data(airway)count<-assay(airway)# 原始整数计数矩阵col_data<-colData(airway)# 样本注释（含 cell 和 dex）# 构建分组信息（也可以直接从 col_data 提取）group<-data.frame(sample=colnames(count),group=col_data$dex,# "trt" / "untrt"cell=col_data$cell# 细胞系信息，用于校正批次)# 确认样本顺序一致（必须为 TRUE）stopifnot(all(colnames(count)==group$sample))# 2.2 过滤低表达基因（RNA-seq 标准操作，提高统计功效）# 保留至少在 3 个样本中 count >= 10 的基因（可根据实际情况调整阈值）keep<-rowSums(count>=10)>=3count_filter<-count[keep,]cat("过滤前基因数：",nrow(count),"\n")cat("过滤后基因数：",nrow(count_filter),"\n")# 2.3 构建 DESeq2 对象（配对设计，校正细胞系效应）dds<-DESeqDataSetFromMatrix(countData=count_filter,# 使用过滤后的矩阵colData=group,design=~cell+group# 加入 cell 作为协变量)# 明确指定对照组（确保 log2FC 方向正确）dds$group<-relevel(dds$group,ref="untrt")# 2.4 运行 DESeq2dds<-DESeq(dds)# 查看可用的系数名（用于后续对比）resultsNames(dds)

可以看到过滤前后基因数为，过滤前基因数： 63677 ，过滤后基因数： 16596 。此步完成了差异分析，核心包为DESeq2，下步进行差异分析的提取。

5. 提取差异基因

提取差异基因，并且分别保存全部的分析结果和显著的差异基因。

# 2.5 提取两组比较的结果# contrast 参数：c(分组列名, 实验组, 对照组)res<-results(dds,contrast=c("group","trt","untrt"))# 按 p 值排序（最显著的在前）res<-res[order(res$pvalue),]# 转为数据框，并添加基因名列res_df<-as.data.frame(res)res_df$gene<-rownames(res_df)# 保存全部基因的差异结果write.csv(res_df,file="DESeq2_all_gene_result.csv",row.names=FALSE)# 2.6 筛选显著差异基因（阈值：|log2FC| > 1 且 padj < 0.05）res_df$change<-"Stable"res_df$change[res_df$log2FoldChange>1&res_df$padj<0.05]<-"Up"res_df$change[res_df$log2FoldChange<-1&res_df$padj<0.05]<-"Down"# 统计上下调基因数量table(res_df$change)cat("显著上调基因数：",sum(res_df$change=="Up"),"\n")cat("显著下调基因数：",sum(res_df$change=="Down"),"\n")# 2.7 （可选）保存显著差异基因子集sig_genes<-res_df[res_df$change!="Stable",]write.csv(sig_genes,file="DESeq2_sig_genes.csv",row.names=FALSE)# 打印前几个显著基因head(res_df[res_df$change!="Stable",])

可以看到 有显著上调基因511，下调基因489。

运行到这里我们完成了：提取差异基因，并且保存全部的分析结果和显著的差异基因结果。

后续绘制火山图、富集分析，就可直接使用这数据啦！

今天的生信知识分享就到这里啦！咱们下期见！

以解牛之法析生信，观微雀之形览科研。