Pipeline Detail

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

时间序列 RNA-seq 分析

面向多时间点 bulk RNA-seq 的动态表达分析流程，覆盖实验设计、DESeq2 LRT、maSigPro/ImpulseDE2、表达模式聚类、时间趋势可视化和功能解释。

创建时间

2026/6/3

分析难度

中级

推荐场景

转录组分析

预计耗时

1-3 天

Metadata

流程元数据

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

Difficulty

中级

Scenario

转录组分析

Estimated Time

1-3 天

Tools

DESeq2STARmaSigProImpulseDE2

Inputs

FASTQcount matrix

Outputs

heatmap

Workflow DAG

流程图

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

STEP 1

时间序列实验设计

→

STEP 2

FASTQ/样本 QC

→

STEP 3

表达定量与 count matrix

→

STEP 4

VST/标准化

→

STEP 5

DESeq2 LRT 时间效应

→

STEP 6

maSigPro/ImpulseDE2 趋势建模

→

STEP 7

动态表达模式聚类

→

STEP 8

时间模块富集

→

STEP 9

趋势图与报告

Protocol

流程文档

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

时间序列 RNA-seq 分析

一、适用场景

时间序列 RNA-seq 用于研究基因表达随时间变化的动态过程，例如发育阶段、药物处理后响应、病原感染过程、胁迫恢复过程或细胞分化诱导过程。它关注的不只是某个时间点的差异，而是完整的表达轨迹。

推荐场景：

场景	目标	推荐方法
单一处理多时间点	找随时间变化的基因	DESeq2 LRT、maSigPro
处理组 vs 对照组多时间点	找时间和处理交互基因	`~ condition + time + condition:time`
非线性表达响应	捕捉脉冲、峰值和恢复	ImpulseDE2、spline model
表达模式归类	找早期/中期/晚期响应模块	k-means、Mfuzz、TCseq
功能解释	每类动态基因代表什么过程	GO/KEGG/GSEA

二、整体流程图

flowchart TD
    A[样本信息: condition/time/replicate] --> B[FASTQ QC 和表达定量]
    B --> C[raw count matrix]
    C --> D[DESeq2 VST / normalized count]
    D --> E[样本 PCA 与时间轨迹检查]
    E --> F[DESeq2 LRT 检测时间效应]
    E --> G[maSigPro / ImpulseDE2 建模趋势]
    F --> H[显著动态基因]
    G --> H
    H --> I[表达模式聚类]
    I --> J[模块趋势图]
    I --> K[GO/KEGG/GSEA 富集]
    J --> L[动态表达报告]
    K --> L

三、实验设计要点

设计项	建议
时间点	至少 3 个，越多越利于趋势建模
生物学重复	每个时间点每组至少 3 个
时间编码	可以作为 factor，也可以作为 numeric
批次	多批建库必须记录 batch
配对设计	同一对象重复采样时记录 subject

样本表：

sample_id,condition,time,batch,replicate
Ctrl_T0_1,Ctrl,0,B1,1
Ctrl_T6_1,Ctrl,6,B1,1
Treat_T0_1,Treat,0,B1,1
Treat_T6_1,Treat,6,B1,1

四、DESeq2 LRT 检测时间效应

LRT 适合回答“加入时间或交互项后模型是否显著更好”。

library(DESeq2)
library(tidyverse)

counts <- read.csv("counts.csv", row.names = 1, check.names = FALSE)
coldata <- read.csv("sample_info.csv", row.names = 1)
counts <- counts[, rownames(coldata)]

coldata$time <- factor(coldata$time, levels = c("0", "6", "12", "24", "48"))
coldata$condition <- factor(coldata$condition)

dds <- DESeqDataSetFromMatrix(
  countData = round(as.matrix(counts)),
  colData = coldata,
  design = ~ condition + time + condition:time
)

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

dds_lrt <- DESeq(
  dds,
  test = "LRT",
  reduced = ~ condition + time
)

res_interaction <- results(dds_lrt)
dynamic_genes <- as.data.frame(res_interaction) |>
  rownames_to_column("gene_id") |>
  filter(padj < 0.05) |>
  arrange(padj)

write.csv(dynamic_genes, "time_condition_interaction_genes.csv", row.names = FALSE)

解释：

full model: condition + time + condition:time
reduced model: condition + time
显著基因说明处理组和对照组的时间变化轨迹不同。

五、趋势可视化

vsd <- vst(dds, blind = FALSE)
expr <- assay(vsd)

plot_gene_trend <- function(gene_id) {
  df <- data.frame(
    expression = expr[gene_id, ],
    coldata
  )

  ggplot(df, aes(as.numeric(as.character(time)), expression, color = condition)) +
    stat_summary(fun = mean, geom = "line", linewidth = 1) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 1) +
    geom_point(size = 2, alpha = 0.8) +
    theme_bw() +
    labs(x = "Time", y = "VST expression", title = gene_id)
}

plot_gene_trend(dynamic_genes$gene_id[1])
ggsave("top_dynamic_gene_trend.pdf", width = 6, height = 4)

六、maSigPro 动态建模

maSigPro 适合多时间点、多组别趋势建模。

library(maSigPro)

expr_matrix <- expr[dynamic_genes$gene_id, ]
design <- make.design.matrix(
  data = data.frame(
    Time = as.numeric(as.character(coldata$time)),
    Replicate = coldata$replicate,
    Group = coldata$condition
  ),
  degree = 3
)

fit <- p.vector(expr_matrix, design, Q = 0.05)
tstep <- T.fit(fit, step.method = "backward")
sig_genes <- get.siggenes(tstep, rsq = 0.6, vars = "all")

七、表达模式聚类

library(pheatmap)

selected <- head(dynamic_genes$gene_id, 1000)
mat <- expr[selected, ]
mat_z <- t(scale(t(mat)))

set.seed(123)
km <- kmeans(mat_z, centers = 6, nstart = 30)

cluster_df <- data.frame(
  gene_id = rownames(mat_z),
  cluster = km$cluster
)

write.csv(cluster_df, "time_series_gene_clusters.csv", row.names = FALSE)

pheatmap(
  mat_z,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  annotation_col = coldata[, c("condition", "time"), drop = FALSE],
  show_rownames = FALSE,
  filename = "time_series_dynamic_heatmap.pdf"
)

图例解释：

模式	含义
Cluster 1 早期上升	可能是快速应答基因
Cluster 2 晚期上升	可能是继发响应或发育后期通路
Cluster 3 先升后降	可能是 transient response
Cluster 4 持续下降	可能是被处理抑制的过程

八、功能富集

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

for (k in sort(unique(cluster_df$cluster))) {
  genes <- cluster_df |>
    filter(cluster == k) |>
    pull(gene_id)

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

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

  write.csv(as.data.frame(ego), paste0("cluster_", k, "_GO_BP.csv"), row.names = FALSE)
}

九、结果解释示例

Cluster 1 在 0-6h 快速上升，随后恢复，富集到 inflammatory response。
这说明处理早期诱导了急性炎症通路。

Cluster 4 在 12-48h 持续下降，富集到 cell cycle。
这可能提示处理后细胞增殖逐步受抑。

十、主要交付物

样本 PCA 与时间轨迹图
LRT 显著动态基因表
每个基因的趋势图
动态基因 heatmap
表达模式聚类表
每个 cluster 的 GO/KEGG/GSEA 富集
时间响应机制解释报告