基于加权基因共表达网络分析法筛选多囊卵巢综合征Hub基因

陈 莉,邹 铸，徐靖原，沈梦溪，李 俊,姚新生

(1.遵义医科大学 免疫学教研室，贵州 遵义 563099；2.遵义医科大学附属医院 医学检验科，贵州 遵义 563099)

多囊卵巢综合征(Polycystic ovarian syndrome，PCOS)是临床最常见的妇科内分泌疾病，也是引起育龄期妇女不孕的最常见原因，至少影响8%～13% 育龄妇女[1]。患者除了卵巢多囊样改变，长期无排卵、血清雄激素水平升高以外[2]，还常常伴发肥胖等代谢障碍[3]。多达 60% 的 PCOS 女性超重或肥胖[4]。此外，身体质量指数(BMI)是PCOS的最强相关因素，校正相关因素后，BMI高的妇女每增加一个BMI点，PCOS的风险增加9%[5]，妊娠率下降4%[6]。然而，目前对肥胖型PCOS的形成原因知之甚少，探究肥胖与PCOS之间的关系，将有助于解析PCOS的发病机制。

近年来，随着高通量测序技术飞速发展，应用基因表达谱技术可以筛选出与疾病分子机制相关的显著差异表达基因(DEG)，并从中筛选出可通过不同信号通路、生物过程或分子机制参与疾病发生的基因。然而，PCOS的发生发展受多基因调控，单个基因具有明显的局限性。WGCNA可鉴别基因表达和临床特征之间复杂的相关性，其根据样品基因表达谱数据对基因进行聚类，构建共表达模块，将模块与临床表型相关联，发掘与临床性状高度相关的重要模块，并筛选出关键基因[7]。WGCNA已广泛应用于基因表达与表型的生物学关系的探索中，例如癌症、小鼠遗传学的研究，疾病的关键通路和中枢基因的鉴定[8]，以及生物标志物的确定[9]。

本研究应用WGCNA方法，识别与肥胖型PCOS患者相关的共表达模块，对关键模块进行GO和KEGG富集分析，并构建关键富集蛋白互作网络，筛选出肥胖型PCOS发生发展过程中的关键基因，为肥胖型PCOS发病机制的研究提供理论支持。

1 资料与方法

1.1 数据获取与处理 用于本研究的数据集(GSE5090)来自GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)，该数据集基于GPL96平台，一共包含15份在减肥手术时的脂肪组织样本。其中有8名肥胖的PCOS患者和7名肥胖的非PCOS患者。通过提取样本信息、构建基因表达矩阵、去除低质量数据和PCA分析，删除一个偏差样本，剩余 14个样本用于后续分析。

1.2 加权基因共表达网络的构建 以标准化后的基因表达矩阵，共14个转录组数据作为输入数据，用R语言genefiliter软件包筛选出各个样本间表达量变异较大的前50%的基因作为备选基因，并将其导入WGCNA 软件包进行后续分析。使用pickSoftThreshold函数计算共表达网络的软阈值，使用powerEstimate函数选择出最优软阈值，此时可得到拟合指数 R2>0.87的近似无尺度网络。使用blockwiseModules自动网络构建函数构建网络，最小模块基因数为30，相似模块的合并阈值设置为0.38(mergeCutHeight=0.38)，其他参数取默认值。

1.3 PCOS特异性模块的鉴定 对每个共表达模块的所有基因进行主成分分析(Principle component analysis,PCA)，将主成分1(PC1)称为此模块的特征向量(Module eigengene,ME)，为了筛选PCOS相关特异性模块，计算各个模块的 ME 值与不同性状之间的相关系数r和P值。r>0 表示正相关，r<0 表示负相关。

1.4 GO和KEGG功能富集分析 对Hub基因集进行功能富集分析，使用GO rest API和KEGG rest API(https://www.kegg.jp/kegg/rest/keggapi.html)获取最新的GO Terms和KEGG Pathway的注解信息作为背景，将候选基因映射到背景集合中，以获得基因的富集结果。

1.5 PCOS特异性模块核心基因的鉴定及基因互作网络构建 将挑选出的Hub基因的表达量与模块 ME 值进行关联，得到模块关系值(Module membership，MM)，其可以体现基因在模块中的重要性，MM 值越高越重要。选取PCOS正相关特异性模块中 MM 值排名前5的基因为模块中的核心基因，并利用Cytoscape软件展示核心基因互作网络。

2 结果

2.1 PCOS加权共表达网络的构建 在对所有15组样本进行聚类时发现一组离群样本，遂去除离群样本再次聚类，数据分布符合要求。将剩余14组取样本按基因表达差异大小从大到小排序，对表达量变异较大的前50%的共计13 515个基因进行加权基因共表达网络构建。当Power值为6时，无尺度网络拟合指数 R2>0.87，平均连通性趋于0(见图1A)，提示该Power值可以获得符合要求的无尺度网络。采用动态剪切法划分模块，合并相似度大于 75%的模块，最终构建了13个共表达模块(见图1B)。各模块含基因数目如图1C所示，其中Darkorange 2模块所含基因数最多，为2 270个，lightpink 4模块中基因数目最少为33个。其他模块的基因数介于33～2 270个之间。Grey(灰色模块)中的基因是一组没有分配到其他模块中的基因集，故该模块被排除在进一步分析之外。通过模块之间基因表达的相关性，发现Ivory模块和Green模块以及Darkorange 2模块和Darkturquoise模块之间相关性较高(见图1D)。

A：软阈值的选择(a：不同软阈值对应的无尺度网络拟合指数R2；b：不同软阈值对应的平均连通性)；B ：基因聚类树及模块划分； C：模块中基因数目分布；D：模块间的相关性分析。

2.2 PCOS相关特异性模块的鉴定及Hub基因筛选 将模块与性状进行相关性分析(见图2A)，挑选出与肥胖型PCOS性状相关度相对较高的4个模块，包括2个正相关模块Darkmagenta和Ivory模块；2个负相关模块Lightsteelblue 1和Darkturquiose模块。利用Limma差异基因表达分析(|log2FC|>0.5及P<0.05)和基因的模块关系值 MM(大于 0.90 )，从这4个模块中一共筛选到189个Hub基因，其中正相关模块中筛选到13个，负相关模块中筛选到176个。在肥胖型PCOS患者中表达量显著上调的有179个，表达量显著下调的有10个，MM-R值为正即为上调，为负即为下调(见表1)。

表1 富集到IL-6和TNF-e信号通路与 PCOS相关的核心基因名称及其功能注释

2.3Hub基因富集分析 对筛选到的关键基因进行GO和KEGG富集分析。KEGG通路注释分析结果显示，矫正后P<0.05的一共涉及21条相关信息通路，包括TNF 信号通路、细胞因子与细胞因子受体相互作用、NF-路、信号通路和IL-17 信号通路等，图2B显示最显著富集的前10条相关通路信息。图2C中GO通路注释分析结果显示，关键基因主要富集在免疫相关进程，包括对刺激反应的调节、免疫系统的过程、对外界刺激的反应、细胞的激活和髓系白细胞介导免疫等。

A：模块与性状的相关性热图； B：PCOS相关特异性模块 KEGG富集分析； C： PCOS特异性模块GO富集分析。

2.4 核心基因的挖掘及互作网络的构建 根据富集分析结果，关键基因主要参与免疫相关信号通路，且已有文献报道PCOS的发生发展与本文富集到的TNF-α信号通路和IL-17信号通路相关。显著富集到IL-17信号通路上的差异表达基因包括PTGS2、IL6、CXCL8、CXCL2、IL1B、S100A9、S100A8和CCL2；显著富集到TNF-α信号通路上的差异表达基因包括PTGS2、IL6、ICAM1、LIF、CXCL2、IL1B、TNFRSF1B、IRF1、CCL2(见表1)。因此将富集到这2个通路中的共计12个基因在String网站下用Cytoscape构建可视化互作网络(见图3)，且这些基因在肥胖PCOS患者中均显著上调。

A：IL-17信号通路的基因互作网络； B：TNF-α信号通路的基因互作网络。

3 讨论

减轻体重可以作为 PCOS女性的治疗目标。有研究显示：相较于非肥胖型PCOS患者，肥胖型PCOS患者往往表现出更为严重的临床表型，包括更严重的月经不调、不孕、高雄激素血症、胰岛素抵抗(IR)以及代谢综合征[10]。本研究通过对肥胖型PCOS患者和肥胖型非PCOS患者进行WGCNA分析，发现关键基因主要富集在免疫相关信号通路，且这些基因大多为炎症介质相关。

PCOS患者通常伴发低度慢性炎症，存在促炎和抗炎细胞因子分泌及功能的不平衡[11-15]。与同年龄和 BMI 匹配的对照组相比，PCOS患者的 C 反应蛋白 (CRP)、白细胞介素18 (IL-18)、肿瘤坏死因子α (TNF-α)、白细胞介素6 (IL-6)、白细胞总数 (WBC)、单核细胞水平均升高，趋化蛋白-1 (MCP-1) 和巨噬细胞炎症蛋白-1α (MIP-1α)也高于健康人群[12]。本研究分析发现 TNF-α信号通路和IL-17信号通路相关基因显著富集于肥胖型 PCOS 组。TNF-α信号通路和IL-17信号通路都与机体多种炎症反应相关。TNF-α属于TNF家族，可以激活Caspase蛋白酶，JNK和NF-κB 3条信号通路。TNF-α激活 NF-κB后，可以介导大量参与细胞存活和炎症反应的蛋白质的转录[16]，参与机体各种炎症反应[17]，此外，PCOS的重要临床表现为不孕，而TNF-α在卵泡膜细胞、卵母细胞和颗粒细胞中均表达，提示本研究发现的与TNF-α通路相关的基因也可能参与肥胖型PCOS患者的不孕。IL-17信号通路参与促进和延长炎症反应[18]。低度慢性炎症使得PCOS患者卵巢出现病增生性理性改变，主要表现为双侧卵巢增大，白膜增厚硬化，皮质纤维化，卵巢间质增生明显，这可能也是肥胖型PCOS患者发生排卵障碍的原因之一。因此，本研究富集到的IL-17信号通路相关基因可能直接参与PCOS患者的低度慢性炎症以及影响卵泡发育。

PCOS患者常伴有IR，而肥胖患者高水平的炎症因子可加剧IR[19]。过多的炎症因子打破体内活性氧自由基(ROS)的平衡，导致ROS增多，促进胰岛β细胞内一氧化氮生成，加速β细胞凋亡；ROS还可通过加速脂肪分解导致血中游离脂肪酸水平升高，引起PCOS患者IR；此外，高水平的ROS还可抑制胰岛素信号传导以及胰岛素介导的葡萄糖转运。提示炎症因子通过破坏机体内的ROS平衡产生IR[20]。除此以外，TNF-α还可以抑制肌肉和脂肪细胞中胰岛素受体和胰岛素受体底物1的酪氨酸磷酸化[21]，而导致IR甚至糖尿病的发生。因此，本研究发现的炎症相关基因可能参与肥胖型PCOS患者的IR。

促炎细胞因子还影响卵巢功能。PCOS女性比正常女性拥有更多卵泡，但它们没有成熟为优势卵泡，并且滤泡周围颗粒细胞层闭锁、退化或凋亡[22]。TNF-熟导致NF-κB被激活以介导抗凋亡因子的转录，实现其调节细胞凋亡的功能[16]。此外，本研究发现的核心基因白细胞介素6 (IL-6) 和白细胞介素(IL-1B)，是主要的促炎症因子。PCOS与许多其他慢性炎症性疾病具有相似的特性，除经典的炎症介质以外，患者的白介素6 (IL-6) 也被发现升高，其参与调控多种生殖过程，包括月经、排卵、着床和分娩[23]。据报道，IL-6分泌导致JAK2/STAT3激活，介导雌激素诱导的HIF-1α和血管内皮生长因子 (VEGF)在猪卵巢颗粒细胞中的表达[24]。VEGF在卵巢，尤其是在卵泡发育和排卵过程中的血管生成中起着重要作用。IL-1B升高可能会破坏卵泡闭锁，削弱细胞凋亡能力，抑制卵母细胞的成熟，最终导致小卵泡的发生；炎症因子的异常变化改变了其原有与AMH水平的负相关关系，从而削弱糖脂的代谢，促进胰岛素抵抗，破坏女性正常的排卵和受精系统，导致以月经稀少和异常排卵为特征的多囊卵巢综合征[25]。因此，本研究发现的肥胖型PCOS患者中的关键基因可能还参与到PCOS的卵泡异常发育进程。

趋化因子可调节炎症反应。本研究发现的核心基因中包括趋化因子相关的CXCL8、CCL2、CCR2和CXCL2。CXCL8是趋化因子受体 CXCR1/2 的配体，属于促炎趋化因子。采用中药系统药理数据库与分析平台(TCMSP)证实CXCL8在PCOS治疗中的关键作用[26]。CCR2是CCL2的受体。在一项前瞻性研究中，发现肥胖型PCOS患者中CCL2在基质成纤维细胞和子宫内膜上皮细胞中表达均上调[27]，与本研究结果一致。CXCL2也称为巨噬细胞炎症蛋白2，在本研究中，其表达量上调。在对非肥胖POCS患者脂肪细胞表达谱的研究中发现，CXCL2的表达量相较于对照组(与PCOS组具有相同身体指数的非PCOS志愿者)下调了近20倍[28]，但肥胖女性卵母细胞中的CXCL2基因表达高于正常体重女性[29]，提示CXCL2在PCOS的肥胖表型和不孕表型中可能扮演者不同的角色。

S100-A8/S100-A9 主要由活化的粒细胞释放，可触发参与炎症的信号通路，并在细胞周期进程、细胞存活、增殖和迁移等许多细胞过程中发挥重要作用[30]。有研究者采用网络分析方法，确定了S100A8控制 PCOS 的靶向蛋白的地位[31]。另有研究发现，PCOS患者外泌体中S100-A9的浓度显着高于对照组，携带 S100 钙结合蛋白 A9 (S100-A9) 的外泌体能够激活NF-κB通路，增加炎症并破坏类固醇生成，这可能与 PCOS 的发生有关。在 PCOS 过程中，富含 S100-A9 的外泌体可能由颗粒细胞、卵巢炎细胞、外周白细胞或其他细胞类型分泌，影响颗粒细胞与局部或远处环境之间的串扰[32]。

综上所述，本研究利用WGCNA研究肥胖型PCOS的枢纽基因和通路，发现参与炎症信号通路的基因在肥胖型PCOS脂肪组织中差异表达。由于PCOS伴随有慢性低度炎症状态，我们推测这些基因的失调是导致肥胖型PCOS发生的重要因素。