DNA甲基化WGBS分析Bismark¶

DNA甲基化检测原理：

正常DNA：  ...C G T A C G A C G...
                ^       ^     ^    ← CpG中的C可能被甲基化

亚硫酸氢盐处理：
  甲基化的C → 不变（还是C）
  未甲基化的C → 变成U → 测序读成T

处理前：   mC G T A  C G A mC G   （m = 甲基化）
处理后：    C G T A  T G A  C G
                     ↑            ← 这个C变成T = 未甲基化
                              ↑   ← 这个C不变 = 甲基化

Bismark的工作：
1. 把基因组变成两个版本：
   - C→T版本（模拟正链处理后的基因组）
   - G→A版本（模拟负链处理后的基因组）
2. 把测序reads比对到这两个基因组
3. 对每个CpG位点统计：
   - 读成C的reads数 = 甲基化reads
   - 读成T的reads数 = 未甲基化reads
   - 甲基化水平 = C/(C+T) × 100%

# ========== 安装Bismark ==========
conda create -n wgbs python=3.10  # 创建WGBS分析环境
conda activate wgbs  # 激活环境

conda install -c bioconda bismark  # 安装Bismark
conda install -c bioconda bowtie2  # 安装Bowtie2比对器
conda install -c bioconda samtools  # 安装SAMtools
conda install -c bioconda trim-galore  # 安装Trim Galore质控
conda install -c bioconda fastqc  # 安装FastQC

# 验证安装
bismark --version  # 查看Bismark版本
bowtie2 --version  # 查看Bowtie2版本

# ========== 准备基因组 ==========
# 下载参考基因组（以人类hg38为例）
mkdir -p ~/genomes/hg38  # 创建基因组目录
cd ~/genomes/hg38
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz  # 下载基因组
gunzip hg38.fa.gz  # 解压

# Bismark基因组准备（建索引，非常重要！）
bismark_genome_preparation \
    --bowtie2 \               # 使用Bowtie2
    --parallel 4 \            # 4线程加速
    ~/genomes/hg38/           # 基因组所在目录
# 这一步会创建C→T和G→A两个转换版本的基因组索引
# 耗时约30-60分钟

# ========== 第一步：质量控制 ==========
# 运行FastQC
fastqc -t 4 -o fastqc_raw/ sample_R1.fq.gz sample_R2.fq.gz  # 对原始数据做质控

# ========== 第二步：Trim Galore质量修剪 ==========
# WGBS数据需要特殊的trim参数
trim_galore \
    --paired \                  # 双端数据
    --quality 20 \              # 质量阈值（Phred≥20）
    --length 30 \               # 最短保留长度
    --fastqc \                  # trim后自动运行FastQC
    --cores 4 \                 # 4核心并行
    -o trimmed/ \               # 输出目录
    sample_R1.fq.gz \           # R1文件
    sample_R2.fq.gz             # R2文件

# Trim Galore对WGBS数据会自动：
# - 去除接头序列
# - 去除低质量碱基
# - 对于PBAT文库还会去除前几个碱基的偏差

# ========== 第三步：Bismark比对 ==========
bismark \
    --genome ~/genomes/hg38/ \     # 基因组目录（含Bismark索引）
    --bowtie2 \                    # 使用Bowtie2比对
    --parallel 4 \                 # 4个并行实例
    --temp_dir /tmp/ \             # 临时文件目录
    -o bismark_output/ \           # 输出目录
    -1 trimmed/sample_R1_val_1.fq.gz \  # trim后的R1
    -2 trimmed/sample_R2_val_2.fq.gz    # trim后的R2

# 关键输出文件：
# sample_R1_val_1_bismark_bt2_pe.bam  ← 比对后的BAM文件
# sample_R1_val_1_bismark_bt2_PE_report.txt  ← 比对报告

# 查看比对率
grep "Mapping efficiency" bismark_output/*_report.txt  # 查看比对效率
# 正常应该在 40-80% 之间（WGBS比对率比普通WGS低）

# ========== 第四步：去除PCR重复 ==========
deduplicate_bismark \
    --paired \                     # 双端数据
    --bam \                        # 输入BAM格式
    -o deduplicated/ \             # 输出目录
    bismark_output/sample_R1_val_1_bismark_bt2_pe.bam  # 输入BAM文件

# 查看去重比例
grep "Total number duplicates" deduplicated/*.deduplication_report.txt
# 重复率通常在 5-30% 之间

# ========== 第五步：提取甲基化信息（核心步骤！） ==========
bismark_methylation_extractor \
    --paired-end \                 # 双端数据
    --comprehensive \              # 输出所有类型的甲基化信息
    --merge_non_CpG \              # 合并非CpG甲基化数据
    --bedGraph \                   # 生成bedGraph格式文件
    --cytosine_report \            # 生成全基因组CpG报告
    --genome_folder ~/genomes/hg38/ \  # 基因组目录
    --parallel 4 \                 # 4线程
    -o methylation/ \              # 输出目录
    deduplicated/sample.deduplicated.bam  # 去重后的BAM

# 关键输出文件说明：
# CpG_context_sample.txt            ← 每个CpG位点的甲基化状态
# sample.bedGraph.gz                ← bedGraph格式甲基化水平
# sample.bismark.cov.gz             ← 覆盖度文件
# sample.CpG_report.txt             ← 全基因组CpG报告（最全面）

# ========== 第六步：生成分析报告 ==========
bismark2report  # 生成HTML可视化报告
bismark2summary  # 生成多样品汇总报告

#!/usr/bin/env python3
"""解析Bismark甲基化提取结果"""

import pandas as pd  # 数据处理
import numpy as np  # 数值计算
import matplotlib.pyplot as plt  # 画图

# ========== 读取bismark.cov文件 ==========
# 格式：chr start end methylation_pct count_methylated count_unmethylated
cov = pd.read_csv(
    "methylation/sample.bismark.cov.gz",
    sep="\t",
    header=None,
    names=["chr", "start", "end", "meth_pct", "count_meth", "count_unmeth"],  # 列名
    compression="gzip"  # 压缩格式
)

print(f"总CpG位点数: {len(cov)}")
print(f"平均甲基化水平: {cov['meth_pct'].mean():.1f}%")

# ========== 过滤低覆盖度位点 ==========
cov["total_reads"] = cov["count_meth"] + cov["count_unmeth"]  # 总reads数
cov_filtered = cov[cov["total_reads"] >= 5]  # 至少5x覆盖度
print(f"≥5x覆盖的CpG位点: {len(cov_filtered)} ({len(cov_filtered)/len(cov)*100:.1f}%)")

# ========== 甲基化水平分布 ==========
plt.figure(figsize=(10, 6))
plt.hist(cov_filtered["meth_pct"], bins=100, color="steelblue", edgecolor="none", alpha=0.8)
plt.xlabel("Methylation Level (%)", fontsize=14)
plt.ylabel("Number of CpG Sites", fontsize=14)
plt.title("CpG Methylation Level Distribution", fontsize=16)
plt.axvline(50, color="red", linestyle="--", alpha=0.5, label="50%")
plt.legend()
plt.tight_layout()
plt.savefig("methylation_distribution.png", dpi=300)
plt.close()
print("甲基化分布图已保存")

# ========== 按染色体统计 ==========
chr_stats = cov_filtered.groupby("chr").agg(
    n_sites=("meth_pct", "count"),  # CpG位点数
    mean_meth=("meth_pct", "mean"),  # 平均甲基化
    median_meth=("meth_pct", "median")  # 中位数甲基化
).reset_index()

# 只看常染色体
chr_order = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY"]
chr_stats = chr_stats[chr_stats["chr"].isin(chr_order)]
chr_stats["chr"] = pd.Categorical(chr_stats["chr"], categories=chr_order, ordered=True)
chr_stats = chr_stats.sort_values("chr")

print("\n各染色体甲基化水平：")
print(chr_stats.to_string(index=False))

========================================
DNA甲基化WGBS分析Bismark 速查表
========================================

【完整流程】
1. 质控              → fastqc + trim_galore
2. 基因组准备         → bismark_genome_preparation
3. 比对              → bismark --bowtie2
4. 去重              → deduplicate_bismark
5. 甲基化提取         → bismark_methylation_extractor
6. 报告              → bismark2report

【关键参数】
比对器               → --bowtie2（默认推荐）
文库类型             → 默认(directional) / --pbat / --non_directional
并行                 → --parallel 4
输出格式             → --bedGraph --cytosine_report

【甲基化类型】
CpG                  → 最重要，哺乳动物中绝大部分甲基化在CpG
CHG                  → 植物中常见（H=A/T/C）
CHH                  → 植物中常见

【质量标准】
比对率               → 40-80%（WGBS正常范围）
重复率               → <30%
覆盖度               → ≥5x才可靠，建议10x+
转换率               → >99%（用lambda DNA或spike-in验证）

【输出文件】
.bam                 → 比对结果
.bismark.cov.gz      → CpG甲基化水平+覆盖度
.bedGraph.gz         → bedGraph格式甲基化
.CpG_report.txt      → 全基因组CpG报告
_splitting_report.txt → 链信息分割报告

【面试考点】
Q: 亚硫酸氢盐测序的原理？
A: 未甲基化C→U→T，甲基化C不变，据此区分甲基化状态

Q: WGBS的比对率为什么比WGS低？
A: 亚硫酸氢盐处理降低了序列复杂度（C→T），增加了比对难度

Q: 怎么验证亚硫酸氢盐转换效率？
A: 加入lambda噬菌体DNA(完全未甲基化)，检测其转换率应>99%
========================================