影響前列腺癌風(fēng)險(xiǎn)的關(guān)鍵基因識(shí)別

李熹陽(yáng) 谷明宇 華琳

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

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

中圖分類(lèi)號(hào)：R737.25? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? DOI：10.3969/j.issn.1006-1959.2020.02.022

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

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

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

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

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

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

2結(jié)果

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

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

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

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

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

參考文獻(xiàn)：

[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]前列腺癌成中國(guó)增速最快的男性惡性腫瘤[J].腫瘤防治研究，2019，46（7）：666.

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

[4]王煒，李傳剛，劉輝，等.前列腺特異性抗原對(duì)前列腺癌診斷價(jià)值的探討[J].中國(guó)醫(yī)科大學(xué)學(xué)報(bào)，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在前列腺癌及前列腺增生中的表達(dá)及其意義[J].中國(guó)現(xiàn)代醫(yī)學(xué)雜志，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表達(dá)臨床意義的生物信息學(xué)分析[J].中國(guó)腫瘤臨床，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

編輯/肖婷婷