五、WGBS数据分析流程

1. 原始数据质控
- TrimGalore：去除接头序列和低质量碱基（对亚硫酸盐数据优化）
- FastQC：检查碱基质量、GC含量、长度分布
- 特别关注：亚硫酸盐转化效率（通过未甲基化对照λ噬菌体DNA评估）
2. 序列比对
亚硫酸盐转化的特殊性：
- C→U转化后，DNA序列中C极为稀少（仅剩甲基化C）
- 两条链不再互补（仅一条链发生C→T转化）
- 需要专门的比对算法
两种主要比对策略：
三字符法（3-letter）：将参考基因组中所有C替换为T，reads也进行C→T转换，然后比对。优点：简单快速；缺点：信息量损失，部分位点无法区分。
通配符法（Wild-card）：将参考基因组中所有C替换为Y（C或T），保留C的两种可能性。优点：信息完整；缺点：比对复杂度增加，可能引入偏差。
常用比对软件：
- Bismark：基于Bowtie/Bowtie2，进行C→T和G→A双重转换，准确率高，使用最广泛
- BSMAP：通配符法，不转换基因组，创建种子列表进行比对
- Bwa-meth：基于BWA-MEM，支持长读长比对
- Walt：针对亚硫酸盐数据的专门优化算法
3. 甲基化水平推断
对于每个CpG位点：
$$\text{甲基化水平} = \frac{\text{覆盖该位点的甲基化reads数}}{\text{覆盖该位点的总reads数}}$$
影响因素和校正：
- DNA片段末端修复引入的非甲基化C
- 接头序列污染
- 3'端测序质量下降
- 不均衡PCR扩增
- 5-hmC的干扰（亚硫酸盐无法区分5-mC和5-hmC）
4. 差异甲基化区域（DMR）鉴定
滑动窗口法：
- 将基因组划分为固定大小的窗口（如1 kb）
- 比较两组样本间每个窗口的平均甲基化水平
- 使用t检验、Wilcoxon检验或β-二项模型计算显著性
- 筛选连续显著窗口，合并为DMR
常用软件：
- methylKit：R包，支持滑动窗口和CpG-wise分析
- BSmooth：基于局部似然平滑，适合WGBS数据
- MOABS：基于分层贝叶斯模型，整合生物学重复信息
- dmrFinder：结合统计检验和区域合并
5. DMR注释和功能分析
- 基因组位置注释：启动子、外显子、内含子、基因间区
- 基因关联：关联到最近的基因或TSS
- 功能富集：GO/KEGG富集分析、motif分析
- 整合分析：结合RNA-seq数据，分析DMR与差异表达基因的关系