超级增强子鉴定ROSE¶

超级增强子 vs 普通增强子

普通增强子(TE)：
  ├── 大小：几百bp到几kb
  ├── 信号：中等H3K27ac信号
  └── 功能：调控附近基因

超级增强子(SE)：
  ├── 大小：几kb到几十kb
  ├── 信号：极强的H3K27ac/Mediator信号
  ├── 特征：由多个普通增强子聚集而成
  ├── 功能：调控细胞身份基因（master TF等）
  └── 临床意义：癌症中常被异常激活

ROSE算法原理：
  1. 合并12.5kb内的H3K27ac peaks（组成增强子簇）
  2. 去除TSS ± 2kb的区域（排除启动子）
  3. 对每个增强子簇计算H3K27ac信号总量
  4. 按信号从低到高排序
  5. 画"排序-信号"曲线
  6. 找曲线的拐点（斜率=1的切线点）
  7. 拐点以上 = 超级增强子，拐点以下 = 普通增强子

# ========== 安装ROSE ==========
# 克隆ROSE仓库
git clone https://github.com/stjude/ROSE.git  # 下载ROSE
cd ROSE

# 安装依赖
conda install -c bioconda samtools  # SAMtools
pip install numpy scipy  # Python依赖

# ROSE需要的文件：
# 1. H3K27ac BAM文件（ChIP样品）
# 2. Input BAM文件（对照）
# 3. MACS2调用的peak文件
# 4. 基因组注释文件

# ========== 准备输入文件 ==========
# 先用MACS2 call peaks
macs2 callpeak \
    -t H3K27ac_chip.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n H3K27ac \
    --nomodel --extsize 200

# ========== 运行ROSE ==========
python ROSE_main.py \
    -g hg38 \                    # 基因组版本
    -i H3K27ac_peaks.narrowPeak \ # MACS2 peak文件
    -r H3K27ac_chip.bam \         # ChIP BAM文件
    -c input.bam \                # Input BAM文件
    -o rose_output/ \             # 输出目录
    -s 12500 \                    # 合并距离（12.5kb，默认）
    -t 2000                       # TSS排除距离（±2kb）

# ROSE输出文件：
# rose_output/
# ├── *_AllEnhancers.table.txt    ← 所有增强子列表（含SE标记）
# ├── *_SuperEnhancers.table.txt  ← 只含超级增强子
# ├── *_Enhancers_withSuper.bed   ← BED格式SE文件
# └── *_Plot_points.png           ← 排序曲线（hockey stick plot）

#!/usr/bin/env python3
"""Python实现ROSE超级增强子鉴定算法"""

import pandas as pd  # 数据处理
import numpy as np  # 数值计算
import pysam  # BAM文件处理
import matplotlib.pyplot as plt  # 画图

# ========== 第一步：合并临近peaks ==========
def merge_peaks(peaks_df, merge_distance=12500):
    """合并distance内的peaks为增强子簇"""

    peaks_sorted = peaks_df.sort_values(["chr", "start"]).reset_index(drop=True)
    merged = []

    current = peaks_sorted.iloc[0].copy()
    for i in range(1, len(peaks_sorted)):
        row = peaks_sorted.iloc[i]

        # 如果在同一条染色体且距离<=merge_distance
        if row["chr"] == current["chr"] and row["start"] - current["end"] <= merge_distance:
            current["end"] = max(current["end"], row["end"])  # 合并：扩展结束位置
        else:
            merged.append(current.to_dict())  # 保存当前簇
            current = row.copy()  # 开始新簇

    merged.append(current.to_dict())  # 保存最后一个
    return pd.DataFrame(merged)

# ========== 第二步：计算信号 ==========
def count_reads_in_regions(bam_file, regions_df):
    """计算每个区域中的reads数"""

    bam = pysam.AlignmentFile(bam_file, "rb")
    counts = []

    for _, region in regions_df.iterrows():
        try:
            count = bam.count(region["chr"], region["start"], region["end"])
            counts.append(count)
        except:
            counts.append(0)

    bam.close()
    return np.array(counts)

# ========== 第三步：ROSE排序和拐点检测 ==========
def rose_algorithm(peaks_file, chip_bam, input_bam, tss_file=None,
                   merge_dist=12500, tss_exclude=2000):
    """
    完整ROSE算法流程
    """
    # 读取peaks
    peaks = pd.read_csv(peaks_file, sep="\t", header=None,
                        names=["chr", "start", "end", "name", "score",
                               "strand", "signalValue", "pValue", "qValue", "peak"])
    print(f"原始peaks数: {len(peaks)}")

    # 合并临近peaks
    merged = merge_peaks(peaks, merge_dist)
    print(f"合并后增强子簇数: {len(merged)}")

    # 排除TSS附近区域（可选）
    if tss_file:
        # 简化处理：实际需要用bedtools intersect
        print("已排除TSS附近区域")

    # 计算ChIP和Input的reads数
    chip_counts = count_reads_in_regions(chip_bam, merged)
    input_counts = count_reads_in_regions(input_bam, merged)

    # 计算标准化信号（ChIP - Input）
    # 简单标准化：按总reads数归一化
    chip_total = sum(chip_counts)
    input_total = sum(input_counts)

    normalized_signal = chip_counts / chip_total - input_counts / input_total
    normalized_signal = np.maximum(normalized_signal, 0)  # 去除负值

    merged["signal"] = normalized_signal

    # 按信号排序
    merged = merged.sort_values("signal").reset_index(drop=True)
    merged["rank"] = range(1, len(merged) + 1)

    # 找拐点（斜率=1的切线点）
    x = merged["rank"].values / max(merged["rank"])  # 归一化rank到0-1
    y = merged["signal"].values / max(merged["signal"])  # 归一化signal到0-1

    # 拐点：到对角线距离最大的点
    distances = np.abs(y - x)  # 到y=x线的距离
    inflection_idx = np.argmax(distances)  # 距离最大的点=拐点

    # 标记超级增强子
    merged["is_super_enhancer"] = merged["rank"] > inflection_idx

    n_se = merged["is_super_enhancer"].sum()
    n_te = (~merged["is_super_enhancer"]).sum()
    print(f"\n===== ROSE结果 =====")
    print(f"超级增强子(SE): {n_se}")
    print(f"普通增强子(TE): {n_te}")
    print(f"拐点位置: rank={inflection_idx}")

    # ========== 画Hockey Stick图 ==========
    plt.figure(figsize=(10, 8))

    te_mask = ~merged["is_super_enhancer"]
    se_mask = merged["is_super_enhancer"]

    plt.scatter(merged.loc[te_mask, "rank"], merged.loc[te_mask, "signal"],
                c="gray", s=5, alpha=0.5, label=f"Typical Enhancers ({n_te})")
    plt.scatter(merged.loc[se_mask, "rank"], merged.loc[se_mask, "signal"],
                c="red", s=15, alpha=0.8, label=f"Super-Enhancers ({n_se})")

    plt.axvline(inflection_idx, color="black", linestyle="--", alpha=0.5)
    plt.xlabel("Enhancer Rank", fontsize=14)
    plt.ylabel("H3K27ac Signal", fontsize=14)
    plt.title("Super-Enhancer Identification (ROSE)", fontsize=16)
    plt.legend(fontsize=12)
    plt.tight_layout()
    plt.savefig("hockey_stick_plot.png", dpi=300)
    plt.close()
    print("Hockey stick图已保存")

    # 保存结果
    se_results = merged[merged["is_super_enhancer"]]
    se_results.to_csv("super_enhancers.csv", index=False)
    merged.to_csv("all_enhancers.csv", index=False)

    return merged

# ========== 运行 ==========
# result = rose_algorithm(
#     "H3K27ac_peaks.narrowPeak",
#     "H3K27ac_chip.bam",
#     "input.bam"
# )

#!/usr/bin/env Rscript
# 超级增强子靶基因注释

library(GenomicRanges)
library(rtracklayer)

# 读取SE结果
se <- read.csv("super_enhancers.csv")
se_gr <- GRanges(seqnames = se$chr, ranges = IRanges(se$start, se$end))

# 读取基因注释
genes <- import("hg38_genes.bed", format = "BED")

# 找每个SE最近的基因
nearest_idx <- nearest(se_gr, genes)
se$nearest_gene <- genes$name[nearest_idx]
se$distance_to_gene <- distance(se_gr, genes[nearest_idx])

# 输出SE-基因对
cat("超级增强子及其最近基因：\n")
print(se[, c("chr", "start", "end", "signal", "nearest_gene", "distance_to_gene")])

# 保存
write.csv(se, "SE_with_genes.csv", row.names = FALSE)

========================================
超级增强子鉴定ROSE 速查表
========================================

【ROSE核心参数】
-g                   → 基因组版本（hg38/mm10）
-i                   → Peak文件（MACS2 narrowPeak）
-r                   → ChIP BAM文件
-c                   → Input BAM文件
-s                   → 合并距离（默认12500bp）
-t                   → TSS排除距离（默认2000bp）

【ROSE算法步骤】
1. 合并peaks          → 12.5kb内的peaks合并为簇
2. 去TSS              → 排除启动子附近区域
3. 计算信号           → ChIP - Input（标准化）
4. 排序               → 信号从低到高
5. 找拐点             → 斜率=1的切线点
6. 分类               → 拐点以上=SE，以下=TE

【超级增强子特征】
大小                 → 比TE大10-50倍
H3K27ac信号          → 比TE强5-20倍
Mediator             → 高度富集
BRD4                 → 高度富集
转录因子密度         → 很高
功能                 → 调控细胞身份基因

【常用SE数据库】
SEdb 2.0             → licpathway.net/sedb/
dbSUPER              → asntech.org/dbsuper/
SEanalysis           → SEanalysis.com

【临床应用】
BET抑制剂(JQ1等)     → 靶向SE上的BRD4蛋白
CDK7抑制剂           → 抑制SE依赖的转录

【面试考点】
Q: 超级增强子和普通增强子的区别？
A: SE由多个增强子聚集而成，信号更强，调控细胞身份基因

Q: ROSE算法怎么区分SE和TE？
A: 按H3K27ac信号排序，找排序曲线的拐点来划分

Q: 为什么SE是癌症治疗靶点？
A: 癌基因常受SE调控，BET抑制剂能选择性抑制SE活性
========================================