BioF3 组学数据分析

表观组学实践教程

导出日期：2026年5月13日

表观组学实践教程

表观组学关注的不是基因本身的序列，而是基因表达能不能被"打开"：染色质开放程度、转录因子结合位置、组蛋白修饰、DNA 甲基化。同一个基因组在不同细胞里呈现不同的表观图谱，这些图谱决定了细胞状态和对环境的反应。

本专栏会围绕最常用的三类实验展开：ATAC-seq（染色质开放区域）、ChIP-seq（蛋白-DNA 结合）、WGBS/RRBS（DNA 甲基化）。思路和 bulk RNA-seq 类似：先搞清楚每步产物是什么，再决定用什么工具把它们跑出来。

一个项目大致是什么样子

ATAC-seq 和 ChIP-seq 的分析主线非常像。从一批样本的 FASTQ 出发，到"差异开放区域"或"差异结合位点"表，大致要经过：

步骤 典型产物 常用工具
接头剪切与质控 清洗后的 FASTQ fastp、FastQC、MultiQC
比对到参考基因组 sorted.bam Bowtie2、BWA-MEM
过滤线粒体 / 重复 清洁 BAM samtools、Picard
peak calling narrowPeak / broadPeak MACS2、MACS3
peak 注释 基因附近 peak 映射表 ChIPseeker、HOMER
motif 分析 显著富集的 motif HOMER、MEME Suite
差异分析 差异 peak 列表 DiffBind、DESeq2 on peak counts
与表达整合 peak ↔ 基因关联 自定义脚本 + ggplot2

WGBS/RRBS 的流程不同：比对要用 bisulfite-aware 工具（Bismark、BWA-Meth），输出是每个 CpG 的甲基化率，差异分析常用 methylKit 或 DSS。

常见工具栈

下面是 BioF3 例子里会优先使用的组合：

阶段 工具 说明
ATAC/ChIP 比对 Bowtie2、BWA-MEM 两者都行，Bowtie2 对 ATAC 友好
BAM 处理 samtools、Picard 去重、过滤 MAPQ、去线粒体
peak calling MACS2 ATAC 和 ChIP 都支持，参数略不同
peak 注释 ChIPseeker R 包，输出表格和图
motif HOMER、MEME HOMER 一条命令出 motif 报告
差异 peak DiffBind、DESeq2 样本数少时 DiffBind 更便捷
可视化 deepTools、IGV、pygenometracks 覆盖度曲线、heatmap、基因浏览器图
甲基化 Bismark、methylKit 标准 WGBS/RRBS 流水线

推荐公开数据集

表观组教学比较依赖公开数据，下面几份是常用参照：

数据集 类型 适合 入口
ENCODE K562 ATAC-seq ATAC-seq，细胞系 ATAC 流程练习 ENCODE
ENCODE H3K27ac ChIP-seq ChIP-seq + input ChIP 流程、peak 注释 ENCODE
10x Genomics PBMC scATAC 10k scATAC 单细胞方向的过渡 10x Genomics
TCGA / GDC 甲基化 450K / EPIC 芯片 差异甲基化入门 GDC

ENCODE 的好处是样本类型全、质量稳定、每个实验都有对应的 input 或 control，新手跟着跑最不容易出问题。

最小可跑的例子

下面用 R 的 ChIPseeker 包做一次最简单的 peak 注释，它自带一个 Nature Neuroscience 论文发布的真实 narrowPeak 文件。数据很小，跑完只需要几秒：

# 一次性安装依赖（如果还没装）
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("ChIPseeker", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

# ChIPseeker 自带的示例 peak 文件
peak_files <- getSampleFiles()
peak_files

# 载入其中一个样本的 peak
peak <- readPeakFile(peak_files[[1]])
peak

# 注释到最近的基因
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
anno <- annotatePeak(peak, TxDb = txdb, tssRegion = c(-3000, 3000))

# 看每个 peak 落在什么类型的区域
head(as.data.frame(anno))

# 基因组区域饼图
plotAnnoPie(anno)

plotAnnoPie 会画出 peak 主要落在哪些基因组区域（启动子、内含子、基因间等），这也是 ChIP/ATAC 文章里最常见的一张配图。把这个例子读懂，再去跑真实数据的 peak calling、motif 分析就会顺很多。

专栏模块规划

模块 主题 状态
01 实验类型与数据格式 已上线
02 Peak 注释与多样本比较 已上线
03 DiffBind 差异结合分析 已上线
04 Peak 可视化与多样本比较 已上线
05 Motif 富集与 HOMER 已上线
06 ATAC-seq 分析要点 已上线
07 DNA 甲基化分析入门 已上线
08 表观组与转录组的整合 已上线

目前 01-08 全部上线。01-04 带可跑脚本，05-08 为理论 + 代码示例。

推荐前置知识

参考资源