Pipeline Detail

RNA-Seq转录组基础与表达变化

Bulk RNA-seq 标准差异表达分析增强版

面向 bulk RNA-seq 项目的标准差异表达增强流程，覆盖实验设计、FASTQ 质控、STAR/Salmon 定量、DESeq2 统计建模、可视化、富集分析与结果交付。

创建时间

2026/6/3

分析难度

入门

推荐场景

转录组分析

预计耗时

0.5-1 天

Metadata

流程元数据

先看应用场景、输入输出和工具依赖，再进入正文命令细节。

Difficulty

入门

Scenario

转录组分析

Estimated Time

0.5-1 天

Tools

DESeq2STARfeatureCountsSalmontximportWGCNArMATSSUPPA2

Inputs

FASTQBAMGTFcount matrixTPM

Outputs

DEG tableheatmapvolcano plotreport

Workflow DAG

流程图

用步骤节点快速理解这个分析从原始数据到结果报告的流转关系。

STEP 1

实验设计与样本表

→

STEP 2

FastQC/MultiQC 原始质控

→

STEP 3

fastp/Trim Galore 清洗

→

STEP 4

STAR+featureCounts 或 Salmon

→

STEP 5

基因 count 矩阵

→

STEP 6

PCA/相关性/离群样本

→

STEP 7

DESeq2 差异表达

→

STEP 8

基因注释与 ID 转换

→

STEP 9

GO/KEGG/GSEA 富集

→

STEP 10

报告与结果交付

Protocol

流程文档

正文保留 Markdown 排版、代码语言标识和表格样式，适合边学边复现。

Bulk RNA-seq 标准差异表达分析增强版

一、适用场景

bulk RNA-seq 以样本为单位测量组织、细胞群或处理条件下的整体转录水平，最常见目标是寻找不同组之间的差异表达基因，并解释这些基因背后的通路和生物学过程。

本增强版适用于：

处理组 vs 对照组的标准差异表达分析。
多处理、多时间点、多组织的表达变化比较。
动植物、微生物、肿瘤、免疫、胁迫响应等常规转录组项目。
需要同时交付 count matrix、TPM、DEG、火山图、热图、PCA、富集分析和可复用报告的项目。

不适合直接套用的场景：

没有生物学重复的项目：只能做探索性分析，不建议做严格差异检验。
单细胞 RNA-seq：应走 scRNA-seq 专用流程。
长读长转录组：应优先考虑 Iso-Seq、FLAIR、TALON、StringTie2 等流程。
可变剪接为核心目标：应增加 rMATS、SUPPA2、MAJIQ 等模块。

二、bulk RNA-seq 分析类型

类型	主要问题	推荐方法	关键输出
基因表达定量	每个样本中每个基因表达量是多少？	STAR + featureCounts 或 Salmon + tximport	raw count、TPM
差异表达分析	处理组相对对照组哪些基因显著变化？	DESeq2、edgeR、limma-voom	DEG table、log2FC、padj
样本质量评估	样本是否聚类合理？是否有离群样本？	PCA、样本相关性、距离热图	PCA plot、correlation heatmap
功能富集分析	DEG 指向哪些功能和通路？	clusterProfiler、fgsea、GSEA	GO/KEGG/GSEA 结果
表达模式聚类	多组/时间点基因如何动态变化？	k-means、Mfuzz、maSigPro	gene clusters
共表达网络	哪些基因模块与性状相关？	WGCNA	module-trait relationship
转录本层分析	isoform 是否变化？	Salmon/Kallisto + tximport、IsoformSwitchAnalyzeR	transcript abundance、isoform switch

三、最佳实践路线

方法选择建议

场景	推荐路线	理由
常规有参考基因组物种	STAR + featureCounts + DESeq2	结果直观，BAM 可复查
样本量很大，想快速完成表达定量	Salmon + tximport + DESeq2	速度快，占用资源少
关注转录本/isoform	Salmon/Kallisto + tximport	保留 transcript-level 信息
非模式物种且注释较弱	HISAT2/STAR + StringTie	可辅助转录本组装
无参考基因组	Trinity + Salmon/RSEM	需要 de novo transcriptome

四、整体流程图

flowchart TD
    A[实验设计与样本信息表] --> B[原始 FASTQ]
    B --> C[FastQC / MultiQC]
    C --> D{质量是否合格?}
    D -- 否 --> E[fastp / Trim Galore 清洗]
    D -- 是 --> F[选择定量路线]
    E --> F
    F --> G1[STAR 比对]
    F --> G2[Salmon 准比对]
    G1 --> H1[featureCounts 基因计数]
    G2 --> H2[tximport 转为 gene-level count]
    H1 --> I[raw count matrix]
    H2 --> I
    I --> J[样本 QC: PCA / 相关性 / 离群样本]
    J --> K[DESeq2 建模与差异检验]
    K --> L[基因注释 / ID 转换]
    L --> M[火山图 / 热图 / MA plot]
    L --> N[GO / KEGG / GSEA]
    M --> O[结果报告]
    N --> O

五、项目目录建议

bulk_rnaseq_project/
├── 00_metadata/
│   ├── sample_info.csv
│   └── contrasts.csv
├── 01_raw_data/
├── 02_qc_raw/
├── 03_clean_data/
├── 04_qc_clean/
├── 05_alignment/
│   ├── star_index/
│   ├── bam/
│   └── logs/
├── 06_counts/
├── 07_deseq2/
│   ├── tables/
│   ├── plots/
│   └── rds/
├── 08_enrichment/
└── report/

样本信息表 sample_info.csv：

sample_id	condition	batch	replicate	fastq_1	fastq_2
Ctrl_1	Ctrl	B1	1	Ctrl_1_R1.fq.gz	Ctrl_1_R2.fq.gz
Ctrl_2	Ctrl	B1	2	Ctrl_2_R1.fq.gz	Ctrl_2_R2.fq.gz
Treat_1	Treat	B1	1	Treat_1_R1.fq.gz	Treat_1_R2.fq.gz
Treat_2	Treat	B1	2	Treat_2_R1.fq.gz	Treat_2_R2.fq.gz

比较组表 contrasts.csv：

contrast_name	numerator	denominator
Treat_vs_Ctrl	Treat	Ctrl

六、实验设计检查

差异表达的统计可靠性主要取决于实验设计。

检查项	推荐
生物学重复	每组至少 3 个，复杂设计建议更多
技术重复	通常先合并或作为 QC，不替代生物学重复
批次信息	必须记录，必要时进入设计公式
测序深度	常规 mRNA-seq 每样本 20-40M reads 可作为起点
配对设计	同一供体前后处理应写入设计公式

DESeq2 设计公式示例：

design = ~ condition
design = ~ batch + condition
design = ~ donor + condition
design = ~ batch + time + condition + time:condition

解释原则：

condition 是你真正关心的生物学变量。
batch、donor、sex 等是需要控制的协变量。
不要把与 condition 完全混杂的 batch 强行加入模型，否则模型不可估。

七、FASTQ 质控与清洗

1. 原始数据质控

mkdir -p 02_qc_raw/fastqc 02_qc_raw/multiqc

fastqc -t 8   -o 02_qc_raw/fastqc   01_raw_data/*.fq.gz

multiqc 02_qc_raw/fastqc   -o 02_qc_raw/multiqc

重点查看：

指标	关注点
Per base sequence quality	末端质量是否明显下降
Adapter Content	是否有接头污染
Per sequence GC content	是否有异常峰
Sequence Duplication Levels	重复率是否过高
Overrepresented sequences	是否有 rRNA、接头或污染序列

2. fastp 清洗

mkdir -p 03_clean_data 04_qc_clean/fastp

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  fastp     -i 01_raw_data/${sample}_R1.fq.gz     -I 01_raw_data/${sample}_R2.fq.gz     -o 03_clean_data/${sample}_R1.clean.fq.gz     -O 03_clean_data/${sample}_R2.clean.fq.gz     --detect_adapter_for_pe     --cut_front     --cut_tail     --cut_window_size 4     --cut_mean_quality 20     --length_required 30     --thread 8     --html 04_qc_clean/fastp/${sample}.html     --json 04_qc_clean/fastp/${sample}.json
done

multiqc 04_qc_clean/fastp -o 04_qc_clean/multiqc

八、路线 A：STAR + featureCounts

1. 构建 STAR 索引

mkdir -p 05_alignment/star_index

STAR   --runThreadN 16   --runMode genomeGenerate   --genomeDir 05_alignment/star_index   --genomeFastaFiles ref/genome.fa   --sjdbGTFfile ref/genes.gtf   --sjdbOverhang 149

sjdbOverhang 通常设为 read length - 1，例如 PE150 设置为 149。

2. STAR 比对

mkdir -p 05_alignment/bam 05_alignment/logs

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  STAR     --runThreadN 16     --genomeDir 05_alignment/star_index     --readFilesIn       03_clean_data/${sample}_R1.clean.fq.gz       03_clean_data/${sample}_R2.clean.fq.gz     --readFilesCommand zcat     --outFileNamePrefix 05_alignment/bam/${sample}_     --outSAMtype BAM SortedByCoordinate     --quantMode GeneCounts     --outFilterMultimapNmax 10     --outFilterMismatchNoverLmax 0.04

  samtools index 05_alignment/bam/${sample}_Aligned.sortedByCoord.out.bam
done

3. 汇总比对质量

mkdir -p 05_alignment/summary

echo "sample,total_reads,unique_mapping_rate,multi_mapping_rate" > 05_alignment/summary/star_mapping_summary.csv

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  log=05_alignment/bam/${sample}_Log.final.out
  total=$(grep "Number of input reads" ${log} | awk '{print $6}')
  unique=$(grep "Uniquely mapped reads %" ${log} | awk '{print $6}' | sed 's/%//')
  multi=$(grep "% of reads mapped to multiple loci" ${log} | awk '{print $9}' | sed 's/%//')
  echo "${sample},${total},${unique},${multi}" >> 05_alignment/summary/star_mapping_summary.csv
done

4. featureCounts 基因计数

mkdir -p 06_counts

featureCounts   -T 12   -p   -B   -C   -t exon   -g gene_id   -a ref/genes.gtf   -o 06_counts/gene_counts.txt   05_alignment/bam/*_Aligned.sortedByCoord.out.bam

提取 count matrix：

cut -f1,7- 06_counts/gene_counts.txt > 06_counts/count_matrix.raw.txt

九、路线 B：Salmon + tximport

如果项目更关注速度和转录本定量，可以使用 Salmon。

mkdir -p ref/salmon_index

salmon index   -t ref/transcripts.fa   -i ref/salmon_index   -p 12

mkdir -p 06_salmon

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  salmon quant     -i ref/salmon_index     -l A     -1 03_clean_data/${sample}_R1.clean.fq.gz     -2 03_clean_data/${sample}_R2.clean.fq.gz     --validateMappings     --gcBias     --seqBias     -p 12     -o 06_salmon/${sample}
done

tximport 导入：

library(tximport)
library(readr)

sample_info <- read.csv("00_metadata/sample_info.csv")
files <- file.path("06_salmon", sample_info$sample_id, "quant.sf")
names(files) <- sample_info$sample_id

tx2gene <- read.delim("ref/tx2gene.tsv", header = TRUE)

txi <- tximport(
  files,
  type = "salmon",
  tx2gene = tx2gene,
  countsFromAbundance = "lengthScaledTPM"
)

十、DESeq2 差异表达分析

1. 读取 count matrix

STAR + featureCounts 路线：

library(DESeq2)
library(tidyverse)
library(pheatmap)
library(ggrepel)

sample_info <- read.csv("00_metadata/sample_info.csv", row.names = 1)
count_data <- read.delim("06_counts/count_matrix.raw.txt", check.names = FALSE)

rownames(count_data) <- count_data$Geneid
count_data <- count_data[, -1, drop = FALSE]
colnames(count_data) <- gsub("_Aligned.sortedByCoord.out.bam", "", basename(colnames(count_data)))
count_data <- count_data[, rownames(sample_info)]

Salmon + tximport 路线：

dds <- DESeqDataSetFromTximport(
  txi,
  colData = sample_info,
  design = ~ condition
)

featureCounts 路线：

dds <- DESeqDataSetFromMatrix(
  countData = round(as.matrix(count_data)),
  colData = sample_info,
  design = ~ condition
)

2. 低表达过滤

keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

过滤原则：

至少在一个组的多数样本中有一定表达。
不要保留大量全零或极低表达基因，否则会降低多重检验效率。
阈值要结合样本数调整，例如每组 3 个重复可用 >=10 且至少 3 个样本。

3. 样本 QC

vsd <- vst(dds, blind = FALSE)

pca_data <- plotPCA(vsd, intgroup = c("condition"), returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))

p_pca <- ggplot(pca_data, aes(PC1, PC2, color = condition, label = name)) +
  geom_point(size = 4) +
  geom_text_repel(max.overlaps = 20) +
  xlab(paste0("PC1: ", percent_var[1], "% variance")) +
  ylab(paste0("PC2: ", percent_var[2], "% variance")) +
  theme_bw()

ggsave("07_deseq2/plots/PCA_samples.pdf", p_pca, width = 7, height = 5)

样本距离热图：

sample_dist <- dist(t(assay(vsd)))
sample_dist_matrix <- as.matrix(sample_dist)

pheatmap(
  sample_dist_matrix,
  annotation_col = sample_info,
  filename = "07_deseq2/plots/sample_distance_heatmap.pdf",
  width = 7,
  height = 6
)

4. 差异检验

dds <- DESeq(dds)

contrast_info <- read.csv("00_metadata/contrasts.csv")
dir.create("07_deseq2/tables", showWarnings = FALSE, recursive = TRUE)
dir.create("07_deseq2/plots", showWarnings = FALSE, recursive = TRUE)

results_list <- list()

for (i in seq_len(nrow(contrast_info))) {
  contrast_name <- contrast_info$contrast_name[i]
  numerator <- contrast_info$numerator[i]
  denominator <- contrast_info$denominator[i]

  res <- results(
    dds,
    contrast = c("condition", numerator, denominator),
    alpha = 0.05
  )

  res_shrink <- lfcShrink(
    dds,
    contrast = c("condition", numerator, denominator),
    res = res,
    type = "ashr"
  )

  res_df <- as.data.frame(res_shrink) |>
    tibble::rownames_to_column("gene_id") |>
    arrange(padj)

  res_df <- res_df |>
    mutate(
      regulation = case_when(
        padj < 0.05 & log2FoldChange >= 1 ~ "Up",
        padj < 0.05 & log2FoldChange <= -1 ~ "Down",
        TRUE ~ "Not significant"
      )
    )

  write.csv(
    res_df,
    file.path("07_deseq2/tables", paste0(contrast_name, "_DESeq2_all.csv")),
    row.names = FALSE
  )

  sig_df <- res_df |>
    filter(padj < 0.05, abs(log2FoldChange) >= 1)

  write.csv(
    sig_df,
    file.path("07_deseq2/tables", paste0(contrast_name, "_DESeq2_significant.csv")),
    row.names = FALSE
  )

  results_list[[contrast_name]] <- res_df
}

5. 火山图

plot_volcano <- function(res_df, contrast_name) {
  plot_df <- res_df |>
    mutate(
      neg_log10_padj = -log10(ifelse(is.na(padj), 1, padj)),
      label = ifelse(padj < 0.05 & abs(log2FoldChange) >= 1, gene_id, NA)
    )

  top_label <- plot_df |>
    filter(!is.na(label)) |>
    arrange(padj) |>
    head(20)

  ggplot(plot_df, aes(log2FoldChange, neg_log10_padj, color = regulation)) +
    geom_point(alpha = 0.7, size = 1.4) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
    geom_text_repel(data = top_label, aes(label = label), size = 3, max.overlaps = 20) +
    scale_color_manual(values = c("Up" = "#d73027", "Down" = "#4575b4", "Not significant" = "#bdbdbd")) +
    theme_bw() +
    labs(title = contrast_name, x = "log2 fold change", y = "-log10 adjusted P")
}

for (contrast_name in names(results_list)) {
  p <- plot_volcano(results_list[[contrast_name]], contrast_name)
  ggsave(file.path("07_deseq2/plots", paste0(contrast_name, "_volcano.pdf")), p, width = 7, height = 6)
}

6. DEG 热图

for (contrast_name in names(results_list)) {
  top_genes <- results_list[[contrast_name]] |>
    filter(padj < 0.05, abs(log2FoldChange) >= 1) |>
    arrange(padj) |>
    head(50) |>
    pull(gene_id)

  if (length(top_genes) >= 2) {
    mat <- assay(vsd)[top_genes, ]
    mat <- mat - rowMeans(mat)

    pheatmap(
      mat,
      annotation_col = sample_info,
      show_rownames = TRUE,
      filename = file.path("07_deseq2/plots", paste0(contrast_name, "_top50_heatmap.pdf")),
      width = 8,
      height = 10
    )
  }
}

十一、TPM 与表达矩阵交付

DESeq2 差异分析使用 raw count，但报告中常需要 TPM 方便展示表达量。

TPM 计算需要基因长度：

counts_to_tpm <- function(counts, gene_length_bp) {
  rate <- counts / (gene_length_bp / 1000)
  t(t(rate) / colSums(rate) * 1e6)
}

gene_length <- read.delim("ref/gene_length.tsv")
gene_length <- gene_length[match(rownames(count_data), gene_length$gene_id), ]

tpm <- counts_to_tpm(as.matrix(count_data), gene_length$length)
write.csv(tpm, "06_counts/gene_tpm.csv")

注意：

TPM 适合样本内和展示层面的表达量解释。
DESeq2/edgeR/limma-voom 的差异检验应使用 count 或模型要求的输入。

十二、基因注释与 ID 转换

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

annotate_genes <- function(res_df) {
  res_df$symbol <- mapIds(
    org.Hs.eg.db,
    keys = res_df$gene_id,
    column = "SYMBOL",
    keytype = "ENSEMBL",
    multiVals = "first"
  )

  res_df$entrez <- mapIds(
    org.Hs.eg.db,
    keys = res_df$gene_id,
    column = "ENTREZID",
    keytype = "ENSEMBL",
    multiVals = "first"
  )

  res_df
}

非人类物种需要替换对应 OrgDb，或使用自定义 GTF 注释表。

十三、GO/KEGG/GSEA 富集分析

1. ORA 富集

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

sig_genes <- results_list$Treat_vs_Ctrl |>
  filter(padj < 0.05, abs(log2FoldChange) >= 1) |>
  pull(gene_id)

sig_entrez <- mapIds(
  org.Hs.eg.db,
  keys = sig_genes,
  column = "ENTREZID",
  keytype = "ENSEMBL",
  multiVals = "first"
) |>
  na.omit()

ego <- enrichGO(
  gene = sig_entrez,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.2
)

write.csv(as.data.frame(ego), "08_enrichment/Treat_vs_Ctrl_GO_BP.csv", row.names = FALSE)

2. GSEA

GSEA 使用全基因排序，不只依赖阈值筛出来的 DEG。

res_df <- results_list$Treat_vs_Ctrl

gene_rank <- res_df$log2FoldChange
names(gene_rank) <- res_df$gene_id
gene_rank <- sort(gene_rank[!is.na(gene_rank)], decreasing = TRUE)

gsea_go <- gseGO(
  geneList = gene_rank,
  OrgDb = org.Hs.eg.db,
  keyType = "ENSEMBL",
  ont = "BP",
  minGSSize = 10,
  maxGSSize = 500,
  pvalueCutoff = 0.05
)

十四、结果解释原则

结果	解读建议
`log2FoldChange`	正值表示 numerator 组更高，负值表示 denominator 组更高
`pvalue`	原始显著性，不建议作为最终筛选标准
`padj`	BH 多重检验校正后的 FDR，常用阈值 `< 0.05`
`baseMean`	平均表达强度，极低表达基因的 fold change 要谨慎解释
`lfcShrink`	收缩低表达和高噪声基因的 log2FC，适合排序和可视化

常用 DEG 阈值：

padj < 0.05
abs(log2FoldChange) >= 1
baseMean >= 10

不要只看差异倍数，也要结合表达量、重复一致性、功能背景和实验设计。

十五、常见问题与排查

问题	可能原因	处理建议
PCA 按批次分开	建库/测序批次效应	设计公式加入 batch，必要时做批次评估
某个样本远离同组	样本污染、标签错误、低质量	回看 FastQC、mapping rate、样本相关性
DEG 数量极少	效应弱、重复少、批次噪声大	检查设计和统计功效，不要强行放宽阈值
DEG 数量异常多	批次混杂、样本组差异过大、污染	检查 PCA 和样本信息
富集结果不稳定	DEG 太少或 ID 转换损失严重	改用 GSEA 或检查注释版本
TPM 与 DESeq2 结果不一致	输入尺度不同	DESeq2 以 count 建模，TPM 用于展示

十六、主要交付物

multiqc_report.html
clean FASTQ 质控报告
STAR BAM、mapping summary 或 Salmon quant.sf
gene raw count matrix
TPM matrix
DESeq2 normalized count
PCA plot
sample distance heatmap
每个比较组的完整差异表
每个比较组的显著 DEG 表
volcano plot
top DEG heatmap
GO/KEGG/GSEA 富集结果
可复用 R 脚本和最终分析报告

十七、参考资料

DESeq2 官方 vignette：差异表达建模、size factor、dispersion、Wald test、LFC shrinkage。
Bioconductor RNA-seq workflow：从 reads 到基因层统计分析的标准实践。
STAR manual：剪接感知 RNA-seq 比对和 GeneCounts 输出。
Subread featureCounts 文档：基因/外显子 read summarization。
MultiQC 文档：多样本质控报告汇总。

Bulk RNA-seq 标准差异表达分析增强版

流程元数据

流程图

流程文档

Bulk RNA-seq 标准差异表达分析增强版

一、适用场景

二、bulk RNA-seq 分析类型

三、最佳实践路线

推荐主路线

方法选择建议

四、整体流程图

五、项目目录建议

六、实验设计检查

七、FASTQ 质控与清洗

1. 原始数据质控

2. fastp 清洗

八、路线 A：STAR + featureCounts

1. 构建 STAR 索引

2. STAR 比对

3. 汇总比对质量

4. featureCounts 基因计数

九、路线 B：Salmon + tximport

十、DESeq2 差异表达分析

1. 读取 count matrix

2. 低表达过滤

3. 样本 QC

4. 差异检验

5. 火山图

6. DEG 热图

十一、TPM 与表达矩阵交付

十二、基因注释与 ID 转换

十三、GO/KEGG/GSEA 富集分析

1. ORA 富集

2. GSEA

十四、结果解释原则

十五、常见问题与排查

十六、主要交付物

十七、参考资料