DNA甲基化检测方法比较¶

# === WGBS数据分析流程 ===
# 使用Bismark做比对，methylKit做差异分析

# --- Bismark比对流程（Shell命令） ---
# 1. 基因组准备
# bismark_genome_preparation --bowtie2 genome/  # 构建bisulfite转换后的索引

# 2. 比对
# bismark --genome genome/ -1 sample_R1.fq.gz -2 sample_R2.fq.gz  # 比对

# 3. 去重复
# deduplicate_bismark sample.bam  # 去PCR重复

# 4. 提取甲基化信息
# bismark_methylation_extractor --cytosine_report sample.deduplicated.bam

# --- methylKit差异甲基化分析（R） ---
library(methylKit)

# 1. 读取Bismark输出的CpG覆盖度文件
file_list <- list("sample1.cov", "sample2.cov", "sample3.cov",
                   "sample4.cov", "sample5.cov", "sample6.cov")
sample_ids <- c("disease1","disease2","disease3","ctrl1","ctrl2","ctrl3")
treatment <- c(1, 1, 1, 0, 0, 0)  # 1=疾病组，0=对照组

myobj <- methRead(
  file_list,
  sample.id = sample_ids,
  treatment = treatment,
  assembly = "hg38",       # 参考基因组版本
  context = "CpG",         # 只看CpG位点
  mincov = 10              # 最少10x覆盖度
)

# 2. 过滤和归一化
filtered <- filterByCoverage(myobj, lo.count=10, hi.perc=99.9)  # 过滤极端覆盖
normalized <- normalizeCoverage(filtered)  # 归一化覆盖度

# 3. 合并样本（只保留所有样本都有覆盖的位点）
meth <- unite(normalized, destrand=FALSE)

# 4. 差异甲基化分析
diff <- calculateDiffMeth(meth)  # 计算差异甲基化

# 5. 筛选差异甲基化位点（DMC）
dmc <- getMethylDiff(diff,
                     difference = 25,  # 甲基化差异>25%
                     qvalue = 0.01)    # FDR<0.01
print(paste("差异甲基化CpG数:", nrow(dmc)))

# 6. 注释到基因组区域
library(genomation)
gene_obj <- readTranscriptFeatures("hg38_genes.bed")
annotateWithGeneParts(as(dmc, "GRanges"), gene_obj)  # 看DMC在启动子/外显子/内含子的分布

# === Illumina EPIC芯片数据分析要点 ===
# 通常用R的minfi或ChAMP包，这里展示基本概念

import pandas as pd
import numpy as np

# 1. 理解β值和M值的转换
beta = 0.75   # β值=0.75表示75%的DNA分子在该位点被甲基化
M_value = np.log2(beta / (1 - beta))  # β→M值转换
print(f"β={beta} → M={M_value:.2f}")  # M=1.58

# β值直觉好理解（0-1），M值统计性质更好（正态分布）
# 差异分析用M值，可视化用β值

# 2. 简单的差异甲基化位点筛选
data = pd.read_csv("methylation_beta_values.csv", index_col=0)  # CpG × 样本
groups = pd.read_csv("sample_groups.csv", index_col=0)

# 计算组间β值差异
disease = data.loc[:, groups["Group"] == "Disease"].mean(axis=1)
control = data.loc[:, groups["Group"] == "Control"].mean(axis=1)
delta_beta = disease - control  # Δβ
print(f"|Δβ| > 0.2 的CpG数: {(abs(delta_beta) > 0.2).sum()}")