超级增强子鉴定与分析¶

# ===== 超级增强子鉴定流程 =====

# 1. H3K27ac ChIP-seq 数据比对（假设已完成）
# 2. Peak calling
macs2 callpeak -t H3K27ac.bam -c Input.bam \
    -f BAM -g hs -n H3K27ac \
    --nomodel --extsize 200 -p 1e-9

# 3. 使用 ROSE 鉴定超级增强子
# 先安装 ROSE（https://github.com/stjude/ROSE）
# ROSE 需要的输入：peaks 文件 + BAM 文件 + 基因组

# 把 MACS2 的 narrowPeak 转成 ROSE 需要的 GFF 格式
awk 'OFS="\t" {print $4, $1, "", $2, $3, "", $6, "", ""}' \
    H3K27ac_peaks.narrowPeak > H3K27ac_peaks.gff

# 运行 ROSE
python ROSE_main.py \
    -g hg38 \
    -i H3K27ac_peaks.gff \
    -r H3K27ac.bam \
    -c Input.bam \
    -o rose_output/ \
    -s 12500 \  # 12.5kb 以内的 peaks 合并
    -t 2500     # 排除 TSS ±2.5kb

# ROSE 输出文件说明：
# *_SuperEnhancers.table.txt → 超级增强子列表
# *_Enhancers_withSuper.bed → 所有增强子（标注了是否为 SE）
# *_Plot_points.png → hockey-stick 曲线图

# ===== R: 超级增强子下游分析 =====
library(GenomicRanges)
library(ggplot2)

# 读取 ROSE 输出的超级增强子
se_table <- read.delim("rose_output/SuperEnhancers.table.txt")
cat("鉴定到的超级增强子数：", nrow(se_table), "\n")
# 通常一个细胞类型有 200-600 个 SE

# 读取所有增强子数据（含信号强度排序）
all_enh <- read.delim("rose_output/AllEnhancers.table.txt")
all_enh <- all_enh[order(all_enh$enhancerRank), ]

# 画 hockey-stick 曲线
all_enh$rank <- seq_len(nrow(all_enh))
all_enh$is_super <- ifelse(all_enh$isSuper == 1, "SE", "TE")

ggplot(all_enh, aes(x = rank, y = signal, color = is_super)) +
    geom_point(size = 0.5) +
    scale_color_manual(values = c("SE" = "red", "TE" = "gray")) +
    labs(x = "增强子排名", y = "H3K27ac 信号",
         title = "Hockey-stick 曲线：超级增强子鉴定") +
    theme_minimal()

# SE 关联的基因分析
# 将 SE 注释到最近基因
se_gr <- GRanges(seqnames = se_table$CHROM,
    ranges = IRanges(start = se_table$START, end = se_table$STOP))

# 常见 SE 驱动的癌基因
oncogenes <- c("MYC", "BCL2", "CCND1", "MDM2", "NOTCH1")
se_genes <- se_table$CLOSEST_GENE

# 检查哪些癌基因被 SE 驱动
hit <- intersect(se_genes, oncogenes)
if (length(hit) > 0) {
    cat("被超级增强子驱动的癌基因：", paste(hit, collapse = ", "), "\n")
}

# GO 富集分析
library(clusterProfiler)
library(org.Hs.eg.db)

# SE 关联基因的功能富集
gene_ids <- bitr(se_genes, fromType = "SYMBOL",
                 toType = "ENTREZID", OrgDb = org.Hs.eg.db)
go_results <- enrichGO(gene = gene_ids$ENTREZID,
    OrgDb = org.Hs.eg.db,
    ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05)