子宮內膜異位癥不同分期患者內膜組織異常表達基因及相關細胞死亡方式分析

趙一鳴 倪喆鑫 劉新敏

子宮內膜異位癥是一種雌激素依賴性疾病，以子宮內膜基質細胞與腺細胞在子宮被覆黏膜及子宮肌層以外部位種植生長為特點[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

圖5 子宮內膜異位癥患者內膜組織異常表達基因與細胞死亡相關基因交集情況Figure 5 Intersection of abnormal expression genes and cell death related genes in endometrium of patients with endometriosis

討 論

子宮內膜異位癥嚴重影響婦女及其家庭的生活質量，給社會帶來的代價與其他慢性疾病(如2型糖尿病、類風濕性關節炎)相似[5]。目前，子宮內膜異位癥治療主要是由臨床癥狀(不孕癥或盆腔疼痛)決定的，包括外科手術、激素治療和止痛藥，但這些治療有許多不良反應，很少能長期緩解疾病帶來的痛苦[6-7]?？梢娮訉m內膜異位癥的治療仍有很長之路要走，而發掘子宮內膜異位癥發生之前的預警信號以及進展過程中的治療靶點具有巨大潛在經濟價值。本研究從在位子宮內膜組織著手，對異常表達的基因進行了分析，并挖掘了部分具有潛在研究價值的異常表達基因和通路。

表2 子宮內膜異位癥與細胞死亡交集相關度前10基因列表Table 2 List of top 10 genes of intersection between endometriosis and cell death

表2(續)

目前，廣泛認可的子宮內膜異位癥發病學說是“經血逆流”學說，即子宮內膜碎片隨著月經血通過輸卵管進入盆腹腔，但是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]。由此，本研究推測子宮內膜異位癥患者子宮內膜細胞死亡方式存在異常，但具體哪種死亡方式的異常對疾病影響較大尚待驗證。

綜上所述，本研究利用在線數據庫已有測序數據，對子宮內膜異位癥患者子宮內膜組織基因表達情況進行了分析。子宮內膜異位癥患者子宮內膜組織基因表達存在顯著異常，涉及多條異常富集通路。微輕度與中重度子宮內膜異位癥患者子宮內膜組織基因表達也存在差異。經過數據交集發現，子宮內膜異位癥患者子宮內膜組織與多種細胞死亡方式均有交集，而與凋亡交集基因最多。希望這些結果在子宮內膜異位癥的診療過程中能起到一定參考借鑒作用。