染色质可及性分析ATAC-seq进阶¶

一句话概述：ATAC-seq通过转座酶Tn5"打碎"开放的染色质区域来检测基因组哪些地方是"打开的"——那些调控基因表达的开关（启动子、增强子）就藏在这些开放区域里。

核心知识点表¶

知识点	白话解释	重要程度
染色质可及性	DNA是"打开"还是"关闭"的状态	⭐⭐⭐⭐⭐
ATAC-seq	用Tn5转座酶检测开放染色质的技术	⭐⭐⭐⭐⭐
Peak calling	找出信号显著富集的区域（开放区域）	⭐⭐⭐⭐⭐
核小体定位	ATAC-seq片段大小反映核小体位置	⭐⭐⭐⭐
转录因子足迹	开放区域中转录因子结合留下的"脚印"	⭐⭐⭐⭐
差异可及性分析	比较不同条件下哪些区域开放程度改变	⭐⭐⭐⭐

一、ATAC-seq原理¶

ATAC-seq = "测试DNA哪里是打开的"

原理类比：
  想象DNA缠绕在蛋白质上，像毛线缠在线轴上
  - 缠紧的地方 = 关闭的染色质 → Tn5进不去 → 没信号
  - 松开的地方 = 开放的染色质 → Tn5能插入 → 有信号

Tn5转座酶的作用：
  1. 带着测序接头，只能插入"松开"的DNA
  2. 插入后把DNA剪断
  3. 测序这些碎片，就知道哪里是"松开"的

片段大小的信息（进阶知识）：
  < 100bp    → 无核小体区域(NFR)，转录因子结合位点
  ~200bp     → 1个核小体
  ~400bp     → 2个核小体
  ~600bp     → 3个核小体

  所以ATAC-seq不仅能看开放性，还能看核小体定位！

二、ATAC-seq分析流程¶

2.1 环境准备¶

# 创建ATAC-seq分析环境
conda create -n atacseq python=3.10  # 创建环境
conda activate atacseq  # 激活

# 安装所有需要的工具
conda install -c bioconda \
    fastqc trim-galore \    # 质控
    bowtie2 samtools \      # 比对
    picard \                # 去重
    macs2 \                 # Peak calling
    deeptools \             # 信号可视化
    bedtools \              # 区间操作
    homer                   # motif分析

# Python包
pip install pysam pyBigWig pybedtools  # Python分析包

2.2 数据预处理¶

# ========== 质控与Trim ==========
trim_galore \
    --paired \
    --quality 20 \
    --length 20 \           # ATAC-seq片段可能很短
    --fastqc \
    --cores 4 \
    -o trimmed/ \
    sample_R1.fq.gz sample_R2.fq.gz

# ========== 比对 ==========
# ATAC-seq推荐用Bowtie2
bowtie2 \
    -x ~/genomes/hg38/bt2_index \  # Bowtie2索引
    -1 trimmed/sample_R1_val_1.fq.gz \
    -2 trimmed/sample_R2_val_2.fq.gz \
    --very-sensitive \       # 高灵敏度模式
    --no-mixed \             # 不允许单端比对
    --no-discordant \        # 不允许不一致比对
    -X 2000 \                # 最大片段长度2000bp
    --threads 8 \            # 8线程
    -S sample.sam 2> sample_bowtie2.log  # 输出SAM + 日志

# SAM转BAM并排序
samtools view -bS -q 30 sample.sam | \  # 过滤低质量比对(MAPQ≥30)
    samtools sort -@ 8 -o sample_sorted.bam  # 排序
samtools index sample_sorted.bam  # 建索引

# ========== 去除线粒体reads ==========
# ATAC-seq中线粒体reads很多（正常现象）
samtools idxstats sample_sorted.bam | grep "chrM"  # 查看线粒体reads比例

# 去除线粒体reads
samtools view -b sample_sorted.bam \
    $(samtools idxstats sample_sorted.bam | grep -v "chrM" | cut -f1 | tr '\n' ' ') \
    > sample_noMT.bam  # 只保留非线粒体reads
samtools index sample_noMT.bam  # 重建索引

# ========== PCR去重 ==========
picard MarkDuplicates \
    I=sample_noMT.bam \
    O=sample_dedup.bam \
    M=dup_metrics.txt \
    REMOVE_DUPLICATES=true  # 直接删除重复reads
samtools index sample_dedup.bam  # 建索引

# ========== 黑名单区域过滤 ==========
# 下载ENCODE黑名单
wget https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz
gunzip hg38-blacklist.v2.bed.gz

# 去除黑名单区域
bedtools intersect -v \
    -a sample_dedup.bam \     # 输入BAM
    -b hg38-blacklist.v2.bed \ # 黑名单BED
    > sample_final.bam         # 最终BAM
samtools index sample_final.bam

2.3 Peak Calling¶

# ========== MACS2 Peak Calling ==========
macs2 callpeak \
    -t sample_final.bam \       # 输入BAM
    -f BAMPE \                   # 双端BAM格式
    -g hs \                      # 基因组大小：hs=人类
    -n sample \                  # 输出文件前缀
    --outdir peaks/ \            # 输出目录
    -q 0.05 \                    # FDR阈值
    --nomodel \                  # 不建模（ATAC-seq推荐）
    --shift -100 \               # 移位-100bp
    --extsize 200 \              # 延伸200bp
    --keep-dup all \             # 保留所有reads（已去重）
    --broad                      # 宽peak模式（可选）

# MACS2输出文件：
# peaks/sample_peaks.narrowPeak  ← 窄peak结果（默认）
# peaks/sample_peaks.broadPeak   ← 宽peak结果（用--broad时）
# peaks/sample_summits.bed       ← peak峰顶位置
# peaks/sample_treat_pileup.bdg  ← 信号堆叠文件

# 统计peak数量
wc -l peaks/sample_peaks.narrowPeak  # 查看peak数量
# 正常范围：20,000 - 150,000 个peaks

2.4 信号可视化¶

# ========== 生成BigWig信号文件 ==========
bamCoverage \
    -b sample_final.bam \
    -o sample.bw \              # 输出BigWig文件
    --binSize 10 \              # 10bp窗口
    --normalizeUsing RPKM \     # RPKM标准化
    --extendReads \             # 延伸reads
    -p 8                        # 8线程

# ========== TSS富集分析（质量评估） ==========
# 计算TSS附近的信号富集
computeMatrix reference-point \
    -S sample.bw \              # BigWig信号文件
    -R hg38_tss.bed \           # TSS位置BED文件
    --referencePoint TSS \      # 参考点：TSS
    -b 2000 \                   # TSS上游2000bp
    -a 2000 \                   # TSS下游2000bp
    --skipZeros \               # 跳过零值
    -o tss_matrix.gz            # 输出矩阵

plotProfile \
    -m tss_matrix.gz \          # 输入矩阵
    -o tss_enrichment.png \     # 输出图片
    --plotTitle "TSS Enrichment" \
    --perGroup                  # 按组画线

# TSS富集分数(TSS enrichment score)：
# >7 = 高质量数据
# 5-7 = 可接受
# <5 = 质量差

三、差异可及性分析¶

#!/usr/bin/env Rscript
# ATAC-seq差异可及性分析（使用DiffBind）

BiocManager::install("DiffBind")
library(DiffBind)  # 加载DiffBind

# ========== 创建样品表 ==========
samples <- data.frame(
    SampleID = c("C1", "C2", "C3", "T1", "T2", "T3"),
    Condition = c("Control", "Control", "Control", "Treatment", "Treatment", "Treatment"),
    bamReads = c("C1.bam", "C2.bam", "C3.bam", "T1.bam", "T2.bam", "T3.bam"),
    Peaks = c("C1_peaks.narrowPeak", "C2_peaks.narrowPeak", "C3_peaks.narrowPeak",
              "T1_peaks.narrowPeak", "T2_peaks.narrowPeak", "T3_peaks.narrowPeak"),
    PeakCaller = rep("narrow", 6)
)

# ========== DiffBind分析流程 ==========
dba <- dba(sampleSheet = samples)  # 创建DiffBind对象
dba <- dba.count(dba)  # 计算reads计数
dba <- dba.contrast(dba, categories = DBA_CONDITION)  # 设置对比
dba <- dba.analyze(dba, method = DBA_DESEQ2)  # DESeq2差异分析

# 提取结果
results <- dba.report(dba)  # 获取差异结果

# 查看差异peaks数量
cat("差异可及性区域:", length(results), "\n")
cat("更开放:", sum(results$Fold > 0), "\n")
cat("更关闭:", sum(results$Fold < 0), "\n")

# 可视化
dba.plotMA(dba)  # MA图
dba.plotVolcano(dba)  # 火山图
dba.plotHeatmap(dba, contrast = 1)  # 热图

四、转录因子足迹分析¶

# ========== HINT-ATAC足迹分析 ==========
# 安装RGT (Regulatory Genomics Toolbox)
pip install rgt  # 安装RGT

# 运行足迹分析
rgt-hint footprinting \
    --atac-seq \                # ATAC-seq模式
    --paired-end \              # 双端数据
    --organism hg38 \           # 基因组版本
    --output-location footprints/ \
    sample_final.bam \          # 输入BAM
    peaks/sample_peaks.narrowPeak  # Peak文件

# 差异足迹分析
rgt-hint differential \
    --organism hg38 \
    --bc \                      # 偏差校正
    --nc 8 \                    # 8线程
    --output-location diff_footprints/ \
    control.bam treatment.bam \
    control_peaks.bed treatment_peaks.bed

常见报错与解决¶

报错信息	原因	解决方法
`Very high mitochondrial rate`	线粒体reads>50%	正常现象，过滤后继续分析
`Too few peaks called`	数据量不够或参数不对	降低q值阈值，检查文库质量
`Fragment size wrong`	不是双端数据	用`-f BAM`代替`-f BAMPE`
`TSS enrichment <5`	数据质量差	检查建库质量，可能需要重做实验
`DiffBind count error`	BAM文件路径错误	使用绝对路径

速查表¶

========================================
ATAC-seq进阶分析 速查表
========================================

【完整流程】
1. 质控+Trim         → trim_galore
2. 比对              → bowtie2 --very-sensitive -X 2000
3. 过滤              → 去线粒体+去重+去黑名单
4. Peak calling      → macs2 callpeak -f BAMPE --nomodel
5. 信号可视化        → deeptools (bamCoverage, computeMatrix)
6. 差异分析          → DiffBind / DESeq2
7. 足迹分析          → HINT-ATAC / TOBIAS

【质量标准】
线粒体reads          → <50%（过滤后）
唯一比对率           → >70%
重复率               → <30%
Peak数               → 20,000-150,000
TSS enrichment       → >7为高质量
FRiP                 → >20%（Fraction of Reads in Peaks）

【MACS2参数（ATAC-seq专用）】
-f BAMPE             → 双端BAM格式
--nomodel            → 不建模（ATAC-seq必须）
--shift -100         → 移位-100bp
--extsize 200        → 延伸200bp
--keep-dup all       → 保留所有（已去重时用）

【片段大小含义】
<100bp               → 无核小体区域(NFR)
~200bp               → 单核小体
~400bp               → 双核小体

【面试考点】
Q: ATAC-seq和DNase-seq的区别？
A: ATAC-seq用Tn5转座酶，需要细胞少，操作简单；DNase-seq用DNase I酶

Q: 为什么ATAC-seq要去除线粒体reads？
A: 线粒体DNA无组蛋白包装，全部开放，Tn5大量插入造成干扰

Q: TSS enrichment score是什么？
A: TSS处信号vs背景的比值，反映数据质量，>7为高质量
========================================

参考资料：Buenrostro et al. Nature Methods 2013 | ENCODE ATAC-seq Pipeline | DiffBind | MACS2