组蛋白变体分析 H2A.Z / H3.3¶

# ===== 组蛋白变体 ChIP-seq / CUT&Tag 分析 =====

# 1. 比对（以 H2A.Z CUT&Tag 为例）
bowtie2 --very-sensitive --no-mixed --no-discordant \
    --maxins 700 -x hg38_index \
    -1 H2AZ_R1.fastq.gz -2 H2AZ_R2.fastq.gz \
    -S H2AZ.sam -p 8

# 2. 过滤、排序、去重
samtools view -bS -q 10 H2AZ.sam | \
    samtools sort -@ 8 -o H2AZ.sorted.bam
samtools markdup -r H2AZ.sorted.bam H2AZ.dedup.bam
samtools index H2AZ.dedup.bam

# 3. 生成 bigWig（用于可视化）
bamCoverage -b H2AZ.dedup.bam -o H2AZ.bw \
    --normalizeUsing RPKM --binSize 10 -p 8

# 4. Peak calling
# H2A.Z 在启动子形成尖峰，用 narrowPeak 模式
macs2 callpeak -t H2AZ.dedup.bam -c IgG.dedup.bam \
    -f BAMPE -g hs -n H2AZ --nomodel --keep-dup all

# 5. TSS 处的 H2A.Z 分布（deepTools）
# 计算 TSS ±3kb 范围内的信号分布
computeMatrix reference-point -S H2AZ.bw \
    -R hg38_tss.bed \
    --referencePoint TSS \
    -a 3000 -b 3000 \
    -o H2AZ_tss_matrix.gz

# 画 TSS 处的信号热图
plotHeatmap -m H2AZ_tss_matrix.gz \
    -o H2AZ_tss_heatmap.pdf \
    --colorMap RdYlBu_r \
    --sortUsing mean

# ===== R: H2A.Z 与转录活性的关联分析 =====
library(GenomicRanges)
library(ChIPpeakAnno)

# 读取 H2A.Z peaks
h2az_peaks <- read.table("H2AZ_peaks.narrowPeak",
    col.names = c("chr","start","end","name","score",
                  "strand","signal","pval","qval","summit"))

# 读取 RNA-seq 表达数据
expr <- read.csv("gene_expression.csv", row.names = 1)
# 按表达量分组
expr$group <- cut(expr$FPKM,
    breaks = c(0, 1, 10, 100, Inf),
    labels = c("silent", "low", "medium", "high"))

# 统计不同表达水平基因启动子处的 H2A.Z 丰度
# 读取 TSS ±1kb 区域的 H2A.Z 信号
library(rtracklayer)
h2az_bw <- import("H2AZ.bw")

# 计算每个基因 TSS ±1kb 的 H2A.Z 平均信号
tss_gr <- GRanges(seqnames = expr$chr,
    ranges = IRanges(start = expr$tss - 1000,
                     end = expr$tss + 1000))

h2az_signal <- sapply(seq_along(tss_gr), function(i) {
    ov <- subsetByOverlaps(h2az_bw, tss_gr[i])
    if (length(ov) > 0) mean(ov$score) else 0
})

expr$h2az_signal <- h2az_signal

# 箱线图：不同表达水平的 H2A.Z 信号
library(ggplot2)
ggplot(expr, aes(x = group, y = h2az_signal, fill = group)) +
    geom_boxplot() +
    labs(x = "基因表达水平", y = "TSS 处 H2A.Z 信号",
         title = "H2A.Z 在活跃基因启动子处富集") +
    theme_minimal()