非梗阻性無精子癥的生物信息學分析

盧振權 侯 健 梁松武 易 翔 段永剛 許 祥 羅兵鋒袁 淵 程小寶 廖蘇才 季瑞東 呂嫻媛 羅光彥

香港大學深圳醫院泌尿外科（廣東 深圳 518000）

世界衛生組織將不孕癥定義為在無保護措施的性交后一年內未能懷孕。而無精癥是男性不育的原因之一，可分為兩大類：輸精管道機械性梗阻（阻塞性無精子癥，OA）和原發性生精衰竭（非阻塞性無精子癥，NOA）。其中，約60%的無精癥患者被診斷為NOA[1]。精原干細胞的維持、分化的準備以及隨后的減數分裂和精子生成是精子發生的關鍵事件[2]。由于精子生成失敗，射精時精子缺失[3]，NOA成為最嚴重的男性不育癥。NOA可繼發于多種原因，常見原因包括性染色體異常、隱睪、易位、Y染色體微缺失和其他導致睪丸功能衰竭的遺傳疾病[1,4,5]。遺憾的是，NOA治療較為棘手，因此了解NOA的發生機制甚為重要。隨著二代測序技術的逐漸成熟，我們可以通過基因層面進一步了解疾病發生發展的過程。2019年的一項回顧性分析顯示[6]，已有78個基因與男性不育相關。但精子發生是一個復雜而又高度有序的持續供應精子的過程。關于NOA的發生機制尚不明確，為診治NOA帶來了一定的困難。幸運的是，本研究通過對比76例NOA患者睪丸組織及10例正常生精男性睪丸組織的表達譜數據，提取兩者差異基因并分析其生物學特性，共篩選出489個差異基因及30個樞紐基因，并且預測其關鍵通路，以便為臨床治療尋找潛在的生物標記物和治療靶點，為男性不育治療靶點進行探索。

資料和方法

一、數據來源

我們從GEO數據庫（https://www.ncbi.nlm.nih.gov/geo/）下載了人類無精子癥的表達譜數據相關數據集：GSE145467，GSE25518，GSE9210。并 使 用R包“sva”combat函數完成了數據集的合并及批次矯正。并利用R包“VennDiagram”對交集基因進行可視化展示combat（如圖1A-C）。

二、提取相關差異表達基因

我們使用R包“Limma”對轉錄組數據進行分析，截取值標準為|log2(FC)|＞1，FDR＜0.05。在此基礎上利用R包“Pheatmap”包和聚類分析的方法，對兩組樣本表達的基因進行聚類分析。同時我們使用R包“ggplot2”繪制火山圖。將顯著表達的基因進行可視化展示，截取值標準為|log2(FC)|＞1，FDR＜0.05。我們對兩組樣本進行同樣的分析。

三、繪制蛋白質-蛋白質互作網絡圖及Hub Gene提取

互作基因數據庫檢索工具(STRING；https://stringdb.org/cgi/input.pl)是分析PPI信息的在線工具。我們構建了一個基于差異基因（DEGs）的PPI網絡圖，（最低要求交互評分為0.7）。然后，我們使用Cytoscape軟件v3.7.1(https:/cytoscape.org/)可視化來自STRING數據庫的PPI網絡。利用Cytoscape中的CytoHubba插件，根據程度計算方法（12）對PPI網絡中的節點進行排序，其中前32個基因被認為是樞紐基因。

四、GO/KEGG富集分析及與NOA相關的Hub gene GO富集分析

分別使用“colorspace、stringi、ggplot2”這三個R包基于可視化和集成發現數據庫(DAVID；https://david.ncifcrf.gov/)來進行基因本體(Gene Ontology，GO)功能和KEGG通路富集分析，前者通過對其生物過程（BP），或細胞成分（CC）來識別確定其基因的功能，后者通過對KEGG途徑進行富集分析，將P值設定為P≤0.05作為截止值。將GO和KEGG途徑富集分析的結果作為TXT文件下載，以供后續分析。使用R軟件版本3.6.2對結果進行可視化。接著使用在線數據庫Metascape（http://metascape.org）同樣將差異基因進行GO和KEGG分析，以描繪基于不同群體之間DEG的獨特生物學意義。京都基因與基因組百科全書（KEGG）數據庫用于確定重要途徑，并與DAVID富集結果進行對比。將P值設定為P≤0.05作為截止值。

五、統計學分析

我們使用了R軟件版本3.6.2進行統計分析。采用Wilcox檢驗比較兩組。P值＜0.05被認為是有統計學意義。

結 果

一、提取無精子癥樣本中的基因集

我們從GEO數據庫下載了3個人類無精子癥的表達譜數據的相關數據集：GSE145467，GSE25518，GSE9210。分別使用在線工具對其提取每個數據集的基因轉錄組表達譜（圖1A-B）。在GSE145467中我們提取到19749個基因，GSE25518數據集中我們提取到了21654個基因，GSE9210芯片中我們提取到了16539個基因，使用在線工具Sangerbox對其進行數據合并和進行批次矯正，然后使用R包“VennDiagram”將三個基因集的交集基因數進行可視化展示（圖1C），并將這個基因集留后續分析。

圖1 NOA數據集下載及差異分析

二、提取NOA相關的DEGS

我們使用R包“Limma”包對這12782個基因進行差異基因的提取（P＜0.001,Log2FC）。我們得到了489個與無精子癥相關差異基因。其中上調的10個，下調的479個。這些結果在圖1D-E中展示。

三、GO富集分析

常見DEG的GO富集分析包括以下三個部分：生物過程（BP），分子功能（MF）和細胞成分（CC）。我們將常見的DEGs使用R包“colorspace、stringi、ggplot2”對489個無精癥相關的差異基因進行基因本體富集分析(Gene Ontology)。從結果我們可以看到，對于BP大部分基因與多細胞生物繁殖中涉及的細胞過程，細胞生物發育，細胞分化，精子發育和精子活動力的調控等方面相關。在CC中大部分基因與精子的組成，細胞外囊泡，細胞外空間，細胞外區域、細胞核，細胞質，中心體，運動性纖毛和頂體囊泡等相關。接下來對使用Cytoscape軟件基于degree值篩選得到的與NOA相關的Hub基因同樣進行GO分析：對于BP大部分基因同樣與多細胞生物繁殖中涉及的細胞過程，細胞生物發育，細胞分化，精子發育和精子活動力的調控等方面相關。在CC中大部分基因同樣基本上與精子的組成，細胞外囊泡，細胞外空間，細胞外區域、細胞核，細胞質，中心體，運動性纖毛和頂體囊泡等相關。在細胞成分（MF）方面與蛋白質的空間結構，微管運動活性，各種生物酶的活性組成相關，這些結果在圖2-A中及表1中展示。

圖2 DEGs的富集分析及相互作用的結果可視化

表1 GO信號通路

四、KEGG富集分析

通過R包基于DAVID數據庫分析常見的DEGs的相關KEGG通路。NOA相關的DEGs主要參與糖代謝，脂肪酸代謝，炭代謝、氨基酸合成、細胞周期及肌肉萎縮這些通路的調控。這些結果在圖2-B中及表1-2中展示。同樣我們也用Metascape在線工具進行DEGs的GO/KEGG富集分析：從結果我們可以看到富集的結果與DAVID富集的結果大致相同，具體結果如圖2C-E展示。

表2 KEGG信號通路

五、PPI網絡圖的構建及Hub基因的分析

根據STRING數據庫，構建了DEGS的PPI網絡，映射了328個節點和604個邊緣，如圖3-A所示。使用Cytohubba插件的Degree算法對前30個Hub基因進行了評估，如圖3-B和表3所示。

圖3 DEGs的PPI網絡

討 論

精子形成是一個復雜的過程，它包括精原細胞增殖（有絲分裂)、精母細胞減數分裂和精子細胞分化[7]。在本研究中，我們整合了GSE145467，GSE25518和GSE9210數據集，并利用生物信息學方法鑒定了489個與NOA相關的DEGs，包括NOA中10個常見的上調基因和479個常見的下調基因。GO富集分析表明，對于生物過程（biological process,BP）中的DEGs主要與細胞生物發育，細胞分化，精子發育和精子活動力的調控相關，在細胞組分（cellular component,CC）中大部分基因同樣基本上與精子的組成，細胞外囊泡，細胞外空間，細胞外區域、細胞核，細胞質，中心體，運動性纖毛和頂體囊泡等相關。在分子功能（molecular function,MF）方面與蛋白質的空間結構，微管運動活性，各種生物酶的活性組成相關，此外NOA相關的DEGs主要參與糖代謝，脂肪酸代謝，炭代謝、氨基酸合成、細胞周期及肌肉萎縮相關通路。已有研究證實上述大多數KEGG通路參與精子發生，細胞周期通路及代謝調控通路，并常常對細胞分化與分裂起著至關重要的作用[8-10]。因此，這些常見的DEGS可能與無精子癥患者的生精障礙相關。

表3 cytoscape篩選的Hub基因列表

接下來我們篩選了30個與NOA相關的Hub基因（BUB1、BIRC5、TTK、HMMR、RACGAP1、NEK2、DYN LL2、PBK、KIF15、UBB、KIF18A、OIP5、PTTG1、SPAG5、CKS2、CEP55、CDKN3、HIST1H2BA、TYMS、KPNA2、SGOL2、KIF2B、H2AFJ、DNALI1、SYCP1、SMC1B、DNAI1、DYNLRB2、FKBP6、SYCP3）。通過檢索文獻，我們發現SYCP3和SYCP1是一種蛋白質編碼基因，兩者與男性不育和無精子癥相關，其作用的相關途徑主要包括減數分裂，細胞周期和有絲分裂。它是突觸復合物的橫向絲的主要成分，形成于減數分裂前期的同源染色體之間。它是突觸復合物所必需的，同時也是減數分裂過程中正常著絲粒配對所必需，也是在卵母細胞和精母細胞發育過程中正常減數分裂染色體突觸以及正常的男性和女性生育力所必需的[11-13]。HMMR是一種透明質酸運動受體，該基因編碼的蛋白質參與細胞運動。有研究表明，它在乳腺組織中表達參與BRCA1和BRCA2復合物的形成。該基因編碼的蛋白質參與細胞運動。它在乳腺組織中表達，與BRCA1和BRCA2形成復合物[14]。一項研究表明，HMMR的下調與精子形態學、活力和數量等指標下降有關[15]。NIMA相關激酶2（Nek2）是一種參與有絲分裂調控的絲氨酸/蘇氨酸激酶，它的功能可能與細胞周期的G2/M相變有關[16]。有研究表明，Nek2在雄性小鼠減數分裂過程中的染色質濃縮起重要作用[17]。Structural Maintenance Of Chromosomes 1B(SMC1B)，是屬于減數分裂和有絲分裂期間染色單體凝聚和DNA重組所需的蛋白質家族，它主要與減數分裂和細胞周期，有絲分裂有關。有研究表明，SMC1B的自發性突變可導致小鼠黏著蛋白功能發生障礙從而導致小鼠不育[18]。TTK又被稱為TTK蛋白激酶，該基因可編碼一種雙特異性蛋白激酶，其產物參與酪氨酸、絲氨酸和蘇氨酸磷酸化的過程。TTK還與細胞增殖相關，這種蛋白質對于有絲分裂期間著絲粒的染色體對齊至關重要，并且是中心體復制所必需的，是在有絲分裂期間用于染色體的精確分離的關鍵有絲分裂檢查點蛋白[19]。有研究表明，在雄性斑馬魚的生殖組織中，TTK表達異常可能會導致雄性生殖細胞染色體異常[19]。PDZ Binding Kinase（PBK），又稱PDZ結合激酶，該基因編碼與雙特異性絲裂原激活的蛋白激酶激酶（MAPKK）家族有關的絲氨酸/蘇氨酸蛋白激酶。有證據表明，該激酶在精子形成過程中發揮重要作用，它主要參與催化有絲分裂的磷酸化酶的活性中心所必需的物質。該激酶可能同樣也參與了淋巴樣細胞的活化并支持睪丸功能，從而影響精子的形成及分化。研究表明，PBK蛋白在成年小鼠的睪丸的生精細胞中，例如精原細胞，精細胞和圓形精子中表達，它的表達異常可能會引起精子的形成異常[20]。Centrosomal Protein 55（CEP55）又稱為中心體蛋白55，在胞質分裂過程中體細胞脫落需要中心體蛋白（CEP55）。在胞質分裂過程中將PDCD6IP和TSG101募集到中體。研究表明，在生殖細胞中，CEP55在生殖細胞減數分裂中起重要的作用。CEP55過表達可能通過持續激活PI3K/Akt信號通路而抑制Foxo1的表達，從而引起雄性小鼠的不育[21]。另一項研究表明，CEP55表達的降低會抑制小鼠精原細胞增生[22]。由此可見，CEP55的表達與精子的增殖緊密相關。此外，DYNLL2、PTTG1、CDKN3、TYMS、DYNLRB2等基因的表達可能與非梗阻無精子癥相關[23]。然而，BUB1及BIRC5在精子發生方面的文獻中還沒有報道。其他11個 與NOA相關 的Hub基 因，即UBB，KIF15、KIF18A、OIP5、SPAG5、CKS2、KPNA2、SGOL2、KIF2B、H2AFJ和DNALI1也沒有報道。這些基因可能是與無精子癥發生發展潛在密切相關的新的標志物，為無精子癥患者的臨床診治提供新的思路。

總之，本研究對NOA和NS病人睪丸組織之間的DEGs進行了生物信息學整合分析，成功篩選到了13個與NOA相關的樞紐基因，即BUB1、BIRC5、UBB、KIF15、KIF18A、OIP5、SPAG5、CKS2、KPNA2、SGOL2、KIF2B、H2AFJ和DNALI1。我們還預測了與NOA相關的HUB基因的一些潛在的調控途徑，包括細胞發育，細胞分化，精子發育和精子活動力的調控等在內的精子發生和成熟過程。我們的結果有助于更詳細的理解NOA發生的分子機制，為其治療提供潛在的藥物靶點。