RRBS 简化代表性亚硫酸盐测序¶

# ===== RRBS 数据分析流程 =====

# 1. 质控
fastqc -t 4 sample_R1.fastq.gz sample_R2.fastq.gz

# 2. 去接头 + RRBS 专用模式
# --rrbs: 自动处理 MspI 切割位点引入的填充碱基
# 这一步会去掉端部人工引入的碱基，避免甲基化水平偏差
trim_galore --paired --rrbs -o trimmed/ \
    sample_R1.fastq.gz sample_R2.fastq.gz

# 3. 准备 Bismark 参考基因组索引（只需做一次）
bismark_genome_preparation --bowtie2 /ref/hg38/

# 4. 比对：Bismark 处理亚硫酸盐转化后的序列
# Bismark 会自动考虑 C→T 转化，三种情况对应三种比对
bismark --genome /ref/hg38/ --bowtie2 -p 4 \
    -1 trimmed/sample_R1_val_1.fq.gz \
    -2 trimmed/sample_R2_val_2.fq.gz \
    -o aligned/

# 5. 去重（RRBS 不建议用常规 PCR 去重）
# 因为 MspI 切割位点固定，同一位置的 reads 不一定是 PCR 重复
# 但如果 PCR 扩增过多，仍需去重
deduplicate_bismark --paired aligned/sample_pe.bam

# 6. 提取甲基化信息
# --comprehensive: 合并 CpG/CHG/CHH 上下文
# --cytosine_report: 生成全基因组胞嘧啶报告
bismark_methylation_extractor --paired-end --comprehensive \
    --cytosine_report --genome_folder /ref/hg38/ \
    -o methylation/ aligned/sample_pe.deduplicated.bam

# 7. 用 methylKit 做差异甲基化分析（R 代码）
# 以下在 R 中运行

# ===== R: methylKit 差异甲基化分析 =====
library(methylKit)

# 读取 Bismark 的 cytosine report 文件
# 每个文件对应一个样本
file_list <- list(
    "sample1_CpG_report.txt",  # 对照组样本1
    "sample2_CpG_report.txt",  # 对照组样本2
    "sample3_CpG_report.txt",  # 处理组样本1
    "sample4_CpG_report.txt"   # 处理组样本2
)

# 读取数据，设定样本 ID、处理组信息、最低覆盖度
myobj <- methRead(file_list,
    sample.id = list("ctrl1", "ctrl2", "treat1", "treat2"),
    treatment = c(0, 0, 1, 1),  # 0=对照，1=处理
    assembly = "hg38",
    context = "CpG",
    mincov = 10)  # 至少 10x 覆盖度

# 过滤极高覆盖（可能是 PCR 偏差）
filtered <- filterByCoverage(myobj,
    lo.count = 10,   # 最低 10 reads
    hi.perc = 99.9)  # 去掉 top 0.1%

# 合并所有样本的共同 CpG 位点
meth <- unite(filtered, destrand = FALSE)

# 差异甲基化检测
diff <- calculateDiffMeth(meth)

# 筛选显著的差异甲基化位点
# 甲基化差异 > 25%，q值 < 0.01
hyper <- getMethylDiff(diff, difference = 25, qvalue = 0.01, type = "hyper")
hypo <- getMethylDiff(diff, difference = 25, qvalue = 0.01, type = "hypo")