影響前列腺癌風險的關(guān)鍵基因識別

李熹陽 谷明宇 華琳

摘要：目的? 基于TCGA數(shù)據(jù)庫篩選影響前列腺癌（PCa）風險水平的關(guān)鍵基因，并建立PCa患者生存風險預測模型。方法? 從TCGA數(shù)據(jù)庫下載PCa患者基因表達數(shù)據(jù)及相關(guān)臨床數(shù)據(jù)，通過前期研究初步篩選基因，并將患者分為高、低風險兩類;對基因進行差異表達分析和GO和KEGG通路富集分析，篩選相關(guān)基因和信號通路;對差異表達基因進行蛋白互作網(wǎng)絡(luò)分析，標記出關(guān)鍵基因;將關(guān)鍵基因的表達數(shù)據(jù)與PCa患者生存時間納入Cox回歸分析，建立生存風險預測模型。結(jié)果? 前期研究得到620個基因，高風險患者234例，低風險患者285例;差異表達分析獲得30個基因，主要分子功能（MF）為：受體結(jié)合和生長因子活動，生物學過程（BP）主要為細胞-細胞信號傳導、細胞增殖的積極調(diào)節(jié)、血管生成的調(diào)節(jié)和細胞表面受體信號通路，細胞組分（CC）主要定位于細胞外區(qū)域，而KEGG信號通路為細胞因子-細胞因子受體相互作用;蛋白互作分析中共7個基因有相互作用，Cytoscape篩選出5個關(guān)鍵基因：PHYHIPL、CNTFR、GFRA1、EDN3和PROK1。結(jié)論? 通過本研究識別的影響PCa預后的關(guān)鍵基因，發(fā)現(xiàn)潛在的PCa風險靶點，可能為PCa的治療和預后提供幫助。

關(guān)鍵詞：前列腺癌;基因差異表達分析;生物信息學;富集分析

中圖分類號：R737.25? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻標識碼：A? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? DOI：10.3969/j.issn.1006-1959.2020.02.022

文章編號：1006-1959（2020）02-0080-06

Abstract：Objective? To screen the key genes affecting prostate cancer （PCa） risk level based on the TCGA database and establish a survival risk prediction model for PCa patients. Methods? Download PCa patient gene expression data and related clinical data from the TCGA database， preliminary screening of genes through early research， and classify patients into high and low risk categories; perform differential expression analysis of genes and enrichment analysis of GO and KEGG pathways to screen relevant gene and signal pathway; Perform protein interaction network analysis on differentially expressed genes to mark key genes; incorporate expression data of key genes and PCa patient survival time into Cox regression analysis to establish a survival risk prediction model. Results? 620 genes were obtained in previous studies， 234 patients were high-risk patients， 285 patients were low-risk patients; 30 genes were obtained by differential expression analysis. The main molecular functions （MF） were： receptor binding and growth factor activity， and the main biological process （BP） For cell-cell signaling， positive regulation of cell proliferation， regulation of angiogenesis， and cell surface receptor signaling pathways， the cell component （CC） is mainly located in the extracellular region， while the KEGG signaling pathway is a cytokine-cytokine receptor interactions： A total of 7 genes interacted in the protein interaction analysis. Cytoscape screened out 5 key genes： PHYHIPL， CNTFR， GFRA1， EDN3， and PROK1.Conclusion? The key genes affecting the prognosis of PCa identified through this study， and the discovery of potential PCa risk targets may help the treatment and prognosis of PCa.

Key words：Prostate cancer;Differential gene expression analysis;Bioinformatics;Enrichment analysis

前列腺癌（prostate cancer，PCa）是一種上皮性惡性腫瘤，是男性最為常見的癌癥類型[1]，國際癌癥研究署（IARC）公開數(shù)據(jù)顯示，2018年中國PCa標化發(fā)病率為13.9/10萬，為低風險地區(qū)，但近10年來，PCa已成為我國增速最快的男性惡性腫瘤[2]，并且由于人口老齡化進程加快，我國PCa發(fā)病率也高居男性惡性腫瘤第6位[3]。目前在PCa臨床篩查與診斷中，廣泛應(yīng)用的血清前列腺特異抗原（PSA）檢測敏感度高而特異性低[4]，亟需尋找新的PCa標志物;癌癥基因組圖譜公共數(shù)據(jù)集（TCGA）數(shù)據(jù)庫存有大規(guī)模的基因測序和患者相關(guān)臨床指標，本研究基于TCGA數(shù)據(jù)庫下載的前列腺癌組織中多種基因的mRNA三代測序數(shù)據(jù)，利用多種生物信息學和生物統(tǒng)計學方法，識別影響高、低風險前列腺癌形成的關(guān)鍵基因，以期建立了前列腺癌患者生存時間統(tǒng)計預測模型，以期為探索PCa分子治療提供一定的參考。

1數(shù)據(jù)與方法

1.1數(shù)據(jù)庫選擇? 選取在TCGA腫瘤數(shù)據(jù)庫收錄的前列腺癌患者基因表達level 3數(shù)據(jù)（截至2018年初）為研究對象。生物學信息注釋數(shù)據(jù)庫（Database for Annotation，Visualization，and Integrated Discovery，DAVID）6.8在線工具（https：//david.ncifcrf.gov）用于對TCGA數(shù)據(jù)庫中篩選出的差異表達基因進行基因本體論（GO）和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）信號通路富集分析，其中GO分析包括細胞組分（CC），生物過程（BP）和分子功能（MF）三個部分。應(yīng)用STRING（Search Tool for the Retrieval of Interacting Genes）數(shù)據(jù)庫（https：//string-db.org/cgi/input.pl），對差異表達基因編碼蛋白構(gòu)建相互作用網(wǎng)絡(luò)，選擇Cytoscape[5] 3.7.2中cytoHubba插件篩選蛋白互作網(wǎng)絡(luò)中的關(guān)鍵基因。

1.2樣本風險聚類? 從TCGA數(shù)據(jù)庫中下載得到20530個基因的表達數(shù)據(jù)，568例PCa患者及對應(yīng)的生存時間與生存狀態(tài)。基因表達數(shù)據(jù)已經(jīng)通過Tophat2比對到染色體hg19上，并生成標準化的FPKM（fragments per kilobase of transcript per million fragments mapped）值。隨后采用Cox風險比例回歸模型提取與PCa患者生存有關(guān)的基因。剔除含有缺失數(shù)據(jù)的樣本，并基于K均值聚類算法，將其分為兩組。根據(jù)死亡事件發(fā)生的頻率，分別命名兩組樣本為高風險組和低風險組[6]。

1.3數(shù)據(jù)分析? 將前期研究獲得的基因表達數(shù)據(jù)整理、錄入R軟件中，并繪制基因表達箱圖，比較樣本表達數(shù)據(jù)的分布情況。采用R軟件limma[7]包，使用無監(jiān)督聚類方法展示519例樣本間的相似性，繪制多維標度分析（multidimensional scaling，MDS）圖，并使用edgeR[8]包估算所有基因的離散度，即生物變異系數(shù)（biological coefficient of variation，BCV）的平方[9]，以展示基因差異程度。利用limma包進行差異表達分析，對低風險PCa患者腫瘤組織和高風險PCa患者腫瘤組織中的差異表達基因進行篩選。篩選條件為：差異表達超過2倍（|log FC|>1），校正后P<0.1且P<0.05。針對上述條件篩選出的差異表達基因，分別使用ggplot2[10]軟件包和gplots軟件包繪制火山圖與熱圖，驗證并展示差異表達基因的結(jié)果。利用DAVID 6.8在線工具對差異表達基因進行GO富集分析和 KEGG信號通路富集分析，以P<0.05為標準篩選出具有顯著性的差異表達基因功能注釋和KEGG通路富集分析結(jié)果。在STRING數(shù)據(jù)庫對上述篩選出的差異表達基因編碼蛋白構(gòu)建相互作用網(wǎng)絡(luò)，條件為：①蛋白互作數(shù)據(jù)來源：Textmining，Experiments， Databases，Co-expression，Neighborhood，Gene Fusion，Co-occurrence;②最低作用評分要求為：0.4。將構(gòu)建的蛋白互作網(wǎng)路數(shù)據(jù)導入Cytoscape中，利用cytoHubba插件查找hub節(jié)點，即相關(guān)信號通路中的關(guān)鍵基因。使用R軟件rms包對篩選出的關(guān)鍵基因與PCa患者生存時間進行Cox回歸分析，將關(guān)鍵基因的表達數(shù)據(jù)作為自變量，PCa患者生存時間作為因變量，并繪制諾莫圖（Nomograph），建立統(tǒng)計預測模型，預測PCa患者的生存風險。

1.4統(tǒng)計學處理? 使用R軟件（3.6.1）進行統(tǒng)計學分析。通過繪制Kaplan-Meier曲線進行生存分析，并且使用對數(shù)秩檢驗法（Log-rank）檢驗顯著性。差異表達分析采用經(jīng)驗貝葉斯先驗趨勢法（empirical bayes prior trend），亦即“l(fā)imma-trend”法，對均值-方差關(guān)系建模。P<0.05表示差異有統(tǒng)計學意義。

2結(jié)果

2.1樣本風險聚類? 經(jīng)過單變量Cox分析，獲得620個與前列腺癌患者生存相關(guān)的基因;K均值聚類算法將前列腺癌患者聚為兩類，其中第一類（高風險）PCa患者234例，含7個死亡病例;第二類（低風險）PCa患者285例，沒有死亡病例。繪制Kaplan-Meier曲線比較高風險（group1）與低風險（group2）PCa患者的生存率，采用Log-rank檢驗法測定生存率曲線間的顯著性（圖1），結(jié)果顯示兩條生存曲線之間存在統(tǒng)計學差異（P=0.004），說明通過聚類方法有效地將PCa患者聚為了高風險、低風險兩類。

2.2差異表達基因篩選? 基因表達箱式圖，見圖2，可見基因表達數(shù)據(jù)較為整齊，可以進行Limma分析; MDS圖（圖3）顯示在前兩個維度構(gòu)成的坐標系中，所有樣本相似度良好;所有基因的估計離散值（圖4），包括經(jīng)驗貝葉斯穩(wěn)健離散值（Tagwise）、經(jīng)驗貝葉斯穩(wěn)健離散值的均值（Common）和經(jīng)驗貝葉斯穩(wěn)健離散值的擬合值（Trended），曲線顯示樣本之間的差異較小。Limma篩選出30個差異表達基因，其中上調(diào)基因6個，下調(diào)基因24個（表1）。繪制上述差異表達基因的火山圖（圖5）。使用gplots軟件包繪制熱圖（圖6），橫坐標為樣本編號，縱坐標為基因種類，并對樣本和基因雙向聚類，展示表達差異分析結(jié)果。從圖6左側(cè)系統(tǒng)聚類樹狀圖可見，30個表達差異基因恰好將519個樣本分為兩類;上調(diào)的基因顯示為黃色，下調(diào)的基因顯示為紅色，6個上調(diào)基因：KRTAP5-AS1，LOC100128842，SKA3，HPN-AS1，SLC44A5和LRRC31表達情況明顯異于其他基因，與差異表達分析得出的結(jié)果相符。

2.3差異表達基因的功能富集分析? 以P<0.05篩選功能富集分析結(jié)果（圖7）：差異表達基因 MF主要為受體結(jié)合（P=0.010），生長因子活動（P=0.019）;BP主要為細胞-細胞信號傳導（P=0.0043），細胞增殖的積極調(diào)節(jié)（P=0.022），血管生成的調(diào)節(jié)（P=0.040），細胞表面受體信號通路（P=0.049）;CC主要為細胞外區(qū)域（P=0.019），KEGG信號通路為細胞因子-細胞因子受體相互作用（P=0.046）。

腺和胎盤）中表達，該基因編碼的蛋白質(zhì)誘導腺體中的毛細血管內(nèi)皮細胞的增殖、遷移。Jilg CA等[19]研究發(fā)現(xiàn)，PRK1（即PROK1）基因有望成為治療雄激素非依賴性轉(zhuǎn)移性PCa的治療靶標。

本研究還基于上述關(guān)鍵基因表達與生存時間的Cox回歸構(gòu)建了統(tǒng)計預測模型。諾莫圖是一個復雜數(shù)學公式的圖形表示[20]，在腫瘤學和醫(yī)學中被廣泛被用作預測手段，是現(xiàn)代醫(yī)學決策制定的重要組成部分，可以通過綜合不同預后和決定性變量，產(chǎn)生個體發(fā)生臨床事件的數(shù)字概率[21]。本研究繪制出的諾莫圖具有良好的區(qū)分度，C-index達到了0.729，不同關(guān)鍵基因表達水平組合的患者可以得到良好的區(qū)分。

參考文獻：

[1]Howard N，Clementino M，Kim D，et al.New developments in mechanisms of prostate cancer progression[J].Semin Cancer Biol，2019（57）：111-116.

[2]前列腺癌成中國增速最快的男性惡性腫瘤[J].腫瘤防治研究，2019，46（7）：666.

[3]王寧，劉碩，楊雷，等.2018全球癌癥統(tǒng)計報告解讀[J].腫瘤綜合治療電子雜志，2019，5（1）：87-97.

[4]王煒，李傳剛，劉輝，等.前列腺特異性抗原對前列腺癌診斷價值的探討[J].中國醫(yī)科大學學報，2016，45（1）：61-65，69.

[5]Shannon，P.Cytoscape：A Software Environment for Integrated Models of Biomolecular Interaction Networks[J].Genome Research，2003，13（11）：2498-2504.

[6]Hua L，Xia H，Xu W，et al.Risk stratification for prostate cancer via the integration of omics data of The Cancer Genome Atlas[J].Translational Cancer Research，2018，7（3）：706-719.

[7]Phipson B，Lee S，Majewski IJ，et al.Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression[J].The Annals of Applied Statistics，2016，10（2）：946-963.

[8]Robinson MD，Smyth GK.Moderated statistical tests for assessing differences in tagabundance[J].Bioinformatics，2008，23（21）：2881-2887.

[9]Mccarthy DJ，Chen Y，Smyth GK.Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation[J].Nucleic Acids Research，2012，40（10）：4288-4297.

[10]Wickham H.ggplot2：elegant graphics for data analysis[J].J R Stat Soc，2011，174（1）：245-246.

[11]Hendriks RJ，Dijkstra S，Smit FP，et al.Epigenetic markers in circulating cell-free DNA as prognostic markers for survival of castration-resistant prostate cancer patients[J]. Prostate，2018，78（5）：336-342.

[12]Yin X，Yu J，Zhou Y，et al.Identification of CDK2 as a novel target in treatment of prostate cancer[J].Future Oncology，2018，14（8）：709-718.

[13]蔣宏毅，趙曉昆，鐘朝暉，等.PTEN/MMAC1/TEP1、TGF-β1在前列腺癌及前列腺增生中的表達及其意義[J].中國現(xiàn)代醫(yī)學雜志，2008，18（9）：1221-1225.

[14]Hearn JWD，Abuali G，Magi-Galluzzi C，et al.HSD3B1 and resistance to androgen deprivation therapy in prostate cancer[J].Journal of Clinical Oncology，2015，33（7_suppl）：156-156.

[15]Fu H，Ge B，Chen D，et al.Phytanoyl-CoA 2-Hydroxylase-Interacting Protein-Like Gene Is a Therapeutic Target Gene for Glioblastoma Multiforme[J].Med Sci Monit，2019（25）：2583-2590.

[16]Cousinery Mary C，LiRuili，VannitambyAmanda，et al.Neurotrophin signaling in a genitofemoral nerve target organ during testicular descent in wmice[J].Journal of Pediatric Surgery，2016，51（8）.

[17]Grasso M，F(xiàn)uso A，Dovere L，et al.Distribution of GFRA1-expressing spermatogonia in adult mouse testis[J].Reproduction，2012，143（3）：325-332.

[18]謝萍芳，趙東怡，周美容，等.乳腺癌中GFRA1表達臨床意義的生物信息學分析[J].中國腫瘤臨床，2018，45（15）：769-773.

[19]Jilg CA，Ketscher A，Metzger E，et al.PRK1/PKN1 controls migration and metastasis of androgen-independent prostate cancer cells[J].Oncotarget，2014，5（24）：12646-12664.

[20]Grimes David A.The nomogram epidemic：resurgence of a medical relic[J].Annals of Internal Medicine，2008，149（4）：273-275.

[21]Vinod P Balachandran，MithatGonen，J Joshua Smith，et al. Nomograms in oncology：more than meets the eye[J].The Lancet Oncology，2015，16（4）：e173-e175.

收稿日期：2019-10-23;修回日期：2019-11-10

編輯/肖婷婷