组蛋白修饰ChIP-seq数据分析¶

组蛋白修饰 = 基因的"交通灯系统"

组蛋白 = DNA缠绕的蛋白质轮子（核小体核心）
修饰 = 在组蛋白"尾巴"上贴标签

常见修饰及含义：
──────────────────────────────────────
修饰          含义                Peak形态
──────────────────────────────────────
H3K4me3      活跃启动子            尖锐窄peak
H3K4me1      活跃/准备态增强子      中等宽peak
H3K27ac      活跃增强子/启动子      尖锐窄peak
H3K27me3     Polycomb抑制          宽域(broad domain)
H3K36me3     转录延伸（基因体）    宽域
H3K9me3      异染色质/沉默          宽域
H3K9ac       活跃启动子            尖锐窄peak
──────────────────────────────────────

染色质状态组合：
  H3K4me3 + H3K27ac = 活跃启动子
  H3K4me1 + H3K27ac = 活跃增强子
  H3K4me1 + H3K27me3 = "准备态"(poised)增强子
  H3K4me3 + H3K27me3 = 双价(bivalent)域（干细胞特有）

# ========== 环境设置 ==========
conda activate chipseq  # 激活ChIP-seq分析环境

# ========== 比对 ==========
# BWA-MEM比对（ChIP-seq常用BWA）
bwa mem \
    -t 8 \                      # 8线程
    ~/genomes/hg38/bwa_index \  # BWA索引
    sample_R1.fq.gz \           # R1
    sample_R2.fq.gz |           # R2
    samtools sort -@ 4 -o sample.bam  # 排序输出BAM

samtools index sample.bam  # 建索引

# ========== 过滤 ==========
# 1. 过滤低质量比对
samtools view -b -q 30 -F 1804 sample.bam > sample_filtered.bam  # MAPQ≥30，去未比对/重复等
# -F 1804: 去除 未比对+mate未比对+非主要比对+补充比对+重复

# 2. 去除重复
picard MarkDuplicates \
    I=sample_filtered.bam \
    O=sample_dedup.bam \
    M=dup_metrics.txt \
    REMOVE_DUPLICATES=true  # 删除重复

# 3. 去除黑名单区域
bedtools intersect -v -a sample_dedup.bam -b hg38-blacklist.v2.bed > sample_final.bam
samtools index sample_final.bam

# ========== 窄peak（H3K4me3, H3K27ac等） ==========
macs2 callpeak \
    -t sample_final.bam \       # 处理样品（ChIP）
    -c input_final.bam \        # 对照样品（Input）
    -f BAMPE \                  # 双端BAM
    -g hs \                     # 人类基因组
    -n sample_H3K4me3 \         # 输出前缀
    --outdir peaks/ \           # 输出目录
    -q 0.05 \                   # FDR阈值
    --call-summits              # 找peak峰顶

# ========== 宽peak（H3K27me3, H3K36me3等） ==========
macs2 callpeak \
    -t sample_final.bam \
    -c input_final.bam \
    -f BAMPE \
    -g hs \
    -n sample_H3K27me3 \
    --outdir peaks/ \
    --broad \                   # 宽peak模式！
    --broad-cutoff 0.1          # 宽peak的q值阈值

# 查看结果
wc -l peaks/sample_H3K4me3_peaks.narrowPeak  # 窄peak数量
wc -l peaks/sample_H3K27me3_peaks.broadPeak  # 宽peak数量

# ========== 生成BigWig ==========
# ChIP信号
bamCoverage -b sample_final.bam -o sample_chip.bw --binSize 10 --normalizeUsing RPKM -p 8

# Input信号
bamCoverage -b input_final.bam -o sample_input.bw --binSize 10 --normalizeUsing RPKM -p 8

# ChIP/Input比值（推荐）
bamCompare \
    -b1 sample_final.bam \      # ChIP BAM
    -b2 input_final.bam \       # Input BAM
    -o sample_log2ratio.bw \    # 输出log2比值
    --binSize 50 \              # 50bp窗口
    --normalizeUsing RPKM \     # RPKM标准化
    --scaleFactorsMethod None \
    --operation log2ratio \     # 计算log2(ChIP/Input)
    -p 8

# ========== TSS/基因体热图 ==========
computeMatrix scale-regions \
    -S sample_chip.bw sample_input.bw \  # BigWig文件
    -R genes.bed \              # 基因注释BED
    -b 3000 \                   # 基因上游3000bp
    -a 3000 \                   # 基因下游3000bp
    --regionBodyLength 5000 \   # 基因体缩放到5000bp
    --skipZeros \
    -o gene_matrix.gz \         # 输出矩阵
    -p 8

plotHeatmap \
    -m gene_matrix.gz \
    -o chipseq_heatmap.png \    # 输出热图
    --colorMap RdBu_r \         # 红蓝色板
    --sortUsing mean \          # 按平均信号排序
    --zMin 0 --zMax 5 \         # 颜色范围
    --plotTitle "H3K4me3 at Gene Body"

#!/usr/bin/env Rscript
# ChIP-seq差异修饰分析（DiffBind）

library(DiffBind)

# ========== 样品信息表 ==========
samples <- data.frame(
    SampleID = c("WT_1", "WT_2", "KO_1", "KO_2"),
    Condition = c("WT", "WT", "KO", "KO"),
    Factor = rep("H3K27ac", 4),  # 修饰类型
    bamReads = c("WT1_chip.bam", "WT2_chip.bam", "KO1_chip.bam", "KO2_chip.bam"),
    ControlID = c("WT_input", "WT_input", "KO_input", "KO_input"),
    bamControl = c("WT_input.bam", "WT_input.bam", "KO_input.bam", "KO_input.bam"),
    Peaks = c("WT1_peaks.narrowPeak", "WT2_peaks.narrowPeak",
              "KO1_peaks.narrowPeak", "KO2_peaks.narrowPeak"),
    PeakCaller = rep("narrow", 4)
)

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

# 获取差异结果
diff_peaks <- dba.report(dba_obj, th = 0.05)  # FDR<0.05的差异peaks
cat("差异修饰区域:", length(diff_peaks), "\n")

# 可视化
dba.plotMA(dba_obj)  # MA图
dba.plotVolcano(dba_obj)  # 火山图

# ========== ChromHMM：整合多个组蛋白修饰 ==========
# ChromHMM将基因组分为不同的"染色质状态"

# 安装ChromHMM
wget https://compbio.mit.edu/ChromHMM/ChromHMM.zip
unzip ChromHMM.zip

# 第一步：二值化数据
java -jar ChromHMM.jar BinarizeBam \
    ~/genomes/hg38/chromsizes.txt \  # 染色体大小文件
    bam_dir/ \                       # BAM文件目录
    cellmarkfiletable.txt \          # 样品-修饰对应表
    binarized_output/                # 输出目录

# 第二步：学习模型
java -jar ChromHMM.jar LearnModel \
    binarized_output/ \  # 二值化数据
    model_output/ \      # 模型输出目录
    15 \                 # 状态数（常用15或18）
    hg38                 # 基因组版本

# 输出：
# - 发射概率矩阵（每个状态对应哪些修饰组合）
# - 转移概率矩阵（状态之间的转换关系）
# - 基因组注释文件（每个区间的状态标注）

========================================
组蛋白修饰ChIP-seq 速查表
========================================

【修饰-功能对照表】
H3K4me3   → 活跃启动子（窄peak）
H3K27ac   → 活跃增强子/启动子（窄peak）
H3K4me1   → 增强子标记（中等peak）
H3K27me3  → Polycomb抑制（宽域）
H3K36me3  → 转录延伸（宽域）
H3K9me3   → 异染色质（宽域）

【MACS2参数选择】
窄peak修饰           → 默认模式（不加--broad）
宽域修饰             → --broad --broad-cutoff 0.1
必须有Input          → -c input.bam

【质量标准】
比对率               → >70%
重复率               → <30%（ChIP-seq正常偏高）
FRiP                 → >1%（最低），>5%（良好）
Peak数               → 因修饰类型而异

【下游分析工具】
DiffBind             → 差异修饰分析
ChromHMM             → 染色质状态注释
deepTools            → 信号可视化+热图
HOMER                → Motif分析+Peak注释

【面试考点】
Q: ChIP-seq需要对照吗？
A: 必须！Input或IgG对照，用于去除背景噪音

Q: H3K4me3和H3K27me3的peak为什么不同？
A: H3K4me3标记启动子（窄区域），H3K27me3覆盖大片段（宽域）

Q: 什么是双价域(bivalent domain)？
A: 同时有H3K4me3和H3K27me3，基因"准备好了但没开启"，常见于干细胞
========================================