Bisulfite测序数据分析¶

# 计算转化率（使用lambda DNA对照）
# lambda DNA无甲基化，所有C应该被转化为T
# 转化率 = (转化的C数) / (总C数) × 100%

# 在Bismark比对后查看lambda DNA的转化效率
grep "lambda" bismark_report.txt
# C methylated in CpG context: 0.3%  → 转化率 = 99.7%

# 安装工具
conda install -c bioconda trim-galore bismark bowtie2 samtools

# 质量检查
fastqc -t 8 sample_R1.fastq.gz sample_R2.fastq.gz

# Trim Galore预处理（自动检测adapter）
# 对WGBS数据：
trim_galore --paired --fastqc \
    --clip_R1 10 --clip_R2 10 \
    --three_prime_clip_R1 10 --three_prime_clip_R2 10 \
    -o trimmed/ \
    sample_R1.fastq.gz sample_R2.fastq.gz

# 对RRBS数据（需要加--rrbs参数）：
trim_galore --paired --fastqc --rrbs \
    -o trimmed/ \
    sample_R1.fastq.gz sample_R2.fastq.gz

# --rrbs: 自动处理MspI位点末端的人工甲基化偏差
# --clip_R1/R2: 去除5'端指定bp数
# --three_prime_clip: 去除3'端指定bp数

# 1. 构建Bismark参考基因组索引
bismark_genome_preparation --bowtie2 /path/to/genome/

# 2. Bismark比对
bismark --genome /path/to/genome/ \
    --bowtie2 \
    --parallel 4 \
    --temp_dir /tmp/ \
    -1 trimmed/sample_R1_val_1.fq.gz \
    -2 trimmed/sample_R2_val_2.fq.gz \
    -o aligned/

# 关键参数：
# --parallel 4: 使用4个并行实例（实际用8个CPU核）
# --non_directional: 如果是非方向性文库
# --pbat: 如果是PBAT文库

# 3. 去重复
deduplicate_bismark --bam aligned/sample_R1_val_1_bismark_bt2_pe.bam \
    --output_dir dedup/

# 4. 查看比对报告
cat aligned/sample_R1_val_1_bismark_bt2_PE_report.txt
# Mapping efficiency: 72.5%
# Unique best hit: 68.3%
# Sequences with no alignment: 27.5%

# 5. 排序和索引
samtools sort dedup/sample_R1_val_1_bismark_bt2_pe.deduplicated.bam \
    -o sorted_sample.bam
samtools index sorted_sample.bam

# 甲基化提取
bismark_methylation_extractor --paired-end --gzip \
    --bedGraph --cytosine_report \
    --genome_folder /path/to/genome/ \
    --parallel 4 \
    --comprehensive \
    --output methylation/ \
    dedup/sample_R1_val_1_bismark_bt2_pe.deduplicated.bam

# 输出文件说明：
# CpG_context_*.txt.gz     - CpG位点甲基化信息
# CHG_context_*.txt.gz     - CHG位点信息
# CHH_context_*.txt.gz     - CHH位点信息  
# *.bedGraph.gz            - bedGraph格式甲基化水平
# *.cov.gz                 - coverage文件（含甲基化/未甲基化读数）
# *_splitting_report.txt   - 报告文件
# *.cytosine_report.txt    - 全基因组cytosine报告

# M-bias检查
# 查看M-bias plot，决定是否需要额外trim
# 如果发现末端有明显偏差，需要在提取时加参数：
# --ignore 5 --ignore_r2 5 --ignore_3prime 5 --ignore_3prime_r2 5

# 生成全基因组报告
bismark2report
bismark2summary

# 使用DSS进行DMR分析
library(DSS)

# 读入Bismark coverage文件
# 格式: chr, pos, N (total), X (methylated)
read_bismark_cov <- function(file, sample_name) {
    df <- read.table(file, header=FALSE,
        col.names=c("chr","start","end","meth_pct","count_M","count_U"))
    data.frame(
        chr = df$chr,
        pos = df$start,
        N = df$count_M + df$count_U,
        X = df$count_M
    )
}

# 读取样本
ctrl1 <- read_bismark_cov("ctrl1.bismark.cov.gz", "ctrl1")
ctrl2 <- read_bismark_cov("ctrl2.bismark.cov.gz", "ctrl2")
treat1 <- read_bismark_cov("treat1.bismark.cov.gz", "treat1")
treat2 <- read_bismark_cov("treat2.bismark.cov.gz", "treat2")

# 创建BSseq对象
BSobj <- makeBSseqData(
    list(ctrl1, ctrl2, treat1, treat2),
    c("ctrl1", "ctrl2", "treat1", "treat2")
)

# DML检测（单位点差异）
dmlTest <- DMLtest(BSobj,
    group1 = c("ctrl1", "ctrl2"),
    group2 = c("treat1", "treat2"),
    smoothing = TRUE)

# DMR检测（区域差异）
dmrs <- callDMR(dmlTest,
    delta = 0.2,     # 最小甲基化差异20%
    p.threshold = 0.001,
    minlen = 50,     # 最小区域长度50bp
    minCG = 3,       # 最少包含3个CpG
    dis.merge = 50,  # 相邻DMR合并距离
    pct.sig = 0.5)   # 至少50%的CpG显著

# 查看结果
head(dmrs)
nrow(dmrs)

# 导出
write.csv(dmrs, "DMR_results.csv", row.names=FALSE)

library(methylKit)

# 读取Bismark coverage文件
file.list <- list(
    "ctrl1.bismark.cov.gz",
    "ctrl2.bismark.cov.gz",
    "treat1.bismark.cov.gz",
    "treat2.bismark.cov.gz"
)

# 创建methylRawList对象
myobj <- methRead(file.list,
    sample.id = list("ctrl1","ctrl2","treat1","treat2"),
    assembly = "hg38",
    treatment = c(0, 0, 1, 1),
    context = "CpG",
    mincov = 10,          # 最小覆盖度10x
    pipeline = "bismarkCoverage")

# 过滤极端覆盖度
filtered <- filterByCoverage(myobj,
    lo.count = 10,        # 最低10x
    lo.perc = NULL,
    hi.count = NULL,
    hi.perc = 99.9)       # 去除top 0.1%异常高覆盖

# 合并样本（取交集位点）
meth <- unite(filtered, destrand=FALSE)

# 计算差异甲基化
myDiff <- calculateDiffMeth(meth,
    overdispersion = "MN",  # beta-binomial模型
    adjust = "BH")          # BH校正

# 筛选DMC（差异甲基化位点）
myDiff25 <- getMethylDiff(myDiff,
    difference = 25,        # 至少25%差异
    qvalue = 0.01)          # FDR < 0.01

# 获取hyper/hypo甲基化
hyper <- getMethylDiff(myDiff, difference=25, qvalue=0.01, type="hyper")
hypo <- getMethylDiff(myDiff, difference=25, qvalue=0.01, type="hypo")

# 区域分析（tiling window）
tiles <- tileMethylCounts(myobj, win.size=1000, step.size=1000)
meth_tile <- unite(tiles, destrand=FALSE)
diff_tile <- calculateDiffMeth(meth_tile)

# 基因注释
library(genomation)
gene.obj <- readTranscriptFeatures("hg38_refseq.bed")
diffAnn <- annotateWithGeneParts(as(myDiff25,"GRanges"), gene.obj)

# 甲基化水平热图
library(ComplexHeatmap)
library(circlize)

# 提取DMR区域的甲基化值
meth_matrix <- percMethylation(meth)
dmr_sites <- as(myDiff25, "GRanges")

# 绘制热图
col_fun <- colorRamp2(c(0, 50, 100), c("blue", "white", "red"))
Heatmap(meth_matrix[1:100,],
    col = col_fun,
    name = "Methylation%",
    cluster_rows = TRUE,
    cluster_columns = TRUE,
    show_row_names = FALSE)

# DMR基因组分布
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

peakAnno <- annotatePeak(dmr_sites,
    TxDb = txdb,
    level = "gene")
plotAnnoPie(peakAnno)
plotDistToTSS(peakAnno)

# 甲基化与表达关联
# 合并甲基化和RNA-seq数据
# 启动子甲基化通常与基因沉默负相关

# === WGBS完整分析流程 ===

#!/bin/bash
set -e
GENOME="/ref/hg38"
THREADS=16

# 批量处理多个样本
for SAMPLE in ctrl1 ctrl2 treat1 treat2; do
    R1="raw/${SAMPLE}_R1.fastq.gz"
    R2="raw/${SAMPLE}_R2.fastq.gz"

    # 质控
    trim_galore --paired --fastqc --cores 4 \
        --clip_R1 10 --clip_R2 10 \
        -o trimmed/ $R1 $R2

    # 比对
    bismark --genome $GENOME --bowtie2 --parallel 4 \
        -1 trimmed/${SAMPLE}_R1_val_1.fq.gz \
        -2 trimmed/${SAMPLE}_R2_val_2.fq.gz \
        -o aligned/

    # 去重
    deduplicate_bismark --bam \
        aligned/${SAMPLE}_R1_val_1_bismark_bt2_pe.bam \
        --output_dir dedup/

    # 甲基化提取
    bismark_methylation_extractor --paired-end --gzip \
        --bedGraph --cytosine_report \
        --genome_folder $GENOME --parallel 4 \
        -o methylation/ \
        dedup/${SAMPLE}*.deduplicated.bam
done

# 汇总报告
bismark2report
bismark2summary

启动子CpG岛高甲基化 → 基因沉默
基因体高甲基化 → 活跃转录（矛盾但真实）
增强子低甲基化 → 增强子活化
重复序列高甲基化 → 转座子沉默
印记基因 → 等位基因特异性甲基化
X染色体失活 → 全染色体高甲基化