全基因组亚硫酸氢盐bisulfite测序——区分5mC与5hmC¶

# 5mC/5hmC区分策略总结
# | 方法 | 读作C的修饰 | 读作T的修饰 |
# |------|-------------|-------------|
# | BS-seq | 5mC + 5hmC | C, 5fC, 5caC |
# | oxBS-seq | 5mC only | C, 5hmC(→5fC), 5fC, 5caC |
# | TAB-seq | 5hmC only | C, 5mC(→5caC), 5fC, 5caC |

# === 数据质控 ===

# 1. FastQC质检
fastqc -t 8 -o qc_raw/ BS_R1.fq.gz BS_R2.fq.gz oxBS_R1.fq.gz oxBS_R2.fq.gz

# 2. Trim Galore 去接头（bisulfite模式）
trim_galore --paired --clip_R1 10 --clip_R2 10 \
    --three_prime_clip_R1 10 --three_prime_clip_R2 10 \
    --cores 8 --fastqc \
    BS_R1.fq.gz BS_R2.fq.gz

trim_galore --paired --clip_R1 10 --clip_R2 10 \
    --three_prime_clip_R1 10 --three_prime_clip_R2 10 \
    --cores 8 --fastqc \
    oxBS_R1.fq.gz oxBS_R2.fq.gz

# 3. Bismark比对
bismark --genome genome_folder/ \
    --bowtie2 --parallel 4 \
    -1 BS_R1_val_1.fq.gz -2 BS_R2_val_2.fq.gz \
    -o bismark_BS/

bismark --genome genome_folder/ \
    --bowtie2 --parallel 4 \
    -1 oxBS_R1_val_1.fq.gz -2 oxBS_R2_val_2.fq.gz \
    -o bismark_oxBS/

# 4. 去重
deduplicate_bismark --bam bismark_BS/BS_pe.bam
deduplicate_bismark --bam bismark_oxBS/oxBS_pe.bam

# 5. 提取甲基化信息
bismark_methylation_extractor --paired-end --comprehensive \
    --CX --bedGraph --cytosine_report \
    --genome_folder genome_folder/ \
    bismark_BS/BS_pe.deduplicated.bam

bismark_methylation_extractor --paired-end --comprehensive \
    --CX --bedGraph --cytosine_report \
    --genome_folder genome_folder/ \
    bismark_oxBS/oxBS_pe.deduplicated.bam

# 6. 转换效率评估（lambda DNA spike-in）
# lambda DNA完全未甲基化，BS-seq后应>99%读作T
# 如果lambda C→T转换率 < 99%，说明BS转化不完全
grep "lambda" BS_CpG_report.txt | \
    awk '{m+=$4; u+=$5} END{print "BS conversion rate:", u/(m+u)*100, "%"}'

# === 5hmC计算 ===
library(methylKit)

# 读取BS-seq和oxBS-seq的CpG甲基化数据
bs_data <- methRead("BS_CpG_report.txt",
                     sample.id = "BS", assembly = "hg38",
                     treatment = 0, context = "CpG",
                     mincov = 10)  # 最小覆盖度10x

oxbs_data <- methRead("oxBS_CpG_report.txt",
                       sample.id = "oxBS", assembly = "hg38",
                       treatment = 0, context = "CpG",
                       mincov = 10)

# 提取共同位点
bs_df <- getData(bs_data)
oxbs_df <- getData(oxbs_data)

# 合并
merged <- merge(bs_df, oxbs_df, by = c("chr", "start", "end"),
                suffixes = c("_BS", "_oxBS"))

# 计算5hmC水平
# 5hmC% = BS_methylation% - oxBS_methylation%
merged$methyl_BS <- merged$numCs_BS / merged$coverage_BS
merged$methyl_oxBS <- merged$numCs_oxBS / merged$coverage_oxBS
merged$hydroxymethyl <- merged$methyl_BS - merged$methyl_oxBS

# 处理负值（统计噪声）
merged$hydroxymethyl[merged$hydroxymethyl < 0] <- 0

# 统计学方法：MLML（最大似然估计）
# 更精确地估计5mC和5hmC水平
library(MLML2R)

# MLML方法联合估计（注意参数名对应关系）
# BS-seq: T.matrix=未转化C(=5mC+5hmC), U.matrix=转化后T(=未甲基化C)
# oxBS-seq: M.matrix=未转化C(=5mC only), L.matrix=转化后T(=C+5hmC被氧化)
mlml_results <- MLML(T.matrix = merged$numCs_BS,   # BS中的C（未转化=5mC+5hmC）
                      U.matrix = merged$numTs_BS,   # BS中的T（转化=未甲基化）
                      M.matrix = merged$numCs_oxBS,  # oxBS中的C（=5mC only）
                      L.matrix = merged$numTs_oxBS)  # oxBS中的T（转化=C+5hmC→5fC）

merged$mC_MLML <- mlml_results$mC     # 5mC
merged$hmC_MLML <- mlml_results$hmC    # 5hmC
merged$C_MLML <- mlml_results$C        # 未修饰C

# === 5hmC基因组分布 ===
library(GenomicRanges)
library(annotatr)

# 定义5hmC位点（hydroxymethyl > 5%的位点）
hmC_sites <- GRanges(
  seqnames = merged$chr,
  ranges = IRanges(merged$start, merged$end)
)
hmC_sites$hmC_level <- merged$hmC_MLML

# 高5hmC位点
high_hmC <- hmC_sites[hmC_sites$hmC_level > 0.05]

# 注释基因组区域
annots <- build_annotations(genome = "hg38",
                             annotations = c("hg38_cpgs", "hg38_basicgenes",
                                              "hg38_genes_intergenic"))
annotated <- annotate_regions(regions = high_hmC, annotations = annots)

# 统计各区域的5hmC分布
region_summary <- summarize_annotations(annotated)
print(region_summary)

# 与已知增强子区域重叠
enhancers <- import("H3K27ac_enhancers.bed")
overlap_enhancer <- countOverlaps(high_hmC, enhancers) > 0
cat(sprintf("5hmC sites at enhancers: %.1f%%\n",
            sum(overlap_enhancer) / length(high_hmC) * 100))

# 5hmC与基因表达的关联
# 计算每个基因body的平均5hmC水平
gene_body_hmC <- aggregate(hmC_level ~ gene, data = annotated_df, FUN = mean)
# 与RNA-seq FPKM合并
merged_expr <- merge(gene_body_hmC, gene_expression, by = "gene")
cor.test(merged_expr$hmC_level, log2(merged_expr$FPKM + 1))
# 期望正相关：基因体5hmC与表达正相关

# === 差异5hmC分析 ===
library(DSS)

# 使用DSS进行差异甲基化/羟甲基化检测
# DSS支持对5hmC的差异分析

# 准备BSseq对象
bs_normal <- BSseq(chr = normal_data$chr,
                    pos = normal_data$start,
                    M = as.matrix(normal_data$numCs),
                    Cov = as.matrix(normal_data$coverage),
                    sampleNames = "Normal")

bs_tumor <- BSseq(chr = tumor_data$chr,
                   pos = tumor_data$start,
                   M = as.matrix(tumor_data$numCs),
                   Cov = as.matrix(tumor_data$coverage),
                   sampleNames = "Tumor")

# 差异5hmC区域（DMR）检测
dml_test <- DMLtest(BSobj = bs_combined, group1 = "Normal", group2 = "Tumor")
dmrs <- callDMR(dml_test, delta = 0.1, p.threshold = 0.001)

# 5hmC loss区域（在肿瘤中常见）
hmC_loss <- dmrs[dmrs$diff.Methy < -0.1, ]  # 肿瘤5hmC降低
cat(sprintf("Regions with 5hmC loss: %d\n", nrow(hmC_loss)))

# === 多组学整合 ===

# 5hmC与组蛋白修饰的共定位
# 5hmC与H3K4me1（增强子标记）高度共定位
# 5hmC与H3K27me3（抑制标记）互斥

# 在增强子区域比较5mC vs 5hmC
enhancer_5mC <- mean(merged$mC_MLML[overlap_enhancer])
enhancer_5hmC <- mean(merged$hmC_MLML[overlap_enhancer])
promoter_5mC <- mean(merged$mC_MLML[overlap_promoter])
promoter_5hmC <- mean(merged$hmC_MLML[overlap_promoter])

cat(sprintf("Enhancer: 5mC=%.3f, 5hmC=%.3f\n", enhancer_5mC, enhancer_5hmC))
cat(sprintf("Promoter: 5mC=%.3f, 5hmC=%.3f\n", promoter_5mC, promoter_5hmC))
# 期望：增强子5hmC高于启动子，启动子5mC高于增强子

# oxBS-seq分析流程
trim_galore --paired R1.fq.gz R2.fq.gz
bismark --genome ref/ -1 R1_val.fq.gz -2 R2_val.fq.gz
deduplicate_bismark --bam aligned.bam
bismark_methylation_extractor --cytosine_report --genome_folder ref/ dedup.bam
# R中计算5hmC
# 5hmC = BS_methylation - oxBS_methylation (MLML2R最大似然估计)