印记基因与等位基因特异性甲基化¶

# ===== 等位基因特异性甲基化分析流程 =====

# 1. 亚硫酸盐测序比对（已完成，使用 Bismark）
# 2. 获取杂合 SNP 信息（从 WGS/WES 数据或公共数据库）

# 3. 用 SNPsplit 按等位基因分割 BAM 文件
# SNPsplit 根据杂合 SNP 将 reads 分配到父源或母源
SNPsplit --snp_file sample_snps.txt \
    --bisulfite \
    sample.deduplicated.bam

# 输出：
# sample.deduplicated.genome1.bam  → 等位基因1（如父源）
# sample.deduplicated.genome2.bam  → 等位基因2（如母源）
# sample.deduplicated.unassigned.bam → 无法分配的 reads

# 4. 分别提取两个等位基因的甲基化信息
for allele in genome1 genome2; do
    bismark_methylation_extractor --paired-end \
        --comprehensive --cytosine_report \
        --genome_folder /ref/hg38/ \
        -o methylation/${allele}/ \
        sample.deduplicated.${allele}.bam
done

# ===== R: 等位基因特异性甲基化检测 =====
library(data.table)

# 读取两个等位基因的 CpG 报告
allele1 <- fread("methylation/genome1/CpG_report.txt",
    col.names = c("chr","pos","strand","meth","unmeth","ctx","tri"))
allele2 <- fread("methylation/genome2/CpG_report.txt",
    col.names = c("chr","pos","strand","meth","unmeth","ctx","tri"))

# 合并两个等位基因的数据
merged <- merge(
    allele1[, .(chr, pos, m1 = meth, u1 = unmeth)],
    allele2[, .(chr, pos, m2 = meth, u2 = unmeth)],
    by = c("chr", "pos"))

# 过滤低覆盖位点（两个等位基因都至少 5x）
merged <- merged[(m1 + u1) >= 5 & (m2 + u2) >= 5]

# 计算各等位基因甲基化率
merged$rate1 <- merged$m1 / (merged$m1 + merged$u1)
merged$rate2 <- merged$m2 / (merged$m2 + merged$u2)

# Fisher 精确检验：检测等位基因间甲基化差异
merged$pval <- apply(merged[, .(m1, u1, m2, u2)], 1, function(row) {
    mat <- matrix(as.numeric(row), nrow = 2, byrow = TRUE)
    fisher.test(mat)$p.value
})

# 校正多重检验
merged$padj <- p.adjust(merged$pval, method = "BH")

# 筛选显著的 ASM 位点
# 两个等位基因甲基化差异 > 30%，校正 p 值 < 0.05
asm_sites <- merged[
    padj < 0.05 &
    abs(rate1 - rate2) > 0.3
]

cat("检测到的 ASM 位点数：", nrow(asm_sites), "\n")

# 与已知印记区域比较
known_icr <- fread("known_imprinted_regions_hg38.bed",
    col.names = c("chr", "start", "end", "gene"))

# 检查 ASM 位点是否落在已知印记区域
library(GenomicRanges)
asm_gr <- GRanges(seqnames = asm_sites$chr,
    ranges = IRanges(start = asm_sites$pos, width = 1))
icr_gr <- GRanges(seqnames = known_icr$chr,
    ranges = IRanges(start = known_icr$start, end = known_icr$end))

overlaps <- findOverlaps(asm_gr, icr_gr)
cat("落在已知印记区域的 ASM 位点：",
    length(unique(queryHits(overlaps))), "\n")