ATAC-seq 开放染色质分析¶

一句话概述¶

ATAC-seq 利用转座酶 Tn5 探测开放染色质区域，只需极少量细胞即可全基因组检测染色质可及性，是表观基因组学最流行的技术之一。

核心知识点表格¶

知识点	说明
ATAC-seq	Assay for Transposase-Accessible Chromatin with sequencing
Tn5 转座酶	超活性转座酶，优先插入开放染色质区域并同时添加测序接头
核小体信号	片段长度呈周期性分布：<150bp 为无核小体片段，~200bp 为单核小体
TSS 富集	转录起始位点周围信号富集，是关键 QC 指标（>6 合格，>10 优秀）
FRiP	Fraction of Reads in Peaks，>30% 为良好
Tn5 偏移校正	正链读段+4bp，负链读段-5bp，使切割位点居中
足迹分析	在峰内检测转录因子的物理结合"印迹"
MACS2/MACS3	最常用峰调用工具
Genrich	专为 ATAC-seq 设计的峰调用工具
HMMRATAC	用 HMM 区分信号和背景，专为 ATAC-seq 设计

白话解释原理¶

想象一下： DNA 缠绕在核小体上就像线缠在线轴上。某些区域的线是松开的（开放染色质），这些位置通常有基因在活跃表达。

ATAC-seq 就像派一个"探测兵"（Tn5 转座酶）去试探： - 遇到松开的线（开放区域）→ 插入成功，留下标记 - 遇到缠得很紧的线（紧缩区域）→ 插不进去

最后看哪些位置被插入了标记 → 那就是开放染色质区域。

巧妙之处： Tn5 插入时会同时带上测序接头，省去了传统建库步骤，所以只需 500-50000 个细胞就够了！

各步骤详解¶

第一步：数据预处理¶

# 1. 质量控制
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc/  # 检查双端读段质量

# 2. 去接头（ATAC-seq 经常有接头残留）
trim_galore \
    --paired \                       # 双端模式
    --quality 20 \                   # Q20质量过滤
    --length 20 \                    # 最短20bp
    --fastqc \                       # 自动FastQC
    sample_R1.fastq.gz \             # R1文件
    sample_R2.fastq.gz               # R2文件

# 3. 比对到参考基因组（Bowtie2）
bowtie2 \
    -x /data/genome/hg38_index \     # 参考基因组索引
    -1 sample_R1_val_1.fq.gz \       # 修剪后R1
    -2 sample_R2_val_2.fq.gz \       # 修剪后R2
    --very-sensitive \                # 高灵敏度
    --maxins 2000 \                   # 最大插入片段长度2000bp（ATAC-seq片段可能很长）
    --no-mixed \                      # 不允许只比对一端
    --no-discordant \                 # 不允许不一致配对
    -p 8 \                            # 8线程
    -S sample.sam                     # 输出SAM

# 4. 后处理
# 排序和转BAM
samtools sort -@ 8 -o sample.sorted.bam sample.sam  # 排序
samtools index sample.sorted.bam                      # 建索引

# 去除线粒体reads（ATAC-seq中线粒体比例通常很高）
samtools view -h sample.sorted.bam | \
    grep -v "chrM" | \               # 移除线粒体比对
    samtools sort -@ 8 -o sample.noMT.bam  # 重新排序

# 去除PCR重复
picard MarkDuplicates \
    I=sample.noMT.bam \              # 输入
    O=sample.dedup.bam \             # 输出
    M=dup_metrics.txt \              # 统计
    REMOVE_DUPLICATES=true           # 移除重复

# 只保留高质量的正常配对比对
samtools view -b \
    -f 2 \                            # 正确配对
    -F 1804 \                         # 去除未比对/重复/次要比对/补充比对
    -q 30 \                           # MAPQ >= 30
    sample.dedup.bam \                # 输入
    > sample.clean.bam                # 输出

samtools index sample.clean.bam       # 建索引

# 去除黑名单区域
bedtools intersect -v \
    -abam sample.clean.bam \          # 输入BAM
    -b hg38-blacklist.v2.bed \        # ENCODE黑名单
    > sample.final.bam                # 最终BAM

第二步：质量控制¶

# 1. 片段长度分布（核心QC）
# 应该看到周期性分布: <150bp峰(无核小体), ~200bp峰(单核小体), ~400bp峰(双核小体)
picard CollectInsertSizeMetrics \
    I=sample.final.bam \              # 输入BAM
    O=insert_size_metrics.txt \       # 输出统计
    H=insert_size_histogram.pdf       # 输出直方图

# 2. TSS 富集分数
# 使用 deeptools 计算
computeMatrix reference-point \
    -S sample.bw \                     # bigWig信号文件
    -R genes_tss.bed \                 # TSS坐标
    --referencePoint center \          # 以TSS为中心
    -a 2000 -b 2000 \                  # 上下游各2kb
    -o tss_matrix.gz                   # 输出矩阵

plotProfile \
    -m tss_matrix.gz \                 # 输入矩阵
    -o tss_enrichment.pdf              # 输出图片
# TSS 富集分数 > 6: 合格; > 10: 优秀

# 3. FRiP 计算
total_reads=$(samtools view -c -F 4 sample.final.bam)  # 总比对reads
reads_in_peaks=$(bedtools intersect -u -a sample.final.bam -b peaks.narrowPeak | \
    samtools view -c)  # 峰内reads
echo "FRiP = $(echo "scale=4; $reads_in_peaks / $total_reads" | bc)"  # 计算FRiP
# FRiP > 0.3 (30%): 良好

第三步：峰调用¶

# === 方法1: MACS2/MACS3（最常用） ===
macs2 callpeak \
    -t sample.final.bam \             # 输入BAM（不需要input对照！）
    -f BAMPE \                         # 双端BAM格式（使用真实片段信息）
    -g hs \                            # 人类基因组
    -n atac_sample \                   # 输出前缀
    --outdir peaks/ \                  # 输出目录
    --nomodel \                        # 不建模型（ATAC-seq不需要）
    --shift -100 \                     # 偏移-100bp
    --extsize 200 \                    # 延伸200bp（模拟Tn5切割位点±100bp）
    --keep-dup all \                   # 保留所有reads（已去重）
    -q 0.05 \                          # FDR阈值
    --bdg                              # 输出bedGraph

# 注意: 用 BAMPE 模式时不需要 --shift 和 --extsize
# BAMPE 会直接使用配对读段推断的真实片段

# === 方法2: Genrich（专为ATAC-seq设计） ===
# 安装: conda install -c bioconda genrich
# Genrich 需要 name-sorted BAM
samtools sort -n -@ 8 \
    sample.final.bam \                 # 输入
    -o sample.namesorted.bam           # 按名排序

Genrich \
    -t sample.namesorted.bam \         # 输入name-sorted BAM
    -o peaks_genrich.narrowPeak \      # 输出峰文件
    -j \                               # ATAC-seq模式（使用Tn5切割位点）
    -e chrM \                          # 排除线粒体
    -r \                               # 移除PCR重复
    -v                                 # 显示详细信息

# === 方法3: HMMRATAC（HMM方法） ===
# 专门为ATAC-seq设计，区分open/nucleosome/background状态
HMMRATAC \
    -b sample.final.bam \              # 输入BAM
    -i sample.final.bam.bai \          # BAM索引
    -g hg38.genome \                   # 基因组大小文件
    -o hmmratac_peaks                  # 输出前缀

第四步：Tn5 偏移校正和足迹分析¶

白话解释： Tn5 插入DNA时会产生9bp的间隔，需要校正后才能精确定位转录因子结合位点。

# Tn5 偏移校正（用于足迹分析）
import pysam  # 导入BAM处理库

# 打开BAM文件
bam_in = pysam.AlignmentFile("sample.final.bam", "rb")   # 读取BAM
bam_out = pysam.AlignmentFile("sample.shifted.bam", "wb",  # 写入新BAM
                               template=bam_in)             # 使用相同的header

for read in bam_in:  # 遍历每条read
    if read.is_reverse:                    # 如果是负链
        read.reference_start -= 5          # 负链向左移5bp
    else:                                   # 如果是正链
        read.reference_start += 4          # 正链向右移4bp
    bam_out.write(read)                    # 写入

bam_in.close()   # 关闭输入
bam_out.close()  # 关闭输出

# 然后用 TOBIAS 做足迹分析
# pip install tobias

# 1. 校正Tn5偏差
# TOBIAS ATACorrect \
#     --bam sample.shifted.bam \      # 偏移校正后的BAM
#     --genome hg38.fa \              # 参考基因组
#     --peaks peaks.bed \             # 峰文件
#     --outdir corrected/             # 输出目录

# 2. 计算足迹分数
# TOBIAS FootprintScores \
#     --signal corrected/sample_corrected.bw \  # 校正后的信号
#     --regions peaks.bed \                      # 峰区域
#     --output footprint_scores.bw               # 足迹分数

# 3. 绑定差异分析
# TOBIAS BINDetect \
#     --motifs JASPAR2024_vertebrates.meme \  # JASPAR motif数据库
#     --signals condition1.bw condition2.bw \  # 两个条件的信号
#     --genome hg38.fa \                       # 参考基因组
#     --peaks peaks.bed \                      # 峰区域
#     --outdir binding_results/                # 输出目录

第五步：差异可及性分析¶

# 使用 DiffBind 或 DESeq2 做差异分析
library(DiffBind)   # 加载差异结合分析包

# 类似 ChIP-seq 分析，但不需要 input 对照
# 1. 创建样本表
samples <- data.frame(
    SampleID = c("WT_1", "WT_2", "KO_1", "KO_2"),   # 样本名
    Condition = c("WT", "WT", "KO", "KO"),             # 条件
    Replicate = c(1, 2, 1, 2),                          # 重复
    bamReads = c("wt1.bam", "wt2.bam", "ko1.bam", "ko2.bam"),  # BAM
    Peaks = c("wt1.narrowPeak", "wt2.narrowPeak",      # 峰
              "ko1.narrowPeak", "ko2.narrowPeak"),
    PeakCaller = rep("narrow", 4)                       # 峰格式
)

# 2. 分析
dba_obj <- dba(sampleSheet = samples)    # 创建对象
dba_obj <- dba.count(dba_obj)            # 计数
dba_obj <- dba.normalize(dba_obj)        # 标准化
dba_obj <- dba.contrast(dba_obj,          # 设对比
                        categories = DBA_CONDITION)
dba_obj <- dba.analyze(dba_obj,           # 差异分析
                        method = DBA_DESEQ2)

# 3. 提取差异可及区域（DARs）
results <- dba.report(dba_obj, th = 0.05)  # FDR < 0.05的DARs
cat("差异可及区域数量:", length(results))    # 打印数量

# 4. 在DARs中做motif富集
# 导出为BED后用HOMER分析
export(results, "DARs.bed")  # 导出BED文件
# findMotifsGenome.pl DARs.bed hg38 motif_output/ -size 200

常见报错与解决¶

报错	原因	解决方案
`>80% mitochondrial reads`	线粒体DNA占比太高	实验优化核提取步骤；分析时过滤chrM
`No nucleosome pattern`	片段大小分布异常	可能 Tn5 浓度不对，或核提取失败
`Low TSS enrichment (<4)`	信号噪声比低	检查实验质量，增加测序深度
`Too many peaks (>300k)`	阈值太松	收紧 q 值阈值(如0.01)或使用更严格模式
`BAMPE error in MACS2`	BAM未正确配对	确保用 `samtools fixmate` 修复配对信息
`Genrich: no alignments`	BAM未按名排序	用 `samtools sort -n` 按名重新排序

速查表¶

┌──────────────────────────────────────────────────────────┐
│                ATAC-seq 分析速查                          │
├──────────────────────────────────────────────────────────┤
│ 流程:                                                    │
│   FASTQ → Trim → Bowtie2 → 去chrM → 去重 → 峰调用       │
│                                                          │
│ QC 标准:                                                 │
│   线粒体比例: <20% (过滤前)                               │
│   TSS 富集: >6 (合格), >10 (优秀)                        │
│   FRiP: >30%                                             │
│   片段分布: 看到核小体周期性                               │
│   NSC/RSC: >1.05 / >0.8                                  │
│                                                          │
│ Tn5 偏移:                                                │
│   正链 +4bp, 负链 -5bp                                   │
│   峰调用和差异分析不需要偏移                               │
│   足迹分析和核小体分析需要偏移                             │
│                                                          │
│ 峰调用工具:                                               │
│   MACS2 -f BAMPE --nomodel -- 最通用                     │
│   Genrich -j -- 专为ATAC设计                             │
│   HMMRATAC -- HMM区分信号/核小体/背景                    │
│                                                          │
│ 足迹分析: TOBIAS / HINT-ATAC                             │
│ Motif分析: HOMER / chromVAR                              │
│ 差异分析: DiffBind / DESeq2                              │
└──────────────────────────────────────────────────────────┘

面试高频问题¶

ATAC-seq 和 ChIP-seq 有什么区别？各自的优缺点？ 答：ATAC-seq 检测的是染色质开放性（所有开放区域），不针对特定蛋白；ChIP-seq 用抗体针对特定蛋白/修饰。ATAC-seq 优点：细胞需求少(500个即可)、操作简单(2小时建库)、无需抗体。缺点：不能区分不同蛋白的结合。ChIP-seq 优点：特异性高。缺点：需大量细胞(通常>10^6)、依赖抗体质量。
ATAC-seq 的 Tn5 偏移校正是什么？什么时候需要做？ 答：Tn5 转座酶插入 DNA 时实际产生9bp的切割间隔(staggered cut)，所以正链读段的起始位置比真实切割位点偏右4bp，负链偏左5bp。校正方法：正链+4bp，负链-5bp。在做足迹分析(footprinting)和核小体定位分析时必须做偏移校正；在做峰调用和差异可及性分析时不需要，因为偏差在bin尺度上可以忽略。
ATAC-seq 数据中为什么线粒体 reads 比例这么高？ 答：线粒体 DNA 没有核小体包裹，是完全"裸露"的，对 Tn5 转座酶来说是最容易插入的目标。正常情况下线粒体 reads 占30-80%。解决方法：实验上优化核提取步骤减少线粒体 DNA 污染；分析时用 samtools view | grep -v chrM 过滤。这些 reads 虽然被丢弃，但不影响核基因组的分析结果。