DNA 甲基化与基因表达关联¶

# ===== R: DNA 甲基化与基因表达关联分析 =====
library(methylKit)
library(DESeq2)
library(GenomicRanges)

# --- 1. 差异甲基化分析 ---
# 读取甲基化数据（Bismark 输出）
methyl_data <- methRead(
    list("ctrl1.cov", "ctrl2.cov", "treat1.cov", "treat2.cov"),
    sample.id = list("C1", "C2", "T1", "T2"),
    treatment = c(0, 0, 1, 1),
    assembly = "hg38", context = "CpG", mincov = 10)

# 合并 + 差异分析
meth_united <- unite(methyl_data)
diff_meth <- calculateDiffMeth(meth_united)

# 获取差异甲基化位点（差异>25%，q<0.01）
dmcs <- getMethylDiff(diff_meth, difference = 25, qvalue = 0.01)

# --- 2. 差异表达分析 ---
# RNA-seq 计数矩阵
rna_counts <- read.csv("gene_counts.csv", row.names = 1)
col_data <- data.frame(
    condition = c("ctrl", "ctrl", "treat", "treat"),
    row.names = colnames(rna_counts))

dds <- DESeqDataSetFromMatrix(rna_counts, col_data, ~condition)
dds <- DESeq(dds)
degs <- results(dds, contrast = c("condition", "treat", "ctrl"))
# 显著差异基因
sig_degs <- subset(degs, padj < 0.05 & abs(log2FoldChange) > 1)

# --- 3. 启动子甲基化与基因表达关联 ---
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

# 获取启动子区域（TSS ±2kb）
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
promoters_gr <- promoters(txdb, upstream = 2000, downstream = 200)

# 将差异甲基化位点注释到启动子
dmcs_gr <- GRanges(seqnames = dmcs$chr,
    ranges = IRanges(start = dmcs$start, end = dmcs$end))

# 找启动子内的差异甲基化位点
overlaps <- findOverlaps(dmcs_gr, promoters_gr)
promoter_dmcs <- dmcs[queryHits(overlaps), ]  # 位于启动子内的 DMC

# --- 4. 关联分析：甲基化 vs 表达 ---
# 期望：启动子高甲基化 → 基因低表达（负相关）

# 将 DMC 注释到最近基因
library(ChIPpeakAnno)
anno <- annotatePeakInBatch(dmcs_gr,
    AnnotationData = TSS.human.GRCh38,
    output = "nearestStart")

# 提取基因名和甲基化变化方向
anno_df <- as.data.frame(anno)

# 合并甲基化变化和表达变化
merged <- merge(
    data.frame(gene = anno_df$feature, meth_diff = dmcs$meth.diff),
    data.frame(gene = rownames(sig_degs),
               expr_fc = sig_degs$log2FoldChange),
    by = "gene")

# 相关性分析
cor_test <- cor.test(merged$meth_diff, merged$expr_fc,
                     method = "spearman")
cat("甲基化变化 vs 表达变化相关系数：",
    round(cor_test$estimate, 3), "\n")
cat("p 值：", cor_test$p.value, "\n")
# 启动子甲基化通常与表达负相关（r < 0）

# 可视化
library(ggplot2)
ggplot(merged, aes(x = meth_diff, y = expr_fc)) +
    geom_point(alpha = 0.5) +             # 散点
    geom_smooth(method = "lm", color = "red") +  # 拟合线
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_vline(xintercept = 0, linetype = "dashed") +
    labs(x = "启动子甲基化变化（%）",
         y = "基因表达变化（log2FC）",
         title = "甲基化-表达负相关") +
    theme_minimal()