《Immunity, Inflammation and Disease》：The Intersection of m6A Methylation and Immune Response in PCOS: A Bioinformatics Perspective
编辑推荐：
本文从生物信息学角度，深入探讨了RNA m6A甲基化修饰在多囊卵巢综合征(PCOS)中的关键作用及其与免疫细胞浸润的关联。研究整合了多个基因表达数据集，鉴定出关键的m6A调控因子（如WTAP、METTL14、ZC3H13等），并构建了一个高准确度的诊断模型。分析揭示了m6A基因与PCOS患者免疫微环境（如T细胞、巨噬细胞等）的显著相关性，提示m6A修饰可能是连接PCOS内分泌失调与免疫异常的重要分子桥梁，为发现新的诊断生物标志物和潜在治疗靶点提供了重要线索。
1 背景
多囊卵巢综合征(PCOS)是一种复杂的女性生殖功能障碍和代谢紊乱疾病，其特征包括持续性无排卵、胰岛素抵抗(IR)以及雄激素分泌过多。全球约有10%–15%的育龄期女性受其影响。PCOS患者长期的代谢紊乱会增加其未来罹患糖尿病、心血管疾病和子宫内膜癌的风险，严重危害女性健康。因此，进一步了解PCOS的病因学对其诊断、治疗和预防至关重要。
PCOS的病因复杂，无种族差异，是遗传、环境和生活方式等多种因素共同作用的结果。PCOS患者的低孕酮、高雄激素状态可导致免疫系统异常活化，引起患者体内产生多种炎症因子和自身抗体，形成慢性炎症状态。近年来，免疫功能失调在PCOS发生中的作用逐渐受到关注。研究表明，PCOS患者的卵巢组织中比非PCOS女性有更多的淋巴细胞和巨噬细胞浸润，并且多种细胞因子（如炎症因子、粒细胞-巨噬细胞集落刺激因子、胰岛素生长因子-1）可导致排卵障碍。
目前，关于PCOS遗传病因学的研究众多，部分结果表明表观遗传学在该病的发生发展中起着一定作用。表观遗传调控是生命过程中普遍存在的调控机制，在干细胞维持、自我更新与分化、个体衰老与发育异常以及疾病发生发展中起决定性作用。
N6-甲基腺苷(m6A)甲基化是真核生物mRNA中最常见的转录后修饰，占所有修饰的80%。这种修饰在哺乳动物细胞中动态可逆，通过调控基因表达、剪接、RNA编辑和RNA稳定性，同时也控制mRNA寿命和降解，并介导竞争性内源RNA(ceRNA)的翻译。尽管m6A甲基化修饰在免疫相关的生物过程(BP)中具有重要性，但由m6A甲基化修饰介导的PCOS免疫细胞浸润特征的整合图景仍不清楚。
本研究旨在利用基因表达综合数据库(GEO)数据集的RNA测序数据，分析m6A RNA甲基化调节因子的差异表达，并进一步分析m6A RNA甲基化与PCOS免疫细胞浸润特征的相关性。
2 方法
2.1 数据与预处理
从GEO数据库获取了基因表达谱和临床元数据。研究选用了直接与人类卵巢颗粒细胞(GCs)相关的数据集GSE137684、GSE80432和GSE114419，这些数据集具有完整的临床分组信息，满足了样本一致性和数据完整性的要求。对数据集进行了批次效应校正和log2归一化处理，并使用箱线图可视化标准化后的表达分布。
2.2 m6A基因全景图构建
从文献中获取了40个m6A相关基因，包括11个写入者(writer)、26个读取者(reader)和3个擦除者(eraser)基因。在与现有表达谱比较后，保留了36个基因。使用R包pheatmap绘制了所有样本中这些基因表达的热图，并使用ggpubr绘制了两组（对照组和患者组）的箱线图。此外，使用RCircos绘制了这36个基因染色体定位的环状图。
2.3 写入者与擦除者基因间的相关性分析
通过计算皮尔逊相关系数分析任意两个基因表达水平之间的相关性。当绝对相关系数超过0.7且p值小于0.01时，认为存在显著相关性。使用R包ggplot2生成符合条件基因对相关性的散点图并拟合相关曲线，使用ggExtra在图外绘制直方图和密度曲线。
2.4 基于m6A基因的诊断模型构建
首先使用R包glmnet对所有m6A基因进行岭回归筛选，选择最优lambda值，保留系数非零的基因。随后，使用逻辑回归进行进一步的基因筛选，并使用R包forestplot将模型中的基因和系数绘制在森林图上。
基于列线图，使用R包rms实现模型临床意义的可视化。使用R包pROC绘制诊断模型的受试者工作特征(ROC)曲线并计算曲线下面积(AUC)。同时，为说明列线图的有效性，使用外部数据集GSE114419绘制了精确率-召回率(PR)曲线和决策曲线分析(DCA)曲线进行验证。
2.5 蛋白质-蛋白质相互作用网络(PPIN)构建
在STRING数据库中，使用默认置信度阈值0.4构建了PPIN。随后，使用Cytoscape软件导出网络进行进一步分析。计算每个节点的网络特征，并利用Cytohubba插件基于节点度识别枢纽节点，将度值排名前10的节点指定为枢纽节点。
基于miRNet数据库，对10个枢纽节点进行了微小RNA(miRNA)和转录因子(TF)的预测研究，预测结果导出后使用Cytoscape绘图。
2.6 样本的无监督聚类
为了解析患者间的异质性，基于枢纽基因的表达谱对样本进行无监督聚类。使用R包factoextra确定最佳聚类数，使用k-means聚类对所有患者进行无监督聚类，最终获得两个聚类。根据两组样本中10个枢纽基因的表达水平创建热图，并使用ggpubr根据样本的聚类绘制组合箱线图和小提琴图。
2.7 基因富集分析
数据预处理后，进行基因差异表达分析，以揭示两组样本间的生物学差异。将校正后p值小于0.05且绝对log2倍数变化(FC)值大于1的基因视为差异表达基因(DEGs)。分析结果通过火山图和热图展示。
使用R包clusterProfiler对所有DEGs进行基因本体(GO)功能注释和京都基因与基因组百科全书(KEGG)通路富集分析，以识别显著富集的生物过程(BP)。使用条形图和气泡图可视化富集结果。富集分析的显著性阈值设定为0.05。
从MSigDB数据库获取参考基因集，使用R包clusterProfiler中包含的基因集富集分析(GSEA)方法进行富集分析和可视化。校正后p值小于0.05被认为具有统计学意义。
2.8 免疫浸润分析
为探索两组样本间免疫细胞浸润水平的异同，使用R包GSVA执行了基于单样本基因集富集分析(ssGSEA)的方法。从文献中获取了28个免疫细胞的标记基因，并以此为背景基因集对每个样本进行ssGSEA。所有免疫细胞的浸润情况通过热图和箱线图展示。
为最大限度地保证研究结果的准确性，使用R包CIBERSORT通过另一种方法评估免疫细胞浸润水平。基于CIBERSORT提供的LM22背景基因集，计算了每个样本中22种免疫细胞的含量以反映浸润水平。结果通过堆叠柱状图和箱线图呈现。
此外，为直接检验枢纽基因与免疫细胞浸润之间的相关性，为具有显著相关性的基因-免疫细胞对生成相关散点图并拟合相关曲线，同时使用R包ggExtra在图外绘制直方图和密度曲线。
3 结果
3.1 m6A基因全景图
通过整合GSE137684和GSE80432数据集的表达谱，构建了所有样本中m6A相关基因的全景图。为解决批次效应并对原始数据进行归一化，对不同来源的数据进行了批次效应校正和对数归一化。结果显示，批次校正和归一化后所有样本的表达分布一致，提升了下游分析的准确性和稳健性。
使用热图和分组箱线图可视化正常和患者样本中所有m6A基因的表达水平。观察到三种m6A功能类型之间存在差异表达。例如，写入者基因CBLL1和擦除者基因FTO在PCOS患者中表达下调，而读取者基因NXF1表达上调。此外，获取了这些基因的染色体定位并绘制了全景图。结果显示，部分基因的染色体位置非常接近，表明它们在基因组水平上密切相关，可能也具有相似的转录组特征。
3.2 写入者与擦除者基因间的相关性
通过计算相关系数、生成散点图和拟合相关曲线，对写入者和擦除者基因进行了相关性分析。四对基因组合（METTL3-ALKBH3、METTL14-ZC3H13、WTAP-CBLL1和CBLL1-PCIF1）表现出统计学显著的相关性。所有组合均显示出强正相关，相关系数超过0.7。这些基因大多是写入者和擦除者，且写入者基因具有相似功能。因此，如CBLL1-PCIF1这样的写入者基因对可能在表达上具有协同效应。写入者和擦除者基因对（如METTL3-ALKBH3）之间的正相关可能表明两者之间存在负反馈调节机制或存在未知的关键致病机制。
3.3 基于m6A基因的诊断模型
鉴于m6A基因的生物学意义，基于其表达谱建立了PCOS的诊断模型。首先，应用岭回归筛选了36个m6A基因，根据最优lambda值保留了所有基因。随后，使用逻辑回归进行第二次筛选，保留了11个基因。然后使用森林图对基于这11个m6A基因的诊断模型进行了可视化。模型识别出影响系数绝对值最高的前五个基因：WTAP (324.34)、METTL14 (320.34)、ZC3H13 (-228.38)、PCIF1 (-218.09) 和 RBM15 (191.7)。
随后，使用列线图进行可视化以验证模型的临床意义。结果显示，大多数基因对模型的影响方式一致，突显了其准确性。接下来，结合外部数据集，使用PR曲线、ROC曲线和DCA曲线进一步验证了模型的预测性能。结果显示，在三种验证方法中，该模型均表现出优异的预测性能和稳健性。ROC和PR曲线显示，该模型的总体AUC为0.89，表明其具有出色的诊断和预测能力。
3.4 m6A基因的PPIN
为研究m6A调节基因间的相互作用，使用STRING数据库构建了PPIN，并使用Cytoscape软件进行可视化。调节相同生物过程(BP)的基因通常表现出更密切的关系，这有助于深入了解它们的功能联系。
网络中包含少量高度节点，表明与其他节点存在强相关性。此类节点通常在网络中起关键作用。基于度分析，识别出前10个枢纽基因。为了更全面地分析其调控网络和遗传背景，使用miRNet数据库预测了微小RNA(miRNA)和转录因子(TF)。这些基因似乎既有专属的也有共享的miRNA或TF，表明它们受到类似的调控并执行相似的生物学功能。
3.5 枢纽基因的无监督聚类
枢纽基因能够区分不同疾病状态的样本。基于10个枢纽基因的表达谱，使用k-means聚类方法进行无监督聚类，在患者样本中识别出两个不同的聚类。随后，通过基于主成分分析的降维方法对聚类结果进行可视化，结果显示两组样本分类清晰，聚类效果良好。
热图以及结合箱线图和小提琴图展示了两组样本聚类中枢纽基因表达水平的相似性和差异性。观察到两组样本在METTL14、HNRNPA2B1、YTHDF3、YTHDF2、YTHDC1和YTHDC2这六个基因上存在显著的表达差异，表明这六个基因可能是关键区分因子。基于这一发现，聚类结果有效且准确。
3.6 样本组间的生物学差异
为识别两组样本间的生物学差异，进行了基因差异表达分析，共鉴定出65个达到统计学显著性阈值的基因。DEGs通过火山图和热图可视化。
随后，进行了GO、KEGG和GSEA富集分析，以研究这些DEGs影响了哪些生物过程(BP)。GO富集分析显示，DEGs显著富集于“细胞外基质组织”、“血小板脱颗粒”、“巨噬细胞来源的泡沫细胞分化调控”等生物过程，以及“血小板α颗粒”等细胞成分和“碳水化合物结合”、“受体配体活性”等分子功能。KEGG通路富集分析显示，DEGs在“ECM-受体相互作用”、“趋化因子信号通路”、“谷氨酸能突触”等通路中富集。GSEA分析进一步揭示了PCOS样本中显著上调的基因集富集于“嗅觉传导”、“神经活性配体-受体相互作用”、“造血细胞谱系”、“细胞因子-细胞因子受体相互作用”等通路；而显著下调的基因集则富集于“Hippo信号通路”、“ECM-受体相互作用”、“缬氨酸、亮氨酸和异亮氨酸降解”、“类固醇生物合成”等通路。
3.7 免疫浸润分析结果
使用ssGSEA和CIBERSORT两种方法评估了PCOS患者和对照组样本间的免疫细胞浸润差异。分析揭示了PCOS患者卵巢微环境中免疫细胞组成的显著改变。
具体而言，分析显示PCOS组与对照组相比，多种免疫细胞的浸润水平存在差异。例如，某些T细胞亚群、巨噬细胞、自然杀伤细胞等的浸润程度发生了变化。相关性分析进一步表明，多个关键的m6A调控基因（如METTL3、WTAP、YTHDF2等）的表达水平与特定免疫细胞（如CD8+T细胞、调节性T细胞(Tregs)、M2型巨噬细胞等）的浸润丰度存在显著的相关性。其中，METTL3与免疫细胞浸润水平的调节关系最为显著。这些结果提示，m6A甲基化修饰可能通过调节这些枢纽基因的表达，进而影响卵巢局部免疫微环境的构成，参与PCOS的病理过程。
4 结论
N6-甲基腺苷(m6A) RNA甲基化调节因子与PCOS的发展密切相关，并可能影响患者的免疫细胞浸润。本研究通过整合生物信息学分析，系统描绘了PCOS中m6A调控因子的表达图谱，并构建了基于关键m6A基因的高性能诊断模型。研究首次深入揭示了m6A修饰与PCOS免疫微环境特征之间的广泛联系，表明m6A甲基化可能是连接PCOS内分泌代谢异常与免疫失调的关键表观遗传机制。这些发现增强了对PCOS分子相互作用的理解，并提出了用于诊断的潜在生物标志物和用于治疗干预的靶点。
