基于R語言的前列腺癌樣本的關(guān)鍵基因數(shù)據(jù)挖掘

孫澤坤 袁錢圖 胡建新

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

關(guān)鍵詞： R語言;數(shù)據(jù)挖掘;前列腺癌;關(guān)鍵基因

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

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

1 材料與方法

1.1 材料

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

1.2 實(shí)驗(yàn)方法

1.2.1 數(shù)據(jù)獲取及數(shù)據(jù)清洗

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

1.2.2 聚類分析和PCA分析

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

1.2.3 獲得表達(dá)差異基因

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

1.2.4 差異表達(dá)基因的KEGG分析和GO分析

使用R語言中的clusterProfiler包對差異表達(dá)基因中的上調(diào)基因進(jìn)行KEGG分析和GO富集分析。找出該基因的功能和富集的KEGG信號通路等信息。

1.2.5 生存分析

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

2 結(jié)果與分析

2.1 數(shù)據(jù)獲取及數(shù)據(jù)清洗

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

2.2 聚類分析和PCA分析

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

2.3 獲得表達(dá)差異基因

通過T檢驗(yàn)，得到包含基因名、log2（FC）以及pvalue的數(shù)據(jù)框，設(shè)定閾值pvalue<0.05，log2（FC）>2以及l(fā)og2（FC）<-2的基因，并規(guī)定log2（FC）>2的基因?yàn)樯险{(diào)基因，log2（FC）<-2的基因?yàn)橄抡{(diào)基因，得到101個下調(diào)基因和41個上調(diào)基因（見圖2（a））。選擇前列腺癌較正常前列腺組織中的上調(diào)基因41個，導(dǎo)出其基因名及pvalue詳見表1。

2.4 KEGG和GO分析

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

2.5 生存分析

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

3 結(jié)束語

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

參考文獻(xiàn)

[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]萬克松. 手術(shù)去勢間斷聯(lián)合抗雄激素藥物治療晚期前列腺癌臨床療效研究[D]. 廣州：南方醫(yī)科大學(xué)， 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語言在生物信息實(shí)踐中的應(yīng)用[J]. 生物技術(shù)世界， 2015（2）：180.

[7]李喜瑩， 李珊珊. 生物芯片技術(shù)及其在臨床檢驗(yàn)醫(yī)學(xué)中的應(yīng)用進(jìn)展[J]. 分子診斷與治療雜志， 2011，3（1）：62-67.

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

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

[10]張正東. Claudin-1和Claudin-10在肝細(xì)胞癌中的表達(dá)及意義[D]. 合肥：安徽醫(yī)科大學(xué)， 2011

[11]左忠林， 陳鵬， 陳小龍，等. Claudin-18在胃癌中的臨床表達(dá)關(guān)系與治療[J]. 中華臨床醫(yī)師雜志（電子版）， 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á)[D]. 鄭州：鄭州大學(xué)， 2011.

[14]姜偉. 復(fù)雜疾病特異的基因網(wǎng)路與microRNA-TF協(xié)同調(diào)控網(wǎng)絡(luò)的構(gòu)建[D]. 哈爾濱：哈爾濱醫(yī)科大學(xué)，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.