表观基因组数据整合方法¶

# ===== ChromHMM 染色质状态注释 =====

# 1. 准备输入数据：将每种组蛋白修饰的 BAM 二值化
# ChromHMM 需要知道每个 200bp bin 是否有信号
java -jar ChromHMM.jar BinarizeBam \
    hg38_chrom_sizes.txt \
    bam_directory/ \
    cell_mark_table.txt \
    binarized_output/

# cell_mark_table.txt 格式（制表符分隔）：
# 细胞名   标记名 BAM文件   对照BAM
# K562  H3K4me3 K562_H3K4me3.bam    K562_Input.bam
# K562  H3K27ac K562_H3K27ac.bam    K562_Input.bam
# K562  H3K27me3    K562_H3K27me3.bam   K562_Input.bam
# K562  H3K36me3    K562_H3K36me3.bam   K562_Input.bam

# 2. 训练模型：学习 15 种染色质状态
java -jar ChromHMM.jar LearnModel \
    binarized_output/ \
    chromhmm_output/ \
    15 \       # 状态数
    hg38       # 基因组版本

# 输出文件：
# emissions_15.png → 每种状态的组蛋白标记模式（热图）
# transitions_15.png → 状态间转换概率
# K562_15_segments.bed → 基因组分割结果
# K562_15_overlap.png → 各状态与基因组注释的重叠

# 3. 查看典型状态含义
# State 1: H3K4me3 高 → 活跃启动子
# State 2: H3K27ac 高 → 活跃增强子
# State 3: H3K36me3 高 → 转录区域
# State 4: H3K27me3 高 → Polycomb 沉默
# State 5: 全低 → 静默区域

# ===== R: 多组学整合分析 =====
library(MOFA2)

# 1. 准备多组学数据矩阵
# 每种组学一个矩阵，行=特征，列=样本
methylation <- readRDS("meth_matrix.rds")   # CpG 甲基化
accessibility <- readRDS("atac_matrix.rds") # ATAC-seq peaks
expression <- readRDS("rna_matrix.rds")     # 基因表达

# 2. 创建 MOFA 对象
mofa_data <- create_mofa(list(
    methylation = methylation,
    accessibility = accessibility,
    expression = expression
))

# 3. 设置模型参数
data_opts <- get_default_data_options(mofa_data)
model_opts <- get_default_model_options(mofa_data)
model_opts$num_factors <- 10  # 提取 10 个因子

train_opts <- get_default_training_options(mofa_data)
train_opts$convergence_mode <- "slow"  # 严格收敛
train_opts$seed <- 42

# 4. 训练模型
mofa_data <- prepare_mofa(mofa_data,
    data_options = data_opts,
    model_options = model_opts,
    training_options = train_opts)

mofa_model <- run_mofa(mofa_data)

# 5. 结果解读
# 查看每个因子解释了各组学多少变异
plot_variance_explained(mofa_model, plot_total = TRUE)

# 查看因子 1 在各组学中的权重
plot_weights(mofa_model, view = "expression",
             factor = 1, nfeatures = 20)

# 因子与临床表型的关联
plot_factor(mofa_model, factor = 1,
            color_by = "disease_status")

# 6. 简单整合：检查各层信号的一致性
# 在差异基因的启动子区域，甲基化是否下降，可及性是否上升？
# 这是最直观的整合验证
library(GenomicRanges)

# 差异表达上调的基因
up_genes <- rownames(subset(degs, log2FC > 1 & padj < 0.05))

# 检查这些基因启动子的甲基化变化
# 期望：表达上调的基因 → 启动子甲基化下降
promoter_meth_change <- sapply(up_genes, function(g) {
    # 计算该基因启动子区域的平均甲基化变化
    mean(diff_meth[promoter_of(g), ])
})

# 检查 ATAC-seq 信号变化
# 期望：表达上调的基因 → 启动子可及性增加
promoter_atac_change <- sapply(up_genes, function(g) {
    mean(diff_atac[promoter_of(g), ])
})

cat("上调基因启动子甲基化变化中位数：",
    median(promoter_meth_change, na.rm = TRUE), "\n")
cat("上调基因启动子可及性变化中位数：",
    median(promoter_atac_change, na.rm = TRUE), "\n")