33. 表观基因组学入门（ChIP-seq / ATAC-seq）¶

# ========== 第一步：构建 Bisulfite 基因组索引 ==========
# Bismark 会把参考基因组做 C→T 和 G→A 两种转换，分别建索引
bismark_genome_preparation \
    --bowtie2 \                  # 使用 Bowtie2 作为比对器（也可用 HISAT2）
    /path/to/genome/             # 参考基因组所在目录（需含 .fa 文件）

# ========== 第二步：比对 ==========
# 把 bisulfite 处理后的 reads 比对到转换后的基因组
bismark \
    --bowtie2 \                  # 比对引擎
    --genome /path/to/genome/ \  # 基因组索引目录
    -1 sample_R1.fq.gz \         # 双端测序 Read1
    -2 sample_R2.fq.gz \         # 双端测序 Read2
    -o bismark_output/           # 输出目录

# ========== 第三步：去重 ==========
# 去除 PCR 重复（bisulfite-seq 的 PCR 偏差比普通 WGS 更严重）
deduplicate_bismark \
    --bam \                      # 输入为 BAM 格式
    --paired \                   # 双端数据
    bismark_output/sample.bam    # 输入文件

# ========== 第四步：甲基化信息提取 ==========
# 从比对结果中提取每个 CpG/CHG/CHH 位点的甲基化状态
bismark_methylation_extractor \
    --paired-end \               # 双端模式
    --comprehensive \            # 输出所有 context（CpG/CHG/CHH）
    --bedGraph \                 # 生成 bedGraph 格式（可用于可视化）
    --cytosine_report \          # 生成全基因组 CpG 报告
    --genome_folder /path/to/genome/ \  # 基因组目录
    -o methylation_output/ \     # 输出目录
    sample.deduplicated.bam      # 去重后的 BAM

# ========== 第五步：生成报告 ==========
bismark2report                   # 自动生成 HTML 质控报告
bismark2summary                  # 多样本汇总报告

固定（交联蛋白-DNA）→ 超声打断 DNA → 用特异性抗体免疫沉淀
→ 回收 DNA → 建库测序 → 生物信息学分析

# ================================================================
# ChIP-seq 分析流水线（以 H3K27ac 为例）
# 需要：处理组 BAM（ChIP）+ 对照组 BAM（Input）
# ================================================================

# ========== 第一步：质控 ==========
fastqc \
    -t 8 \                       # 8 线程
    -o fastqc_results/ \         # 输出目录
    chip_sample.fq.gz \          # ChIP 样本
    input_control.fq.gz          # Input 对照

# 如需，用 Trim Galore 去接头
trim_galore \
    --paired \                   # 双端模式
    --fastqc \                   # 修剪后自动跑 FastQC
    -o trimmed/ \                # 输出目录
    chip_R1.fq.gz chip_R2.fq.gz  # 输入文件

# ========== 第二步：比对（Alignment）==========
# 用 Bowtie2 将 reads 比对到参考基因组
bowtie2 \
    -p 8 \                       # 8 线程
    -x /path/to/bowtie2_index \  # 基因组索引前缀
    -1 chip_R1_trimmed.fq.gz \   # Read1
    -2 chip_R2_trimmed.fq.gz \   # Read2
    --very-sensitive \           # 高灵敏度模式
    --no-mixed \                 # 不允许单端比对
    --no-discordant \            # 不允许不一致的配对
    -S chip.sam                  # 输出 SAM 文件

# SAM 转 BAM，排序并建索引
samtools view -@ 8 -bS -q 10 chip.sam | \   # -q 10：过滤低质量比对
    samtools sort -@ 8 -o chip_sorted.bam    # 按坐标排序
samtools index chip_sorted.bam               # 建立索引（.bai 文件）

# 对 Input 对照做同样的比对操作（此处省略，命令相同）

# ========== 第三步：去重（Remove Duplicates）==========
# PCR 重复会造成假阳性，必须去除
# 用 Picard MarkDuplicates（也可用 sambamba）
picard MarkDuplicates \
    INPUT=chip_sorted.bam \           # 输入 BAM
    OUTPUT=chip_dedup.bam \           # 去重后 BAM
    METRICS_FILE=chip_dup_metrics.txt \  # 重复率统计
    REMOVE_DUPLICATES=true \          # 直接移除重复 reads
    VALIDATION_STRINGENCY=LENIENT     # 宽松验证

samtools index chip_dedup.bam         # 重新建索引

# ========== 第四步：去除黑名单区域（可选但推荐）==========
# ENCODE 黑名单区域是已知的假信号来源
bedtools intersect \
    -v \                              # 保留不在黑名单中的 reads
    -abam chip_dedup.bam \            # 输入 BAM
    -b hg38-blacklist.v2.bed \        # ENCODE 黑名单 BED 文件
    > chip_clean.bam                  # 清理后的 BAM
samtools index chip_clean.bam

# ========== 第五步：Peak Calling（MACS2）==========
# MACS2 是最主流的 peak caller
# 核心思想：找到 ChIP 信号显著高于 Input 背景的区域

# --- 窄峰（narrow peak）模式 ---
# 适用于：转录因子、H3K4me3、H3K27ac
macs2 callpeak \
    -t chip_clean.bam \               # 处理组（ChIP）
    -c input_clean.bam \              # 对照组（Input）
    -f BAMPE \                        # 双端 BAM 格式（单端用 BAM）
    -g hs \                           # 有效基因组大小（hs=人类，mm=小鼠）
    -n H3K27ac_sample \               # 输出文件名前缀
    --outdir macs2_output/ \          # 输出目录
    -q 0.05 \                         # FDR 阈值（q-value < 0.05）
    --keep-dup all                    # 保留所有 reads（已手动去重）

# --- 宽峰（broad peak）模式 ---
# 适用于：H3K27me3、H3K36me3、H3K9me3 等覆盖大区域的修饰
macs2 callpeak \
    -t chip_clean.bam \
    -c input_clean.bam \
    -f BAMPE \
    -g hs \
    -n H3K27me3_sample \
    --outdir macs2_broad/ \
    --broad \                         # 开启宽峰模式
    --broad-cutoff 0.1 \              # 宽峰 q-value 阈值
    --keep-dup all

# MACS2 输出文件说明：
# _peaks.narrowPeak / _peaks.broadPeak  → Peak 坐标（BED 格式）
# _peaks.xls                            → Peak 详细统计
# _summits.bed                          → Peak 顶点坐标（仅 narrow 模式）
# _model.r                              → 片段大小模型（可用 R 绘图）

# ========== 第六步：差异分析（Differential Binding）==========
# 比较两个条件下 peak 的变化，常用 R 包 DiffBind
# 在 R 中运行：
# library(DiffBind)
# samples <- read.csv("sample_sheet.csv")  # 样本信息表
# dba_obj <- dba(sampleSheet = samples)    # 创建 DBA 对象
# dba_obj <- dba.count(dba_obj)            # 计数
# dba_obj <- dba.contrast(dba_obj)         # 设置对比
# dba_obj <- dba.analyze(dba_obj)          # 差异分析（默认用 DESeq2）
# dba.report(dba_obj)                      # 输出差异 peak

# ========== 第七步：Peak 注释 ==========
# 用 HOMER 的 annotatePeaks.pl 或 R 包 ChIPseeker
# HOMER 方式：
annotatePeaks.pl \
    macs2_output/H3K27ac_sample_peaks.narrowPeak \  # Peak 文件
    hg38 \                            # 基因组版本
    -go go_output/ \                  # GO 富集分析输出
    -annStats ann_stats.txt \         # 注释统计
    > annotated_peaks.txt             # 注释结果

# ========== 第八步：Motif 分析 ==========
# 在 peak 区域寻找富集的转录因子结合基序
# HOMER findMotifsGenome.pl
findMotifsGenome.pl \
    macs2_output/H3K27ac_sample_summits.bed \  # Peak 顶点（±100bp 效果更好）
    hg38 \                            # 基因组版本
    motif_output/ \                   # 输出目录
    -size 200 \                       # 取 summit 两侧各 100bp
    -mask \                           # 屏蔽重复序列
    -p 8                              # 8 线程

# 也可以用 MEME-ChIP 做 de novo motif 发现：
# 先提取 peak 序列
bedtools getfasta \
    -fi hg38.fa \                     # 参考基因组
    -bed peaks_top500.bed \           # 取前 500 个最显著 peak
    -fo peaks_top500.fa               # 输出 FASTA
# 然后上传到 MEME-ChIP 在线服务器分析

# ================================================================
# ATAC-seq 分析流水线
# 注意：与 ChIP-seq 的主要区别用 [ATAC特有] 标记
# ================================================================

# ========== 第一步：质控 + 去接头 ==========
# 与 ChIP-seq 相同
trim_galore \
    --paired \
    --fastqc \
    --nextera \                      # [ATAC特有] ATAC-seq 用 Nextera 接头
    -o trimmed/ \
    atac_R1.fq.gz atac_R2.fq.gz

# ========== 第二步：比对 ==========
# 与 ChIP-seq 基本相同
bowtie2 \
    -p 8 \
    -x /path/to/bowtie2_index \
    -1 atac_R1_trimmed.fq.gz \
    -2 atac_R2_trimmed.fq.gz \
    --very-sensitive \
    -X 2000 \                        # [ATAC特有] 最大插入片段设为 2000bp
    --no-mixed \
    --no-discordant \
    -S atac.sam

# ========== 第三步：过滤 ==========
# 3a. 去除线粒体 reads [ATAC特有]
# ATAC-seq 中线粒体 DNA 占比很高（20%-80%），必须去除
samtools view -@ 8 -bS -q 10 atac.sam | \
    samtools sort -@ 8 -o atac_sorted.bam
samtools index atac_sorted.bam

# 去除 chrM（线粒体）上的 reads
samtools idxstats atac_sorted.bam | \
    cut -f 1 | grep -v chrM | \
    xargs samtools view -b atac_sorted.bam > atac_noMT.bam
samtools index atac_noMT.bam

# 3b. 去重
picard MarkDuplicates \
    INPUT=atac_noMT.bam \
    OUTPUT=atac_dedup.bam \
    METRICS_FILE=atac_dup_metrics.txt \
    REMOVE_DUPLICATES=true
samtools index atac_dedup.bam

# 3c. 去除黑名单区域
bedtools intersect -v \
    -abam atac_dedup.bam \
    -b hg38-blacklist.v2.bed \
    > atac_clean.bam
samtools index atac_clean.bam

# ========== 第四步：Tn5 偏移校正 [ATAC特有] ==========
# Tn5 转座酶插入时会造成 9bp 的偏移：
# 正链 reads 向 +4bp 方向移动，负链 reads 向 -5bp 方向移动
# 这在某些下游分析（如足迹分析）中必须校正
# 方法：用 alignmentSieve（deepTools）或手动用 awk 调整
alignmentSieve \
    -b atac_clean.bam \
    --ATACshift \                    # 自动校正 Tn5 偏移
    -o atac_shifted.bam
samtools sort -o atac_final.bam atac_shifted.bam
samtools index atac_final.bam

# ========== 第五步：Peak Calling（MACS2）==========
# ATAC-seq 不需要 Input 对照
macs2 callpeak \
    -t atac_final.bam \              # 只有处理组
    -f BAMPE \                       # 双端 BAM
    -g hs \                          # 有效基因组大小
    -n ATAC_sample \                 # 输出前缀
    --outdir macs2_atac/ \
    --nomodel \                      # [ATAC特有] 不建模（ATAC 的片段大小分布与 ChIP 不同）
    --shift -75 \                    # [ATAC特有] 配合 --extsize 使用
    --extsize 150 \                  # [ATAC特有] 以切割位点为中心，取 150bp 区域
    --keep-dup all \
    -q 0.05                          # FDR 阈值

# 注意：也有文献推荐 ATAC-seq 直接使用 --nomodel --shift -100 --extsize 200
# 或者用 BAMPE 模式时不设 shift/extsize（MACS2 自动用实际片段大小）

# ========== 第六步：下游分析 ==========
# 与 ChIP-seq 类似：Peak 注释、差异分析、Motif 分析
# 额外可做：
# - 核小体位置分析（按片段大小分类：<150bp 为无核小体区域，150-300bp 为单核小体）
# - TF 足迹分析（Footprinting）：用 TOBIAS 或 HINT-ATAC

# 安装
pip install macs2        # 推荐 Python 3.9+
# 或
conda install -c bioconda macs2

# 核心功能
macs2 callpeak          # Peak calling（最常用）
macs2 bdgcmp            # 比较两个 bedGraph 信号
macs2 bdgdiff           # 差异 peak 检测
macs2 filterdup         # 过滤重复 reads
macs2 predictd          # 预测片段大小（建模）

# 安装
conda install -c bioconda deeptools

# 常用命令
# (1) BAM 转 bigWig（用于 IGV 可视化和热图）
bamCoverage \
    -b chip_clean.bam \              # 输入 BAM
    -o chip.bw \                     # 输出 bigWig
    --binSize 10 \                   # 分辨率（bp）
    --normalizeUsing RPKM \          # 标准化方法（也可用 CPM/BPM/RPGC）
    --effectiveGenomeSize 2700000000 \  # 有效基因组大小（人类）
    -p 8                             # 线程数

# (2) 计算基因/Peak 周围信号矩阵
computeMatrix reference-point \
    -S chip.bw input.bw \            # 输入 bigWig 信号文件
    -R genes.bed \                   # 参考区域（基因/peak 坐标）
    --referencePoint TSS \           # 以 TSS 为参考点（也可用 TES/center）
    -a 3000 -b 3000 \               # TSS 上下游各 3000bp
    -o matrix.gz                     # 输出矩阵

# (3) 绘制热图
plotHeatmap \
    -m matrix.gz \                   # 输入矩阵
    -o heatmap.png \                 # 输出图片
    --colorMap RdBu_r \              # 配色方案
    --whatToShow 'heatmap and colorbar' \
    --zMin 0 --zMax 10               # 信号范围

# (4) TSS Profile 图
plotProfile \
    -m matrix.gz \
    -o profile.png \
    --perGroup                       # 按组绘制

# (5) 样本相关性检查（QC）
multiBamSummary bins \
    --bamfiles sample1.bam sample2.bam sample3.bam \
    -o results.npz \
    -p 8
plotCorrelation \
    -in results.npz \
    --corMethod spearman \           # Spearman 相关系数
    --whatToPlot heatmap \
    -o correlation_heatmap.png

# (6) 指纹图（Fingerprint plot，评估 ChIP 富集质量）
plotFingerprint \
    -b chip.bam input.bam \
    -o fingerprint.png \
    --labels ChIP Input

# 安装
conda install -c bioconda homer

# 核心功能
# (1) 创建标签目录（HOMER 的核心数据结构）
makeTagDirectory tag_dir/ chip_clean.bam

# (2) Peak Calling（HOMER 自带的 peak caller）
findPeaks tag_dir/ \
    -style factor \                  # factor=窄峰, histone=宽峰
    -i input_tag_dir/ \              # Input 对照
    -o peaks.txt                     # 输出

# (3) Peak 注释
annotatePeaks.pl peaks.txt hg38 > annotated.txt

# (4) Motif 发现（HOMER 最强功能）
findMotifsGenome.pl peaks.txt hg38 motif_results/ -size 200 -mask

# (5) 差异分析
getDifferentialPeaks peaks1.txt tag_dir1/ tag_dir2/ > diff_peaks.txt

ChIP-seq 流程：
FASTQ → FastQC → Trim → Bowtie2比对 → SAM→BAM排序 → 去重(Picard)
→ 去黑名单 → MACS2 Peak Calling(需Input) → Peak注释 → Motif分析 → 差异分析

ATAC-seq 流程：
FASTQ → FastQC → Trim(Nextera接头) → Bowtie2比对 → SAM→BAM排序
→ 去线粒体reads → 去重(Picard) → 去黑名单 → Tn5偏移校正
→ MACS2 Peak Calling(无Input, --nomodel) → Peak注释 → Motif分析 → 差异分析

Bisulfite-seq 流程：
FASTQ → FastQC → Trim → Bismark比对(转换基因组) → 去重
→ 甲基化提取 → 差异甲基化分析(methylKit/DSS)