基于R語言的前列腺癌樣本的關鍵基因數據挖掘

孫澤坤 袁錢圖 胡建新

摘 要：為尋找前列腺癌組織與正常前列腺組織的關鍵基因，從Gene Expression Omnibus（GEO）數據庫下載前列腺癌樣本基因表達譜數據集GSE69223。進行芯片標準化處理后設置閾值|log2（FC）|>2且pvalue<0.05篩選出差異表達的基因，選擇其中高表達的41個基因進行GO和KEGG分析，得出8個關鍵基因：FFAR2、THBS4、TRPM4、CLDN3、CLDN8、HPN、PLA2G2A和FOLH1基因。再經UALCAN生存分析的到3個上調后患者生存可能性降低的基因：FFAR2、HPN和FOLH1。得出的8個關鍵基因主要富集在細胞趨化性、細胞-細胞連接、脂肪酸代謝等通路，這些通路與前列腺癌的發生發展有著密切聯系。除文獻已經報道的與前列腺癌有密切聯系的基因外，研究推測：CLDN3、CLDN8和FFAR2基因可能與前列腺癌特別是處于T2、T3分期的前列腺癌有著潛在的聯系。

關鍵詞： R語言;數據挖掘;前列腺癌;關鍵基因

【Abstract】 To search for key genes in Prostate Cancer tissues and normal prostate tissues， the gene expression profile data set GSE69223 of Prostate Cancer samples is downloaded from the Gene Expression Omnibus （GEO） database. After the chip standardization treatment， the threshold ∣log2（FC）∣>2 and pvalue<0.05 are used to screen out the differentially expressed genes， and 41 genes with high expression are selected for GO and KEGG analysis， and 8 key genes are obtained： FFAR2， THBS4， TRPM4， CLDN3， CLDN8， HPN， PLA2G2A and FOLH1 genes. After UALCAN survival analysis， the genes with reduced survival probability after three up-regulations are： FFAR2， HPN and FOLH1. The eight key genes are mainly enriched in cell chemotaxis， cell-cell junction， and fatty acid metabolism. These pathways are closely related to the development of Prostate Cancer. In addition to the genes already reported in the literature that are closely related to Prostate Cancer， it is hypothesized that the CLDN3， CLDN8， and FFAR2 genes may be potentially associated with Prostate Cancer， particularly Prostate Cancer at T2 and T3 stages.

【Key words】 ?R language; data mining; Prostate Cancer; key genes

0 引 言

前列腺癌（Prostate Cancer，PCa）是男性常見惡性腫瘤之一，多發于老年男性，同時具有高轉移性，且早期沒有明顯癥狀，發現可能已經是晚期[1]。據美國癌癥協會估計，2018年美國有大約164 690例新發PCa病例。同年大約有29 430例死于該病，這使其在世界致癌誘因統計榜單中已排至第二位[2]。與大多數其它癌癥一樣，PCa病情的發展取決于其擴散，因此局部疾病患者的5年生存率幾乎為100%，癌癥轉移患者的生存率將下降至28%[3]。中國前列腺癌發病率雖遠低于歐美國家，但隨著中國社會老齡化程度的逐漸提高、飲食結構及生活習慣的不斷改變、診療水平及生產工藝的亟待改進等因素，中國前列腺癌的發病率也有逐年上升的趨勢[4]。

研究可知，R語言是由Ihaka和 Gentleman教授聯合開發的一種計算機語言[5]，現已經主要應用于數據處理、統計計算、數學建模、數據可視化等多個領域，是一款開源、免費、自由的面向對象的編程軟件，并已擁有Linux、（Mac）OS X、Windows等多個版本。R語言使用的拓展包（packages）可根據用戶需要自由開發，同時還可供使用者免費下載[6]。隨著計算機技術及高通量測序技術的發展，生物芯片已然成為臨床樣本分析的一種有效方法，為疾病預測、分子診斷、新藥開發發揮著強有力的助益作用[7-8]。本研究采用了基于R語言的芯片分析方法來研究前列腺癌與正常前列腺組織之間的基因差異，從GEO數據（https：//www.ncbi.nlm.nih.gov/geo/）下載基因表達譜數據集GSE69223后對樣本進行質量檢測，數據清洗后設定閾值∣log2（FC）∣>2， pvalue<0.05，篩選出差異表達基因（FC：fold change 基因倍數變化），對其中的上調基因進行KEGG和GO分析以及UALACN（http：//ualcan.path.uab.edu/）生存分析，從而發現了一些前列腺癌、特別是處于T2、T3分期的前列腺癌的關鍵基因，對研究前列腺癌的分子診斷、抗前列腺癌藥物候選靶點提供了有益參考。

1 材料與方法

1.1 材料

芯片數據集GSE69223及芯片平臺數據GPL570從GEO數據庫（https：//www.ncbi.nlm.nih.gov/geo/）下載得到，R語言版本為R3.6。除內置程序包外，其余拓展包下載自https：//cran.r-project.org/及http：//bioconductor.org/packages。

1.2 實驗方法

1.2.1 數據獲取及數據清洗

GSE69223基因表達譜芯片數據由美國Affymetrix公司制作，使用芯片平臺為GPL570。數據集GSE69223包括15個正常前列腺組織樣本以及15個前列腺癌組織樣本。下載txt格式的原始數據，使用R語言獲取表達矩陣、分組信息、表型數據，過濾掉沒有基因名對應的探針以及對應某個基因名的多個探針。

1.2.2 聚類分析和PCA分析

使用R語言中的dist和hclust函數對30個樣品進行聚類分析，初步判斷15個正常樣本與15個前列腺癌樣本的差異，用以檢測該數據集是否具有數據挖掘的潛力。再對樣本進行主成分分析（PCA），用以判斷是否有潛在因子影響兩者之間的差異性。

1.2.3 獲得表達差異基因

用T檢驗獲得包含基因名、LogFC、pvalue等信息的數據框，以∣log2（FC）∣>2，pvalue<0.05為閾值篩選出差異基因，并規定LogFC>2為上調，LogFC<-2為下調。

1.2.4 差異表達基因的KEGG分析和GO分析

使用R語言中的clusterProfiler包對差異表達基因中的上調基因進行KEGG分析和GO富集分析。找出該基因的功能和富集的KEGG信號通路等信息。

1.2.5 生存分析

將經KEGG分析和GO分析的上調差異基因上傳到UALACN（http：//ualcan.path.uab.edu/），選擇prostate adenocarcinoma（前列腺腺癌）進行生存分析，獲得差異基因與生存時間之間的關系。

2 結果與分析

2.1 數據獲取及數據清洗

數據集包括15個正常前列腺樣本以及15個前列腺癌樣本的、共54 675個基因。通過數據清洗及標準化過程，可得與探針具有一一對應關系的基因有23 521個。為檢驗基因表達量的準確性，研究繪制了管家基因（GAPDH）以及β-actin的箱型圖（見圖1（a）），發現兩者的表達量平均值都在0附近，這表明此數據集中的基因表達未出現異常，在誤差允許范圍內可進行后續分析。將30個樣本納入分析范圍，圖1（b）展示了各樣本中基因的表達情況。

2.2 聚類分析和PCA分析

為初步判斷30個樣本中的前列腺正常樣本（normal）與前列腺癌樣本的差異，研究對樣本進行了聚類分析和PCA分析。分析結果表明，樣本中的某些基因的差異表達，可作為前列腺癌的診斷依據。在此次聚類分析中，有10個正常樣本與前列腺癌樣本分開，準確度達到66.7%，但此數據集中樣本總量為30個，分組數據較少，用聚類分析只能初步揭示正常樣品與前列腺癌樣品具有差異性（見圖1（c））。進一步地，對樣品進行PCA分析。結果表明，主成分1對樣本差異性的貢獻率為11.44%，主成分2對樣本差異性的貢獻率為9.87%，通過主成分1（PC1）和主成分2（PC2）可以將前列腺正常樣本與前列腺癌樣本較好的進行區分（見圖1（d））。

2.3 獲得表達差異基因

通過T檢驗，得到包含基因名、log2（FC）以及pvalue的數據框，設定閾值pvalue<0.05，log2（FC）>2以及log2（FC）<-2的基因，并規定log2（FC）>2的基因為上調基因，log2（FC）<-2的基因為下調基因，得到101個下調基因和41個上調基因（見圖2（a））。選擇前列腺癌較正常前列腺組織中的上調基因41個，導出其基因名及pvalue詳見表1。

2.4 KEGG和GO分析

將得到的41個基因利用超幾何分布原理在KEGG和GO數據庫中進行比對，得到富集結果見表2、表3。GO分析結果表明，前列腺癌細胞較前列腺正常細胞上調的差異基因主要富集的細胞活動過程有：白細胞遷移、細胞趨化性、細胞粘附、外肽酶活性、細胞 - 細胞連接。FFAR2、THBS4、TRPM4、CLDN3、CLDN8以及HPN被富集到多條通路。FFAR2、HBS4和TRPM4三個基因與白細胞遷移與細胞趨化性有關，白細胞遷移，可能導致前列腺癌組織中的白細胞增多，白細胞產生白介素，調控多種生理生化反應。該樣本中前列腺癌樣本集中于T2、T3分期，該時期的前列腺癌存在轉移潛能，因此可能與細胞趨化性有關。CLDN3和CLDN8是Claudin家族基因，該基因編碼的蛋白由Shoichiro Tsukita及其同事在1998年發現，是細胞緊密連接的重要分子，已有報道稱Claudin-1在結腸癌、Claudin-10在肝細胞癌、Claudin-18在胃癌中具有一定的臨床價值[9-11]。HPN基因又叫Hepsin基因，該基因編碼一種II型跨膜絲氨酸蛋白酶，該蛋白酶可能參與多種細胞功能，包括凝血和維持細胞形態。編碼蛋白的表達與癌癥，尤其是前列腺癌的生長和發展有關[12]。KEGG富集分析結果表明，差異表達的基因主要集中在緊密連接信號通路、多種生物分子代謝信號通路（在此列舉一條α-亞麻酸代謝信號通路）、細胞粘附分子（CAMs）信號通路、黏著力信號通路、維生素消化吸收信號通路。除GO分析結果涉及的基因外，KEGG分析中還出現了PLA2G2A、FOLH1兩個基因。PLA2G2A基因編碼的蛋白是磷脂酶A2家族（PLA 2）的成員。該基因產物屬于Ⅱ類，含有分泌型PLA 2，這是一種低分子質量的胞外酶，需要鈣離子進行催化。也可催化磷酸甘油中sn-2脂肪酸酰基酯鍵的水解，釋放游離脂肪酸和溶血磷脂，并參與生物膜磷脂代謝的調控[13]。同時，通過富集的結果來看，該基因還參與其他生物大分子如亞油酸代謝、脂肪消化吸收、醚脂代謝、花生四烯酸代謝、甘油磷脂代謝。而脂肪酸的氧化代謝過程已被證實與前列腺癌的發生和發展有著密切聯系[14]。FOLH1基因編碼屬于M28肽酶家族的Ⅱ型跨膜糖蛋白。該蛋白以谷氨酸羧肽酶的形式存在于不同的替代底物上，包括營養葉酸和神經肽N-乙酰-1-天冬氨酰-1-谷氨酸，在前列腺、中樞神經、外周神經系統和腎臟等多種組織中均有表達。在前列腺中，該基因編碼的蛋白質（PSMA）在癌細胞中被上調，并被用作前列腺癌的有效診斷和預后指標[15]。

2.5 生存分析

在UALCAN得到的生存分析結果中，研究發現，在候選的8個基因中，有5個基因的高表達組的生存可能較高，而FFAR2、FOLH1、HPN高表達組的生存可能性較低（見圖2（b）～（d））。其中，已經有文獻報道HPV編碼的蛋白與前列腺癌有關[12]，FOLH1基因編碼的蛋白已成為前列腺癌的腫瘤標志物[11]，而在相同的數據庫相同樣本的情況下，FFAR2組的P值最小，差異最為顯著，因此研究推斷，FFAR2基因與前列腺癌有較大關聯性。

3 結束語

R語言作為一種操作簡單、免費、開源的編程語言，適用于多種操作系統，為使用者提供了極大的方便。此次研究從GSE69223基因表達譜數據集中獲取了54 675個基因，設定閾值pvalue<0.05，log2（FC）>2，篩選出其中的41個上調基因，并對這些基因進行了KEGG分析和GO分析，獲得8個關鍵基因FFAR2、THBS4、TRPM4、CLDN3、CLDN8、HPN、PLA2G2A以及FOLH1。其中，FFAR2、THBS4、TRPM4三個基因與細胞趨化性相關，查看該數據集的表型數據發現，腫瘤樣本全部處于T2、T3時期，此3個基因的上調，印證了該分期的前列腺癌繼續發展可能轉移的事實。CLDN3和CLDN8屬于Claudin家族基因，該基因編碼的蛋白是細胞緊密連接的重要分子，已有報道稱Claudin-1在結腸癌、Claudin-10在肝細胞癌、Claudin-18在胃癌中具有一定的臨床價值，因此研究推測CLDN3和CLDN8兩個基因可能與前列腺癌有潛在聯系。PLA2G2A參與多種脂類大分子代謝，而脂肪酸的氧化代謝過程已被證實與前列腺癌的發生發展有著密切聯系。HPN和FOLH1已被文獻報道與前列腺癌有著密切聯系，并且FOLH1編碼的蛋白（PSMA）還被用作前列腺癌的腫瘤標志物，在前列腺癌的診斷和預后中起著不可替代的作用。通過生存分析，研究還發現這8個關鍵基因中，FFAR2、HPN以及FOLH1三個基因的高表達會減低患者生存可能性，除文獻已經報道的HPN核FOLH1基因外，本文再次經過分析推測后指出，FFAR2基因與前列腺癌的發生及發展有著潛在的關聯性。但要明確其具體機制，卻還需展開進一步研究。

參考文獻

[1]SHI Wei， DONG Li， BAO Junsheng. Progress in the studies of prostate cancer related molecules[J]. National Journal of Andrology， 2015， 21（4）：357-362.

[2]America Cancer Society. Cancer Information， Answers， and Hope[EB/OL]. https：//www.cancer.org/cancer/prostate-cancer/about/key statistics.html.

[3]MILLER K D， SIEGEL R L， LIN C C， et al. Cancer treatment and survivorship statistics，2016[J]. CA Cancer J Clin. 2016，66（4）：271-289.

[4]萬克松. 手術去勢間斷聯合抗雄激素藥物治療晚期前列腺癌臨床療效研究[D]. 廣州：南方醫科大學， 2012.

[5]IHAKA R， GENTLEMAN R. R： A language for data analysis and graphics[J].Journal of Computational and Graphical Statistics，1996， 5（3） ：299-314.

[6]韓俊偉， 智慧， 王宏， 等. R語言在生物信息實踐中的應用[J]. 生物技術世界， 2015（2）：180.

[7]李喜瑩， 李珊珊. 生物芯片技術及其在臨床檢驗醫學中的應用進展[J]. 分子診斷與治療雜志， 2011，3（1）：62-67.

[8]于穎彥. 生物芯片在胃癌藥物病理學研究中的先導作用[J]. 上海交通大學學報（醫學版）， 2007，27（5）：491-493.

[9]李東升， 王巍， 李晴， 等. 結腸癌組織中Her-2和Claudin-1的表達及意義[J]. 廣東醫學， 2012， 33（2）：237-239.

[10]張正東. Claudin-1和Claudin-10在肝細胞癌中的表達及意義[D]. 合肥：安徽醫科大學， 2011

[11]左忠林， 陳鵬， 陳小龍，等. Claudin-18在胃癌中的臨床表達關系與治療[J]. 中華臨床醫師雜志（電子版）， 2018， 12（3）：173-176.

[12]KIM H J ， HAN J H ， CHANG I H ， et al. Variants in the HEPSIN gene are associated with susceptibility to prostate cancer[J]. Prostate Cancer and Prostatic Diseases， 2012， 15（4）：353-358.

[13]洪雙雙. PLAG1和PLA2G2A在肝癌中的異常表達[D]. 鄭州：鄭州大學， 2011.

[14]姜偉. 復雜疾病特異的基因網路與microRNA-TF協同調控網絡的構建[D]. 哈爾濱：哈爾濱醫科大學，2008.

[15]MARAJ B H ， MARKHAM A F. Prostate-specific membrane antigen （FOLH1）： recent advances in characterising this putative prostate cancer gene[J]. Prostate Cancer and Prostatic Diseases， 1999， 2（4）：180-185.