影響前列腺癌風險的關鍵基因識別

李熹陽 谷明宇 華琳

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

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

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

1數據與方法

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

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

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

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

2結果

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

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

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

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

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

參考文獻：

[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全球癌癥統計報告解讀[J].腫瘤綜合治療電子雜志，2019，5（1）：87-97.

[4]王煒，李傳剛，劉輝，等.前列腺特異性抗原對前列腺癌診斷價值的探討[J].中國醫科大學學報，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].中國現代醫學雜志，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，Fuso 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

編輯/肖婷婷