Bulk RNA Seq Analysis

Author

Anqi Dong

Published

February 11, 2026

本文档使用 DESeq2 包对组织的转录组测序进行下游分析，从 raw counts 数据到差异表达基因（Differentially Expressed Genes, DEGs）。后续可以进一步绘制火山图、做通路富集分析、基因集变异分析等等。

数据读取和清洗

library(dplyr)
library(stringr)
library(ggplot2)
library(DESeq2)

读取counts矩阵，注意此处不能是reads，也不能是 normalize 或者 log 后的基因表达矩阵。 清洗数据，使每行为一个基因，每列为一个样本，并设置行名为基因名，列名为样本名。

f <- file.path('_')
countData <- readxl::read_xlsx(f)
countData

设置分组信息。是一个只有一列的数据框，行名为样本名。此处演示前8个样本为实验组，后10个样本为对照组。

colData <- data.frame(
  row.names = colnames(countData),
  condition = c(rep('treatment', 8), rep('control', 10))
)

DESeq2流程

# 构建 DESeq2 数据对象
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ condition)

# 过滤低表达基因，此步骤可有可无，对结果几乎无影响。
dds <- dds[rowSums(counts(dds)) > 10, ]

# 差异表达分析
dds <- DESeq(dds)

# 获取结果。
# 注意此处contrast参数的顺序，表示分组信息储存在condition中，计算treatment相对于 control的变化，也就是 > 0 表示在实验组上调。
res <- results(dds, contrast = c("condition", "treatment", "control"))

head(res)

明确上调和下调的判断标准，一般log2FoldChange设置为1.5或1.2或1。

# calculate statistical difference ----
res_sig <- res |>
  tibble::rownames_to_column(var = 'gene') |> 
  dplyr::mutate(signif = case_when(
    log2FoldChange > 1 & pvalue < 0.05 ~ 'up',
    log2FoldChange < -1 & pvalue < 0.05 ~ 'down',
    .default = 'ns'
  ))
count(res_sig, signif)

# 保存数据
readr::write_csv(res_sig, file = '_.csv')

GSVA

基因集变异分析，指计算某基因集合在样品中的表达量，从而来评估不同的生物学活动在不同样品间是否富集。与 GSEA 的区别是，GSVA 只需要一个样品即可计算，而 GSEA 至少需要两个（组）样本计算差异倍数再进行分析。另有单样本 GSEA（ssGSEA）可实现类似功能，但计算原理不同。

GSVA 运算前需先制作一个专门的 GSVA 参数对象，其中 exprData 是基因表达 matrix 变量，一般为 logCPM、logTPM 或者 logFPKM，geneSets 是包含数个自定义基因集的 list，kcdf 根据提供的基因表达矩阵分情况设置，如果是上述 log 值时可以设置为 'Gaussian'。理论上 exprData 也可以是 raw counts 值，此时 kcdf 设置为 'Poisson'。但是既然 raw counts 有了那么计算 log 值也不是难事，所以还是推荐使用 log 值和 'Gaussian'。关于 TPM、CPM、FPKM 的优劣已是较为经典且古老的讨论话题，本文不展开讨论。

此处以 DESeq2 计算得到的 dds 为例，先用 DESeq2::counts(dds, normalized = T) 计算 normalized counts，再进行 GSVA 运算。

normalized_counts <- counts(dds, normalized = T)

# GSVA ----
log_tpm <- as.matrix(log2(normalized_counts + 1))
gene_set <- list(GENE_123 = c('GENE_1', 'GENE_2', 'GENE_3'),
                 GENE_456 = c('GENE_4', 'GENE_5', 'GENE_6'))

param <- gsvaParam(
  exprData = log_tpm,
  geneSets = gene_set,
  kcdf = "Gaussian"
)
gsva_mat <- gsva(param)

## make a data.frame for statistic analysis and plot
res <- data.frame(score = gsva_mat[GENE_123,], group = colData$condition)
res
wilcox.test(res$score ~ res$group)

# plot ----
p <- res |> ggplot(aes(group, score)) + geom_boxplot() + theme_bw()
p

基因芯片 Microarray

GSE67501 的数据文件标题是 non-normalized data，但其数据不是整数。GEO 网站注明 Experiment type 为 Expression profiling by array，平台为 Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip 也表明这是基因芯片数据，不是测序数据。使用 limma 包计算：先取 log 再 normalize，如果要计算 logTPM 则需再 log 一次。和测序得到的 raw counts 先 normalize 再 log 不一样。

# limma ----
library(limma)

## Step 1: log2
if (max(data) > 100) {
  expr <- log2(data + 1)
}

## Step 2: between-array normalization
expr.norm <- normalizeBetweenArrays(expr, method = "quantile")

## Step 3: limma DE
design <- model.matrix(~ colData$condition)
fit <- lmFit(expr.norm, design)
fit <- eBayes(fit)
res <- topTable(fit, number = Inf)

# GSVA ----
log_norm <- as.matrix(log2(expr.norm + 1))
gene_set <- list(GENE_123 = c('GENE_1', 'GENE_2', 'GENE_3'),
                 GENE_456 = c('GENE_4', 'GENE_5', 'GENE_6'))

param <- gsvaParam(
  exprData = log_norm,
  geneSets = gene_set,
  kcdf = "Gaussian"
)

gsva_mat <- gsva(param)

res <- data.frame(score = gsva_mat[GENE_123,], group = colData$condition)
res
with(res, wilcox.test(score ~ group))

# plot ----
p <- res |> ggplot(aes(group, score)) + geom_boxplot() + theme_bw()
p