前列腺癌差異基因篩選、功能預測及臨床意義分析

楊 智 米 軍

蘭州大學第二醫院泌尿外科研究所/甘肅省泌尿系疾病重點實驗室/甘肅省泌尿系疾病臨床醫學中心(甘肅 蘭州 73000)

前列腺癌的發病率在全球范圍內不斷升高。 在美國,前列腺癌位于惡性腫瘤及癌癥相關住院率的第二位[1]。 有研究表明,多種mRNA 參與前列腺癌發生、發展過程。 本研究基于美國基因表達綜合數據庫NCBI Gene Expression Omnibus(GEO)和美國癌癥腫瘤基因圖譜(the cancer genome atlas, TCGA)數據庫對前列腺癌組織及良性組織中的差異表達基因進行篩選,并探討其相關分子機制和臨床預后意義。

材料與方法

一、數據提取

從美國基因表達綜合數據庫NCBI Gene Expression Omnibus(GEO)下載目標前列腺癌組織及非癌組織mRNA 轉錄組數據,形成三個數據集。 其中數據集GSE46602 包含36 例前列腺癌樣本,14 例前列腺良性組織樣本;數據集GSE3325 包含13 例前列腺癌樣本,6例前列腺良性增生組織樣本;數據集GSE104749 包含4例前列腺癌組織樣本,4 例前列腺良性增生組織樣本。

二、差異表達分析

在R 語言環境中,利用edgeR package 對GSE46602,GSE3325 和GSE104749 三個選擇的前列腺癌數據集進行正常組織與癌癥組織的差異表達分析,篩選出差異表達基因。 edge R 工具包讀取文件,包括基因ID、log2轉化后的差異倍數(Fold Change)值、顯著性p值以及校正后p值(默認偽發現率,FDR 校正)。 然后,設定篩選閾值為:FC 值＜1 且P＜0.05,用于比較篩選出前列腺癌組織和非癌組織的差異表達的基因。 定義|log2FC|＞1 以及FDR＜0.01 為差異表達基因。 其中log2FC＞1 設定為“up”,即差異表達基因中的上調基因;log2FC＜-1設定為“down”,即差異表達基因中的下調基因。 獲得識別了顯著差異表達的基因后,載入ggplot2,繪制數據集的火山圖。 載入pheatmap package,將獲得的基因差異表達結果,完成基因ID 轉換后,獲得相應的基因表達熱圖。

三、關鍵基因篩選

利用韋恩圖,選定上調的差異表達基因中前六位的基因認定為關鍵基因,進一步進行關鍵基因的功能分析和臨床預后預測。

四、GO 及KEGG 信號通路分析

在R 語言環境下,載入openxlsx R 包,用于讀取之前獲得的差異表達篩選的文件。 載入stringr R 包用于完成基因的ID 轉換;載入enrichplot R 包和clusterProfiler R 包用于對功能分析的基因進行富集分析及結果可視化,最后利用barplot R 包和dotplot R 包完成圖片繪制。

五、關鍵基因的蛋白組學表達情況

在HPA(Human Protein Atlas,https:/ /www.proteinatlas.org/)直接檢索目標基因可獲得目標基因的組織特性、組織及單細胞表達聚類性質及蛋白質功能等一系列信息。

六、基因表達水平及生存分析

比較關鍵基因在前列腺癌組織及非癌組織中的表達水平。 利用GEPIA(http:/ /gepia. cancer-pku. cn/),結合TCGA 和GTEx 相關數據,可直接分析出單一基因在腫瘤組織和非腫瘤組織的表達差異和相關性結果。利用UCSC Xena(https:/ /xena. ucsc. edu),將目標基因選擇“Kaplan Meier plot”,輸出目標基因表達量和無進展間隔期(progression-free interval,PFI)的生存曲線結果。

七、統計學方法

本研究采用SPSS 22、R 語言(R 4.04)及R studio(2021.09. 1 Build 372)處理數據。 Cytoscape(version 3.5.1 java 1.8.2_291)、GraphPad prism 8 繪制圖像。 計量資料若符合正態分布,以均數±標準差(±s) 表示,比較用t檢驗,若非正態性分布則采用非參數檢驗;計數資料以率(%)表示,比較用卡方檢驗。 UCSC Xena采用Kaplan-Meier 法繪制目標基因表達量的生存曲線,分為兩組或三組,比較采用Log rankχ2檢驗。P＜0.05為差異有統計學意義。

結 果

一、差異表達基因

根據差異基因的篩選條件,GSE46602 共讀取229個差異表達基因,包含68 個上調的差異表達基因和161 個下調的差異表達基因;GSE3325 共讀取468 個差異表達基因,包含232 個上調的差異表達基因和236 個下調的差異表達基因;GSE104749 共讀取1007 個差異表達基因,包含402 個上調的差異表達基因和605 個下調的差異表達基因。 篩選出上調基因中,靠前六位的基因作為差異表達基因的關鍵基因,關鍵基因為:PCA3(前列腺癌抗原3), SIM2(Single-Mind 2), AMACR(異構梅2-甲酰輔酶A 消旋梅), HOXC6(同源異型盒-6 基因), ERG(ETS 轉錄調節因子)及TOP2A(拓撲異構梅ⅢA 型)。 (見圖1A-圖1H)

圖1 篩選差異表達基因的熱圖、火山圖及前六位關鍵基因的韋恩圖

二、GO 及KEGG 富集分析結果

GO 結果顯示,差異表達基因的生物功能主要與減數分裂核分裂、減數分裂細胞周期、三價鐵結合和氧化還原酶活性相關。 KEGG 數據庫結果顯示,差異表達基因主要富集的代謝途徑包括:原代膽汁酸生物合成、嘧啶代謝、谷胱甘肽代謝、鉑類耐藥性及p53 信號通路(見圖2)。

圖2 GO+KEGG 可視化結果

三、蛋白互作網絡結果

將數據導入STRING 數據庫獲得蛋白質-蛋白質相互作用網絡(見圖3),PPI 共有17 點節點,35 個面;從類聚結果分析,主要類聚為三個部分,分別以TOP2A、TACC3、AMACR、TP53、SIM2、ERG、HOXC6 等為中心。

圖3 蛋白質-蛋白質相互作用網絡

四、關鍵基因表達水平驗證

借助人類蛋白質表達圖集(The Human Protein Atlas,HPA),能查看和分析關鍵基因編碼蛋白在不同腫瘤組織/器官中的表達情況(見圖4,由于PCA3為非編碼基因,故未納入PCA3 在組織中的染色標本)。

圖4 HPA 數據庫中關鍵基因的組織學染色

五、關鍵基因生存分析

結合TCGA 數據庫中前列腺癌的臨床相關信息,將篩選的關鍵基因做臨床預后分析,結果見圖5-7。 對于PCA3,正常非癌組織中的表達量和T2、T3 和T4 期前列腺癌組織比較,其表達量均有統計學差異。 PCA3 在前列腺癌N 分期中,隨著疾病的進展,PCA3 的表達量階梯式上升,但是PCA3 和前列腺癌的遠處轉移無明顯統計學差異。 SIM2 和AMACR 在前列腺癌中表達水平和T 分期基本無關系,非癌組織中的表達量明顯低于前列腺癌組織;N 分期和M 分期的結果類似,SIM2 和AMACR 的表達量隨著疾病的進展并未明顯變化。HOXC6 和ERG 表達量在TNM 分期的各組間基本無明顯的統計學差異,多為良性前列腺組織的表達量低于前列腺癌組織。 TOP2A 的表達量與PCA3 的趨勢相似,隨著癌癥的進展,TOP2A 的表達量在TNM 分期中呈現階梯式上升。 進一步分析關鍵基因與前列腺癌無進展間隔期(progression-free interval, PFI)的關系。 PCA3 和TOP2A的表達情況對前列腺癌的無進展間隔期(progression-free interval,PFI)時間存在統計學差異(p＜0.05),PCA3 的高表達量組和低表達量組的PFI 存在統計學差異(p=0.005), HR 為0.55(0.36-0.84);TOP2A 的高表達量組和低表達量組的PFI 也存在統計學差異(p＜0.001),HR 為2. 33(1. 51-3. 62),其余關鍵基因(SIM2、AMACR、HOXC6 和ERG)未顯著影響前列腺癌PFI。

圖5 關鍵基因表達量對TNM 分期的影響

圖6 關鍵基因表達量對TNM 分期的影響

圖7 關鍵基因表達量對PFI 的生存分析

討 論

前列腺癌作為美國男性最常見的惡性腫瘤,全球約有21% 的癌癥病例為前列腺癌[1]。 本研究結合GEO、TCGA 等相關數據庫,利用生物信息學技術,篩選出與前列腺癌相關的六個關鍵基因,分別為PCA3、SIM2、AMACR、HOXC6、ERG 和TOP2A。 篩選的關鍵基因中,PCA3 和TOP2A 在前列腺癌中表達較顯著,與前列腺癌TNM 分期,腫瘤的無進展間隔期顯著現關,PCA3 和TOP2A 有望成為預測前列腺癌預后的相關指標及新的腫瘤標志物。

PCA3 是一種長非編碼RNA(lncRNA),最初通過差異顯示實驗法( differential display experimental approach)被命名為DD3[2]。 PCA3 通過調節雄激素受體(AR)信號傳導進而影響前列腺癌(PCa)細胞功能,同時控制多個雄激素應答和癌癥相關基因的表達,包括上皮-間質轉化(EMT)過程。 此外,尿液中的PCA3已成功應用于PCa 的臨床診斷[3]。 Qin Z 等[4]對8139名病例和14116 名對照病例進行薈萃分析,以評估PCA3 的診斷價值,PCA3 的總診斷優勢比(DOR)較高,95%置信區間為5.44(4.53-6.53)。 Ye LF等[5]發明了一種新方法,即LBXexo 評分法,用于測量尿液中的外泌體PCA3/PRAC 表達水平,結果證實LBXexo 可以改進中國人群前列腺癌的診斷。 對于臨床應用,Alshalafa M 等[6]發現PCA3 在初始活檢中預測高級別前列腺癌(Gleason Score,GS≥8)的表現較差,假陰性率較高。Kotova ES 等[7]在PCa 患者尿液樣本中檢測PCA3 和AMACR 的mRNA 水平。 在其研究中,PCA3 評分AUC為0.632(95%置信區間:0.511-0.752),AMACR 評分AUC 為0.711(95%置信區間:0.617-0806),較高的預測價值使其成為潛在的前列腺癌生物標記物。

SIM2(single-Mind 2)基因是具有基本螺旋-環-螺旋/果蠅Per 蛋白-芳香烴受體核轉運蛋白-果蠅Sim 蛋白(helix-loop-helix/per-Arnt-Sim,bHLH/PAS)結構域的轉錄因子家族的成員,影響實體瘤的發生過程。 在乳腺癌中,SIM2s 直接下調SNAI2(Snail Family Transcriptional Repressor 2)的表達,抑制EMT,從而抑制腫瘤生長和侵襲[8]。 Wyatt GL 等[9]發現SIM2s 和NFκB(核因子κB,nuclear factor kappa-B)之間存在串擾,對乳腺腫瘤的生長和轉移具有抑制作用。

AMACR 最為人所知的生理作用是催化支鏈脂肪酸的α-甲基質子在線粒體和過氧化物酶體中的β-氧化[10]。 前列腺癌細針穿刺活檢中,AMACR 表達對PCa檢測具有97%的敏感性和100%的特異性[11]。 Xie H等[12]使用miR200c 抑制前列腺癌中的AMACR 活性,研究證實抑制AMACR 可以阻礙細胞增殖和遷移。

Homeobox C6,HOXC6 屬于同源盒家族(homeobox,HOX),其成員編碼高度保守的轉錄因子,其家族有39個HOX 基因位于4 個染色體位點。 HOXC6 已被證實為侵襲性前列腺癌中可靠的生物標記物,可提高早期診斷率并預測治療后癌癥復發[14]-[17]。 Luo Z 等[13]發現,大多數HOXC4 結合位點與HOXC6 結合位點重疊。通過全基因組分析,發現HOXC6 與侵襲性前列腺癌的臨床預后指標相關。

ERG 是屬于ETS 轉錄因子家族,是紅細胞轉化特異性(ETS)家族成員,與胚胎發育、細胞增殖、分化、血管生成、炎癥和凋亡相關[18]。 ERG 家族基因與致癌基因融合有關,是前列腺癌在內的多種癌癥的典型特征[19]-[21]。 Eerola SK 等[22]通過分析mRNA 表達和染色質免疫沉淀測序(ChIP-Seq)數據集,發現PIM(Proviras Integration of Maloney Kinase)與MYC 原癌基因和ERG癌蛋白在前列腺癌的發生和進展中具有協同作用。

TOP2A 作為拓撲異構酶Ⅱ(TOP2)家族的兩個成員之一,僅在循環細胞中表達,并負責DNA 鏈復制過程中的酶解偶聯[23]。 據報道,在結腸癌組織中TOP2A 過度表達[24]。 然而,TOP2A 在組織發育和腫瘤進展中的作用和機制仍然未知[25]。 在根治性前列腺切除術組織中,TOP2A 與細胞周期進程(cell cycle progression,mCCP)的相關性最顯著(r=0.7)[26]。 Cattrini C 等[27]還發現TOP2A mRNA 過度表達與患者預后不良相關。

綜上,本研究通過對前列腺癌差異表達基因的篩選,建立蛋白質-蛋白質相互作用網絡,探究關鍵基因參與的生物學過程和代謝途徑,結合前列腺癌的臨床數據,評價PCA3 和TOP2A 對前列腺癌臨床預后的關系。 由于基因測序結果、篩選閾值、實驗方法、數據庫選擇等多種因素的影響,得出的結論存在偏差,需要進一步細胞、動物以及大樣本的人體實驗的驗證。