← 返回分析流程中心

Pipeline Detail

Multi-omics机制解释与多组学调控

多组学联合分析：RNA-seq + ATAC-seq/CUT&Tag

整合 RNA-seq 差异表达、ATAC-seq 开放染色质和 CUT&Tag 组蛋白修饰/转录因子结合信号，解释基因表达变化背后的调控机制。

创建时间
2026/6/3
分析难度
高级
推荐场景
多组学机制解释
预计耗时
3-5 天

Metadata

流程元数据

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

Difficulty

高级

Scenario

多组学机制解释

Estimated Time

3-5 天

Tools

STAR

Outputs

DEG tableheatmapreport

Workflow DAG

流程图

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

STEP 1

建立多组学项目目录

→
STEP 2

样本配对与示例数据

→
STEP 3

RNA-seq DEG

→
STEP 4

ATAC/CUT&Tag peak

→
STEP 5

peak 注释到基因

→
STEP 6

DEG 与 peak 联合分型

→
STEP 7

motif/TF 富集

→
STEP 8

基因座浏览图

→
STEP 9

调控机制报告

Protocol

流程文档

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

多组学联合分析：RNA-seq + ATAC-seq/CUT&Tag

一、项目目标

RNA-seq 告诉我们“哪些基因表达改变”，ATAC-seq/CUT&Tag 告诉我们“哪些调控区域状态改变”。多组学联合分析的目标是把两者串起来，解释基因表达变化背后的调控机制。

典型问题：

问题数据证据
某个 DEG 为什么上调？启动子/增强子 ATAC signal 上升，H3K27ac 上升
哪些 TF 可能驱动变化？上调 peak 中 motif 富集，TF 自身表达上调
哪些 enhancer 可能调控目标基因？peak-to-gene 距离、相关性、同方向变化
CUT&Tag mark 是否支持表达变化？H3K4me3/H3K27ac 与表达正相关，H3K27me3 与表达负相关

二、建立项目目录

mkdir -p multiomics_project/{00_metadata,01_rnaseq,02_atac_cuttag,03_peak_annotation,04_integration,05_motif,06_tracks,report}
mkdir -p multiomics_project/01_rnaseq/{counts,deg,plots}
mkdir -p multiomics_project/02_atac_cuttag/{peaks,bigwig,diffbind,plots}

推荐目录：

multiomics_project/
├── 00_metadata/
│   ├── sample_info.csv
│   └── contrast_info.csv
├── 01_rnaseq/
│   ├── counts/
│   ├── deg/
│   └── plots/
├── 02_atac_cuttag/
│   ├── peaks/
│   ├── bigwig/
│   └── diffbind/
├── 03_peak_annotation/
├── 04_integration/
├── 05_motif/
├── 06_tracks/
└── report/

三、示例数据

00_metadata/sample_info.csv：

sample_id,condition,assay,batch
Ctrl_RNA_1,Ctrl,RNA-seq,B1
Ctrl_RNA_2,Ctrl,RNA-seq,B1
Treat_RNA_1,Treat,RNA-seq,B1
Treat_RNA_2,Treat,RNA-seq,B1
Ctrl_ATAC_1,Ctrl,ATAC-seq,B1
Ctrl_ATAC_2,Ctrl,ATAC-seq,B1
Treat_ATAC_1,Treat,ATAC-seq,B1
Treat_ATAC_2,Treat,ATAC-seq,B1

01_rnaseq/deg/Treat_vs_Ctrl_DEG.csv：

gene_id,symbol,log2FoldChange,padj
ENSG00000141510,TP53,1.8,0.0001
ENSG00000171862,PTEN,-1.2,0.003
ENSG00000136997,MYC,2.1,0.00001

02_atac_cuttag/diffbind/Treat_vs_Ctrl_peaks.bed：

chr8	127735000	127736200	peak_MYC_enhancer	4.5
chr17	7668000	7669500	peak_TP53_promoter	3.2
chr10	87860000	87861500	peak_PTEN_enhancer	-2.7

四、整体流程图

flowchart TD
    A[RNA-seq DEG table] --> C[gene-level regulation]
    B[ATAC/CUT&Tag differential peaks] --> D[regulatory-region regulation]
    D --> E[annotate peaks to nearest genes]
    C --> F[merge DEG and peak annotations]
    E --> F
    F --> G[concordant genes: expression and chromatin same direction]
    F --> H[discordant genes: complex regulation]
    G --> I[motif enrichment and TF prioritization]
    I --> J[locus tracks and mechanism report]

五、peak 注释到基因

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

peak_file <- "02_atac_cuttag/diffbind/Treat_vs_Ctrl_peaks.bed"
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

peak_anno <- annotatePeak(
  peak_file,
  TxDb = txdb,
  tssRegion = c(-3000, 3000),
  annoDb = "org.Hs.eg.db"
)

peak_df <- as.data.frame(peak_anno)
write.csv(peak_df, "03_peak_annotation/Treat_vs_Ctrl_peak_annotation.csv", row.names = FALSE)

六、RNA 和 peak 联合分型

library(tidyverse)

deg <- read.csv("01_rnaseq/deg/Treat_vs_Ctrl_DEG.csv")
peak <- read.csv("03_peak_annotation/Treat_vs_Ctrl_peak_annotation.csv")

peak_gene <- peak |>
  transmute(
    gene_id = geneId,
    peak_id = paste(seqnames, start, end, sep = "_"),
    annotation,
    distanceToTSS,
    peak_log2FC = score
  )

integrated <- deg |>
  left_join(peak_gene, by = "gene_id") |>
  mutate(
    pattern = case_when(
      log2FoldChange > 0 & peak_log2FC > 0 ~ "RNA_up_peak_up",
      log2FoldChange < 0 & peak_log2FC < 0 ~ "RNA_down_peak_down",
      log2FoldChange > 0 & peak_log2FC < 0 ~ "RNA_up_peak_down",
      log2FoldChange < 0 & peak_log2FC > 0 ~ "RNA_down_peak_up",
      TRUE ~ "other"
    )
  )

write.csv(integrated, "04_integration/RNA_ATAC_integrated_table.csv", row.names = FALSE)

图例解释：

模式可能解释
RNA_up_peak_up开放染色质/活性修饰增强，支持基因上调
RNA_down_peak_down调控区域关闭，支持基因下调
RNA_up_peak_down可能受远端 enhancer、TF、RNA 稳定性调控
RNA_down_peak_up可能有抑制性 TF 或复杂 chromatin context

七、motif 富集和 TF 优先级

findMotifsGenome.pl   04_integration/RNA_up_peak_up.bed   hg38   05_motif/RNA_up_peak_up_motif   -size given

整合 TF：

tf_expr <- deg |>
  filter(symbol %in% c("MYC", "JUN", "FOS", "STAT1", "IRF1")) |>
  select(symbol, log2FoldChange, padj)

motif <- read.delim("05_motif/RNA_up_peak_up_motif/knownResults.txt")

tf_priority <- motif |>
  select(Motif.Name, P.value, q.value) |>
  left_join(tf_expr, by = c("Motif.Name" = "symbol")) |>
  arrange(q.value, desc(log2FoldChange))

八、基因座图示例

computeMatrix reference-point   -S 02_atac_cuttag/bigwig/Ctrl.bw 02_atac_cuttag/bigwig/Treat.bw   -R 04_integration/RNA_up_peak_up.bed   --referencePoint center   -b 3000 -a 3000   -o 06_tracks/ATAC_matrix.gz

plotHeatmap   -m 06_tracks/ATAC_matrix.gz   -out 06_tracks/ATAC_peak_heatmap.pdf

九、结果解释例子

MYC 在 RNA-seq 中显著上调，附近 enhancer ATAC signal 和 H3K27ac signal 同时增强。
该 enhancer peak 中富集 AP-1 motif，且 JUN/FOS 表达上调。
因此可以提出假设：处理激活 AP-1 相关调控元件，进一步促进 MYC 表达。

十、交付物

  • DEG 表和 differential peak 表
  • peak 注释表
  • RNA + chromatin 联合分型表
  • motif 富集结果
  • TF 优先级列表
  • deepTools heatmap/profile
  • IGV/locus tracks
  • 候选 enhancer-gene-TF 调控轴报告