表观遗传学综合性使用范例

表观遗传学分析在实际研究中有着广泛的应用场景，从癌症甲基化标志物的发现到免疫细胞活化过程中的染色质重塑，从三维基因组互作解析到早期胚胎发育的表观重编程，不同的研究问题需要组合使用不同的表观遗传学技术。本节通过五个具有代表性的实际案例，展示如何将前面介绍的表观遗传学分析方法整合应用于具体的科学问题。这些案例涵盖了DNA甲基化分析、染色质可及性分析、染色质互作分析、组蛋白修饰分析和表观遗传时钟等多个方面，帮助你理解如何根据研究目标设计合适的分析流程并解读结果。

范例一：癌症中DNA甲基化异常与预后标志物（WGBS + 差异甲基化 + EWAS）

DNA甲基化异常是癌症中最常见的表观遗传学改变之一，肿瘤抑制基因启动子的高甲基化导致的基因沉默是肿瘤发生发展的重要机制。一个典型的应用场景是收集100例肝细胞癌组织与配对的癌旁正常组织，进行全基因组重亚硫酸盐测序（WGBS），目的是鉴定与预后相关的甲基化标志物并解释肿瘤抑制基因沉默的机制。这种配对设计可以有效控制个体间的遗传背景差异，使得检测到的甲基化变化更可能与肿瘤相关而非个体差异。

整个分析流程从WGBS数据的质控和比对开始。使用fastp对原始测序数据进行质量过滤和接头修剪后，使用Bismark将重亚硫酸盐处理后的reads比对到人类参考基因组上。Bismark是目前最广泛使用的重亚硫酸盐比对工具，它将参考基因组进行原位重亚硫酸盐转化后建立索引，然后将测序reads与转化后的基因组进行比对，从而正确处理重亚硫酸盐测序中C到T的转化。比对完成后，使用deduplicate_bismark去除PCR重复，然后使用bismark_methylation_extractor提取每个CpG位点的甲基化状态信息。

甲基化水平的计算采用beta值，即beta = mC / (mC + umC)，其中mC是甲基化读段数，umC是未甲基化读段数。beta值的范围在0到1之间，0表示完全未甲基化，1表示完全甲基化。对于下游的统计检验，M值（M = log2(beta / (1 - beta))）由于更符合正态分布假设，在差异分析中更为常用。在计算甲基化水平之前，需要过滤低覆盖度的CpG位点，通常要求每个位点的覆盖度至少为5-10个reads，以避免由于采样误差导致的甲基化水平估计不准确。

差异甲基化分析使用DSS（Dispersion Shrinkage for Sequencing data）包进行。DSS采用基于beta-binomial回归模型的方法，考虑了甲基化数据中常见的过度离散问题，并使用经验Bayes方法对离散参数进行收缩估计以提高统计功效。在这个案例中，设置筛选标准为绝对beta值差异大于0.2且FDR小于0.01，最终鉴定出约2000个差异甲基化位点（DMP）和100个差异甲基化区域（DMR）。这些DMR中，启动子区域的高甲基化尤其值得关注，因为启动子高甲基化通常与基因表达沉默相关。

注释与可视化是理解DMR生物学意义的关键步骤。使用annotatr或ChIPseeker等R包将DMR注释到基因组功能区域，包括启动子、基因体、增强子和重复元件等。在这个肝癌案例中，启动子区域的DMR附近经常关联表达下调的肿瘤抑制基因，如CDKN2A（p16）和RASSF1A。CDKN2A是细胞周期调控的关键基因，其启动子高甲基化导致的沉默在多种癌症中反复出现；RASSF1A则参与Ras信号通路的负调控，其甲基化沉默在肝癌中也是高频事件。使用IGV或UCSC Genome Browser可以在基因组浏览器中直观查看这些区域的甲基化状态变化。

甲基化与基因表达的整合分析需要同一批样本的RNA-seq数据。对100对肝癌和正常组织同时进行RNA-seq测序，可以获得匹配的甲基化和表达数据。通过计算启动子甲基化水平与对应基因表达值之间的相关性，可以验证甲基化对基因表达的调控作用。在这个案例中，超甲基化启动子与基因表达呈显著负相关，特别是那些启动子DMR对应的基因，其表达下调的幅度更为显著。这种整合分析不仅验证了甲基化的功能效应，还帮助筛选出最可能具有生物学意义的甲基化变化。

生存分析是将甲基化标志物与临床预后联系起来的核心步骤。基于关键CpG位点的甲基化值对患者进行层次聚类，可以将患者分为高风险和低风险两组。然后使用Cox比例风险模型评估甲基化水平与生存时间的关联，计算风险比（hazard ratio）和置信区间。在这个案例中，高风险组的患者甲基化模式与更差的总体生存率显著相关，提示这些甲基化标志物具有预后预测价值。进一步的多因素Cox回归分析可以在校正年龄、分期等临床因素后，评估甲基化标志物的独立预后价值。

EWAS分析将甲基化与临床特征进行关联，扩展了甲基化标志物的应用范围。对年龄、肿瘤分期、血清标志物（如AFP）等临床特征做关联回归分析，可以发现与特定临床特征相关的甲基化位点。例如，某些CpG位点的甲基化水平与肿瘤分期显著相关，可能反映了肿瘤进展过程中的表观遗传演变。年龄相关的甲基化变化则需要特别关注，因为年龄是甲基化最强的混杂因素之一，在分析中必须作为协变量进行校正。

# WGBS数据处理流程 (Bismark)
# 1. 构建重亚硫酸盐基因组索引
bismark_genome_preparation /path/to/hg38/

# 2. 比对 (双端测序)
bismark --genome /path/to/hg38/ \
    -1 tumor_R1.fq.gz -2 tumor_R2.fq.gz \
    --parallel 8 --output_dir bismark_tumor/

bismark --genome /path/to/hg38/ \
    -1 normal_R1.fq.gz -2 normal_R2.fq.gz \
    --parallel 8 --output_dir bismark_normal/

# 3. 去重
deduplicate_bismark --paired bismark_tumor/tumor_bismark_bt2_pe.bam
deduplicate_bismark --paired bismark_normal/normal_bismark_bt2_pe.bam

# 4. 甲基化提取
bismark_methylation_extractor \
    --paired-end --gzip --comprehensive \
    --bedGraph --CX_context \
    bismark_tumor/tumor_bismark_bt2_pe.deduplicated.bam

bismark_methylation_extractor \
    --paired-end --gzip --comprehensive \
    --bedGraph --CX_context \
    bismark_normal/normal_bismark_bt2_pe.deduplicated.bam

# 差异甲基化分析 (DSS)
library(DSS)
library(bsseq)

# 构建BSseq对象 (示例: 3对肿瘤/正常样本)
chr <- rep("chr9", 5)
pos <- c(21994789, 21994801, 21994810, 21994825, 21994833)

# 肿瘤样本甲基化计数
M_tumor <- cbind(c(18, 20, 15, 19, 17),
                 c(16, 19, 14, 20, 18),
                 c(17, 21, 16, 18, 16))
# 正常样本甲基化计数
M_normal <- cbind(c(3, 5, 2, 4, 3),
                  c(4, 6, 3, 5, 4),
                  c(2, 4, 1, 3, 2))
# 总覆盖度
Cov_tumor <- cbind(rep(20, 5), rep(20, 5), rep(20, 5))
Cov_normal <- cbind(rep(20, 5), rep(20, 5), rep(20, 5))

BSobj <- BSseq(
    chr = chr, pos = pos,
    M = cbind(M_tumor, M_normal),
    Cov = cbind(Cov_tumor, Cov_normal),
    sampleNames = c("T1", "T2", "T3", "N1", "N2", "N3")
)

# 差异甲基化检验
dmlTest <- DMLtest(BSobj,
                   group1 = c("T1", "T2", "T3"),
                   group2 = c("N1", "N2", "N3"))
head(dmlTest)

# 调用DML (差异甲基化位点)
dmls <- callDML(dmlTest, p.threshold = 0.001)
head(dmls)

# 调用DMR (差异甲基化区域)
dmrs <- callDMR(dmlTest, p.threshold = 0.01, delta = 0.2)
head(dmrs)

# 甲基化与表达整合 + 生存分析
library(survival)
library(survminer)

# 甲基化与表达相关性分析
meth_values <- c(0.85, 0.78, 0.12, 0.91, 0.65, 0.08, 0.72, 0.45)
expr_values <- c(1.2, 2.1, 8.5, 0.8, 3.2, 9.1, 2.5, 5.8)
cor_test <- cor.test(meth_values, expr_values, method = "spearman")
cat(sprintf("启动子甲基化与基因表达相关性: rho = %.3f, p = %.4e\n",
            cor_test$estimate, cor_test$p.value))

# 基于甲基化聚类分组
meth_matrix <- matrix(rnorm(100 * 50, mean = 0.5, sd = 0.2), nrow = 100)
hc <- hclust(dist(t(meth_matrix)), method = "ward.D2")
risk_group <- cutree(hc, k = 2)

# Cox回归生存分析
survival_data <- data.frame(
    time = rnorm(50, mean = 60, sd = 30),
    status = sample(0:1, 50, replace = TRUE, prob = c(0.3, 0.7)),
    group = factor(risk_group)
)
cox_model <- coxph(Surv(time, status) ~ group, data = survival_data)
summary(cox_model)

# Kaplan-Meier曲线
fit <- survfit(Surv(time, status) ~ group, data = survival_data)
ggsurvplot(fit, data = survival_data,
           pval = TRUE, risk.table = TRUE,
           xlab = "时间 (月)", ylab = "生存概率",
           title = "基于甲基化分组的生存曲线")

这个综合分析案例的意义在于它将DNA甲基化从单纯的描述性分析提升到了临床转化的层面。通过鉴定与预后相关的甲基化标志物，可以为肝癌患者的风险分层提供分子层面的依据，补充传统的组织病理学评估。更重要的是，启动子高甲基化导致的肿瘤抑制基因沉默机制为甲基化靶向治疗提供了理论基础——如果能够使用DNMT抑制剂逆转这些异常的甲基化，可能恢复肿瘤抑制基因的表达并抑制肿瘤生长。目前，地西他滨和阿扎胞苷等DNMT抑制剂已经在血液肿瘤中获得临床应用，这类EWAS研究为将它们扩展到实体瘤的治疗提供了候选靶点。

范例二：利用ATAC-seq和RNA-seq揭示免疫细胞活化时的染色质重塑

T细胞活化是适应性免疫应答的核心事件，理解这一过程中染色质可及性的动态变化对于揭示转录调控的分子机制至关重要。一个典型的应用场景是从人外周血分离初始CD4+T细胞，使用抗CD3/CD28抗体进行激活，在0小时、6小时和24小时三个时间点分别收集细胞进行ATAC-seq和RNA-seq测序。这种时间序列设计可以捕捉T细胞从静息状态到完全活化过程中染色质重塑的逐步变化，理解转录因子如何按时间顺序结合到基因组上并驱动基因表达的程序性变化。

ATAC-seq数据的处理从比对开始。使用Bowtie2将测序reads比对到人类参考基因组hg38上，由于ATAC-seq使用Tn5转座酶切割开放染色质区域，比对参数需要设置较大的插入片段上限（通常-X 2000）以容纳核小体包裹的较长片段。比对后需要去除线粒体reads，这是ATAC-seq数据处理中特别重要的一步，因为线粒体基因组没有核小体包裹，Tn5转座酶可以自由切割，导致线粒体reads在ATAC-seq数据中占比极高（可达30-50%），如果不去除会严重影响后续分析的准确性。

Tn5转座酶偏移校正是ATAC-seq数据分析中一个独特而关键的步骤。Tn5转转酶以二聚体形式结合DNA，在两个位置同时切割，导致插入位点的偏移。对于正链reads需要向3'端偏移+4bp，负链reads需要向5'端偏移-5bp，这种+4/-5的偏移校正可以准确定位Tn5插入位点，从而正确反映核小体游离区的边界。如果不进行偏移校正，峰值调用结果可能会出现偏移，影响后续的基序分析和足迹分析的准确性。

片段长度分布分析是ATAC-seq数据质量评估的重要指标。在理想的ATAC-seq数据中，片段长度分布应呈现明显的周期性模式：核小体游离片段（Nucleosome-Free Fragments, NFR）长度小于100bp，代表转录因子结合的开放区域；单核小体片段约200bp，代表被一个核小体包裹的DNA；双核小体片段约400bp，依此类推。这种周期性模式是ATAC-seq数据质量的重要标志，如果NFR峰不明显，可能提示实验质量不佳或Tn5转座酶活性不足。

峰值调用使用MACS2进行，采用--nomodel --shift -100 --extsize 200的参数设置，这是ATAC-seq数据的推荐参数。在每个时间点分别调用峰值后，可以获得各时间点的开放染色质区域集合。0小时样本的峰值代表初始T细胞中已经开放的染色质区域，6小时和24小时样本中新增的峰值则代表T细胞活化过程中新开放的区域。通过比较不同时间点的峰值集合，可以构建染色质可及性的动态变化图谱。

差异可及性分析使用DiffBind包进行。DiffBind以峰值区域的read计数为输入，使用DESeq2或edgeR的统计框架来识别在不同时间点之间显著差异的开放区域。在这个案例中，6小时与0小时比较鉴定出早期响应的差异可及性区域，24小时与0小时比较鉴定出晚期响应的区域。这些差异可及性区域反映了T细胞活化过程中染色质重塑的时序性，早期响应区域可能包含即时早期基因的调控元件，而晚期响应区域可能包含效应基因的调控元件。

基序富集分析使用HOMER进行，它可以识别差异可及性区域中富集的转录因子结合基序。在这个案例中，早期（0-6小时）差异开放区域富集AP-1家族转录因子（Fos、Jun）的结合基序，AP-1是T细胞受体信号通路的即时早期响应因子，在T细胞活化的启动阶段发挥关键作用。晚期（24小时）差异开放区域则富集NFAT和NF-κB的结合基序，这些转录因子分别响应钙信号和炎症信号，驱动T细胞效应功能的获得。这种时序性的转录因子活性模式揭示了T细胞活化过程中信号传导的层级关系。

足迹分析使用TOBIAS工具进行，它可以在核小体游离区内检测转录因子结合留下的"足迹"——即由于转录因子结合保护DNA免受Tn5切割而导致的可及性下降区域。与简单的基序富集分析不同，足迹分析可以推断转录因子的实际结合活性，而不仅仅是结合位点的存在。TOBIAS通过比较不同条件下足迹的深度变化，可以量化转录因子结合活性的动态变化。在这个案例中，AP-1家族转录因子在6小时即显示出强烈的足迹信号，而NFAT的足迹在24小时才变得显著，这与基序富集分析的结果相互印证。

RNA-seq差异表达与ATAC-seq数据的整合是理解染色质重塑功能效应的关键。对同一批样本进行RNA-seq测序，使用DESeq2鉴定差异表达基因后，将差异表达基因的启动子或增强子区域与差异可及性区域进行重叠分析。在这个案例中，共表达基因的启动子或增强子的开放动态与表达变化方向一致——可及性增加的区域对应的基因表达上调，反之亦然。这种一致性验证了染色质可及性变化对基因表达的调控作用。

可视化方面，IGV可以在基因组浏览器中展示特定基因座在T细胞激活前后的染色质可及性变化。例如，IL2基因座在0小时时启动子区域几乎没有ATAC-seq信号，而在6小时和24小时信号逐渐增强，这与IL2作为T细胞关键效应因子在激活后大量表达的模式完全吻合。使用deepTools的computeMatrix和plotHeatmap可以批量生成多个基因区域的可及性热图，直观展示不同时间点的全局染色质重塑模式。

# ATAC-seq数据处理流程
# 1. 比对 (Bowtie2)
bowtie2 -x /path/to/hg38_index \
    -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
    --very-sensitive -X 2000 \
    -p 12 -S sample.sam 2> align_log.txt

# 2. 转换、排序、索引
samtools view -bS sample.sam | samtools sort -o sample.sorted.bam
samtools index sample.sorted.bam

# 3. 去除线粒体reads
samtools idxstats sample.sorted.bam | cut -f 1 | grep -v "chrM" | \
    xargs samtools view -b sample.sorted.bam > sample.noMT.bam

# 4. 去除PCR重复
picard MarkDuplicates INPUT=sample.noMT.bam OUTPUT=sample.dedup.bam \
    METRICS_FILE=sample.metrics.txt REMOVE_DUPLICATES=true
samtools index sample.dedup.bam

# 5. Tn5偏移校正 + 峰值调用 (MACS2)
# 正链偏移+4, 负链偏移-5
samtools view -H sample.dedup.bam > header.sam
samtools view sample.dedup.bam | \
    awk '($2==0 || $2==16)' | \
    awk '{if($2==0) $4=$4+4; else $4=$4-5; print}' | \
    cat header.sam - | samtools view -b - > sample.shifted.bam

macs2 callpeak -t sample.shifted.bam -f BAMPE \
    --nomodel --shift -100 --extsize 200 \
    -g hs -n sample_peaks --keep-dup all -q 0.05

# 差异可及性分析 (DiffBind)
library(DiffBind)

# 创建样本信息表
samples <- dba.sampleSheet(
    SampleID = c("T0_1", "T0_2", "T0_3", "T6_1", "T6_2", "T6_3",
                  "T24_1", "T24_2", "T24_3"),
    Condition = rep(c("T0h", "T6h", "T24h"), each = 3),
    bamReads = list.files("bam/", pattern = "*.dedup.bam$", full.names = TRUE),
    Peaks = list.files("peaks/", pattern = "*.narrowPeak$", full.names = TRUE),
    PeakCaller = "narrow"
)

# 创建DBA对象
dba_obj <- dba(sampleSheet = samples)

# 计算一致性峰值 (至少在2个样本中出现)
dba_obj <- dba.count(dba_obj, minOverlap = 2)

# 差异分析: T6h vs T0h
dba_contrast <- dba.contrast(dba_obj, categories = DBA_CONDITION,
                             minMembers = 2)
dba_analyze <- dba.analyze(dba_contrast)

# 提取差异可及性区域
diff_peaks <- dba.report(dba_analyze, th = 0.05, fold = 1)
cat(sprintf("差异可及性区域数: %d\n", length(diff_peaks)))

# 基序富集分析 (HOMER)
# 早期响应区域 (T6h vs T0h) 的基序富集
findMotifsGenome.pl diff_peaks_early.bed hg38 homer_early_out/ -size 200 -mask

# 晚期响应区域 (T24h vs T0h) 的基序富集
findMotifsGenome.pl diff_peaks_late.bed hg38 homer_late_out/ -size 200 -mask

# 足迹分析 (TOBIAS)
# 1. 纠正Tn5偏移
TOBIAS ATACorrect --bam sample.dedup.bam \
    --genome hg38.fa --out sample_corrected.bw

# 2. 足迹检测
TOBIAS FootprintScores --bam sample.dedup.bam \
    --regions peaks.narrowPeak \
    --genome hg38.fa --out footprints.bw

# 3. 比较不同条件的足迹
TOBIAS BINDetect --motifs JASPAR2022.motif \
    --footprints footprints_T0.bw footprints_T6.bw footprints_T24.bw \
    --regions peaks.narrowPeak \
    --genome hg38.fa --out BINDetect_output/

这个综合分析案例的意义在于它为理解信号诱导的转录调控网络提供了动力学蓝图。通过整合ATAC-seq和RNA-seq的时间序列数据，我们不仅知道了哪些基因在T细胞活化中被激活，还揭示了这些基因被激活的染色质层面的时序机制——AP-1家族转录因子首先打开染色质，随后NFAT和NF-κB等因子结合到新开放的区域驱动效应基因的表达。这种层级性的调控模式对于理解免疫应答的精确控制具有重要意义，也为免疫相关疾病的干预提供了新的思路——如果能够特异性地阻断早期转录因子的染色质开放功能，可能选择性地抑制病理性T细胞活化而不影响已经建立的正常免疫记忆。

范例三：Hi-C鉴定疾病相关增强子-启动子环路（Capture Hi-C + GWAS）

全基因组关联研究（GWAS）已经发现了大量与复杂疾病相关的遗传变异，但其中约90%位于非编码区域，这些变异如何影响疾病风险一直是遗传学中的核心难题。一个典型的应用场景是冠状动脉疾病的GWAS研究发现多个非编码风险SNP位于增强子域内，但这些增强子的靶基因往往不是基因组线性距离上最近的基因，而是通过三维染色质互作连接的远端基因。为了鉴定这些非编码变异的靶基因，对血管平滑肌细胞进行H3K27ac HiChIP实验，利用H3K27ac标记捕获与活性增强子互作的启动子区域，从而将GWAS信号与靶基因联系起来。

HiChIP数据的处理使用HiC-Pro流程进行。HiC-Pro是一个综合性的Hi-C数据处理工具，支持从原始fastq文件到标准化接触矩阵的完整分析流程。对于HiChIP数据，HiC-Pro首先将reads比对到参考基因组，然后识别有效的配对读段，过滤自连接、未连接和PCR重复等噪音读段，最终构建不同分辨率的接触矩阵。HiChIP相比传统Hi-C的优势在于它通过免疫沉淀富集了与特定组蛋白修饰相关的染色质互作，因此在相同的测序深度下可以获得更高信噪比的互作信号，特别适合检测增强子-启动子互作。

峰检测首先对H3K27ac信号进行峰值调用，确定增强子的基因组位置。使用MACS2对HiChIP数据中的ChIP信号部分进行峰值调用，鉴定H3K27ac富集的区域作为候选增强子。这些增强子区域是后续互作分析的锚点，只有在增强子区域内的互作才会被纳入增强子-启动子环路的分析。在冠状动脉疾病的案例中，平滑肌细胞中鉴定出数千个H3K27ac峰值，这些峰值代表了细胞类型特异的活性增强子集合。

互作调用使用FitHiC2进行，它基于接触矩阵中的互作频率和基因组距离的期望分布来识别显著的染色质互作。FitHiC2考虑了基因组距离对互作频率的影响——距离越近的位点互作频率越高，这是Hi-C数据的固有特征。通过拟合距离-频率曲线并计算观测互作频率偏离期望的显著性，FitHiC2可以在不同距离范围内检测显著的互作事件。在这个案例中，FitHiC2产生了数千对显著的增强子-启动子互作，这些互作将增强子与其调控的靶基因连接起来。

将GWAS SNP与增强子区域进行重叠分析是连接遗传变异与靶基因的关键步骤。使用bedtools intersect将GWAS显著SNP（p < 5e-8）与H3K27ac峰值区域进行重叠，找到风险等位基因影响的增强子。在这个案例中，多个冠状动脉疾病风险SNP落在平滑肌细胞的增强子区域内，这些SNP可能通过改变转录因子结合位点来影响增强子的活性。例如，某些风险等位基因可能削弱AP-1或SRF等平滑肌关键转录因子的结合，从而降低增强子活性并影响靶基因的表达。

靶基因识别通过互作网络将增强子与启动子关联的基因联系起来。对于每个包含GWAS SNP的增强子，查找其通过HiChIP互作连接的启动子，启动子对应的基因即为候选靶基因。这种基于三维基因组互作的靶基因映射比简单的线性距离方法更准确，因为它考虑了染色质的空间折叠。在这个案例中，SORT1和CILP2等基因被鉴定为冠状动脉疾病风险SNP的候选靶基因。SORT1编码sortilin蛋白，参与脂质代谢，其表达变化与低密度脂蛋白胆固醇水平密切相关，这为理解冠状动脉疾病的发病机制提供了分子层面的解释。

实验验证是确认增强子-启动子互作和风险等位基因功能的关键步骤。荧光素酶报告实验可以验证增强子的调控活性和风险等位基因的影响。将包含风险等位基因或保护等位基因的增强子序列克隆到荧光素酶报告载体中，转染平滑肌细胞后测量荧光素酶活性。在这个案例中，风险等位基因对应的增强子活性显著低于保护等位基因，表明风险SNP降低了增强子的转录调控能力。表达数量性状位点（eQTL）分析则提供了遗传变异与基因表达关联的人群证据，利用GTEx等公共数据库可以验证互作基因在携带风险等位基因的个体中是否存在表达差异。

# HiChIP数据处理流程 (HiC-Pro)
# 1. 配置HiC-Pro参数
# config_hichip.txt中设置参考基因组、限制性酶切位点等

# 2. 运行HiC-Pro
/path/to/HiC-Pro/bin/HiC-Pro \
    -i fastq_dir/ \
    -o hichip_output/ \
    -c config_hichip.txt \
    -p 12

# 3. H3K27ac峰值调用 (MACS2)
macs2 callpeak -t hichip_ChIP.bam -c hichip_input.bam \
    -f BAMPE -g hs -n H3K27ac_peaks \
    --keep-dup all -q 0.01 --broad

# 4. 互作调用 (FitHiC2)
FitHiC2.py -b hichip_output/iced_matrix/5000/ \
    -o fithic2_output/ \
    -l 5000 -p 12 \
    --lib 50000000

# GWAS SNP与增强子重叠 + 靶基因映射
library(GenomicRanges)
library(rtracklayer)

# 读取GWAS SNP
gwas_snps <- read.csv("CAD_GWAS_significant.csv")
gr_snps <- GRanges(seqnames = gwas_snps$chr,
                   ranges = IRanges(start = gwas_snps$pos, width = 1),
                   rsid = gwas_snps$rsid, pvalue = gwas_snps$pvalue)

# 读取H3K27ac增强子峰值
enhancers <- import("H3K27ac_peaks.broadPeak")
gr_enhancers <- GRanges(seqnames = seqnames(enhancers),
                        ranges = ranges(enhancers),
                        name = enhancers$name)

# 读取互作对
interactions <- read.csv("fithic2_output/significant_interactions.csv")
gr_enh <- GRanges(seqnames = interactions$chr1,
                  ranges = IRanges(start = interactions$start1,
                                   end = interactions$end1))
gr_prom <- GRanges(seqnames = interactions$chr2,
                   ranges = IRanges(start = interactions$start2,
                                    end = interactions$end2),
                   gene = interactions$gene_name)

# SNP与增强子重叠
snp_in_enhancer <- findOverlaps(gr_snps, gr_enhancers)
cat(sprintf("落在增强子内的GWAS SNP数: %d\n", length(snp_in_enhancer)))

# 查找增强子互作的靶基因
for (i in seq_len(length(snp_in_enhancer))) {
    snp_idx <- queryHits(snp_in_enhancer)[i]
    enh_idx <- subjectHits(snp_in_enhancer)[i]
    enh_region <- gr_enhancers[enh_idx]
    target_genes <- gr_prom$gene[overlapsAny(gr_prom, enh_region)]
    if (length(target_genes) > 0) {
        cat(sprintf("SNP %s -> 增强子 -> 靶基因: %s\n",
                    gr_snps$rsid[snp_idx],
                    paste(target_genes, collapse = ", ")))
    }
}

这个综合分析案例的意义在于它阐明了非编码风险变异如何通过三维染色质互作远程调控疾病相关基因的表达。传统的GWAS分析通常将非编码SNP注释到线性距离最近的基因，但这种注释方式在许多情况下是不准确的，因为增强子可以通过染色质环与数百kb甚至数Mb之外的启动子互作。HiChIP等染色质互作技术提供了从变异到基因的物理连接证据，使得GWAS结果的解读更加准确。SORT1基因的发现就是一个经典案例——虽然风险SNP位于SORT1基因约40kb之外，但通过染色质互作分析确认了它们之间的调控关系，后续的功能实验验证了SORT1在脂质代谢中的关键作用，为降脂治疗提供了新的靶点。

范例四：组蛋白修饰图谱描绘早期胚胎发育的表观重编程（ChIP-seq）

早期胚胎发育是表观遗传重编程最剧烈的阶段之一，受精后的合子需要擦除亲本的表观遗传记忆并建立新的表观遗传程序以支持胚胎发育。一个典型的应用场景是研究小鼠卵细胞受精后不同阶段（合子基因激活前后）的组蛋白修饰重塑过程，使用少量细胞的ChIP-seq（Low-input ChIP-seq）分析H3K4me3和H3K27me3两种修饰的动态变化。H3K4me3标记活性启动子，H3K27me3标记多梳抑制区域，这两种修饰的动态重塑反映了胚胎发育过程中基因激活和沉默的程序性变化。

少量细胞ChIP-seq的技术挑战在于起始材料有限，传统的ChIP-seq通常需要数百万个细胞，而早期胚胎每个阶段可能只有数十到数百个细胞可用。Low-input ChIP-seq通过优化交联条件、使用carrier DNA和改进文库构建方法来适应低输入量。例如，使用Micrococcal Nuclease（MNase）代替超声片段化可以更温和地切割染色质，减少材料损失；使用线性扩增（如T7聚合酶扩增）代替PCR扩增可以减少扩增偏倚。这些技术优化使得在少量胚胎细胞中进行组蛋白修饰分析成为可能。

ChIP-seq数据的比对和峰值调用遵循标准流程，但需要注意不同组蛋白修饰的峰值特征差异。使用Bowtie2将reads比对到小鼠参考基因组mm10后，使用MACS2进行峰值调用。H3K4me3产生窄而尖锐的峰值，集中在转录起始位点附近，使用MACS2的默认窄峰模式调用。H3K27me3则产生宽而弥散的峰值，覆盖大范围的基因组区域，需要使用MACS2的宽峰模式（--broad参数）或SICER工具进行调用。SICER特别适合检测宽峰修饰，它通过空间聚类识别富集区域，对H3K27me3等广泛分布的修饰具有更高的灵敏度。

染色质状态分段使用ChromHMM进行，它基于多组组蛋白修饰的组合模式将基因组分割为不同的功能状态。ChromHMM使用隐马尔可夫模型（HMM）学习组蛋白修饰的组合模式，每个状态对应一种特定的染色质特征。在这个案例中，使用H3K4me3和H3K27me3两种修饰作为输入，ChromHMM可以识别出活性启动子状态（H3K4me3阳性）、多梳抑制状态（H3K27me3阳性）和二价状态（H3K4me3和H3K27me3同时阳性）等。二价染色质状态在胚胎干细胞和早期胚胎中特别常见，它将激活和抑制信号同时置于一个基因的启动子上，使基因处于"待命"状态，在后续分化中根据信号决定激活还是沉默。

动态模式分析揭示了受精后组蛋白修饰的非对称性重塑。在受精后的合子中，父源基因组上H3K4me3被广泛建立，而H3K27me3则被非对称性消除。这种非对称性反映了父本和母本基因组在受精后经历不同的表观遗传重编程过程——父本基因组在受精后经历大规模的组蛋白替换，精子的组蛋白被卵细胞的组蛋白变体替换，新组装的染色质倾向于建立活性标记H3K4me3；母本基因组则保留了卵细胞中的组蛋白修饰模式，包括部分H3K27me3标记。这种非对称的组蛋白修饰重塑为理解早期胚胎中亲本基因组的差异性调控提供了重要线索。

基因表达关联分析将组蛋白修饰变化与母源转录本降解和合子基因激活联系起来。在合子基因激活之前，胚胎发育依赖于母源沉积的mRNA和蛋白质，这些母源因子在合子基因激活后逐渐被降解。组蛋白修饰的变化与这一过程密切相关——母源转录本降解伴随其基因座上H3K4me3的减少和H3K27me3的增加，而合子新激活基因的启动子则获得H3K4me3标记。这种修饰变化与表达变化的时序对应关系验证了组蛋白修饰在基因调控中的功能作用。

差异分析使用MAnorm比较不同发育阶段之间重编程区域的峰强度变化。MAnorm最初设计用于比较两个ChIP-seq样本之间的差异结合，它通过将共享峰作为参考来归一化两个样本之间的全局信号差异，然后计算每个峰区域的log2倍数变化。在这个案例中，MAnorm比较合子基因激活前后的H3K4me3和H3K27me3信号变化，鉴定出在重编程过程中获得或丢失修饰的基因组区域。

可视化热图可以直观展示组蛋白修饰的全局重塑模式。使用deepTools的computeMatrix计算所有基因启动子区域的H3K4me3和H3K27me3信号分布，然后使用plotHeatmap生成热图。热图特别适合展示非启动子区的H3K4me3建立——在早期胚胎中，H3K4me3不仅出现在传统的启动子区域，还在大量非启动子区域（如基因间区和内含子中）被建立，这种非经典的H3K4me3分布是早期胚胎特有的表观遗传特征，可能与合子基因激活的准备有关。

# Low-input ChIP-seq数据处理流程
# 1. 比对 (Bowtie2)
bowtie2 -x /path/to/mm10_index \
    -U sample.fq.gz \
    --very-sensitive -p 8 \
    -S sample.sam 2> align_log.txt

# 2. 过滤和排序
samtools view -bS -q 30 sample.sam | samtools sort -o sample.sorted.bam
samtools index sample.sorted.bam

# 3. 去除PCR重复
picard MarkDuplicates INPUT=sample.sorted.bam OUTPUT=sample.dedup.bam \
    METRICS_FILE=sample.metrics.txt REMOVE_DUPLICATES=true

# 4. 峰值调用
# H3K4me3 (窄峰)
macs2 callpeak -t H3K4me3.dedup.bam -c input.dedup.bam \
    -f BAM -g mm -n H3K4me3_peaks -q 0.01

# H3K27me3 (宽峰)
macs2 callpeak -t H3K27me3.dedup.bam -c input.dedup.bam \
    -f BAM -g mm -n H3K27me3_peaks -q 0.01 --broad

# 或使用SICER检测H3K27me3宽峰
# SICER.sh input_dir/ H3K27me3.dedup.bam input.dedup.bam mm10 200 150 0.86 600 0.01

# ChromHMM染色质状态分段
# ChromHMM需要Java运行环境，使用命令行调用
# 1. 准备输入文件 (BAM转换为bed格式)
# java -jar ChromHMM.jar BinarizeBam mm10_chromsizes.txt bam_dir/ cell_mark.txt binarized/

# 2. 学习模型并分割基因组
# java -jar ChromHMM.jar LearnModel -p 4 binarized/ ChromHMM_output/ 4 mm10

# MAnorm差异分析
library(MAnorm)

# 读取峰值文件
peaks1 <- read.csv("zygote_H3K4me3_peaks.narrowPeak", sep = "\t", header = FALSE)
peaks2 <- read.csv("2cell_H3K4me3_peaks.narrowPeak", sep = "\t", header = FALSE)

# 运行MAnorm (需要命令行版本)
# python MAnorm.py zygote_H3K4me3_peaks.narrowPeak 2cell_H3K4me3_peaks.narrowPeak \
#     zygote_H3K4me3.reads.bed 2cell_H3K4me3.reads.bed \
#     -o MAnorm_output/

# deepTools可视化热图
# 1. 生成矩阵
computeMatrix reference-point \
    --regions gene_promoters.bed \
    --scorefiles H3K4me3_zygote.bw H3K4me3_2cell.bw \
                 H3K27me3_zygote.bw H3K27me3_2cell.bw \
    --referencePoint TSS \
    -b 3000 -a 3000 \
    --binSize 50 --missingDataAsZero \
    -o matrix_promoters.gz

# 2. 绘制热图
plotHeatmap -m matrix_promoters.gz \
    --colorMap RdBu_r \
    --whatToShow "heatmap and colorbar" \
    --zMin -2 --zMax 2 \
    -o heatmap_promoters.png

# 3. 非启动子区H3K4me3分析
computeMatrix reference-point \
    --regions non_promoter_H3K4me3.bed \
    --scorefiles H3K4me3_zygote.bw H3K4me3_2cell.bw \
    --referencePoint center \
    -b 3000 -a 3000 \
    -o matrix_non_promoter.gz

plotHeatmap -m matrix_non_promoter.gz \
    --colorMap RdBu_r \
    -o heatmap_non_promoter_H3K4me3.png

这个综合分析案例的意义在于它揭示了早期发育中表观遗传信息的重编程机制。受精后的组蛋白修饰非对称性重塑表明，父本和母本基因组在合子中经历了不同的表观遗传重建过程，这种非对称性可能对亲本印记的维持和合子基因激活的时序控制具有重要意义。非启动子区H3K4me3的建立是早期胚胎特有的表观遗传特征，它可能代表了一种"预标记"机制，为后续发育阶段的基因激活做准备。二价染色质状态的鉴定则解释了胚胎细胞如何同时保持多能性和分化潜能——通过将关键发育基因置于二价状态，细胞可以在接收到分化信号时快速激活或永久沉默这些基因。这些发现不仅深化了我们对发育生物学的理解，也为干细胞研究和再生医学提供了重要的理论基础。

范例五：表观遗传时钟预测衰老速度（EWAS + 公共甲基化数据）

表观遗传时钟是近年来衰老研究领域最重要的突破之一，它通过测量特定CpG位点的甲基化水平来预测个体的生物学年龄，为量化衰老速度和评估抗衰老干预效果提供了分子层面的工具。一个典型的应用场景是从GEO数据库下载多个队列（年龄范围20-90岁）血液样本的Illumina 450K甲基化芯片数据，构建一个预测年龄的模型并评估健康生活方式对表观遗传年龄加速（AgeAccel）的影响。这种利用公共数据的研究策略可以在不产生新测序成本的情况下验证表观遗传时钟的适用性并探索影响衰老速度的因素。

数据质量控制与归一化是甲基化芯片数据分析的首要步骤。使用minfi包读取IDAT格式的原始芯片数据，首先进行样本级别的质控，去除检测p值不合格的探针（通常阈值0.01）、性别不匹配的样本和异常样本。探针级别的质控需要去除检测失败的探针、位于性染色体上的探针（避免性别混淆）、已知含有SNP的探针以及交叉反应性探针。归一化使用ChAMP包或minfi包中的functional normalization方法，后者通过将探针信号与对照探针进行回归来去除技术变异，是目前最推荐的归一化方法。归一化后的beta值矩阵用于后续的所有分析。

筛选年龄相关CpG位点是构建表观遗传时钟的基础。使用随机森林或弹性网络回归等方法，在全基因组范围内筛选与年龄显著相关的CpG位点。弹性网络回归是构建表观遗传时钟最常用的方法，它结合了L1正则化（Lasso，促进稀疏性）和L2正则化（Ridge，处理共线性），可以在数十万个CpG位点中选择一个信息量最大的子集。通过交叉验证确定最优的正则化参数后，弹性网络回归最终选出一组CpG位点及其权重系数，这些位点的甲基化水平组合可以最优地预测年龄。

训练Horvath范本时钟使用选定的CpG位点的beta值拟合年龄。Horvath时钟的原始模型使用353个CpG位点，通过弹性网络回归从超过20000个候选位点中筛选得到。表观遗传年龄的计算公式为：表观年龄 = sum(weight_i × beta_i) + intercept，其中weight_i是第i个CpG位点的回归系数，beta_i是其甲基化水平。由于年龄与甲基化之间的关系可能是非线性的，Horvath时钟还使用了一个反变换函数将预测值映射回年龄尺度。在这个案例中，使用500个CpG位点拟合的模型在训练集上的预测误差（中位绝对误差，MAD）约为3-4年，与原始Horvath时钟的性能相当。

表观遗传年龄加速（AgeAccel）是衡量个体衰老速度偏离人群平均水平的指标。AgeAccel = 表观遗传年龄 - 实际年龄，正值表示表观遗传年龄大于实际年龄（加速衰老），负值表示延缓衰老。另一种计算方式是残差法，即将表观遗传年龄对实际年龄进行线性回归，取残差作为AgeAccel。残差法的优势在于它自动校正了表观遗传年龄与实际年龄之间的系统性偏差，而简单的差值法可能受到回归均值效应的影响。在这个案例中，使用残差法计算AgeAccel，以确保不同年龄段的AgeAccel具有可比性。

生活方式问卷数据分析将AgeAccel与生活方式因素进行关联，评估健康行为对衰老速度的影响。使用线性回归模型将AgeAccel作为因变量，运动频率、饮食质量、吸烟状态等生活方式因素作为自变量，同时校正年龄、性别、细胞比例等协变量。在这个案例中，高运动组（每周运动超过150分钟）的AgeAccel比低运动组减少约3年，表明规律运动可以延缓表观遗传衰老。吸烟者的AgeAccel比非吸烟者增加约2年，与GrimAge时钟中吸烟相关CpG位点的发现一致。这些结果为生活方式干预的抗衰老效果提供了分子层面的证据。

外部验证是评估表观遗传时钟泛化能力的关键步骤。使用一个独立的验证数据集（不参与模型训练）来评估时钟的预测精度。在验证集中，计算每个样本的表观遗传年龄并与实际年龄进行比较，计算相关系数（R²）和中位绝对误差（MAD）。一个可靠的表观遗传时钟在验证集上的R²应大于0.8，MAD应小于5年。在这个案例中，使用GEO中的另一个独立队列进行验证，验证集上的R²为0.85，MAD为3.7年，表明模型具有良好的泛化能力。

# 表观遗传时钟分析流程
library(minfi)
library(ChAMP)
library(glmnet)
library(survival)

# 1. 读取和质控450K芯片数据
# 读取IDAT文件
rgSet <- read.metharray.exp(base = "idat_dir/")
sample_names <- sampleNames(rgSet)

# 质量控制
qc <- getQC(rgSet)
plotQC(qc)

# 检测p值过滤
detP <- detectionP(rgSet)
keep <- colMeans(detP) < 0.01
rgSet <- rgSet[, keep]

# 2. 归一化 (Functional Normalization)
mSet <- preprocessFunnorm(rgSet)
beta_matrix <- getBeta(mSet)

# 去除低质量探针
# (性染色体探针、SNP探针、交叉反应探针等)
annotation <- getAnnotation(mSet)
keep_probes <- !(annotation$chr %in% c("chrX", "chrY")) &
               !grepl("rs", annotation$Name)
beta_matrix <- beta_matrix[keep_probes, ]

# 3. 弹性网络回归筛选年龄相关CpG位点
age <- pheno_data$age

# 标准化beta值
beta_t <- t(beta_matrix)
beta_t[is.na(beta_t)] <- 0

# 弹性网络回归 (alpha = 0.5为弹性网络)
cv_fit <- cv.glmnet(beta_t, age, alpha = 0.5, nfolds = 10)
best_lambda <- cv_fit$lambda.min
cat(sprintf("最优lambda: %.4f\n", best_lambda))

# 提取非零系数的CpG位点
coef <- coef(cv_fit, s = "lambda.min")
selected_cpgs <- rownames(coef)[which(coef != 0)][-1]
cat(sprintf("选中的CpG位点数: %d\n", length(selected_cpgs)))

# 4. 计算表观遗传年龄
weights <- coef[selected_cpgs, ]
intercept <- coef[1, ]
pred_age <- beta_t[, selected_cpgs] %*% weights + intercept

# 5. 计算AgeAccel (残差法)
lm_fit <- lm(pred_age ~ age)
age_accel <- residuals(lm_fit)
cat(sprintf("AgeAccel范围: %.2f 到 %.2f\n",
            min(age_accel), max(age_accel)))

# 6. 生活方式与AgeAccel关联分析
lifestyle_data <- data.frame(
    age_accel = as.numeric(age_accel),
    exercise = factor(c(rep("high", 50), rep("low", 50))),
    smoking = factor(c(rep("smoker", 30), rep("non-smoker", 70))),
    sex = factor(c(rep("M", 45), rep("F", 55))),
    age = age
)

# 线性回归: 运动与AgeAccel
fit_exercise <- lm(age_accel ~ exercise + age + sex, data = lifestyle_data)
summary(fit_exercise)

# 线性回归: 吸烟与AgeAccel
fit_smoking <- lm(age_accel ~ smoking + age + sex, data = lifestyle_data)
summary(fit_smoking)

# 7. 外部验证
# 在验证集上计算表观遗传年龄
# val_beta: 验证集的beta值矩阵
# val_age: 验证集的实际年龄
val_pred_age <- val_beta_t[, selected_cpgs] %*% weights + intercept
cor_val <- cor(val_pred_age, val_age)
mad_val <- median(abs(val_pred_age - val_age))
cat(sprintf("验证集: R² = %.3f, MAD = %.1f年\n", cor_val^2, mad_val))

这个综合分析案例的意义在于它提供了一个分子层面的"衰老时钟"，可以用于抗衰老干预效果的客观评估。传统的衰老研究依赖于临床指标和功能测试，这些指标往往变化缓慢且受多种因素影响。表观遗传时钟则提供了一个整合性的分子指标，它综合了数百个CpG位点的甲基化信息，可以更敏感地检测衰老速度的变化。当评估一种新的抗衰老干预（如热量限制、运动处方或药物干预）时，表观遗传时钟可以在相对较短的时间内（数月到数年）检测到AgeAccel的变化，而不需要等待数十年才能观察到临床终点的差异。这种快速评估的能力大大加速了抗衰老研究的进展，目前多项临床试验正在使用表观遗传时钟作为主要终点来评估候选抗衰老药物的效果。