01 表观组学实践教程

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

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

表观组要回答的核心问题

转录组告诉你"这个基因表达了多少"，表观组告诉你"为什么它能表达 / 不能表达"：

表观组适合的项目：发育、分化、肿瘤进展（耐药机制）、神经精神疾病、细胞身份建立。不适合：单纯比较两个稳态条件的差异（直接做转录组就行，表观组只是补充证据）。

一个项目大致是什么样子

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

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

常见工具栈

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

推荐公开数据集

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

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-08 全部上线。01-04 带可跑脚本，05-08 为理论 + 代码示例。

推荐前置知识

• 编程基础：R、Python、Bash 学习路径
• R 数据整理与 ggplot2 可视化
• 10 scATAC-seq 分析（已有的单细胞版 ATAC 分析）

常见误区

新手做表观组项目最容易的几个判断错误：

误区 1：peak 数量等于实验质量

样本 A 出 50000 个 peak，样本 B 只出 30000 个 — 不一定是 A 测得"更深更好"。peak 数量受测序深度、过滤阈值、call peak 参数严苛程度影响很大。看 FRiP（reads in peaks 比例）、TSS enrichment 这些标准化指标比 peak 数量更可靠。

误区 2：只做 peak 数量比较，不做 peak 强度分析

A 样本和 B 样本 peak 列表 90% 重叠，但同一个 peak 在 A 里可能 reads 数 100，在 B 里只有 10 — 这才是真差异。用 DiffBind 做基于 reads 强度的差异分析，光比较 peak 集合 overlap 信息量不足。

误区 3：把 ChIP-seq 当 RNA-seq 思路做差异

ChIP / ATAC 不能用全转录组的 DESeq2 直接套。peak 是按位置定义的，不同样本的 peak 不一致 — 要先建 consensus peak set，再在 consensus 上数 reads，DiffBind 把这套封装好了。

误区 4：把启动子 peak 当成调控 peak

很多 ChIP / ATAC 的 peak 落在启动子，但真正决定细胞身份的往往是远端增强子（TSS 几十 kb 之外）。注释 peak 时要区分 promoter / enhancer / gene body，不要全归到 "最近基因" 就完事。

误区 5：忽略 input / IgG 对照

ChIP-seq 没有 input 是没法做正经分析的 — peak caller 没办法区分"真信号"和"基因组开放区域的非特异背景"。项目设计阶段必须每个 condition 至少配一个 input。ATAC-seq 不强制要 input，但 IgG 对照能区分非特异 Tn5 切割。

参考资源

• ENCODE Portal
• MACS2 文档
• ChIPseeker Bioconductor 页面
• deepTools 文档
• Bismark 手册
• methylKit 用户指南

02 实验类型与数据格式

表观组学实验类型多，但分析思路可以归为两大类：开放区域检测（ATAC-seq、DNase-seq）和蛋白-DNA 结合检测（ChIP-seq）。两者的分析流程几乎一样，区别在 peak calling 参数和质量指标。

三种主要实验

本专栏前 4 个模块聚焦 ATAC-seq 和 ChIP-seq（它们共享 peak-based 分析框架）。甲基化后续单独开。

从 FASTQ 到 peak 文件

不管是 ATAC 还是 ChIP，从原始数据到可分析的 peak 文件，标准流程是：

FASTQ → fastp (trim) → Bowtie2/BWA (align) → samtools (sort/filter)
     → Picard (dedup) → MACS2 (peak calling) → narrowPeak / broadPeak

BioF3 的表观组教程从 peak 文件开始。如果你需要从 FASTQ 跑起，参考 nf-core/chipseq 或 nf-core/atacseq 流水线。

peak 文件格式

MACS2 输出的 .narrowPeak 是 BED6+4 格式：

chr1  9356548  9356648  peak_1  100  .  5.0  10.5  7.2  50

R 里用 ChIPseeker::readPeakFile() 或 rtracklayer::import() 读入，得到 GRanges 对象。

怎么判断这份数据值不值得分析

拿到 peak 文件先做的不是注释或差异分析，而是判断这份数据本身值不值得花时间分析。三件事看一下：

1. peak 总数在合理范围

数量过少（< 1000）说明 IP 效率低或测序深度不足；数量异常多（> 100 万）多半是没正确去 input。

2. FRiP（fraction of reads in peaks）

FRiP = peak 区域 reads / 总比对 reads。好数据 FRiP > 0.2（ENCODE 严格要求 > 0.3），低于 0.05 基本不可用。这是辨识"真信号 vs 全是噪声"最直接的指标。

3. peak 在 TSS 附近富集

把 peak 中心和 TSS 距离画分布图，应该在 TSS 附近有明显尖峰。ATAC 的 TSS enrichment > 6 是常见门槛。如果 TSS 周围没有富集，说明实验失败或 chromatin 没有被有效打开。

这三件事 5 分钟就能查完，比闷头跑两天后发现数据没法用强多了。

常见坑

坑 1：参考基因组版本不一致

peak 文件、注释 GTF、TxDb 三者必须用同一个基因组版本。hg19 的 peak 配 hg38 的 TxDb 不会报错但坐标全错。看 peak 文件第一行染色体名（"chr1" vs "1"）和 TxDb 文档对一下。

坑 2：narrowPeak 当 broadPeak 用

TF / DNase / ATAC 用 narrowPeak（带 summit 列）；组蛋白修饰（H3K27me3、H3K4me3 etc）用 broadPeak（没有 summit 列）。强行把 broadPeak 当 narrowPeak 读会丢失信息，downstream motif 分析无法定位准确位置。

坑 3：peak 文件没有按染色体大小过滤

线粒体、未定位的 contig（chrUn_）、ENCODE blacklist 区域常出现假阳性 peak。call peak 后用 ENCODE blacklist BED 文件做一次 bedtools intersect -v 过滤掉，能去掉 5-15% 的假信号。

坑 4：把 BAM 里 MAPQ 低的 reads 也用了

ATAC / ChIP 的标准要求 MAPQ ≥ 30（Bowtie2 default）或 ≥ 10。MAPQ 低的 reads 可能是多处比对（重复区域），保留会让重复区域出现假 peak。

坑 5：paired-end ATAC 用 single-end 模式 call peak

-f BAM 是 single-end 模式，但 ATAC 是 paired-end，会丢一半信息。paired-end ATAC 必须 -f BAMPE，这一点和 ChIP-seq 不一样。

下一步

接着深入：

• 02 Peak 注释与多样本比较 — 把 peak 文件里那一坨基因组坐标变成"哪些基因"
• 06 ATAC-seq 分析要点 — 如果你的数据是 ATAC，先看这一篇了解关键差异

横向延伸：

• 02 原始数据处理与 Cell Ranger — bulk ATAC 的 BAM 处理思路和单细胞 Cell Ranger 类似
• nf-core/chipseq 流水线 — 完整的 FASTQ → peak 自动化方案

参考资源

• MACS2 文档
• nf-core/chipseq
• nf-core/atacseq
• ENCODE 实验标准

03 Peak 注释与多样本比较

拿到 peak 文件之后，第一个问题是"这些 peak 落在基因组的什么位置"。ChIPseeker 把 peak 注释到最近的基因、标注它在启动子 / 内含子 / 基因间等区域，一步出图。

本章用 ChIPseeker 自带的 AR（雄激素受体）ChIP-seq 数据演示：3 个剂量（0M / 1nM / 100nM）的 peak 文件，看剂量增加时 peak 数量、位置分布和关联基因如何变化。

核心流程

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

peak <- readPeakFile("peaks.narrowPeak")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

anno <- annotatePeak(peak, TxDb = txdb, tssRegion = c(-3000, 3000))
plotAnnoPie(anno)
plotDistToTSS(anno)

tssRegion = c(-3000, 3000) 定义"启动子"的范围：TSS 上下游 3kb 以内的 peak 算启动子区域。

真实示例

配套脚本 epi02_chipseeker_sci.R 在 ChIPseeker 内置的 AR ChIP-seq 数据上跑完整流程：

Rscript scripts/epigenomics/epi02_chipseeker_sci.R

每张图看什么

图 1：AR 100nM 的 peak 落在哪些基因组区域。启动子占比越高，说明这个 TF 越倾向于结合在基因起始位点附近。AR 是经典的启动子 + 增强子结合 TF，所以启动子和远端基因间区域都有不少。

图 2：三个剂量的 peak 到最近 TSS 的距离分布。剂量越高，靠近 TSS 的 peak 比例越大 —— 说明高剂量下 AR 更多地占据启动子。

图 3：三个剂量的基因组区域分布对比。堆叠条形图一眼看出"启动子占比随剂量增加"。

图 4：三个剂量的 peak 关联基因重叠。"三个剂量都有"的基因是 AR 的核心靶基因；"只在 100nM 出现"的是高剂量特异的。

图 5：AR 100nM peak 关联基因的 GO BP 富集。应该能看到雄激素响应、细胞增殖调控等通路。

图 6：每条染色体上的 peak 数量。分布大致和染色体大小成正比，但某些染色体可能因为基因密度高而偏多。

套到自己数据上

把 getSampleFiles() 换成自己的 .narrowPeak 文件路径即可。注意：

• 参考基因组版本要和 peak 文件一致（hg19 / hg38 / mm10）
• 对应的 TxDb 包：TxDb.Hsapiens.UCSC.hg38.knownGene（hg38）或 TxDb.Mmusculus.UCSC.mm10.knownGene（小鼠）
• 如果是 broadPeak（组蛋白修饰），annotatePeak 的参数不用改，但解读时"启动子"的含义要注意

常见坑

坑 1：tssRegion 范围拍脑袋定

默认 tssRegion = c(-3000, 3000)，但组蛋白活化标记 H3K4me3 通常在 ±1kb 内最强，K27ac 可以扩到 ±5kb。范围拍错了会让"启动子" peak 占比虚高或虚低。根据自己实验类型调，TF 一般 ±2kb 够用。

坑 2：peak 注释到"最近基因"当成调控关系

peak 离哪个基因最近 ≠ peak 调控哪个基因。增强子可以跨数十 kb 调控其他基因，严格的调控关系需要 Hi-C / 4C / ABC model 等额外数据。把"最近基因"当作初筛候选可以，写在文章里要明确说"基于 proximity 的关联"。

坑 3：多样本对比时 peak 数量直接比

样本 A 在更深的测序下 call 出更多 peak，不代表 A 比 B"信号更强"。用 consensus peak set 上的 reads 数（DiffBind 思路）做基于强度的对比，光比较 peak 数量是误导。

坑 4：GO 富集背景错

ChIPseeker 的 enrichGO 默认背景是物种全部基因，但你的 ChIP/ATAC 只在某些基因上有信号，应该用"所有可被检测到 peak 的基因"做 universe。enrichGO(universe = all_peak_genes) 显式传背景。

坑 5：染色体名前缀不一致

peak 文件可能用 "chr1"，但有些 GTF 用 "1"。ChIPseeker 不会自动统一，要么用 seqlevelsStyle() 改一致，要么 manually seqlevels(peak) <- paste0("chr", seqlevels(peak))。否则 annotatePeak 会找不到任何基因。

下一步

接着深入：

• 03 DiffBind 差异结合分析 — 在 consensus peak 上做差异强度分析
• 05 Motif 富集与 HOMER — 知道 peak 在哪后下一步：哪些 TF motif 在这些 peak 富集

横向延伸：

• 04 Peak 可视化与多样本比较 — 注释结果做成可发表的图
• bulk RNA-seq 03 GO/KEGG/GSEA — peak 关联基因的富集分析思路完全一致

下载资源

参考资源

• ChIPseeker Bioconductor 文档
• Yu et al. 2015, ChIPseeker 论文
• TxDb 包列表

04 DiffBind 差异结合分析

ChIP-seq 的差异分析和 RNA-seq 思路一样：把每个 peak 在每个样本里的 reads 数当作"表达量"，用 DESeq2 或 edgeR 做统计检验。DiffBind 把这条流程封装成几个函数，从 peak 文件 + BAM 到差异结合位点一步到位。

本章用 DiffBind 自带的 tamoxifen 数据演示：11 个乳腺癌细胞系的 ER（雌激素受体）ChIP-seq，分为 Responsive（对他莫昔芬敏感）和 Resistant（耐药）两组。

真实示例

配套脚本 epi03_diffbind_sci.R 在 tamoxifen 数据上跑完整的差异结合分析：

Rscript scripts/epigenomics/epi03_diffbind_sci.R

每张图看什么

图 1：样本间 binding affinity 的相关性热图。同一条件的样本应该聚在一起。

图 2：PCA 按条件着色。Responsive 和 Resistant 在 PC1 上分开。

图 3：MA plot。红点是 FDR < 0.05 的差异结合位点。

图 4：火山图。横轴 log2FC，纵轴 -log10(FDR)。

图 5：差异结合位点的 binding affinity 热图。每行一个 peak，每列一个样本。

图 6：差异位点按方向分：Gained（Responsive 里更强）vs Lost（Resistant 里更强）。

核心代码

library(DiffBind)
data(tamoxifen_counts)

# 设置对比
tamoxifen <- dba.contrast(tamoxifen, categories = DBA_CONDITION)

# 差异分析（DESeq2 后端）
tamoxifen <- dba.analyze(tamoxifen, method = DBA_DESEQ2)

# 提取结果
db_report <- dba.report(tamoxifen)

为什么不直接用 DESeq2 而要走 DiffBind

DESeq2 输入要"行 = 基因"的固定矩阵，而 ChIP / ATAC 不同样本的 peak 集合是不一样的。要做差异分析必须先建一个所有样本共用的 consensus peak set，再在每个 peak 上数 reads。DiffBind 把这一步封装了：

1. dba() 读 sample sheet → 知道每个样本的 BAM + peak
2. dba.count() 在 consensus peaks 上数 reads，得到 peak × sample 矩阵
3. dba.analyze() 把这个矩阵交给 DESeq2 / edgeR 做差异

后端就是 DESeq2 / edgeR — 你也可以自己跑，但 DiffBind 处理了 consensus peak 构建、normalization（reads-in-peak vs library-size）、blacklist 过滤等细节。新项目直接用 DiffBind，自己拼这套很容易出错。

常见坑

坑 1：sample sheet 的 BAM 路径写错

DiffBind 的 sample sheet 里要求 bamReads 列指向已经去重 / 过滤好的 BAM。不是原始 BAM，否则 reads-in-peak 计数被 PCR duplicate / 低 MAPQ reads 污染。先做完 Picard MarkDuplicates + samtools view -q 30 再喂给 DiffBind。

坑 2：normalization 选错

DiffBind 默认按 library size 归一化，但 ChIP-seq 一些情况下应该按 reads-in-peak 归一化（spike-in 实验、或样本间总信号变化大时）。dba.normalize(method = DBA_NORM_RLE / DBA_NORM_TMM / DBA_NORM_LIB) 三个选项各有适用场景，看 vignette 不要直接用默认。

坑 3：consensus peaks 用了所有样本的并集

默认 minOverlap = 2 意味着一个 peak 至少在 2 个样本里出现才进 consensus。如果你只想保留高置信 peak，提高到 minOverlap = (n_samples / 2) + 1，结果更稳。

坑 4：FDR 阈值不区分 ChIP / ATAC

ChIP-seq DB 通常 FDR < 0.05 + |FC| > 1 就够；ATAC-seq 信号更稀疏，可以放到 FDR < 0.1。死磕 0.05 在 ATAC 上经常 0 个 DB peak，不是分析做错了。

坑 5：Greylist 没用

Greylist（区域级别的"看上去像 peak 但其实是噪声"区域）和 Blacklist 不同。DiffBind 默认会调用 GreyListChIP，需要联网下载参考。国内有时下载慢，可以预先下载好 cached。

下一步

接着深入：

• 04 Peak 可视化与多样本比较 — DiffBind 输出做成可发表的图
• 05 Motif 富集与 HOMER — 差异 peak 富集了什么 TF motif

横向延伸：

• bulk RNA-seq 02 DESeq2 — DiffBind 后端就是 DESeq2，看一遍 RNA-seq 那边对理解 DiffBind 大有帮助
• 08 表观组与转录组的整合 — 差异 peak 的最终意义在于和差异基因的关系

下载资源

参考资源

• DiffBind vignette
• Ross-Innes et al. 2012, tamoxifen 数据源

05 Peak 可视化与多样本比较

本章把 ChIPseeker + DiffBind 的结果变成能直接用于论文的图。

配套脚本 epi04_visualization_sci.R 输出 6 张图：

Rscript scripts/epigenomics/epi04_visualization_sci.R

每张图看什么

图 1：每个样本的 peak 宽度分布。TF ChIP-seq 的 narrow peak 通常 100-500bp；组蛋白修饰的 broad peak 可以到几 kb。

图 2：5 个样本的基因组区域分布对比。AR（TF）和 CBX6/CBX7（chromatin reader）的分布模式不同。

图 3：AR 100nM 的 peak 在全基因组上的覆盖密度。某些染色体区域 peak 特别密集，可能对应 super-enhancer 或基因密集区。

图 4：AR 三个剂量的 consensus peak 重叠。"三个剂量都有"的是核心结合位点；"只在 100nM 出现"的是剂量依赖的新增位点。

图 5：DiffBind 差异结合位点的基因组区域注释。差异位点主要落在哪里（启动子 vs 增强子）决定了它们的功能解读方向。

图 6：差异结合位点的 fold change vs 到最近 TSS 的距离。靠近 TSS 的位点（启动子区域）和远离 TSS 的位点（增强子区域）可能有不同的 fold change 分布。

下载资源

常见坑

坑 1：peak 宽度做 boxplot 不取 log

peak 宽度分布常跨 2-3 个数量级（200bp 到 100kb），直接 boxplot 大部分箱子被压扁。scale_y_log10() 是默认操作。

坑 2：基因组区域饼图比例没归一化

10 个样本各自的 peak 数差很多，直接画占比饼图看上去都差不多，真实差异被绝对数量稀释。改成堆叠条形图 + 百分比标签。

坑 3：覆盖度 plot 用 raw counts

不同样本测序深度差几倍，raw count 覆盖度图没有可比性。用 deepTools 的 bamCoverage --normalizeUsing CPM/RPGC 做归一化后的 bigWig，再画。

坑 4：fold change vs distance to TSS 不区分方向

差异 peak 的 FC 有正有负（gained vs lost），混在一张图里看趋势会被相互抵消。用 facet 或颜色把 gained / lost 分开。

坑 5：导出 PNG 时分辨率默认

ggsave("fig.png") 默认 dpi=300、单位 inch — 对 panel 图够用，但单图会偏小。论文用 dpi = 600, width = 6, height = 4, units = "in"。

下一步

接着深入：

• 05 Motif 富集与 HOMER — 把 peak 上的序列拿去做 motif 富集
• 08 表观组与转录组的整合 — 这些可视化最终都要回到表达数据上

横向延伸：

• FigCode peak-anno — peak 注释饼图的可重用模板
• FigCode diffbind-ma — DiffBind MA 图的通用代码

参考资源

• ChIPseeker 文档
• DiffBind 文档
• deepTools（Python 端覆盖度可视化）

06 Motif 富集与 HOMER

Peak 本身只是"这里有信号"，motif 分析回答的是"什么转录因子可能在这里结合"。思路是：从 peak 序列里统计过表达的短序列模式（motif），和已知 TF motif 数据库比对。

常用工具

HOMER 典型用法

# 安装 HOMER（一次性）
# http://homer.ucsd.edu/homer/introduction/install.html

# 对 narrowPeak 做 motif 富集
findMotifsGenome.pl peaks.narrowPeak hg38 motif_output/ \
  -size 200 -mask -p 8

-size 200 表示取 peak summit 两侧各 100bp 做分析。输出目录里会有：

• knownResults.html — 已知 motif 的富集排名
• homerResults.html — de novo 发现的 motif
• 每个 motif 的 logo 图

R 里做 motif scanning

library(motifmatchr)
library(TFBSTools)
library(JASPAR2020)

# 获取 JASPAR 的人类 TF motif
pfm_list <- getMatrixSet(JASPAR2020, opts = list(species = "Homo sapiens"))

# 在 peak 序列里扫描 motif
library(BSgenome.Hsapiens.UCSC.hg38)
motif_hits <- matchMotifs(pfm_list, peaks_gr, genome = BSgenome.Hsapiens.UCSC.hg38)

# 统计每个 motif 在 DE peaks vs non-DE peaks 里的富集

这种方式的好处是能和 DiffBind 的差异结合结果直接联动：只看"在耐药细胞里 gained 的 peak 富集了什么 motif"。

解读要点

• 已知 motif 富集：如果你做的是 AR ChIP-seq，排第一的应该是 ARE（androgen response element）。如果不是，说明实验或分析有问题。
• de novo motif：HOMER 会尝试从头发现新 motif。如果 de novo 结果和已知 motif 高度相似，说明信号很强。
• 背景选择：默认用随机基因组区域做背景。如果你的 peak 集中在启动子，用"所有启动子"做背景会更严格。

motif 富集的常见误判

motif 分析容易出"看上去显著但其实没意义"的结果。三件事注意一下：

1. 富集 motif 不一定是"驱动" TF

motif 富集只是说"这些 peak 序列上有这种 motif"，不等于该 TF 真的在这些 peak 上结合。一个 motif 可能被多个 TF 家族识别（e.g. AP-1 家族），具体哪个亚型起作用要靠 ChIP-seq 验证。

2. 富集结果受背景集严重影响

用全基因组随机区域当背景 → 启动子相关 motif 都会显著（启动子本身就富集 motif）；用同等数量的非 DB peak 当背景 → 才能看出"DB 特异的 motif"。设计 motif 分析前先想清楚要回答的问题再选背景。

3. CG 含量需要匹配

Peak 的 CG 含量高，富集出的 motif 也会偏向高 CG（CpG 岛特征）。HOMER 的 -mset 或 -bg 选项可以传匹配 CG 含量的背景，避免 GC bias 假阳性。

常见坑

坑 1：HOMER 装在 conda 但参考基因组没下

findMotifsGenome.pl peaks.bed hg38 ... 需要 HOMER 自己下载 hg38 reference (configureHomer.pl -install hg38)。首次跑会卡很久，国内有时下载失败。预先 configureHomer.pl -list 看可用版本。

坑 2：peak 数量太少（< 200）出不来 motif

motif 算法需要足够多的 peak 序列做统计。< 100 个 peak 几乎跑不出有意义的结果。peak 太少时改用 known motif scanning（每个 peak 单独对照已知 motif）。

坑 3：把"motif 富集"和"motif 强度"混淆

"显著富集" = 出现频率高于背景；"motif 强度"= 序列和共识 motif 的相似度。HOMER Motif.Score 是后者，rank by % of Targets 是前者。报告时分清楚。

坑 4：motif logo 图嵌进 PPT 时模糊

HOMER 默认输出的 logo PNG 分辨率低。用 MEME / Tomtom 出 SVG，或直接 R 里 seqLogo 包重画。

坑 5：JASPAR / HOCOMOCO / CIS-BP 的 motif 名字不一致

跨数据库比较时 STAT1 在 JASPAR 是 MA0137.3，在 HOMER 是 STAT1(IPR007726)。做整合前先 unify 成统一 ID 系统（一般用 TF symbol）。

下一步

接着深入：

• 06 ATAC-seq 分析要点 — ATAC 的 peak motif 分析有自己的细节
• 08 表观组与转录组的整合 — motif 富集的 TF 是否在 RNA-seq 里也有差异

横向延伸：

• 03 DiffBind 差异结合分析 — motif 分析最有价值是和 differential peak 联动
• HOMER 官方文档 — 完整的命令行参数手册

参考资源

• HOMER 官方文档
• MEME Suite
• JASPAR 数据库
• motifmatchr Bioconductor

07 ATAC-seq 分析要点

ATAC-seq 和 ChIP-seq 共享大部分分析流程（比对 → peak calling → 注释 → 差异），但有几个关键差异需要单独说明。

和 ChIP-seq 的区别

Fragment size 分布

ATAC-seq 最重要的 QC 图是 fragment size 分布：

~100-150 bp  → nucleosome-free fragments（信号最好的部分）
~200 bp      → mono-nucleosome
~400 bp      → di-nucleosome
~600 bp      → tri-nucleosome

好的 ATAC-seq 数据应该在 < 150bp 处有一个明显的 peak（nucleosome-free），然后在 ~200bp 处有第二个 peak。如果第一个 peak 不明显，说明 Tn5 转座效率低或者细胞核裂解不充分。

MACS2 参数

macs2 callpeak \
  -t sample.bam \
  -f BAMPE \
  --nomodel \
  --shift -100 --extsize 200 \
  -g hs \
  -n sample_atac \
  --keep-dup all \
  -q 0.05

关键参数：

• -f BAMPE：paired-end 模式，用实际 fragment 长度
• --nomodel --shift -100 --extsize 200：不做 fragment size 建模，直接用 Tn5 切割位点两侧 100bp
• --keep-dup all：ATAC-seq 的 PCR duplicate 已经在前面用 Picard 去过了

和单细胞 scATAC 的关系

单细胞实践 10 scATAC-seq 里用的 Signac 流程，底层思路和 bulk ATAC 一样（TF-IDF + LSI），只是把"每个样本"换成了"每个细胞"。bulk ATAC 的 peak 可以直接作为 scATAC 的参考 peak set。

推荐流水线

如果不想手动跑每一步，推荐用 nf-core/atacseq：

nextflow run nf-core/atacseq \
  --input samplesheet.csv \
  --genome GRCh38 \
  --outdir results/

它会自动完成 trim → align → dedup → shift → peak call → QC report 全流程。

常见坑

坑 1：忘了线粒体 reads 过滤

ATAC-seq 数据里 20-50% 的 reads 来自线粒体（mtDNA 是裸露的，Tn5 切得猛）。如果不过滤会让 peak calling 严重偏向线粒体区域。比对后 samtools view -h ... | grep -v chrM 过掉。

坑 2：Tn5 偏移没校正

Tn5 转座插入位置不是 read 起点 — 它在 R1 后偏移 +4bp，R2 前偏移 -5bp。精细 footprint 分析必须做 shift：deepTools alignmentSieve --ATACshift，或 MACS2 --shift -100 --extsize 200。

坑 3：低质量样本 fragment size 双峰不明显

健康 ATAC 在 fragment size 直方图上有 < 150bp 主峰 + ~200bp 次峰。双峰平坦或丢失说明实验失败（细胞核没裂解、Tn5 用量不对）。继续分析也是浪费时间，这一步要果断。

坑 4：把 input 当 ATAC 的 control

ATAC 没有 input 这个概念。某些教程让你把"基因组 DNA"当 control，但 ATAC 本身就是基因组开放区域的富集，用 input 做 control 反而抹掉信号。直接 macs2 callpeak -t sample.bam 不要 -c。

坑 5：用 RNA-seq 标准做 ATAC 差异分析

差异 ATAC 不应该用 DESeq2 直接套基因。要用 DiffBind 或 csaw，输入是 consensus peak set 上的 reads count，不是基因 count。

下一步

接着深入：

• 03 DiffBind 差异结合分析 — ATAC 的差异分析也走 DiffBind
• 05 Motif 富集与 HOMER — ATAC peak 的 motif 富集能找到驱动 TF

横向延伸：

• 单细胞 10 scATAC-seq 分析 — 把 bulk ATAC 思路扩展到单细胞
• 08 表观组与转录组的整合 — ATAC + RNA-seq 联合分析

参考资源

• ENCODE ATAC-seq pipeline
• nf-core/atacseq
• Buenrostro et al. 2013, ATAC-seq 原始论文
• Yan et al. 2020, ATAC-seq 分析最佳实践

08 DNA 甲基化分析入门

DNA 甲基化（主要是 CpG 位点的 5-甲基胞嘧啶）是最稳定的表观修饰之一。和 ChIP/ATAC 不同，甲基化分析不做 peak calling，而是直接量化每个 CpG 位点的甲基化率（0~100%）。

实验类型

BioF3 这一章聚焦 WGBS/RRBS 的 bisulfite-seq 分析。芯片数据用 minfi 包处理，思路不同。

分析流程

FASTQ → Bismark (bisulfite-aware alignment) → methylation extraction
     → per-CpG methylation table → DMR calling (methylKit / DSS / dmrseq)

Bismark 比对

# 建立 bisulfite 索引（一次性）
bismark_genome_preparation --bowtie2 reference/

# 比对
bismark --genome reference/ -1 sample_R1.fq.gz -2 sample_R2.fq.gz

# 去重
deduplicate_bismark sample_pe.bam

# 提取甲基化信息
bismark_methylation_extractor --paired-end --comprehensive --cytosine_report \
  --genome_folder reference/ sample_pe.deduplicated.bam

输出的 CpG_context 文件每行一个 CpG 位点，包含染色体、位置、甲基化 reads 数、非甲基化 reads 数。

R 里做差异甲基化

library(methylKit)

# 读入 Bismark 的 CpG 报告
file_list <- list("sample1.CpG_report.txt", "sample2.CpG_report.txt",
                  "sample3.CpG_report.txt", "sample4.CpG_report.txt")
obj <- methRead(file_list,
                sample.id = list("ctrl1","ctrl2","trt1","trt2"),
                assembly = "hg38",
                treatment = c(0, 0, 1, 1),
                context = "CpG",
                mincov = 10)

# 合并所有样本的 CpG 位点
meth <- unite(obj, destrand = FALSE)

# 差异甲基化位点
diff <- calculateDiffMeth(meth)
diff_25 <- getMethylDiff(diff, difference = 25, qvalue = 0.01)

# 差异甲基化区域（DMR）
# 用 tileMethylCounts 把基因组分成 1kb 窗口再做差异
tiles <- tileMethylCounts(obj, win.size = 1000, step.size = 1000)
meth_tiles <- unite(tiles)
diff_tiles <- calculateDiffMeth(meth_tiles)

关键概念

• 甲基化率：某个 CpG 位点被甲基化的 reads 占总 reads 的比例
• DMC（Differentially Methylated Cytosine）：单个 CpG 位点的差异
• DMR（Differentially Methylated Region）：连续多个 CpG 位点一起变化的区域
• 覆盖度过滤：覆盖度 < 10x 的位点噪声太大，通常过滤掉

和表达的关系

启动子区域的高甲基化通常和基因沉默相关；基因体内的甲基化和活跃转录正相关。把 DMR 和 RNA-seq 的差异基因做交叉，能找到"甲基化变化驱动表达变化"的候选基因。

但要注意：甲基化和表达的关系不是简单一对一。

实务上：先把 DMR 注释到 promoter / gene body / enhancer / intergenic 各自分类，再分别和 RNA-seq 对比。整体混在一起看会得到无意义的"弱相关"。

常见坑

坑 1：没做 bisulfite 转化率 QC

WGBS 实验的 bisulfite 转化率应该 > 99%。未甲基化的 lambda phage spike-in 是金标准，转化率低意味着大量未甲基化 C 被当成 5mC，所有结果系统偏移。Bismark 输出会自动算这个值，看一眼。

坑 2：覆盖度不足直接做差异

WGBS 单 CpG 位点的可信差异需要每位点至少 10× 覆盖度。一份 30× WGBS 看上去深，但 CpG 只占基因组 1%，单位点覆盖度会被稀释。先 methRead(..., mincov = 10) 过滤。

坑 3：用 DMC 当 DMR 报告

单个 CpG（DMC）噪声大，真正可靠的是连续多个 CpG 一起变化的 DMR。methylKit tileMethylCounts 或 dmrseq 找区域，而不是把 DMC 直接报告。

坑 4：跨平台数据强行合并

WGBS / RRBS / 450K / EPIC 的位点集合不重合（450K 只覆盖 ~2% CpG）。跨平台 meta 分析需要先取交集 + 再 normalize，直接 cbind 是错的。

坑 5：differential = 25% 的阈值死搬

methylKit 默认 difference = 25 意为甲基化率差 ≥ 25%。对小效应（药物处理、衰老）来说 25% 太严，差不到 5% 但显著的位点经常有意义。从分布看实际差异范围再选阈值，不要照搬教程。

下一步

接着深入：

• 08 表观组与转录组的整合 — DMR 最有价值的解读是和差异表达基因的关系
• 多组学整合 module01 — 甲基化是整合分析里最常用的 omics 之一

横向延伸：

• Bismark 完整文档 — bisulfite 比对的细节
• methylKit 官方教程

参考资源

• Bismark 手册
• methylKit 教程
• DSS 差异甲基化
• dmrseq
• minfi（芯片数据）

09 表观组与转录组的整合

表观修饰（开放染色质、TF 结合、甲基化）最终要通过影响基因表达来发挥功能。把表观组数据和转录组数据放在一起看，能回答"哪些表观变化真正驱动了表达变化"。

常见整合策略

Peak-gene 关联（最简单）

library(ChIPseeker)

# 差异 peak 注释到最近基因
db_anno <- annotatePeak(db_peaks, TxDb = txdb)
db_genes <- unique(as.data.frame(db_anno)$SYMBOL)

# 差异基因（来自 DESeq2）
de_genes <- res_df$SYMBOL[res_df$padj < 0.05]

# 交集
overlap <- intersect(db_genes, de_genes)
cat("Overlap:", length(overlap), "genes\n")

# Fisher 检验看是否显著富集
fisher.test(matrix(c(
  length(overlap),
  length(setdiff(db_genes, de_genes)),
  length(setdiff(de_genes, db_genes)),
  total_genes - length(union(db_genes, de_genes))
), nrow = 2))

如果 overlap 显著大于随机期望，说明表观变化和表达变化确实有关联。

方向一致性检查

更严格的验证：不仅看"有没有重叠"，还看"方向是否一致"：

# 合并 peak FC 和 gene FC
merged <- inner_join(
  data.frame(gene = db_genes_df$SYMBOL, peak_fc = db_genes_df$Fold),
  data.frame(gene = de_df$SYMBOL, rna_fc = de_df$log2FoldChange)
)

# 散点图：peak FC vs RNA FC
ggplot(merged, aes(x = peak_fc, y = rna_fc)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = "Peak log2FC (ChIP/ATAC)", y = "RNA log2FC")

正相关 = 开放区域增加的基因表达也增加（符合预期）。如果是甲基化数据，启动子区域应该是负相关（甲基化增加 → 表达下降）。

SCENIC：从 scATAC + scRNA 推断调控网络

如果有配对的单细胞数据（10x Multiome 或分别测的 scRNA + scATAC），SCENIC+ 能推断出"哪个 TF 通过哪个增强子调控哪个基因"：

# Python (pySCENIC+)
import scenicplus

# 输入：scRNA AnnData + scATAC AnnData + motif 数据库
# 输出：TF → enhancer → gene 的三元组网络

这是目前单细胞表观组整合的最前沿方向，计算量大但信息量也最大。

实用建议

1. 先做简单的 peak-gene overlap，确认方向一致性
2. 如果 overlap 显著且方向一致，再做更复杂的网络推断
3. 多组学整合的结果要用独立实验验证（比如 CRISPRi 敲掉某个增强子看表达是否下降）
4. 不要过度解读"相关性 = 因果性"

整合分析的预期值要校准

新手做 ATAC + RNA 整合常常失望："只有 30% 的差异 ATAC peak 关联的基因在 RNA-seq 里也是差异的"。这其实是正常的，甚至是好的：

ATAC 变化先于 RNA 变化、增强子调控不一定立即翻译成表达 — 这些都是真实的生物学。期望"100% 一致"是统计学上的过度自信。

常见坑

坑 1：peak-gene 关联用 nearest-gene 当唯一标准

最近基因 ≠ 真调控基因。严格的关联需要 Hi-C / ABC model / CRISPRi 验证。教程级用 nearest 可以，发表级要补 chromatin loop 数据或至少加距离过滤（比如 < 100kb）。

坑 2：不区分启动子 vs 增强子 peak

启动子 peak 和增强子 peak 的 RNA 关联机制不同（启动子直接，增强子要先经过 looping）。整合时应该按 peak 类型分开做关联，混着看会稀释信号。

坑 3：方向一致性看的是 raw correlation

raw 表达矩阵和 raw peak signal 的 correlation 受 batch / library size 影响大。要先 normalize（vst / TMM）再算 correlation。

坑 4：把"显著富集"误读成"重要驱动"

Fisher test 显著只说明重叠不是随机的，不说明这个 overlap 集合里的基因就是关键调控者。后续要做 motif 富集 + KO 验证才能锁定 driver。

坑 5：SCENIC / MOFA2 跑出结果不验证就上文章

这类高级方法对参数敏感，跑两次结果会不一样。任何整合分析的结论都要用独立数据集（公开数据 + own 验证）confirm，不要依赖单次跑出的网络。

下一步

接着深入：

• 多组学整合 overview — 系统的多组学方法学专栏
• 多组学整合 module01 — 项目设计与样本匹配

横向延伸：

• 03 DiffBind 差异结合分析 — 整合的输入是各 omics 的差异表
• bulk RNA-seq 02 DESeq2 — 表达层的差异分析

参考资源

• SCENIC+
• MOFA2
• LOLA（区域富集）
• Corces et al. 2018, ATAC + RNA 整合