子宫内膜异位症是一种雌激素依赖性疾病，以子宫内膜基质细胞与腺细胞在子宫被覆黏膜及子宫肌层以外部位种植生长为特点[1]。虽然子宫内膜异位症是一种良性疾病，但它所导致的盆腔疼痛、月经异常、不孕等影响着10%左右的育龄女性，并给全球各个国家造成不小的经济负担[2]。目前，有关子宫内膜异位症发病原因，“经血逆流”学说是被最广泛接受的学说， 90%的育龄女性存在经血逆流，但最终只有10%左右女性发展为子宫内膜异位症[3]。经血逆流的子宫内膜组织碎片在子宫内膜异位症发展过程中可能起着重要作用。正常育龄期女性子宫内膜组织会随着体内激素水平变化而呈现周期性更换，而子宫内膜组织在生长、脱落过程中伴随着大量细胞的增殖、死亡[4]。本研究通过已有数据分析，挖掘子宫内膜异位症患者内膜组织细胞死亡途径和异常表达基因及其相关通路，并探讨子宫内膜异位症不同分期的潜在生物学标志物，为子宫内膜异位症的诊治提供参考。
资料与方法
一、基因芯片选取及差异基因筛选
基因芯片数据选取通过美国国立生物技术信息中心(NCBI)的GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)进行检索与筛选完成。检索式为：endometriosis[Title] AND "Homo sapiens"[porgn] AND ("Expression profiling by array"[Filter] AND "attribute name tissue"[Filter])。搜索时间为2021年1月3日。在搜索结果中两人背对背逐条阅读，确认入选研究。利用在线R语言分析工具GEO2R进行分析，设调整后P值(Padj)<0.05的基因为显著差异基因。
二、差异基因生物学功能及通路富集分析
将得到的差异基因去重，并删除无对应Gene symbol的基因后，利用生物学信息注释数据库DAVID(https://david.ncifcrf.gov/)对所有差异基因进行分析，并利用Functional Annotation Tool进行GO(gene ontology)生物学过程富集分析(取P<0.05为显著富集)。利用在线气泡图工具(http://www.bioinformatics.com.cn/login/)绘制通路富集气泡图。
三、子宫内膜异位症不同分期样本差异基因分析
利用在线韦恩图工具(https://bioinfogp.cnb.csic.es/tools/venny/index.html)，获取不同分期子宫内膜异位症患者内膜组织差异表达的基因。利用GEO2R获取差异基因在两组间表达的具体情况，并利用String数据库(https://string-db.org/)获得相对应蛋白关联情况；利用Uniprot数据库(https://www.uniprot.org/uniprot/)获取每个基因的详细名称信息。
四、子宫内膜异位症患者内膜组织细胞死亡相关基因分析
利用GeneCards数据库(https://www.genecards.org/)，获取铁死亡(ferroptosis)、凋亡(apoptosis)、自噬(autophagy)、程序性坏死(necroptosis)和焦亡(pyroptosis)这五种细胞死亡方式相关基因。通过在线韦恩图工具(https://bioinfogp.cnb.csic.es/tools/venny/index.html)，将筛选出的子宫内膜异位症相关基因与细胞死亡相关基因取交集，观察子宫内膜异位症患者内膜组织细胞相关死亡方式。
结 果
一、样本间比较与差异基因筛选
本研究共检索到16篇相关文献。通过数据筛选后，最终纳入研究1项(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE51981)，涉及GEO平台为GPL570， 登记号为GSE51981。该研究包含微轻度子宫内膜异位症患者内膜组织样本28例(设为A组)，中重度子宫内膜异位症患者内膜组织样本49例(设为B组)，正常对照内膜组织样本34例(设为C组)，样本间详细分布见图1A。GEO2R分析发现，A、B组间共有12 410种差异基因ID(Padj<0.05)，A、C组间共有26 695种差异基因ID(Padj<0.05)，B、C组间共有21 230种差异基因ID(Padj<0.05)。C组 vs A组与B 组vs C组交集而不包括A组 vs B组部分，代表子宫内膜异位症患者内膜组织样本与正常对照组相比，共有11 695种基因ID存在显著差异(Padj<0.05)；而这些差异基因ID很可能就是区别健康人群与子宫内膜异位症患者的潜在基因谱(图1B)。
与此同时，本研究利用GEO2R自带的R limma(plotDensities)分析了样本相关基因的表达密度情况，显示样本值分布均一性较好(图1C)；利用limma(qqt)评估了limma测试结果的质量，所有样本点组成了近似直线，提示质量可(图1D)；利用hist包绘制调整后的P值直方图，提示大部分探针检测结果Padj小于0.05，差异基因较多(图1E)；利用R limma(plotSA, vooma)检查了表达式数据的均方差关系，提示整体数据变化较小，较为稳定(图1F)；利用R(boxplot)进一步分析选定样本值的分布情况，提示不同组间基因表达值是被标准化且具有可比性的(图1G)。
组间火山图与平均差图进一步可视化差异基因数量情况(图2)，提示相比对照组，子宫内膜异位症患者子宫内膜样本基因表达存在显著异常；微轻度(组A)与中重度(组B)子宫内膜异位症患者内膜组织基因表达存在显著差异。
二、子宫内膜异位症患者内膜组织异常表达基因富集分析
本研究获取了子宫内膜异位症与对照组之间11 695种显著差异基因ID，并剔除了无对应Gene symbol以及重复的基因ID。由于差异基因数量过多，本研究进一步筛选出|log2FC|>2.5的431种具有Gene symbol基因进行GO功能富集分析发现(图3)，这些差异基因在RNA聚合酶II启动子转录的正调控(GO:0045944)、细胞增殖(GO:0008283)、细胞分裂(GO:0051301)、有丝分裂核分裂(GO:0007067)、凋亡过程正调控(GO:0043065)和细胞内信号转导(GO:0035556)等60条生物功能通路上显著富集(P<0.05)。
(A) UMAP diagram. Uniform Manifold Approximation and Projection (UMAP) is a dimension reduction technology used to visualize how samples are related. The figure shows the number of nearest neighbors in the calculation. (B) Venn diagram. It is used to compare the similarities and differences of different gene IDs between different combination groups. (C) The expression density diagram. It is used to view the distribution of the selected sample values. This diagram complements the boxplot that checks the normalization of the data before differential expression analysis (Figure 1G). (D) Moderated t-statistical quantile-quantile (q-q) plot. It plots the quantiles of data samples according to the theoretical quantiles of Student′s t distribution. This chart helps to evaluate the quality of limma test results. (E) The adjusted p-value histogram. It is used to view the distribution of p-values in the analysis. (F) Mean-variance trend chart. It is used to check the mean variance relationship of expression data after fitting the linear model. Each point represents a gene; the red line is an approximation of the mean variance trend; the blue line is a constant variance approximation. (G) Boxplot chart. It is used to view the distribution of selected sample values.
图1 三组样本分布情况
Figure 1 Distribution of three groups of samples
The volcano maps of differential genes between (A) A vs B groups, (B) C vs A groups, and (C) B vs C groups. They are used to show the relationship between statistical significance (-log10p value) and multiple of change (log2FC), which are helpful to show the differentially expressed genes. The mean difference (MD) maps between (D) A vs B groups, (E) C vs A groups, and (F) B vs C groups. They are used to show the comparison of log2FC and average log2 expression value, which are convenient to visualize the differentially expressed genes. The adjusted p value (Padj) with less than 0.05 is considered to have significant difference.
图2 两组间差异火山图与平均差图
Figure 2 Difference volcano maps and average difference maps between the two groups
图3 子宫内膜异位症相关基因GO富集分析气泡图
Figure 3 Bubble diagram of GO enrichment analysis of endometriosis related genes
三、微轻度与中重度子宫内膜异位症患者内膜组织基因表达差异分析
在图1B的基础上，本研究提取A、B两组间独有的717种显著差异基因ID，剔除了无对应Gene symbol以及重复的基因，获取了最终14种具有Gene symbol的差异基因(图4)。本研究进一步利用GEO2R工具分析了这些基因在AB两组患者内膜组织中的表达差异情况，发现A、B组间差异倍数显著(表1)，提示这14种基因表达情况可能对子宫内膜异位症分期具有潜在作用。
四、子宫内膜异位症患者内膜组织细胞死亡相关基因分析
通过GeneCards数据库，本研究获得了铁死亡(103)、凋亡(13 206)、自噬(5 786)、程序性坏死(653)和焦亡(110)五种细胞死亡相关基因(图5)。在上文的基础上，本研究获取了子宫内膜异位症疾病相关且具有Gene symbol的1 468个差异基因(|log2FC|>2)。与上述五种细胞死亡方式相关基因取交集后发现，子宫内膜异位症相关基因与凋亡基因交集数量最多，有1 058个，占子宫内膜异位症相关基因80.1%，见图5；交集基因在铁死亡、凋亡、程序性坏死、自噬及焦亡基因中分别占12.6%(13/103)，8.01%(1 058/13 206)，14.7%(96/653)，10.5%(605/5 786)和8.18%(9/110)。本研究列举了五大细胞死亡方式中各自排名前十的与子宫内膜异位症相关基因，见表2。
图4 对子宫内膜异位症疾病分期具有潜在作用基因
Figure 4 Genes with potential role in the stage of endometriosis
表1 对子宫内膜异位症疾病分期具有潜在作用基因情况表
Table 1 Genes with potential effect on the stage of endometriosis
IDGene symbolProtein nameslog2FC204415_atIFI6Interferon alpha inducible protein 62.07204669_s_atRNF24Ring finger protein 242.08204846_atCPCeruloplasmin (ferroxidase)-2.06204885_s_atMSLNMesothelin2.19205780_atBIKBCL2 interacting killer2.28210425_x_atGOLGA8BGolgin A8 family member B-2.42219730_atMED18Mediator of RNA polymerase II transcription subunit 18-like2.15228030_atRBM6RNA binding motif protein 6-2.37228272_atDNLZDNL-type zinc finger2.11228550_atRTN4RReticulon 4 receptor2.09229334_atRUFY3RUN and FYVE domain containing 3-2.01229667_s_atHOXB8Homeobox B82.49231769_atFBXO6F-box protein 62.14238762_atMTHFD2LMethylenetetrahydrofolate dehydrogenase (NADP+ dependent) 2-like2.20
图5 子宫内膜异位症患者内膜组织异常表达
基因与细胞死亡相关基因交集情况
Figure 5 Intersection of abnormal expression genes
and cell death related genes in endometrium of
patients with endometriosis
讨 论
子宫内膜异位症严重影响妇女及其家庭的生活质量，给社会带来的代价与其他慢性疾病(如2型糖尿病、类风湿性关节炎)相似[5]。目前，子宫内膜异位症治疗主要是由临床症状(不孕症或盆腔疼痛)决定的，包括外科手术、激素治疗和止痛药，但这些治疗有许多不良反应，很少能长期缓解疾病带来的痛苦[6-7]。可见子宫内膜异位症的治疗仍有很长之路要走，而发掘子宫内膜异位症发生之前的预警信号以及进展过程中的治疗靶点具有巨大潜在经济价值。本研究从在位子宫内膜组织着手，对异常表达的基因进行了分析，并挖掘了部分具有潜在研究价值的异常表达基因和通路。
表2 子宫内膜异位症与细胞死亡交集相关度前10基因列表
Table 2 List of top 10 genes of intersection between endometriosis and cell death
GenesymbolDescriptionRelevance scoreIntersection gene of endometriosis andferroptosis ATG5Autophagy Related 54.59 TFRCTransferrin Receptor2.72 SLC11A2Solute Carrier Family 11 Member 22.72 ACSL1Acyl-CoA Synthetase Long Chain Family Member 12.72 NFE2L2Nuclear Factor, Erythroid 2 Like 22.02 ITGA6Integrin Subunit Alpha 61.97 FHFumarate Hydratase1.82 HELLSHelicase, Lymphoid Specific1.82 OTUB1OTUDeubiquitinase, Ubiquitin Aldehyde Binding 11.82 HMGB1High Mobility Group Box 11.76Intersection gene of endometriosis and apoptosis FASFas Cell Surface Death Receptor65.1 TNFSF10TNF Superfamily Member 1044.4 CASP6Caspase 633.3 PAWRPro-Apoptotic WT1 Regulator31.7 TRADDTNFRSF1A Associated Via Death Domain16.1 HMGB1High Mobility Group Box 115.8 SQSTM1Sequestosome 112.5 STAT1Signal TransducerAnd Activator Of Transcription 112.2 DNM1LDynamin 1 Like12 HSPA4Heat Shock Protein Family A (Hsp70) Member 411.8Intersection gene of endometriosis andnecroptosis SPATA2Spermatogenesis Associated 27.1 TRPM7Transient Receptor Potential Cation Channel Subfamily M Member 76.98 FASFas Cell Surface Death Receptor4.7 TNFSF10TNF Superfamily Member 103.7 SQSTM1Sequestosome 13.61 HMGB1High Mobility Group Box 13.61 MAP3K7Mitogen-Activated Protein KinaseKinase Kinase 73.52 TRADDTNFRSF1A Associated Via Death Domain2.57 JAK1Janus Kinase 12.43 HSP90AB1Heat Shock Protein 90 Alpha Family Class B Member 12.43
表2(续)
GenesymbolDescriptionRelevance scoreIntersection gene of endometriosis and autophagy ATG5Autophagy Related 548.3 SQSTM1Sequestosome 131.4 ATG2BAutophagy Related 2B23.9 ULK2Unc-51 Like Autophagy Activating Kinase 222.9 HMGB1High Mobility Group Box 122.2 TFEBTranscription Factor EB16 NBR1NBR1 Autophagy Cargo Receptor13.5 SOGA1Suppressor Of Glucose, Autophagy Associated 112.7 NRBF2Nuclear Receptor Binding Factor 212.6 STX17Syntaxin 1711.7Intersection gene ofendometriosis and pyroptosis NAIPNLR Family Apoptosis Inhibitory Protein6.17 DHX9DExH-Box Helicase 95.48 DDX3XDEAD-Box Helicase 3 X-Linked3.43 HMGB1High Mobility Group Box 12.51 MALAT1Metastasis Associated Lung Adenocarcinoma Transcript 12.33 GBP1Guanylate Binding Protein 12.06 SQSTM1Sequestosome 11.99 ATF6Activating Transcription Factor 60.34 IFI16Interferon Gamma Inducible Protein 160.34
目前，广泛认可的子宫内膜异位症发病学说是“经血逆流”学说，即子宫内膜碎片随着月经血通过输卵管进入盆腹腔，但是90%的输卵管未阻塞女性存在经血逆流，其数量远大于子宫内膜异位症发病人数[8]。所以，逆流的子宫内膜细胞本身特性在子宫内膜异位症病灶形成过程中扮演了重要作用，这也是子宫内膜异位症研究的方向之一。而在子宫内膜异位症发生之前，子宫内膜细胞是否已经异常值得探讨，目前尚缺乏大规模前瞻性研究验证，但是不少病例对照研究发现子宫内膜异位症患者在位内膜细胞基因表达存在异常[9-11]。尽管子宫内膜瘤和深部浸润型子宫内膜异位症可以通过成像技术(超声或MRI)检测到，但子宫内膜异位症的最终诊断仍需要通过手术(最常见为腹腔镜)和病理检测来完成[12]。目前，子宫内膜异位症诊断仍然缺乏特异性的分子标志物，无法帮助临床工作人员对子宫内膜异位症进行有效无创诊断，从而减少患者因经济负担及手术带来的痛苦[13]。
本研究挖掘了431种|log2FC|>2.5的异常表达基因，而这些基因在RNA聚合酶II启动子转录的正调控、细胞增殖、细胞分裂、有丝分裂核分裂、凋亡过程正调控和细胞内信号转导等60条生物功能通路上显著富集。这些结果提示子宫内膜异位症患者内膜组织相比正常内膜组织存在异常，而在这些异常表达基因中很可能存在潜在生物学标志物，可用于疾病的诊断及跟踪，但是进一步的分子验证需要更加深入和广泛的实验研究验证。
此外，子宫内膜异位症分期主要以美国生育学会(AFS)和美国生殖医学学会(ASRM)制定的标准为主[14]，但是症状的严重程度与分期系统之间不存在相关性[15]。可见，仍然缺少可用于子宫内膜异位症分期参考的特异性生物标志物。此外，本研究发现微轻度与中重度子宫内膜异位症患者内膜组织中IFI6、RNF24、CP等14种基因表达存在显著差异。这些差异基因挖掘工作只是一个开始，而在下一步的研究中，将针对这些差异基因开展深入实验并结合临床标本加以验证。
在众多细胞死亡途径中，本研究发现子宫内膜异位症与凋亡途径交集基因最多。细胞凋亡是一种主动行为，是由基因控制的细胞自主的有序死亡方式，具有维持机体内环境稳定的作用[16]。在子宫内膜组织周期性脱落的过程中，必然伴随着细胞适应性和主动性死亡。但不可忽视交集基因占比较高的程序性坏死和铁死亡途径。前者会导致细胞膜破裂，释放大量炎症介质，甚至引发机体对其产生过度反应，造成局部免疫环境异常[17]；后者可能引起细胞内脂质过氧化物过度堆积，并以波形传播方式诱发临近细胞相继功能异常和死亡，产生大量炎症因子[18]。由此，本研究推测子宫内膜异位症患者子宫内膜细胞死亡方式存在异常，但具体哪种死亡方式的异常对疾病影响较大尚待验证。
综上所述，本研究利用在线数据库已有测序数据，对子宫内膜异位症患者子宫内膜组织基因表达情况进行了分析。子宫内膜异位症患者子宫内膜组织基因表达存在显著异常，涉及多条异常富集通路。微轻度与中重度子宫内膜异位症患者子宫内膜组织基因表达也存在差异。经过数据交集发现，子宫内膜异位症患者子宫内膜组织与多种细胞死亡方式均有交集，而与凋亡交集基因最多。希望这些结果在子宫内膜异位症的诊疗过程中能起到一定参考借鉴作用。
1 Zondervan KT,Becker CM,Missmer SA.Endometriosis.N Engl J Med,2020,382:1244-1256.
2 Koninckx PR,Ussia A,Adamyan L,et al.Pathogenesis of endometriosis:the genetic/epigenetic theory.Fertil Steril,2019,111:327-340.
3 Czyzyk A,Podfigurna A,Szeliga A,et al.Update on endometriosis pathogenesis.Minerva Ginecol,2017,69:447-461.
4 Keeler E,Fantasia HC,Morse BL.Interventions and Practice Implications for the Management of Endometriosis.Nurs Womens Health,2020,24:460-467.
5 Zondervan KT,Becker CM,Koga K,et al.Endometriosis.Nat Rev Dis Primers,2018,4:9.
6 Falcone T,Flyckt R.Clinical Management of Endometriosis.Obstet Gynecol,2018,131:557-571.
7 Chapron C,Marcellin L,Borghese B,et al.Rethinking mechanisms,diagnosis and management of endometriosis.Nat Rev Endocrinol,2019,15:666-682.
8 郎景和.对子宫内膜异位症认识的历史、现状与发展.中国实用妇科与产科杂志,2020,36:193-196.
9 Marquardt RM,Kim TH,Shin JH,et al.Progesterone and Estrogen Signaling in the Endometrium:What Goes Wrong in Endometriosis? Int J Mol Sci,2019,20:3822.
10 Tamaru S,Kajihara T,Mizuno Y,et al.Endometrial microRNAs and their aberrant expression patterns.Med Mol Morphol,2020,53:131-140.
11 Holzer I,Machado Weber A,et al.GRN,NOTCH3,FN1,and PINK1 expression in eutopic endometrium - potential biomarkers in the detection of endometriosis-a pilot study.J Assist Reprod Genet,2020,37:2723-2732.
12 Scioscia M,Virgilio BA,Laganà AS,et al.Differential Diagnosis of Endometriosis by Ultrasound:A Rising Challenge.Diagnostics (Basel),2020,10:848.
13 Hudson QJ,Perricos A,Wenzl R,et al.Challenges in uncovering non-invasive biomarkers of endometriosis.Exp Biol Med (Maywood),2020,245:437-447.
14 American Society for Reproductive Medicine.Revised American Society for Reproductive Medicine classification of endometriosis:1996. Fertil Steril,1997,67:817-821.
15 Vercellini P,Fedele L,Aimi G,et al.Association between endometriosis stage,lesion type,patient characteristics and severity of pelvic pain symptoms:a multivariate analysis of over 1000 patients.Hum Reprod,2007,22:266-271.
16 Kaczanowski S.Apoptosis:its origin,history,maintenance and the medical implications for cancer and aging.Phys Biol,2016,13:031001.
17 Khoury MK,Gupta K,Franco SR,et al.Necroptosis in the Pathophysiology of Disease.Am J Pathol,2020,190:272-285.
18 Stockwell BR,Friedmann Angeli JP,Bayir H,et al.D.Ferroptosis:A Regulated Cell Death Nexus Linking Metabolism,Redox Biology,and Disease.Cell,2017,171:273-285.