150_Nanopore甲基化检测¶

# 概念演示：甲基化如何影响电流
# 实际中由Dorado模型自动处理
import numpy as np

# 模拟一个5-mer通过纳米孔的电流
# 未甲基化的 ACGTA
signal_unmethylated = np.array([105.2, 98.7, 112.3, 89.5, 95.1])
# 甲基化的 A(5mC)GTA - CpG位置的C被甲基化
signal_methylated = np.array([105.2, 103.4, 115.8, 89.5, 95.1])
# 甲基化导致电流特征改变，AI模型学习这种差异

# === Dorado basecalling + 5mC检测 ===
# 下载合适的模型
dorado download --model dna_r10.4.1_e8.2_400bps_sup@v4.3.0

# Basecalling + 5mCpG检测 (一步完成)
dorado basecaller \
  dna_r10.4.1_e8.2_400bps_sup@v4.3.0 \
  pod5_dir/ \
  --modified-bases 5mCG \
  --reference hg38.fa \
  --device "cuda:0" \
  > aligned_mod.bam

# 检测5mC和5hmC (双修饰)
dorado basecaller \
  dna_r10.4.1_e8.2_400bps_sup@v4.3.0 \
  pod5_dir/ \
  --modified-bases 5mCG_5hmCG \
  --reference hg38.fa \
  > aligned_mod_5mc_5hmc.bam

# 排序和索引
samtools sort -o sorted_mod.bam aligned_mod.bam
samtools index sorted_mod.bam

# === 检查modBAM的MM/ML tag ===
samtools view sorted_mod.bam | head -1 | tr '\t' '\n' | grep -E "^(MM|ML)"
# MM:Z:C+m,0,1,3;     # 甲基化位置(相对于比对位置)
# ML:B:C,230,15,200    # 甲基化概率(0-255)

# === 安装modkit ===
# 从GitHub releases下载: https://github.com/nanoporetech/modkit
# 或 conda install -c nanoporetech modkit

# === 基本用法：生成bedMethyl ===
modkit pileup \
  sorted_mod.bam \
  output.bedmethyl \
  --ref hg38.fa \
  --cpg \
  --combine-strands \
  --filter-threshold 0.75 \
  --threads 16

# 参数说明:
# --cpg               : 只输出CpG位点
# --combine-strands   : 合并正负链的CpG信息
# --filter-threshold  : 概率阈值(>=0.75判为甲基化)
# --threads           : 线程数

# === 输出bedMethyl格式说明 ===
# chr  start  end  modified_base  score  strand  ...  N_valid  fraction_modified  N_mod  N_canon  N_other  N_delete  N_fail  N_diff  N_nocall
# chr1 10468  10469  m             953    +      ...  21       0.952              20     1        0       0         0       0       0

# === 按区间统计 ===
# 统计promoter区域的甲基化
modkit pileup \
  sorted_mod.bam \
  promoter_methyl.bedmethyl \
  --ref hg38.fa \
  --cpg \
  --combine-strands \
  --include-bed promoters.bed \
  --threads 16

# === 甲基化摘要统计 ===
modkit summary sorted_mod.bam
# 输出全局甲基化比例、检测到的修饰类型等

# === 提取特定区域 ===
modkit pileup \
  sorted_mod.bam \
  region_methyl.bedmethyl \
  --ref hg38.fa \
  --region chr1:1000000-2000000 \
  --cpg

library(ggplot2)
library(GenomicRanges)

# 1. 读取bedMethyl文件
bedmethyl <- read.delim("output.bedmethyl", header = FALSE,
  col.names = c("chr", "start", "end", "mod_base", "score", "strand",
                "start2", "end2", "color", "N_valid", "fraction_mod",
                "N_mod", "N_canon", "N_other", "N_delete", "N_fail",
                "N_diff", "N_nocall"))

# 2. 全基因组甲基化分布
ggplot(bedmethyl, aes(x = fraction_mod)) +
  geom_histogram(bins = 50, fill = "steelblue", alpha = 0.7) +
  labs(x = "Methylation Fraction", y = "Number of CpG sites",
       title = "Genome-wide CpG Methylation Distribution") +
  theme_minimal()

# 3. 染色体水平甲基化
chr_methyl <- aggregate(fraction_mod ~ chr, data = bedmethyl, FUN = mean)
ggplot(chr_methyl, aes(x = chr, y = fraction_mod)) +
  geom_bar(stat = "identity", fill = "coral") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Average Methylation by Chromosome")

# 4. 特定基因区域的甲基化
gene_region <- bedmethyl[bedmethyl$chr == "chr17" &
                          bedmethyl$start >= 43044295 &
                          bedmethyl$end <= 43170245, ]  # BRCA1区域
ggplot(gene_region, aes(x = start, y = fraction_mod)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_smooth(method = "loess", span = 0.1) +
  labs(x = "Genomic Position", y = "Methylation",
       title = "BRCA1 Locus Methylation") +
  theme_minimal()

# 5. IGV可视化
# 将modBAM加载到IGV，开启"Color by modification"模式
# 红色=甲基化, 蓝色=未甲基化

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

# 1. 准备输入数据 (从bedMethyl转换)
# DSS需要: chr, pos, N(总覆盖), X(甲基化reads)
prepare_dss_input <- function(bedmethyl_file) {
  bed <- read.delim(bedmethyl_file, header = FALSE)
  data.frame(
    chr = bed$V1,
    pos = bed$V2,
    N = bed$V10,        # N_valid
    X = bed$V12         # N_mod (甲基化reads数)
  )
}

# 读取tumor和normal样本
tumor1 <- prepare_dss_input("tumor1.bedmethyl")
tumor2 <- prepare_dss_input("tumor2.bedmethyl")
normal1 <- prepare_dss_input("normal1.bedmethyl")
normal2 <- prepare_dss_input("normal2.bedmethyl")

# 2. 创建BSseq对象
bs_obj <- makeBSseqData(
  dat = list(tumor1, tumor2, normal1, normal2),
  sampleNames = c("Tumor1", "Tumor2", "Normal1", "Normal2")
)

# 3. DML检测 (差异甲基化位点)
dml_test <- DMLtest(
  bs_obj,
  group1 = c("Tumor1", "Tumor2"),
  group2 = c("Normal1", "Normal2"),
  smoothing = TRUE       # 平滑化(推荐)
)

# 4. DMR调用
dmrs <- callDMR(
  dml_test,
  delta = 0.2,           # 甲基化差异阈值(20%)
  p.threshold = 0.001,   # P值阈值
  minlen = 50,           # 最小DMR长度
  minCG = 3,             # 最小CpG数
  dis.merge = 50,        # 合并距离
  pct.sig = 0.5          # 区域内显著CpG比例
)

head(dmrs)
# chr  start    end      length  nCG  meanMethy1  meanMethy2  diff.Methy  areaStat

# 5. 注释DMR
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
dmr_gr <- GRanges(seqnames = dmrs$chr,
                   ranges = IRanges(start = dmrs$start, end = dmrs$end))
dmr_anno <- annotatePeak(dmr_gr, TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene)
plotAnnoPie(dmr_anno)

# === 单倍型标记 + 甲基化 ===
# 1. 变异检测 (用于分型)
clair3 --bam_fn sorted_mod.bam \
  --ref_fn hg38.fa \
  --output clair3_output/ \
  --platform ont \
  --model_path ont_model/

# 2. 使用WhatsHap进行单倍型分型
whatshap phase \
  clair3_output/merge_output.vcf.gz \
  sorted_mod.bam \
  --reference hg38.fa \
  --output phased.vcf.gz

# 3. 单倍型标记reads
whatshap haplotag \
  phased.vcf.gz \
  sorted_mod.bam \
  --reference hg38.fa \
  --output haplotagged_mod.bam

samtools index haplotagged_mod.bam

# 4. 按单倍型分离甲基化
# Haplotype 1
samtools view -b -d HP:1 haplotagged_mod.bam > hap1_mod.bam
samtools index hap1_mod.bam
modkit pileup hap1_mod.bam hap1.bedmethyl --ref hg38.fa --cpg --combine-strands

# Haplotype 2
samtools view -b -d HP:2 haplotagged_mod.bam > hap2_mod.bam
samtools index hap2_mod.bam
modkit pileup hap2_mod.bam hap2.bedmethyl --ref hg38.fa --cpg --combine-strands

# 5. 或使用modkit的单倍型感知模式
modkit pileup \
  haplotagged_mod.bam \
  haplotype_methyl.bedmethyl \
  --ref hg38.fa \
  --cpg \
  --partition-tag HP \
  --prefix hap
# 输出: hap_1.bedmethyl, hap_2.bedmethyl, hap_ungrouped.bedmethyl

# === 完整流程: POD5 → 甲基化分析 ===

# 1. Basecalling + 甲基化检测
dorado basecaller sup pod5_dir/ \
  --modified-bases 5mCG \
  --reference hg38.fa \
  --device "cuda:all" \
  > calls.bam

# 2. 排序索引
samtools sort -@ 16 -o sorted.bam calls.bam
samtools index sorted.bam

# 3. 比对统计
samtools flagstat sorted.bam
mosdepth --threads 4 --no-per-base coverage sorted.bam

# 4. 生成bedMethyl
modkit pileup sorted.bam methylation.bedmethyl \
  --ref hg38.fa \
  --cpg \
  --combine-strands \
  --filter-threshold 0.75 \
  --threads 16

# 5. 甲基化统计摘要
modkit summary sorted.bam

# 6. 质控
modkit sample-probs sorted.bam sample_probs.tsv
# 查看修饰概率分布是否合理

# 7. DMR分析 (在R中)
Rscript run_dmr_analysis.R \
  --tumor tumor.bedmethyl \
  --normal normal.bedmethyl \
  --output dmr_results/ \
  --delta 0.2 \
  --pvalue 0.001