基于盆腔器官脱垂相关衰老基因的GEO 数据库和LASSO 回归算法的生物信息学分析

宁敏琦, 何 勇, 李秉枢, 黄国涛, 左晓虎, 赵芷晗, 韩武岳, 洪 莉

（武汉大学人民医院妇产科 湖北省盆底疾病临床医学研究中心，湖北 武汉 430060）

盆腔器官脱垂（pelvic organ prolapse， POP）是由于盆底肌肉和筋膜组织异常造成盆腔器官下移而引发的器官位置异常及功能障碍，表现为阴道口肿物脱出，可伴排尿、排便和性功能障碍［1］。盆底组织肌源性损伤、神经源性损伤、结缔组织变化和家族遗传因素均可导致POP［2］，衰老是重要的危险因素之一。DUARTE 等［3］研究显示：无论盆底脱垂状态如何，骨盆支撑功能均会随着年龄增长而降低。既往研究［4-5］表明：生物信息学分析在探究POP 相关差异表达基因（differenntly expressed genes，DEGs）的功能和挖掘免疫系统作用方面具有重要意义，但基于POP 结合衰老筛选关键基因的研究较少。本研究旨在应用基因表达汇编（Gene Expression Omnibus， GEO） 数据库、Aging Atlas 数据库、CellAge 数据库和人类衰老基因组资源（Human Ageing Genomic Resources，HAGR） 数据库筛选出有诊断价值的POP 相关衰老DEGs，并结合LASSO 回归算法筛选关键基因，探讨其在POP 发生发展过程中的潜在机制，为POP 的诊断和治疗提供参考。

1 资料与方法

1.1 数据收集和标准化①POP 相关基因：以“pelvic organ prolapse”为检索词，从GEO 数据库（https：//www. ncbi. nlm. nih. gov/geo/） 下载已公布的人类POP 数据集GSE53868 和GSE151188。其中GSE53868 数据集共24 例样本，对照组12 例，POP 组12 例；GSE151188 数据集共15 例样本，对照组2 例，POP 组13 例。②衰老相关基因：分别从Aging Atlas 数据库（https：//ngdc.cncb.ac.cn/aging/index） 、 CellAge 数 据 库 （https：//genomics.senescence.info/cells/）和HAGR 数据库（https：//genomics.senescence.info/）筛选出502、279 和307 个衰老相关基因，将3 组基因合并后去重得724 个衰老相关基因。

1.2 数据处理和筛选DEGs利用GEO2R 分析工具（https：//www.ncbi.nlm.nih.gov/geo/geo2r/）结合R 4.2.1 软件对GSE53868 和GSE151188 数据集进行分析。其中，GSE151188 数据库通过对数转换［log2（原数值+1）］减小数据间差异，采用“hugene10srrranscriptcluster. db”工具包和人类基因组注释包“org.Hs.eg.db”转换基因名后去除表达量低于 50% 的基因。 将 GSE151188 和GSE53868 数据集标准化后合并分组得到14 例对照组样本，25例POP组样本，共39例样本，29 864 个POP 相关基因。将724 个衰老相关基因与前者取交集提取到含624 个基因的POP 相关DEGs，采用limma 工具包［6］以|log2FC|＞0.5 且校正后P＜0.05为标准进行差异分析，筛选出29 个衰老相关POP基因，其中27 个基因表达下调，2 个基因表达上调，使用ggplot2 和pheatmap 程序包进行可视化分析，生成火山图和热图。

1.3 基因本体论（Gene Ontology，GO）功能富集分析和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）信号通路富集分析使用标准化后的数据表达矩阵进行基因集富集分析（gene set enrichment analysis，GSEA），以KEGG pathway 基因集为分类标准，置换检验的置换次数设置为1 000，当P值和假发现率（false discovery rate，FDR）均≤0.05 表示富集明显。 通过DAVID （https：//david/ncifcrf.gov/） 在线分析工具对DEGs 进行GO 功能和KEGG 信号通路富集分析，以P＜0.05 为差异有统计学意义。其中，GO 功能包括生物学过程（biological process， BP）、 细胞组成 （cellular component，CC）和分子功能（molecular function，MF）。

1.4 蛋白-蛋白互作（protein-protein interaction，PPI）网络构建和核心基因（Hub 基因）筛选使用STRING（https：//string-db.org/）在线分析工具预测并构建DEGs 的PPI 网络，交互得分阈值以0.4 作为临界标准。采用Cytoscape 3.9.1 软件进行可视化，并使用插件cytoHubba 筛选排名前10 位（Top10）的Hub 基因，筛选条件按度值（degree）由高向低排序取前10 位。

1.5 获取免疫浸润矩阵采用CIBERSORT 反卷积算法基于线性支持向量回归评估组织中22 种免疫细胞所占百分率［7］。使用R 软件CIBERSORT反卷积法对中性粒细胞、嗜酸性粒细胞、单核细胞、巨噬细胞、B 淋巴细胞、T 淋巴细胞、自然杀伤细胞和树突状细胞等22 种免疫细胞转录特征的矩阵模拟计算，计算次数设置为100 次，以P＜0.05 为差异有统计学意义。

1.6 免疫细胞与DEGs 相关性分析和免疫浸润差异分析采用R 软件计算22 种免疫细胞浸润矩阵数据和DEGs 间的相关系数（r）并可视化。

1.7 LASSO 回归算法筛选关键基因使用“glmnet”包进行LASSO 回归算法分析，进一步从29 个DEGs 中筛选出10 个关键基因。

1.8 统计学分析采用R Studio 软件制图并进行统计学分析。绘制受试者工作特征曲线（receiver operating characteristic curve，ROC），评估LASSO回归算法筛选的10 个关键基因与POP 的关系，计算ROC 曲线下面积（area under curve，AUC） 评估关键基因的诊断效能［8］。以P＜0.05 为差异有统计学意义。

2 结 果

2.1 DEGs 筛 选 结 果对 GSE53868 和GSE151188 数据集合并后的数据进行分析，与3 个衰老数据库中筛选的724 个衰老基因取交集得624 个POP 相关衰老基因，以校正后P＜0.05 和|log2FC|＞0.5 为标准筛选出29 个POP 相关衰老DEGs，其中表达上调基因2 个，表达下调基因27 个。见图1。

图1 DEGs 的火山图 (A) 和热图 (B)Fig. 1 Volcano map (A) and heat map (B) of DEGs

2.2 GO 功能富集分析和KEGG 信号通路富集分析结果以KEGG 信号通路基因集为分类标准，富集到5 条上调通路和5 条下调通路。上调通路包括糖尿病并发症中晚期糖基化终末产物-晚期糖基化终末产物受体（advanced glycation end productsreceptor for advanced glycation end products，AGERAGE）信号通路、细胞衰老、恰加斯病、核苷酸结合寡聚化结构域样受体（nucleotide binding oligomerization domain like receptors，NLR） 信号传导途径和朊病毒疾病。下调通路包括阿尔茨海默病、缺氧诱导因子1 （hypoxia inducible factor-1，HIF-1）信号通路、白细胞介素17（interleukin-17，IL-17）信号通路、卡波西肉瘤相关的疱疹病毒感染和多种疾病的神经变性通路。见图2。

图2 KEGG 信号通路基因集的GSEAFig. 2 GSEA of KEGG signaling pathway gene sets

GO 功能富集分析显示：DEGs 富集于细胞对脂多糖的反应、炎症反应、细胞增殖的负向调节和DNA 诱导的转录正向调节等BP；富集于细胞核等CC；富集于转录因子结合、激酶结合和转录调控区序列特异性DNA 结合等MF。见图3A。KEGG信号通路富集分析显示：DEGs 参与IL-17、军团病、肿瘤坏死因子（tumor necrosis factor，TNF）和核因子κB（nuclear factor-kappa B，NF-κB） 等信号通路。见图3B。

图3 DEGs 的GO 功能富集分析（A）和KEGG 信号通路富集分析（B）Fig. 3 GO functional enrichment analysis (A) and KEGG signaling pathway enrichment analysis on DEGs (B)

2.3 PPI 网络和Hub 基因采用STRING 软件绘制29 个DEGs 的PPI 网络，得到1 个包含28 个结点和127 条边的PPI 网络图。 见图4A。 使用cytoHubba 插件筛选出Top10 Hub 基因，分别为白细胞介素6（interleukin-6，IL-6）、白细胞介素1B（interleukin-1B，IL-1B）、MYC 基因、前列腺素氧化环化酶2（prostaglandin endoperoxide synthase 2，PTGS2）、NF-κB 抑制因子α（nuclear factor-kappa B inhibitory factor α，NFKBIA）、早期生长反应蛋白1 （early growth response protein 1， EGR1）、CXC 趋化因子配体1 （C-X-C motif ligand 1，CXCL1）、CCAAT 增强子结合蛋白β （CCAAT enhancer binding protein β，CEBPB）、细胞周期依赖性激酶抑制因子1A （cyclin dependent kinase inhibitor 1A，CDKN1A） 和CXC 趋化因子配体2（C-X-C motif ligand 2，CXCL2）。见图4B。

图4 DEGs 的PPI 网络图 （A）和Top10 Hub 基因 （B）Fig. 4 PPI network diagram of DEGs (A) and Top10 Hub genes (B)

图5 免疫细胞浸润百分率箱图Fig. 5 Box plot of percentages of immune cell infiltration

2.4 免疫浸润矩阵结果基于CIBERSORT 反卷积法以P＜0.05 为标准的筛选结果显示：中性粒细胞和活化的肥大细胞在POP 组患者中浸润百分率较高。见图 5。

2.5 免疫细胞和DEGs 相关性分析巨噬细胞与IL-1B 呈明显正相关关系（r＞0.6），活化的肥大细胞与超氧化物歧化酶2 （superoxide dismutase 2，SOD2）、 双调蛋白 （amphiregulin， AREG）、MYC、Jun D 原癌基因（Jun D proto-oncogene，JUND）、Snail 同源物1（Snail homolog 1，SNAI1）和金属硫蛋白1E （metallothionein 1E， MT1E）呈正相关关系（r＞0.5）。静息的肥大细胞与核纤层蛋白A/C（lamin A/C，LMNA）表达呈负相关关系（r＜-0.5）。见图6。

图6 免疫细胞与DEGs 相关性分析Fig. 6 Correlation analysis on immune cells and DEGs

2.6 LASSO 回归算法筛选关键基因使用LASSO 回归算法进一步筛选DEGs 得到10 个关键基因：AREG、胰岛素样生长因子2（insulin like growth factor 2， IGF2）、 白细胞介素2 受体β（interleukin-2 receptor beta， IL-2RB）、 JUND、LMNA、 MT1E、 同源盒基因SIX1 （homebox SIX1， SIX1）、 SNAI1、 SOD2 和锌指蛋白36（zinc finger protein 36，ZFP36）。见图7。

图7 LASSO 回归算法筛选关键基因Fig. 7 Key genes screened by LASSO regression algorithm

图8 10 个关键基因免疫浸润差异小提琴图Fig. 8 Violin plot of immune infiltration differences of 10 key genes

图9 10 个关键基因的效能验证Fig. 9 Efficacy validations of 10 key genes

2.7 关键基因的免疫浸润分析与对照组比较，POP 组患者JUND、 SNAI1、 AREG、 LMNA、SOD2、MT1E、IL-2RB、IGF2 和SIX1 免疫浸润差异有统计学意义（P＜0.05）。见图 8。上述基因与“2.5”中的结果大部分相对应，进一步说明关键基因与免疫细胞之间存在密切关联。

2.8 POP 诊断效能检验结果绘制10 个关键基因的ROC 曲线判断诊断效能， JUND、 SNAI1、AREG、LMNA 和SOD2 的AUC 分别为0.880、0.840、0.831、0.814 和0.803，均大于0.750，即5 个关键基因均有较高的检验效能。见图 9。

3 讨 论

POP 严重影响患者生活质量，中国70 岁以上人群患病率高达26.11%［9-10］。盆底结缔组织的稳定性依赖于完整且具有功能的胶原纤维支撑。研究［11］显示：POP 患者盆底结缔组织胶原蛋白和弹性蛋白含量减少，胶原蛋白亚型比例改变。炎性因子是调节胶原蛋白分泌的重要因素之一。更年期前后女性阴道萎缩，阴道壁变薄，易发生炎症，发生盆底相关障碍性疾病可能性增加［12］。研究［13］显示： POP 患者主韧带组织中白细胞介素1（interleukin-1，IL-1）、IL-6 和TNF-α 等炎性因子水平均高于正常人，表明衰老作为POP 发生发展的重要因素，炎性因子的释放和胶原蛋白的代谢在POP 发生发展中发挥了重要作用。本研究采用生物信息学方法筛选出29 个POP 相关衰老DEGs，GO 功能富集分析结果显示：DEGs 主要参与细胞对脂多糖的反应、炎症反应和细胞增殖的负向调节过程，MF 方面集中表现为调控蛋白的表达，如转录因子结合和激酶结合等；KEGG 信号通路富集分析显示：DEGs 主要富集于促炎相关通路，如IL-17、TNF-α 和NF-κB 信号通路。本研究结果提示：衰老过程中POP 的发生发展与炎症反应有密切关联。

采用cytoHubba 插件从29 个DEGs 中选出Top10 Hub 基因，其中IL-6、IL-1B、PTGS2 和NFKBIA 均为炎症相关因子。 本研究中，CIBERSORT 反卷积法分析结果显示：与对照组比较，POP 组患者中性粒细胞和活化的肥大细胞浸润相对较多；DEGs 与免疫细胞的相关性分析结果显示：活化的肥大细胞与大部分DEGs 均呈正相关关系，推测在POP 病理生理机制中，活化的肥大细胞发挥着重要作用。本研究采用LASSO 回归算法进一步筛选出10 个关键基因进行免疫浸润分析，结果显示： JUND、 SNAI1、 AREG、 LMNA、SOD2 和MT1E 基因免疫浸润差异有统计学意义。本研究中ROC 曲线AUC 结果显示： JUND、SNAI1、AREG、LMNA 和SOD2 对于POP 患者有较高的检验效能，提示上述基因在衰老过程中不仅影响免疫细胞浸润，还对POP 发生发展有重要作用。

JUND 是转录因子激活物蛋白1 （activator protein-1，AP-1） 的二聚体活化形式，AP-1 可响应各种刺激结合DNA，启动转录调节基因表达。肾小球细胞外基质的重要成分层黏蛋白和Ⅰ型胶原启动子区域均有AP-1 的结合位点。研究［14］显示：AP-1 参与细胞外基质相关蛋白表达的生物学过程。人体衰老过程中，JUND 表达下调，细胞外基质相关蛋白和Ⅰ型胶原蛋白表达和分泌降低，盆底结缔组织支撑功能受影响，从而介导POP 的发生发展。SNAI1 是一种锌指蛋白，可促进上皮-间质转化，下调上皮细胞特征性标志物，上调基质细胞特征性标志物［15］，上皮-间质转化过程可减少细胞凋亡，抑制免疫反应［16］。相关研究［17-18］表明：NF-κB 通过结合上皮-间质转化基因启动子区域上调SNAI1基因表达， SNAI1 直接激活微小 RNA-21（microRNA-21，miR-21） 的转录，间接抑制M1巨噬细胞标志物的表达，增加M2 巨噬细胞标志物的表达。根据本研究中cytoHubba 插件筛选出NFKBIA 基因推测： 衰老过程中可能表达NFKBIA，抑制NF-κB 通路，进而抑制上皮-间质转化过程，增加细胞凋亡，引起免疫细胞浸润。AREG 是表皮生长因子家族成员之一，IL-17 信号通路中的炎症和纤维化相关基因表达于上皮细胞和间充质细胞，可调节成纤维细胞和免疫细胞的增殖、凋亡和迁移［19］。AREG 由多种衰老细胞分泌，研究［20-22］表明：AREG 是感染和组织修复中免疫反应的重要介质，肌肉损伤时被诱导，在肌肉再生过程中发挥重要作用。PERUGORRIA 等［23］研究发现：AREG 可刺激细胞增殖并可有效抑制分离的成纤维细胞凋亡，敲除AREG 基因会导致平滑肌肌动蛋白和胶原蛋白的沉积明显减少。LMNA 表达产物是人类核纤层A 型蛋白，作为核纤层的主要成分，在维持核结构、基因转录和正常细胞的生理性老化中发挥着重要作用［24-25］。研究［26］显示：成熟LMNA 蛋白的减少可提高成纤维细胞对活性氧的敏感性，增加DNA 损伤。LMNA 异构体聚集激活NF-κB 触发炎症相关转录程序，增加IL-6 和TNF-α 的分泌，导致细胞炎性衰老［27］。SOD2 是超氧化物歧化酶家族的主要成员之一，其编码蛋白主要分布于线粒体，可清除线粒体产生的超氧阴离子自由基，主要作用是抗氧化［28］。SOD2 还参与细胞增殖和凋亡［29］等生理过程。研究［30］显示：SOD2 信号激活可增加人体内硫酸化糖胺聚糖和胶原蛋白的含量。在心脏老化导致心力衰竭的分子机制研究中，细胞活力的增加和SOD2 的表达明显增加伴随着衰老标志物p53 的表达抑制及氧化应激的减少［31］。结合本研究结果推测：衰老过程中，SOD2 表达降低，机体自由基清除障碍，细胞外基质成分及胶原蛋白水平降低，细胞衰老加速，形成恶性循环，从而导致结缔组织支撑功能下降，加速POP 的发展进程。

综上所述，JUND、SNAI1、AREG、LMNA和SOD2 基因通过多种途径参与POP 的发生发展。SNAI1 和LMNA 可能通过参与NF-κB 通路触发炎症相关转录程序，AREG 可能通过IL-17 信号通路影响盆底平滑肌肌动蛋白修复，JUND、AREG 和SOD2 均可能参与盆底结缔组织胶原蛋白分泌和降解的关键环节。本研究通过生物信息学的方法为衰老相关POP 发生发展的分子机制研究提供了理论依据，关键基因JUND、SNAI1、AREG、LMNA和SOD2 有可能成为预测老年人POP 发生发展的分子标志物，同时有望成为POP 新的临床干预及治疗靶点。

利益冲突声明：

所有作者声明不存在利益冲突。

作者贡献声明：

宁敏琦参与研究设计、论文构思、数据分析和论文撰写，何勇参与数据采集和数据分析，李秉枢、黄国涛、左晓虎、赵芷晗和韩武岳参与论文审阅和编辑，洪莉参与研究项目的监督和指导。