《International Journal of Molecular Sciences》：A Batch-Aware, Subject-Grouped Pipeline to Prioritize HIIT-Associated DNA Methylation Candidates with Cross-Cohort Evaluation in Skeletal Muscle
编辑推荐：
骨骼肌DNA甲基化与运动训练引起的分子适应密切相关，但该领域的纵向研究常受限于重复肌肉活检采样的困难及延长随访的可行性。由此产生的数据集通常样本量较小，且包含个体内重复测量及复杂的亚队列或批次结构。为应对上述局限，研究人员开发了一个甲基化候选位点发现流程，该流
骨骼肌DNA甲基化与运动训练引起的分子适应密切相关，但该领域的纵向研究常受限于重复肌肉活检采样的困难及延长随访的可行性。由此产生的数据集通常样本量较小，且包含个体内重复测量及复杂的亚队列或批次结构。为应对上述局限，研究人员开发了一个甲基化候选位点发现流程，该流程整合了批量感知统计筛选、受试者分组嵌套交叉验证及二元/多类/时间三轴一致性筛选。研究人员将该流程应用于公共骨骼肌HIIT队列GSE171140（来自52名独立参与者的195个样本；EPIC v1.0；股外侧肌），其中8周和12周时间点源自单一亚队列。该流程在每个亚队列内估计逐位点证据，使用Fisher方法合并这些估计值，实施按individual_id分组的5×3嵌套交叉验证，并要求候选CpG位点通过全部三个任务轴。在内部评估中，LR-L2实现了AUC=0.754（受试者水平bootstrap 95% CI [0.689, 0.832]；置换p≤0.002），并辅以合成信号恢复和方法比较。一致性筛选将95,477个高方差CpG位点缩减至21个L6候选位点及一个L5哨兵位点（靠近SIM1的cg06639166）。在参与者感知混合模型再分析中，21个候选位点中有20个仍与训练状态相关，经错误发现率校正后21个中有19个达到全基因组显著水平，且敏感性分析中几乎没有亚队列异质性证据。在扩展MHC（6p21.33）区域内，TNXB五-CpG簇构成了一个内部一致的区域信号，并得到25-CpG DMR（Stouffer p=3.15×10???）的支持，尽管基于仅有甲基化数据无法排除遗传、祖先相关及拷贝数变异的影响。研究人员将外部评估视为有界检验。具体而言，GSE268211提供了一个小型合作实验室、同平台、同组织一致性检验（n=10；五对配对受试者；转移AUROC 0.96–1.00，但置信区间不稳定）。GSE60655显示同组织耐力训练方向一致性（15/15可检测CpG位点方向一致，4/15原始p<0.05，0/15 BH-FDR，来自17名配对参与者的34个分析样本）。最后，GSE213363提供了一项跨组织/跨疾病压力测试（AUC=0.673），其中全部21个候选位点均可检测，但配对分析中点水平一致性有限。在配对DNAm年龄分析中（15名参与者），仅Horvath2013在多检验校正后仍显著（平均变化+1.40岁，95% CI +0.56至+2.23；BH-FDR=0.015）；其他时钟均未达到显著。本研究提供了一个可复用的分析框架及一组保守的探索性甲基化候选位点，供更大规模、基因型感知的同组织运动队列进行评估。
**研究背景与现存问题**
DNA甲基化作为哺乳动物基因组上稳定且可逆的表观遗传标记，通过CpG二核苷酸上的胞嘧啶甲基化参与细胞身份、谱系分化及基因表达的精细调控。骨骼肌作为典型的具有高度可塑性的有丝分裂后组织，持续暴露于机械负荷、能量代谢及炎症应激中，是研究运动诱导组织水平表观遗传适应的理想模型。高强度间歇训练（high-intensity interval training, HIIT）以短时、近最大强度运动与短暂恢复交替为特征，在较短训练时间内可获得与中等强度持续训练相当或更优的生理适应效益，包括心肺适能提升、骨骼肌线粒体生物发生及胰岛素敏感性改善。然而，针对HIIT的骨骼肌表观遗传研究仍较少，现有报道多集中于单次运动的急性反应或仅持续数周的短期干预，包含多个训练时间点纵向采样的研究极为稀缺。此外，小样本运动队列常包含多次招募波次或独立亚队列，其结构与性别比例及干预时长相互混杂；若直接对合并样本进行合并t检验或整体limma分析，亚队列效应可能被错误识别为运动应答；当批次标签与生物学变量完全混杂时，ComBat类预处理批次校正还可能削弱真实生物学信号。同时，仅基于单一二元对比的筛选会丢弃多时间点剂量-反应关系和受试者内时间轨迹信息；跨队列独立评估亦十分罕见，多数研究依赖单一发现队列，难以排除队列特异性假阳性信号。为应对上述局限，本研究提出了整合批量感知统计筛选、受试者分组嵌套交叉验证及二/多类/时间三轴一致性筛选的甲基化候选位点发现流程，并应用于公共HIIT骨骼肌甲基化队列GSE171140，同时纳入三个独立外部队列进行跨队列稳健性评估。
**研究方法概述**
研究人员开发并应用了一套可复用的候选位点发现流程，其核心设计包括三个可迁移组件：在合并样本前先于每个亚队列内估计位点水平证据并以Fisher方法合并（批量感知筛选）；以individual_id为分组键实施5×3受试者分组嵌套交叉验证（GroupKFold）；以及二元、多类、时间三轴一致性筛选，仅当CpG位点在全部三个任务轴同时达到阈值时才纳入最终候选列表。发现队列为GSE171140，共195个样本（男性149例、女性46例），来自52名独立参与者（男性32名、女性20名），采样部位为股外侧肌，平台为Illumina Infinium MethylationEPIC v1.0 BeadChip。三个外部评估队列包括：GSE268211（合作实验室、同平台、同组织，n=10，5对配对受试者）；GSE213363（独立实验室、多囊卵巢综合征女性患者、全血、抗阻加有氧训练，n=112，56对配对）；GSE60655（前一代450K平台、股外侧肌、三个月耐力训练，来自17名配对参与者的34个分析样本）。评估体系涵盖合成注入基准、随机种子稳定性、受试者水平bootstrap与置换检验、合并折外模型比较、与limma和minfi的方法学比较及多个独立外部队列评估。
**研究结果**
*2.1 Pipeline Evaluation*基于二元HIIT vs. 对照任务，LR-L2分类器经5×3受试者分组嵌套交叉验证在L5和L6阈值下均获得平均外折AUC=0.754；52名独立个体的1000次受试者水平bootstrap重采样提供95%置信区间[0.689, 0.832]；500次置换零分布显示观测AUC超过最大零值（0.676），单侧p≤0.002。合成注入实验显示，单轴二元注入的TPR随效应量增大而上升（d=0.5时平均TPR=0.87），而协调三轴注入在d=0.8时协调TPR仅0.018，交集集合缩减约325倍，表明三轴交集是一种保守的优先排序程序。三个随机种子产生的L5和L6候选列表完全一致（Jaccard指数=1.0），证明流程实现层面无随机性。K曲线比较显示，在K=100–500区间内本研究流程优于limma和minfi，在K=21时三种方法接近（本研究0.658、limma 0.635、minfi 0.698）。
*2.2 Binary HIIT Classification Performance*五种模型系列的内部交叉验证AUC为：LR-L1=0.679±0.066、LR-L2=0.754±0.042（最优）、RF=0.727±0.040、ElasticNet=0.755±0.039、LightGBM=0.732±0.063。合并折外预测的参与者聚类bootstrap比较显示五个模型置信区间大幅重叠，仅两个逻辑回归模型显著优于LightGBM。五模型top-100重要特征的平均配对Jaccard指数为0.13，表明二元信号分布在众多CpG位点上而非集中于紧凑特征集。多类五向任务与时间四类有序任务的macro-OVR平均AUC分别为0.760±0.029与0.736±0.026。
*2.3 The Corrected 21-CpG Shortlist*流程漏斗从865,859个原始CpG开始，经预处理保留368,596个，方差预筛选缩减至95,477个，各轴L6筛选分别识别出12,978个（二元）、705个（多类）和30个（时间）CpG，三轴交集产生21个L6候选位点。应用更严格的L5阈值进一步将短列表缩减至单一CpG（cg06639166，位于SIM1 TSS上游约1kb）。21个CpG中17个位于或靠近已知基因体，包括TNXB内五个CpG构成165bp岛体区域簇，以及与JARID2、FERMT3、MYO16等基因相关的位点。在参与者感知线性混合模型再分析中（受试者随机截距，校正亚队列和性别），21个候选位点中的20个在β和M尺度上均保持与训练状态相关，19个在全基因组368,596个CpG的BH-FDR校正后仍显著。DMRcate识别出3540个差异甲基化区域（DMR），TNXB区域（chr6:32,054,561–32,056,249，25个CpG，Stouffer p=3.15×10???）及两个SIM1基因体DMR提供了区域水平佐证。
*2.4 Candidate Case Studies*TNXB五-CpG簇（chr6:32,055,370–32,055,534）全部呈低甲基化，平均Cohen's d介于?0.65至?0.72，五个位点强力共甲基化（配对Pearson r=0.92–0.96）；GSE60655中5/5可检测位点方向一致，2/5达名义显著。cg06639166在21个L6候选位点中效应量最大（Cohen's d=?0.921），位于SIM1（单基因肥胖关键基因）TSS上游约1kb，DMRcate识别出两个SIM1附近DMR；GSE60655中方向一致但未达显著（raw p=0.159）。cg00075402（PIWIL4 TSS1500）为既往未报道位点，在GSE213363配对分析中无关联（Wilcoxon p=0.18），方向反转的名义关联未通过校正。三个候选位点均需谨慎解读：TNXB簇位于多态性丰富、mQTL富集且易发生拷贝数变异的扩展MHC/RCCX区域；cg06639166在肌肉中为多梳抑制状态而非活性增强子；cg00075402显示组织与疾病背景依赖性方向差异。
*2.5 Cross-Cohort Evaluation*GSE268211五个模型的转移AUROC为0.96–1.00，但由于仅10个样本（5对配对受试者），置信区间宽或不稳定，被视为小型合作实验室一致性检验而非充分效力的外部验证。GSE213363中五折受试者分组交叉验证平均AUC=0.673±0.065，21个候选位点均可检测，仅一个CpG（cg20064462）在发现队列观察方向上低于BH-FDR阈值，TNXB五-CpG簇无变化（|Δβ|≤0.002）。GSE60655中15个可检测候选位点全部与发现队列方向一致，但无一在BH-FDR校正后显著，且配对方向背景分析（matched-background p=0.21）显示该方向一致性并不优于队列总体背景模式；面板水平的参与者内符号翻转检验p=0.042提示存在协调训练应答。
*2.6 Sensitivity*男性单独重跑（n=149）获得L5 AUC=0.769±0.065和L6 AUC=0.774±0.067，与全队列基线相当或略高；加入性别协变量的扩展模型产生L5 AUC=0.747±0.040和L6 AUC=0.745±0.037，性别项系数接近零。PERMANOVA分析显示个体身份解释最大多变量方差比例（R2=0.624），亚队列次之（R2=0.195），五类时间标签仅占较小比例（R2=0.041）。21个候选位点在仅含Inter2的纵向模型和跨三亚队列的基线至4周共同模型中分别有15/21和19/21显著，16/21在三个亚队列中方向一致，未发现候选信号反映亚队列结构的证据。
*2.7 DNAm Age Clocks*交叉比较中四个时序时钟均显示HIIT样本年龄加速较高（Δ≈+0.8至+1.2岁），DunedinPACE无变化；由于HIIT与对照样本大多来自相同参与者的不同时间点，该比较仅作描述性报告。以15名Inter2配对受试者（30个样本形成15对基线至12周配对）进行的受试者内配对分析中，仅Horvath2013在五时钟BH-FDR校正后仍显著（Δ=+1.40岁，95% CI +0.56至+2.23；BH-FDR=0.015），Horvath2018 SkinBlood和PhenoAge仅名义显著，Hannum2013不显著，DunedinPACE为零结果。
**讨论与结论总结**
讨论部分指出，本研究的主要方法学贡献在于提供了一个可复用的候选发现流程，其迁移性源于设计原则而非单一算法：批量感知筛选、受试者分组交叉验证和多轴一致性筛选均可适用于其他小样本、高维度、配对或多时间点候选发现任务。TNXB区域作为数据驱动优先排序中最连贯的区域信号，与运动相关骨骼肌细胞外基质重塑存在生物学关联，但其位于扩展MHC区域，缺少匹配基因型、祖先和拷贝数信息，未来需基因型感知研究加以区分。cg06639166/SIM1和cg00075402/PIWIL4作为产生假说的信号，其机制解读仍需匹配的转录组、小RNA测序和功能验证。DNAm年龄时钟结果最恰当的解释是骨骼肌中组织特异性训练相关甲基化应答，而非系统生物学衰老加速的证据，且与血液中报道的运动相关DNAmAge减速现象存在根本性组织差异。跨队列评估显示同组织条件下一致性最强，跨组织和跨疾病背景下迁移性减弱。局限性包括发现队列样本量有限、性别分布偏倚、8周和12周时间点仅存在于Inter2亚队列的残余混杂、外部队列的样本量与探针覆盖限制、缺乏匹配基因型/祖先/CNV信息以及bulk组织中的细胞异质性。未来工作应优先推进匹配多组学分析、长期训练与去训练队列、关键区域功能研究及参与者特征更均衡的大型独立队列验证。
研究结论：本研究开发并应用了一个面向纵向骨骼肌DNA甲基化数据的批量感知、参与者感知候选位点优先排序框架，整合统计筛选、受试者分组嵌套交叉验证和多任务一致性分析。应用于HIIT队列GSE171140后，该框架优先排序出21个L6 CpG位点，其中20个在参与者感知混合模型分析中保持与训练状态相关，19个在全部368,596个分析CpG的BH-FDR校正后仍显著。TNXB簇作为最连贯的区域信号出现，得到25-CpG DMR和同组织方向性证据支持，而SIM1和PIWIL4附近的位点提供了额外的探索性分子追踪线索。跨队列评估显示匹配组织条件下一致性最强，跨组织和跨疾病背景下迁移性减弱。总体而言，这些发现建立了一个保守且可复用的框架，用于在小样本纵向队列中优先排序运动相关甲基化候选位点，并为后续基因型感知、转录组学和功能研究提供了一组聚焦位点。
