前列腺間質(zhì)狀態(tài)在良性前列腺增生的發(fā)生發(fā)展中的轉(zhuǎn)錄組學分析

徐向來， 楊正青， 王佳駿， 朱延軍*

1. 復旦大學附屬中山醫(yī)院泌尿外科，上海 200032 2. 上海楊思醫(yī)院泌尿外科，上海 200126

良性前列腺增生癥(BPH)是造成老年男性膀胱出口梗阻的常見原因，其引起的下尿路癥狀(LUTS)大大降低了患者生活質(zhì)量[1]。BPH的患病率隨著年齡增長逐漸增加，50%的50歲以上男性存在組織學上的前列腺良性增生[2]。目前，BPH的治療藥物主要包括5α還原酶抑制劑、α受體阻滯劑和5型磷酸二酯酶抑制劑。然而，仍有大量患者患有難以通過藥物改善的LUTS癥狀及并發(fā)癥，出現(xiàn)反復的尿潴留[3-4]。這些患者需手術(shù)干預(yù)，例如經(jīng)尿道前列腺切除術(shù)(TURP)和激光剜除前列腺(HoLEP)[5]。

BPH的組織學改變是基質(zhì)細胞和上皮細胞的過度生長，主要發(fā)生于前列腺的移行帶。BPH的發(fā)病與雄激素信號通路的改變、間質(zhì)活化和炎癥免疫反應(yīng)[1-2]相關(guān)。雄激素信號通路激活前列腺基質(zhì)相關(guān)的信號通路，造成增生相關(guān)生長因子的表達，進一步造成間質(zhì)細胞的增殖和分化[6-7]。前列腺間質(zhì)改變在BPH發(fā)生發(fā)展中起重要作用，而局部缺血缺氧可能是這種間質(zhì)病變的誘因。BPH相關(guān)的前列腺炎癥也被認為是BPH的重要因素[8-9]。盡管已有研究團隊對前列腺間質(zhì)增生、變化在BPH中的機制進行過探索，但是尚無轉(zhuǎn)錄組水平的深入研究。

轉(zhuǎn)錄組學技術(shù)已廣泛用于人類疾病的研究。幾十年前，就有團隊使用芯片進行了BPH表達譜的初探[10-11]。近年來，基因組學、轉(zhuǎn)錄助學和表觀基因組學分析等，已被用于BPH研究[12-14]。這些研究在一定程度上揭示了BPH發(fā)生發(fā)展的生物學過程。研究人員發(fā)現(xiàn)了2種分子分型的BPH亞組，代表不同的生物學行為[12]。另外，Middleton等[13]發(fā)現(xiàn)骨形態(tài)發(fā)生蛋白5(BMP5)等基因的表達也與BPH疾病進展有關(guān)。本研究擬從轉(zhuǎn)錄組學入手，觀察BPH患者前列腺間質(zhì)組織成分變化，并分析正常前列腺和BPH患者前列腺標本的間質(zhì)相關(guān)基因差異及特征。

1 材料與方法

1.1 研究人群和數(shù)據(jù)采集 研究工作流程如圖1。本研究搜索了Gene Expression Omnibus(GEO)以獲取BPH的轉(zhuǎn)錄組數(shù)據(jù)，共收集了34個來自GEO數(shù)據(jù)集(GSE119195、GSE7307、GSE101486)的BPH樣品和13個正常前列腺樣品的轉(zhuǎn)錄組測序數(shù)據(jù)(表 1)。其中，GSE119195的BPH標本來自全膀胱切除標本中的前列腺組織，正常前列腺組織來自腦死亡患者捐贈的標本。GSE7307的BPH和正常前列腺組織，均來自腦死亡患者捐贈的標本。GSE101486的BPH標本來自TURP標本，與前列腺剜除標本。所有BPH標本均為前列腺移行帶標本，均經(jīng)過組化病理分析排除前列腺惡性腫瘤；所有正常對照標本，均經(jīng)過病理分析排除BPH及前列腺惡性腫瘤。另一個數(shù)據(jù)集(GSE132714)含18個BPH樣本(來自TURP標本)和4個正常前列腺樣本(來自T2a前列腺腫瘤標本的對側(cè)正常移行帶)用作驗證數(shù)據(jù)集。R軟件(版本4.0.2；https://www.r-project.org/)用于數(shù)據(jù)分析。sva 程序包用于去除批次效應(yīng)[15]。

表1 GEO數(shù)據(jù)庫中調(diào)取的BPH及正常組織的個數(shù)

圖1 分析流程圖

1.2 樣本微環(huán)境分析 定量基質(zhì)評分采用ESTIMATE算法[16]，每個樣本的微環(huán)境細胞組成計算采用加州大學舊金山分校開發(fā)的xCell算法[17]。

1.3 加權(quán)基因共表達網(wǎng)絡(luò)分析(WGCNA) WGCNA是利用基因表達數(shù)據(jù)構(gòu)建無標度網(wǎng)絡(luò)的系統(tǒng)生物學方法[18]。本研究通過WGCNA分析基因組和表型之間的關(guān)聯(lián)，發(fā)現(xiàn)可能與BPH發(fā)育有關(guān)的潛在生物標志物基因。將基因表達矩陣和表型數(shù)據(jù)構(gòu)建為輸入數(shù)據(jù)。使用R軟件中的“WGCNA”軟件包執(zhí)行分析。

1.4 功能富集分析 通過“Limma”包計算BPH組織和正常前列腺組織的差異表達基因(differentially expressed gene, DEG)，采用t檢驗進行統(tǒng)計分析，將P<0.001和表達倍數(shù)絕對值>2的基因視為差異表達基因[19]。使用Metascape進行基因本體論(gene onotology, GO)分析和Kyoto encyclopedia of genes and genomes(KEGG)富集分析。

1.5 蛋白質(zhì)相互作用網(wǎng)絡(luò)分析(PPI) 使用search tool for the retrieval of interacting genes/proteins (STRING)工具進行蛋白互作網(wǎng)絡(luò)分析。此外，對PPI數(shù)據(jù)，進一步使用Cytoscape 3.7.2版中的cytoHubba插件進行關(guān)鍵基因(hub gene)的篩選。使用以下4種算法計算hub基因：maximum neighborhood component (MNC)、maximal clique centrality (MCC)、maximum neighborhood component (DMNC)、degree。根據(jù)分值的高低篩選出排名前15的關(guān)鍵基因。

1.6 主成分分析(PCA) PCA是利用降維的方法，將多個變量轉(zhuǎn)化降為成為幾個綜合的變量(主成分)，每個成分之間互不相關(guān)。這些主成分能在信息不重疊的情況下，較好地反映原始變量的絕大部分信息。因此，PCA是一種能夠在最大化降低保持數(shù)據(jù)信息損失的原則下，對高維變量空間進行降維處理的線性映射方法。本研究使用PCA評估篩選的PCA基因是否能夠良好地區(qū)分BPH標本及正常前列腺標本。

1.7 LASSO回歸方法 LASSO回歸的原理是在最小二乘的基礎(chǔ)上，增加一個懲罰項來對估計系數(shù)進行壓縮。當系數(shù)在懲罰項的作用下無限趨于0時，就可選擇出對因變量影響較大的自變量和對應(yīng)的系數(shù)估計值。LASSO回歸常用于處理存在多重共線性的樣本數(shù)據(jù)。參數(shù)λ決定了LASSO回歸懲罰項的壓縮程度，λ越大則懲罰項的影響力越大。通過LASSO回歸，可使最終納入模型的變量均為與因變量顯著相關(guān)(P<0.05)且考慮變量共線性問題的自變量集合。LASSO回歸可避免參數(shù)過度擬合，適合可給出一個具備優(yōu)良性能但是自變量個數(shù)相對較少的診斷模型。

2 結(jié) 果

2.1 BPH組織和正常前列腺組織中間質(zhì)細胞占比不同 本研究從GEO數(shù)據(jù)庫下載了3個數(shù)據(jù)集，去除批次效應(yīng)，合并成為包含了34個BPH和13個正常前列腺樣品的轉(zhuǎn)錄組數(shù)據(jù)和表型數(shù)據(jù)的探索數(shù)據(jù)集(圖1)。通過ESTIMATE程序包計算了在BPH和正常前列腺組織中的基質(zhì)評分發(fā)現(xiàn)，BPH組織中間質(zhì)評分顯著高于正常前列腺組織(P=0.025，圖2)。結(jié)果表明，BPH組織中的間質(zhì)成分占比顯著多于正常前列腺組織。為了進一步研究BPH發(fā)生發(fā)展過程中的間質(zhì)變化，本研究使用xCell算法對BPH和正常前列腺組織樣本中的間質(zhì)細胞進行詳細計算(圖3A、圖3B)。目前，尚無研究團隊對BPH組織的間質(zhì)組成進行分型分析，特別是不同類型或表型的間質(zhì)細胞。通過xCell算法發(fā)現(xiàn)，BPH樣品中的成纖維細胞、周細胞顯著增加(P均<0.05)。相反，內(nèi)皮細胞、肌細胞、成骨細胞和平滑肌細胞組織減少(P均<0.05，圖3A)。本研究還發(fā)現(xiàn)在BPH組織中成纖維細胞/平滑肌細胞的比例顯著升高(P=0.001 8，圖3C)。從轉(zhuǎn)錄組數(shù)據(jù)可推斷，成纖維細胞或平滑肌細胞的升高，也是BPH組織的特征之一。

圖2 BPH組織中間質(zhì)評分顯著高于正常前列腺組織

圖3 xCell分析BPH和正常前列腺組織樣本中的間質(zhì)細胞A-B：BPH樣品中的成纖維細胞、周細胞顯著增加。相反，內(nèi)皮細胞、肌細胞、成骨細胞和平滑肌細胞組織減少；C： BPH組織中成纖維細胞/平滑肌細胞的比例顯著升高(P=0.001 8)。

2.2 WGCNA模塊分析與BPH及間質(zhì)相關(guān)的關(guān)鍵模塊 WGCNA構(gòu)建的輸入表達譜矩陣和表型矩陣包含常見的14 24個基因和34個樣品的表型數(shù)據(jù)。輸入的表型數(shù)據(jù)包括間質(zhì)評分、成纖維細胞、平滑肌細胞、內(nèi)皮細胞。根據(jù)表達譜數(shù)據(jù)，計算每個樣品中每個基因的方差。選擇標準偏差大于1.2的基因并進一步聚類所有樣品，并標注樣本特征(圖4)。所有樣品均符合WGCNA的篩選標準,未排除樣品。使用GCNA包進行構(gòu)建權(quán)重共表達網(wǎng)絡(luò)，使用分析包自動選擇的軟閾值計算得到軟閾值β=18作為軟閾值進行后續(xù)分析。確定軟閾值后，通過動態(tài)剪切樹法進行模塊初步識別并合并相似模塊，設(shè)置每個基因網(wǎng)絡(luò)模塊最少的基因數(shù)目為30。總共獲得18個模塊，其中灰色模塊是不能聚合到其他模塊中的基因的集合(圖5A)。

圖4 WGCNA聚類所有樣品，標注樣本特征，切割合并模塊A：WGCNA算法對樣品進行聚類，未發(fā)現(xiàn)需要剔除的樣本；B：以0.2相異性系數(shù)為閾值合并相關(guān)的基因模塊。

根據(jù)每個模塊的特征向量，計算這些模塊和每個輸入表型(免疫評分、成纖維細胞、平滑肌細胞、內(nèi)皮細胞)之間的相關(guān)性。對富集到的模塊進行GS分布預(yù)測。每個模塊中每種表型的GS分布可顯示表型和基因之間的相關(guān)性(圖5B、圖6)。結(jié)果提示，volet、darkgrey、lightyellow、darkorange、purple、darkolivegreen、pink、greenyellow、lightcyan1模塊和間質(zhì)評分、成纖維細胞、平滑肌細胞相關(guān)。因此，提取這9模塊所包含的基因作為進一步分析。

圖5 WGCNA合并模塊，對表型信息進行相關(guān)性分析A-B：通過動態(tài)剪切方法確定模塊后，依次計算每個模塊的特征向量，然后聚類模塊并將更近的模塊合并到新模塊中，共獲得18個模塊，其中灰色模塊是不能聚合到其他模塊中的基因的集合。

圖6 WGCNA基因畫拓撲重疊熱圖和特征向量基因A：隨機選擇400個基因畫拓撲重疊熱圖，行和列都表示單個基因，深黃色表示高度的拓撲重疊；B：計算每個模塊的特征向量基因，繪制特征向量基因臨近熱圖，代表該模塊內(nèi)基因表達的整體水平。

2.3 前列腺間質(zhì)增生相關(guān)的關(guān)鍵基因 使用“l(fā)imma”程序包，對BPH和正常前列腺樣本的表達譜數(shù)據(jù)進行了差異表達分析。|logFC|≥1且P<0.05為閾值，共篩選出98個上調(diào)基因和204個下調(diào)基因，DEGs見火山圖(圖7A)。將DEGs與上一小節(jié)中的9模塊的基因取交集。結(jié)果發(fā)現(xiàn)violet、darkgrey、lightyellow、darkorange、lightcyan1模塊與差異基因并無交集，因此剔除了這5個模塊。從其余4模塊purple、greenyellow、pink、darkolivegreen中，共提取到122個基因(圖7B)。進一步的GO富集分析發(fā)現(xiàn)多個纖維增生，肌細胞分化、增殖，細胞外基質(zhì)等功能模塊富集(圖8A)。KEGG富集分析發(fā)現(xiàn)，篩選的基因集在Hippo、絲裂原活化蛋白激酶(MAPK)、哺乳動物雷帕霉素靶點(mTOR)等通路富集(圖8B)。

圖7 聯(lián)合差異基因進一步篩選WGCNA前列腺間質(zhì)增生相關(guān)基因A：BPH組織和正常前列腺組織的差異表達基因；B，差異表達基因及WGCNA的purple、greenyellow、pink、darkolivegreen模塊聯(lián)合分析。

圖8 WGCNA前列腺間質(zhì)增生相關(guān)差異基因的富集分析A： GO富集分析發(fā)現(xiàn)，多個纖維增生，肌細胞分化、增殖，細胞外基質(zhì)等功能模塊富集；B：KEGG富集分析發(fā)現(xiàn)，篩選的基因集在Hippo、MAPK、mTOR、TGF-β等通路富集。

GO和KEGG富集分析的結(jié)果進一步證明，BPH組織中的成纖維細胞相關(guān)通路激活，誘導成纖維細胞增生、分化，這一過程可能參與了BPH的發(fā)生和發(fā)展。

本研究將篩選所得的122個基因輸入STRING數(shù)據(jù)庫，構(gòu)建PPI網(wǎng)絡(luò)(圖9A)。并使用Cytoscape的cytoHubba插件對hub基因進行篩選。使用了MCC、DMNC、MNC、Degree4算法計算top15的hub基因，取交集獲得與前列腺間質(zhì)增生的相關(guān)hub基因(圖9B)。通過PCA發(fā)現(xiàn)，12個hub基因能夠很好的區(qū)分前列腺組織和BPH組織，驗證了這12個基因在前列腺增生發(fā)生中可能發(fā)揮了重要作用(圖10A～10C)。

圖9 蛋白互作網(wǎng)絡(luò)分析及關(guān)鍵基因篩選A：本研究篩選所得的122個基因構(gòu)建PPI網(wǎng)絡(luò)；B：使用Cytoscape的cytoHubba插件對hub基因進行篩選。MCC、DMNC、MNC、Degree 4種算法計算top15的關(guān)鍵基因，取交集獲得與前列腺間質(zhì)增生的相關(guān)關(guān)鍵基因。

圖10 hub基因的表達和主成分分析A：12個hub基因在前列腺組織和BPH組織的表達熱圖；B：12個hub基因的蛋白互作網(wǎng)絡(luò)關(guān)系；C：12個hub基因的主成分分析能夠很好的區(qū)分前列腺組織和BPH組織。

2.4 hub基因與前列腺的間質(zhì)狀態(tài)相關(guān) 獲得了hub基因后，本研究進一步分析了這12個差異基因與增生前列腺間質(zhì)狀態(tài)的關(guān)系(圖11)。肌球蛋白輕鏈1(MYL1)、肌球蛋白輕鏈2(MYL2)、伴肌動蛋白(NEB)、Kirsten大鼠肉瘤病毒癌基因同源物(KRAS)、Snail家族轉(zhuǎn)錄抑制子2(SNAI2)等與間質(zhì)狀態(tài)正相關(guān)。其余基因未有統(tǒng)計學意義可能與較少的樣本量有關(guān)。相反，人肌球蛋白重鏈3(MYH3)、人肌球蛋白重鏈2(MYH2)等與間質(zhì)狀態(tài)負相關(guān)。

圖11 關(guān)鍵基因與增生前列腺間質(zhì)狀態(tài)的相關(guān)性分析

2.5 BPH間質(zhì)相關(guān)標記基因和模型的建立和驗證 本研究分析了12個hub基因在BPH組織和正常前列腺組織中的表達量(圖10A)發(fā)現(xiàn)，WNT2、胰島素樣生長因子1(IGF1)、SNAI2在BPH組織中表達較高。目前，在BPH研究領(lǐng)域，如何鑒定BPH標本或正常前列腺標本是一個突出的難題。前列腺體積較大，肉眼觀察難以直接分辨病變組織與正常組織，在收集患者的手術(shù)標本以后，研究者往往難以明確標本的具體病例性質(zhì)。因此，本研究嘗試建立BPH的基因診斷模型，利用LASSO回歸從12個hub基因中進行篩選，最終篩選出4個變量(圖12A、12B)。分別為NEB、DMD、IGF1、KRAS回歸系數(shù)分別為-0.980 750、-4.091 432、4.339 909、-2.959 698。利用建立的模型在驗證集中進行驗證(圖12C)顯示，該模型能夠很好地對BPH組織進行基因診斷分類(AUC=0.889)。

圖12 BPH基因模型的建立和驗證A：LASSO回歸從12個hub基因中進行篩選，最終篩選出NEB、DMD、IGF1、KRAS建立模型，回歸系數(shù)分別為-0.980 750、-4.09 143 2、4.339 909、-2.959 698；B：利用建立的模型在驗證集中進行驗證。

3 討 論

BPH是老年人膀胱出口梗阻的主要原因，能夠引起LUTS癥狀，影響生活質(zhì)量，但BPH發(fā)生的分子生物學機制尚不明確。本研究應(yīng)用了生物信息學技術(shù)，通過分析BPH組織和正常前列腺組織的轉(zhuǎn)錄組數(shù)據(jù)，確認了BPH組織的間質(zhì)狀態(tài)，得到了一個由4個DEG組成的標記基因，同時在驗證數(shù)據(jù)集中得到了驗證(AUC=0.889)。這可能為BPH的分子學診斷提供了新的思路，也為進一步研究BPH的發(fā)生發(fā)展的分子生物學機制提供了新的方向。

研究[20]表明，BPH患者存在不同程度的間質(zhì)成分改變，并和LUTS癥狀的相關(guān)。但尚無研究團隊對BPH間質(zhì)的細胞組成進行精準分析，并對其進行基因分型。BPH患者前列腺間質(zhì)組織中膠原成分與肌纖維比例低于正常前列腺組織，該比例與LUTS癥狀正相關(guān)。而膠原蛋白正是由成纖維細胞分泌，并受TGF-β等通路調(diào)節(jié)。這從側(cè)面驗證了本研究的結(jié)論。同時，慢性炎癥也是BPH的重要驅(qū)動因素[21]。有研究[22]分離出BPH患者前列腺的原代間質(zhì)細胞與永生化的前列腺上皮細胞株BPH-1進行體外共培養(yǎng)，能夠加快BPH-1的增殖，而從正常前列腺分離出的間質(zhì)細胞與BPH-1細胞系體外共培養(yǎng)則不能刺激BPH-1的增殖，這說明BPH的間質(zhì)細胞能夠外分泌特定的生長因子刺激前列腺上皮增生。另一項研究[23]顯示，前列腺間質(zhì)細胞能夠分泌多種生長類細胞因子如VEGF、成纖維細胞生長因子(FGF)、TGF-β等細胞因子，刺激前列腺間質(zhì)細胞增生，特別是成纖維細胞，誘導成纖維細胞分泌膠原蛋白。而與之相反，前列腺上皮細胞不能分泌這類細胞因子。因此，前列腺間質(zhì)狀態(tài)改變在BPH發(fā)生發(fā)展中起重要作用。

WGCNA是利用基因表達數(shù)據(jù)和表型數(shù)據(jù)構(gòu)建無標度網(wǎng)絡(luò)。它能尋找協(xié)同作用的基因模塊，探索基因網(wǎng)絡(luò)和感興趣的疾病表型之間的關(guān)系，及網(wǎng)絡(luò)中的中樞基因。WGCNA分析方法已成為疾病研究領(lǐng)域的一種重要生物信息學分析手段，它側(cè)重于基因共表達模塊和表型特征之間的關(guān)聯(lián)，因而具有相對較高的可靠性和生物學意義。本研究通過WGGCNA方法構(gòu)建共表達模塊，分析了與BPH組織中間質(zhì)評分、成纖維細胞、平滑肌細胞相關(guān)的表達模塊，確定了4個與BPH間質(zhì)評分、成纖維細胞、平滑肌細胞顯著相關(guān)的共表達模塊，用于檢測BPH轉(zhuǎn)錄組與BPH組織間質(zhì)狀態(tài)之間的關(guān)系。進一步對4個表達模塊和差異基因的交集分析，功能富集分析，和PPI分析發(fā)現(xiàn)，NEB、DMD、IGF1、KRAS可作為BPH和正常前列腺組織基因分型的hub基因。其中IGF1經(jīng)由雄激素受體通路誘導，主要通過旁分泌刺激前列腺間質(zhì)細胞和前列腺上皮細胞的病理性增殖，特別是成纖維細胞的增殖[24]。間質(zhì)細胞自身表達IGF1受體，因此能夠被IGF1激活，經(jīng)由IGF1通路，加快有絲分裂。BPH組織中KRAS基因活性抑制也同時得到了其他研究團隊的驗證，其在前列腺中的激活狀態(tài)也收到雄激素受體通路的調(diào)控[25]。而NEB和DMD基因在BPH疾病中的功能和作用目前尚無文獻報道。

綜上所述，本研究對BPH和正常前列腺組織的轉(zhuǎn)錄組數(shù)據(jù)進行生物信息學分析發(fā)現(xiàn)，BPH組織中間質(zhì)細胞成分顯著高于正常前列腺組織。BPH組織存在較多的成纖維細胞和較少的平滑肌細胞，這可能和前列腺的發(fā)生發(fā)展相關(guān)。本研究也建立了4個基因的BPH診斷模型，這些信息對BPH的分子生物學機制研究及藥物研發(fā)可能具有一定的參考價值。