《Heredity》：A comprehensive evaluation of computational models for identifying methylation variability in DNA methylation arrays
编辑推荐：
DNA甲基化（DNAm）是一种关键的表观遗传修饰，通过调控基因表达参与发育、疾病进展及对环境因素应答等生物学过程。尽管以平均甲基化差异为焦点的差异甲基化分析已被广泛用于表观基因组关联研究（epigenome-wide association study, EW
DNA甲基化（DNAm）是一种关键的表观遗传修饰，通过调控基因表达参与发育、疾病进展及对环境因素应答等生物学过程。尽管以平均甲基化差异为焦点的差异甲基化分析已被广泛用于表观基因组关联研究（epigenome-wide association study, EWAS），甲基化变异性（variable methylation, VM）分析因其揭示多种表型结局背后动态分子过程的能力而日益受到关注。VM分析聚焦于甲基化模式的变异性，在检测可能预示疾病早期阶段或复杂机制的细微表观遗传改变方面尤为重要，例如在癌症中所观察到的情形，这些改变反映了基因调控的转变，而这种转变无法通过平均甲基化水平加以捕获。用于VM分析的计算工具采用多种策略，大致可归为基于特征选择的方法、联合位置-尺度检验（joint location-scale test）方法与基于方差的方法，这使得工具选择颇具挑战性。在本研究中，研究人员对八种可用于VM分析的工具进行了综合评估，以期为研究人员选择最合适的工具提供经验证据与指导。通过广泛的模拟研究，研究人员基于多种模拟情景与评价指标评估了这些工具的性能。研究发现表明，除DiffVar在特定情景下外，大多数工具均能较好地控制错误发现率（false discovery rate, FDR）。在较小样本量与中等效应量设置下，所有工具的统计效能均有限，且该效能随样本量与效应量的增加而变化。在所有模拟情景与评价指标中，与其它模型相比，JLSsc在识别DNAm微阵列数据中的VM方面展现出更优的性能。
**甲基化变异性分析的统计模型**
本文方法部分综述了可用于甲基化变异性（variable methylation, VM）分析的多种统计方法，包括经典方法Bartlett检验（Bartlett's test）、Levene检验（Levene's test）、Brown-Forsythe检验（Brown-Forsythe test）与变离差广义线性模型（DGLM），以及针对DNAm数据开发的DiffVar、EVORA、iEVORA、pETM、JLSp与JLSsc等模型。Bartlett检验为参数化方法，假定各组数据服从正态分布，基于方差的对数变换构造服从卡方分布的检验统计量。Levene检验通过计算各观测相对其组均值的绝对偏差并进行方差分析，对正态性偏离的敏感性较低。Brown-Forsythe检验以中位数或截尾均值替代均值，在数据偏离正态假设时更为稳健，但无法直接调整协变量。DGLM在广义线性模型框架下同时以两个线性预测子分别对均值与离差（方差）参数进行建模。DiffVar借鉴Levene检验的思想，采用经验贝叶斯调节t统计量以控制跨数千个CpG位点检验时的假阳性。EVORA利用Bartlett检验识别异常甲基化位点，并结合癌症异常谱分析（COPA）计算风险评分，但因其仅输出排序列表而未提供明确的显著性阈值，未被纳入本研究的基准比较；其改进版本iEVORA则被纳入评估，该模型对Bartlett检验筛选出的位点进行再评估，赋予平均甲基化差异较大的位点更高权重。pETM通过基于网络的惩罚化指数倾斜模型联合建模均值与方差，利用网络正则化考量相邻CpG位点间的相关性，并推荐以大于60%的选择概率筛选位点。JLSp采用Fisher方法组合位置检验与尺度检验的p值，并允许在Breusch-Pagan检验与Brown-Forsythe检验之间进行选择。JLSsc为联合位置与尺度得分检验，同时考虑位置与尺度统计量间的相关性，支持连续与分类变量及协变量调整，并可选择四种方差检验方法。
**异常值检测**
研究采用Cook距离量化每个观测对拟合模型的影响，将Cook距离超过均值三倍的观测判定为潜在异常值，并指出现有VM检测模型大多未能同时妥善处理方差差异与异常值问题。
**模拟研究**
模拟研究采用两种框架。基于beta分布的框架利用真实数据集（EOPET 450K对照样本）的经验均值和方差估计每个CpG位点的α与β参数，通过将经验方差乘以预设方差效应量（σ
e=3和5）构造差异方差位点，并在部分无方差效应位点中按0%或10%的样本比例引入服从均匀分布U(q,1)的异常值。基于正态分布的框架参照Phipson和Oshlack的层级模型，从缩放逆卡方分布抽样方差并据此生成M值，设定方差效应量为1.5和2.5；异常值引入方式为在部分位点中随机选择一个样本的M值替换为全部模拟CpG的最大值。此外还开展了基因组背景分层、AR(1)相关结构及细胞组成模拟等补充模拟。每个模拟情景重复100次。
**计算资源比较**
所有分析均在印第安纳大学Quartz集群的单一节点上完成，记录各模型的总运行时间与峰值内存使用量；各方法均采用默认实现，无内置并行支持。
**结果**
**DNA甲基化微阵列数据中甲基化变异性存在的证据**
对Illumina 450K阵列数据应用Levene检验后，在多个CpG位点观察到疾病组与对照组间甲基化方差的显著差异，证实甲基化变异性并非随机现象，而是甲基化微阵列数据的常见特征。
**模拟设置**
模拟评估涵盖样本量40至150、均衡（50%对照/50%病例）与非均衡（70%对照/30%病例）分组配置、无异常值与含10%异常值及不同方差效应量等情景；变异甲基化位点（VMC）比例在零模型中设为0%，在效能模型中设为10%。
**模型识别变异性甲基化的I型错误率**
在无异常值的均衡设计下，除iEVORA、DiffVar与JLSp外，其余模型均能将I型错误率控制在5%的名义水平以下，其中Brown-Forsythe与JLSsc最贴近名义水平，DGLM与pETM则较为保守。引入10%异常值后，Brown-Forsythe与JLSsc的I型错误率超过5%，DiffVar与JLSp仍无法控制错误率，iEVORA则在含异常值情景下恢复了对I型错误的控制。
**模型识别变异性甲基化的FDR**
在10% VMC、方差效应为3且无异常值的条件下，所有模型均能将FDR控制在5%阈值内，其中JLSp最接近名义水平。引入异常值后，仅Bartlett检验与pETM仍能维持有效的FDR控制。在方差效应增强至5时，除pETM在非均衡分组含异常值情景外，各模型均能有效控制FDR。
**模型识别变异性甲基化的统计效能**
在均衡设计、无异常值且方差效应为3的条件下，所有模型的统计效能均随样本量增加而逐步提升。非均衡分组导致除JLSsc与JLSp外的大多数模型效能下降。含异常值时DGLM的效能出现升高，但其FDR亦相应失控，提示该增益反映的是反保守行为而非稳健性的改善。
**模型识别变异性甲基化的AuROC**
在相同条件下，pETM的AuROC约为0.5，与随机猜测相当且不随样本量改善；其余模型的AuROC随样本量增加而提高。非均衡分组导致DGLM与iEVORA的AuROC明显下降，其余模型受影响较小。
**模型识别变异性甲基化的MCC与F1分数**
MCC与F1分数的变化趋势与统计效能一致，均随样本量增加而稳步上升；非均衡设计下JLSsc与JLSp略有改善，其余模型下降；含异常值情景使所有模型的MCC与F1分数普遍下降。
**群体水平DNAm数据集的基因组范围变异性甲基化分析**
将模拟中效能达到40%以上的四种模型（DiffVar、JLSsc、JLSp与iEVORA）应用于EOPET数据集（GSE44667，早发型子痫前期胎盘，40个样本，450K芯片，预处理后430,685个探针）。结果显示JLSp鉴定VMC数量最多（65,575个，其中31,681个为独有），JLSsc鉴定31,402个（1,041个独有），iEVORA鉴定6,506个（1,723个独有），DiffVar仅鉴定26个（1个独有）；四种模型共同识别的VMC仅3个，提示模型间一致性极低。乳腺癌数据集（GSE237036）的分析结果及计算资源比较详见补充材料。
**讨论**
综合各模拟情景与评价指标，JLSsc在灵敏度与错误控制之间取得了最佳平衡，Brown-Forsythe在均衡设计且无异常值条件下表现良好，DiffVar常无法充分控制FDR，pETM则在多数条件下分类表现接近随机。真实数据分析中模型间极低的重叠性凸显了当前VM分析缺乏统一统计框架的困境，亟需标准化报告格式，涵盖效应量、标准误、检验统计量、置信区间及多重检验校正后p值等关键量，以支持元分析与跨研究验证。研究还强调样本量与均衡设计对方差层面信号检测至关重要，并建议未来发展能够联合建模位点与区域水平均值及方差差异的统一分析流程，以促进VM分析在表观基因组研究中的可靠应用。
