跳到正文
百沐一下百沐一下
登录
本节目录

RNA-seq 分析全流程：从 count 矩阵到差异基因与富集

跟完这篇，您能拿一份原始 count 矩阵和一张分组表，从头跑到底：对齐样本、质控看 PCA、用 DESeq2 算差异、卡阈值取出上调 / 下调基因、导出结果表，再做富集把机制讲出来——每一步都有能照着敲的代码，深水区链到各专篇。

一条常规 bulk RNA-seq 主线，就是下面这五步。每步的输入都是上一步的输出，中途出问题往往能在对应节点定位。

步骤做什么常用工具细节专篇
1. 读入数据整理 count 矩阵 + 分组表，对齐样本名read.csv矩阵准备
2. 质控看 PCA判断组间分不分得开、有没有离群样本DESeq2 · vstPCA 怎么看
3. 跑差异算每个基因的 log2FC 与 padjDESeq2DESeq2 实操
4. 定阈值取基因卡出「上调 / 下调」差异基因并导出——本页
5. 富集讲机制把基因列表变成通路和功能clusterProfiler富集分析

示例以人类数据、DESeq2 较新版本为准；代码用于说明每步「做什么、调什么」，不追求整段复制即能跑。

第 1 步：装包、读入数据、对齐样本名

做什么：把 count 矩阵和分组表读进来，确认两者的样本名严丝合缝地对上。为什么：下游所有统计都建立在「样本 → 分组」这层对应上，样本名一旦错位，后面全盘皆错——这也是最高频的报错来源。

先备齐两样东西，缺一不可：

  1. count 矩阵：行是基因、列是样本，每个格子是原始 read count（非负整数）。
  2. 分组表 coldata：行是样本、列是分组等变量（如 condition、batch），行名必须和矩阵的列名一致。
library(DESeq2)

# 读入原始 count 矩阵，第一列作行名（基因），check.names 保留原始样本名
counts <- read.csv("counts.csv", row.names = 1, check.names = FALSE)
counts <- as.matrix(counts)

# 读入分组表，第一列作行名（样本）
coldata <- read.csv("coldata.csv", row.names = 1)
coldata$condition <- factor(coldata$condition)   # 分组变量一定要转成因子

# 让 coldata 的行严格按 counts 的列排序，再确认完全对齐
coldata <- coldata[colnames(counts), , drop = FALSE]
all(colnames(counts) == rownames(coldata))        # 必须为 TRUE

关键提醒：

  • 用原始 count，别用 TPM / FPKM。DESeq2 的负二项模型需要整数原始计数在内部自己估计文库大小因子；喂进归一化或取过对数的值，结果不再可信。三者差在哪、该喂哪种，见 count / TPM / FPKM。
  • 样本名的大小写、下划线、前后缀都算数，sample_01 和 Sample01 是两个名字。
  • 用 Salmon / kallisto 得到的是估计 count（可能带小数），不要手动四舍五入，改用 tximport 导入。矩阵从上游结果怎么整理成这个形状，见 矩阵准备。

第 2 步：建对象、质控看 PCA

做什么：把两份数据包成一个 DESeq2 对象，先做方差稳定化变换、画一张 PCA，别急着算差异。为什么：PCA 一眼回答两个决定成败的问题——组间分得开吗？有没有明显离群的坏样本？这两点决定后面的差异结果值不值得信。

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = coldata,
  design    = ~ condition        # 想比较的分组变量
)

# 明确指定参照组，否则按字母顺序默认，会决定后面倍数变化的方向
dds$condition <- relevel(dds$condition, ref = "normal")

# 方差稳定化后再做 PCA，避免高表达基因主导
vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = "condition")

怎么读这张图：

  • 同组样本聚在一起、不同组分开 —— 理想情况，可以放心往下走。
  • 某个样本远远飞出去 —— 疑似离群，先查是不是建库 / 测序批次问题，必要时记录原因再决定去留，别随手删。
  • 分组完全混在一起 —— 主要变异不来自您关心的分组，先回头查分组标注和潜在批次效应，见 批次效应怎么处理。

PCA、样本距离热图这些质控图怎么系统地读，见 PCA 怎么看。

第 3 步：跑差异分析（DESeq2）

做什么：一行 DESeq() 完成归一化、离散度估计和假设检验，再用 contrast 取出「谁比谁」的结果表。为什么：DESeq2 是被广泛接受的 bulk RNA-seq 差异分析标准方法之一，对小样本、低计数基因稳健，输出的 padj 已做多重检验校正。

# 可选：预过滤几乎不表达的基因，加速并去噪，阈值按数据调
dds <- dds[rowSums(counts(dds)) >= 10, ]

dds <- DESeq(dds)

# contrast 明确方向：tumor 相对 normal
res <- results(dds, contrast = c("condition", "tumor", "normal"))
summary(res)

# 收缩 log2FC，让低表达基因的倍数不被高估，画火山图、排序更稳
res <- lfcShrink(dds, coef = "condition_tumor_vs_normal", type = "apeglm")

提示：coef 的名字要和 resultsNames(dds) 里的一致；上面这个名字，正是因为第 2 步把 normal 设成了参照组才成立。收缩只影响 log2FoldChange 的展示，不改变 padj 与显著性判断。

结果表里每列的含义：

列含义怎么用
baseMean所有样本归一化后表达量的平均太低的基因倍数变化通常不稳
log2FoldChangelog2 尺度的倍数变化，正上调、负下调=1 即 2 倍，=−1 即降到一半；方向由 contrast 决定
lfcSElog2FoldChange 的标准误越小说明倍数估计越稳
pvalue原始 p 值别直接当阈值，没做多重检验校正
padjBH 法校正后的 p 值（即 FDR）筛差异基因就看这一列

padj 出现 NA 是正常的：DESeq2 会对表达过低、几乎不可能显著的基因做独立过滤，直接记为 NA，筛选时当作「未通过」即可。参数细节、多因素设计（如校正批次 ~ batch + condition）见 DESeq2 实操。

第 4 步：定阈值、取差异基因、导出

做什么：在结果表上卡两条线——统计显著性（padj）和效应量（log2FC）——取出「值得讲」的差异基因，分成上调 / 下调，并导出结果表。为什么：只卡 padj 会混进一堆变化极小的基因；同时卡效应量，留下的才既显著、变化又够大。

# 常用惯例：padj < 0.05 且 |log2FC| > 1（约 2 倍变化）
sig  <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

up   <- subset(sig, log2FoldChange > 0)   # 上调
down <- subset(sig, log2FoldChange < 0)   # 下调

nrow(up); nrow(down)                        # 看看各取到多少个

# 导出完整结果表，按 padj 排序，方便后续查阅与复现
res_df <- as.data.frame(res)
res_df <- res_df[order(res_df$padj), ]
write.csv(res_df, "deseq2_results.csv")

关于阈值，务必心里有数：

  • padj < 0.05 和 |log2FC| > 1 是领域内常用惯例，不是铁律。样本量大、噪声小时可收紧到 padj < 0.01；差异基因太少时可放宽到 |log2FC| > 0.585（约 1.5 倍）。
  • 阈值应在看数据之前大致定好，而不是反复调到「结果好看」为止，后者会人为夸大显著性。
  • 把最终阈值写进方法学部分，保证可复现。

取到的差异基因常配一张火山图来展示，画法见 火山图。

第 5 步：富集，把基因列表讲成机制

做什么：把差异基因列表丢给富集分析，看它们集中在哪些通路和功能上。为什么：一串基因名讲不出生物学；富集把它们翻译成「这批基因主要参与某某过程」，才是机制层面的结论。

多数富集函数认 ENTREZ ID，而矩阵行名常是 SYMBOL 或 ENSEMBL，先转换：

library(clusterProfiler)
library(org.Hs.eg.db)

# 差异基因行名若是 SYMBOL，转成 ENTREZ ID；ENSEMBL 就把 fromType 改成 "ENSEMBL"
id_map <- bitr(rownames(sig), fromType = "SYMBOL",
               toType = "ENTREZID", OrgDb = org.Hs.eg.db)

两条常规路线，任选或并用：

  1. ORA（过表达分析） —— 输入上一步筛出的差异基因，问「它们在哪些通路里显著扎堆」。

    ego <- enrichGO(
  gene          = id_map$ENTREZID,
  OrgDb         = org.Hs.eg.db,
  ont           = "BP",            # 生物学过程，另有 MF、CC
  pAdjustMethod = "BH",
  qvalueCutoff  = 0.05,
  readable      = TRUE             # 把结果里的 ID 转回基因名，方便读
)

dotplot(ego, showCategory = 15)    # 点图展示最显著的若干条通路

  2. GSEA（基因集富集） —— 不预先卡阈值，输入全部基因按 log2FC 排序的列表，能捕捉「整体轻微但一致偏移」的通路，是 ORA 的补充。做法见 GSEA 怎么做。

选 GO 还是 KEGG、ID 怎么转换、图怎么调，见 富集分析 与 基因 ID 转换。

把五步串成一个能写进论文的故事

做什么：把前面的结果拼成一条经得起审稿的逻辑线。为什么：审稿人看的不是单张图，而是「数据可信 → 差异明确 → 机制合理」这条链是否环环相扣。

一条常见的叙述骨架：

  1. 数据可信：PCA 显示分组分得开、无异常样本 —— 交代质控，为后续结论托底；
  2. 差异明确：在既定阈值下得到 N 个差异基因（上调 X、下调 Y）—— 给出规模，火山图 / 热图佐证；
  3. 机制成形：富集显示这些基因集中在某某通路 —— 从基因列表上升到生物学过程；
  4. 回扣问题：把富集到的通路与您的科学问题对上，给出一句话主结论。

每一步都留好可复现的代码、参数和导出的结果表，这条链才立得住。

在平台上跑一遍

上面这条主线，在百沐一下不写代码也能走完：上传 count 矩阵和分组表，用一句话说清目的（如「对 tumor vs normal 做差异表达，再做 GO / KEGG 富集」），组学智能体 会自动规划 DESeq2 差异、阈值筛选与富集，并在关键处（分组、参照、阈值）跟您确认，最后把结论、图形、结果表和可复现的 R 代码一并交付。

打开 https://app.baimuyixia.com 的组学智能体，选好数据、描述目标即可开跑。

相关