RET/RAS/BRAF基因驅動甲狀腺乳頭狀癌共表達關鍵基因的生物信息學分析

仇以利，陳一峰，楊金玲，張宏圖

甲狀腺乳頭狀癌(papillary thyroid carcinoma, PTC)是甲狀腺最常見的惡性腫瘤，部分地區的發病率10年增長了近10倍[1]，胞外信號→RET/RAS/BRAF→MAPK轉導通路是PTC最常見的致癌機制，RET/RAS/BRAF這三個分子事件間互相排斥。當前腫瘤分子演進機制探索已從單基因單通路發展到多基因多通路。本實驗利用生物信息學方法分析RET/RAS/BRAF驅動PTC的共表達基因，從中尋找關鍵基因，旨在獲得更多PTC發生、發展的相關分子標志物，為PTC的精準危險度分層進一步拓寬有關基因變異，有利于針對放射性碘難治性或手術不能切除的PTC患者靶向藥物的研發、篩選及預期效果的評估。

1 材料與方法

1.1 數據集及樣本收集在基因表達綜合數據庫(Gene Expression Omnibus database, GEO)的網站(www.ncbi.nlm.nih.gov)搜索“papillary thyroid carcinoma AND normal AND Homo sapiens”，根據目的選取數據集GSE58545，芯片源于GPL96[HG-U133A] Affymetrix Human Genome U133A Array，其包括正常甲狀腺組織18例，PTC組織27例，后者包含BRAF基因突變者(BRAF+)18例、RET基因突變者(RET+)8例、RAS基因突變者(RAS+)1例。

1.2 差異表達基因篩選在GEO網站運用GEO2R對話框搜索所有樣品后，設定分組，包括BRAF+對比正常甲狀腺組織、RET+對比正常甲狀腺組織、RAS+對比正常甲狀腺組織三組，然后借助R及Limma包進行線上基因表達差異分析，導出所有結果。將2為底數取值差異倍數(fold change, FC)，以校正后的P值<0.01，|log2FC|≥2為標準篩選差異表達基因，把3組篩選所得基因名稱導入Draw Venn Diagrams在線工具(http://bioinformatics.psb.ugent.be/webtools/Venn/)，交集三組的共同表達基因。

1.3 基因本體(Gene Ontology, GO)功能分析和通路富集分析打開Metascape主頁(http://metascape.org)，在對話框提交所篩選的交集基因，在物種選項欄中選擇H.sapiens，點擊custom analysis進行自定義分析，包含GO功能、京都基因與基因組百科全書(Kyoto encyclopedia of genes and genomes, KEGG)信號通路、反應組學基因集(Reactome gene sets)和經典通路(Canonical Pathways)等富集分析，把Min Overlap、P Value Cutoff、Min Enrichment值分別設置為3、0.01、1.5，系統隨之劃分基因類別、構建網絡。

1.4 Cytoscape軟件和STRING網絡分析關鍵基因打開STRING(http：//www.string-db.org/)官方網站，在“SEARCH”對話框輸入所篩選的基因名稱，數據庫(Version 11.0)將可視化輸出富集網絡圖，保存TSV格式的結果文件。同時下載Cytoscape軟件及Cytohubba插件，把TSV文件導入軟件中，挖掘關鍵基因。

1.5 Oncomine數據庫分析關鍵基因的表達登陸Oncomine(www.oncomine.org)網站，在“search”對話框中輸入關鍵基因名稱，在“Primary Filters”目錄下根據需要依次點擊選擇具體類別，設置分析參數為：P-value<0.05，Fold change≥2，gene rank=top 10%，分別檢索上述每個基因在單個或多個PTC數據集的整合表達。

1.6 TCGA數據庫中驗證關鍵基因在PTC中的表達趨勢及與臨床病理的聯系登陸Ualcan TCGA analysis網頁(http：//ualcan. path. uab. edu/)，訪問甲狀腺癌基因表達譜數據，其中包括正常甲狀腺組織59例，PTC 505例(經典型358例、高細胞亞型36例、濾泡亞型102例、其他亞型9例)，進一步驗證關鍵基因在TCGA數據庫PTC中的差異表達情況，獲取關鍵基因表達的臨床病理因素。

1.7 HPA數據庫分析PTC中關鍵基因編碼蛋白的表達登陸HPA網站(https//www.proteinatlas.org/)，在搜索框中輸入關鍵基因的名稱，點擊腫瘤欄，查看在PTC中的免疫組化酶標表達情況，亦可顯示在細胞系的定位。

2 結果

2.1 差異表達基因篩選結果應用GEO2R工具，點擊value distribution，先確認芯片數據質量可靠，共有45例樣本，再以預先設定的閾值，BRAF+組18例、RET+組8例、RAS+組1例對比正常甲狀腺組織18例，分別篩選出差異表達基因310、824、280個，通過韋恩圖交集得出共表達基因34個(圖1)：KCNJ2、GDF15、MET、SOX4、COL1A1、DPP4、MUC1、RYR1、CITED1、TENM1、RUNX1、TYMS、GALNT7、PTP4A3、SCEL、NELL2、HMGA2、PDZK1IP1、RXRG、TUSC3、MDK、SPINT1、BHLHE40、TGFA、PHLDA2、ENTPD1、NRCAM、TBC1D2、RAB27A、LRP4、CLDN1、TESC、SCG5、PSD3。

圖1 差異表達基因交集的韋恩圖

2.2 差異表達基因的富集分析結果應用Metascape對34個共表達差異基因進行功能及通路分析，發現被顯著富集的相關基因與16條GO通路、1條KEGG通路、1條Reactome通路和1條Canonical通路相關，其中GO通路涉及生物過程調控功能14項、分子功能2項，按照P值從小到大，生物過程依次包含從皮膚的發育到細胞對非生物刺激的反應(圖2)。分子功能包含肝素及鈣離子的結合功能。KEGG、Canonical、Reactome通路分析分別顯示共表達基因參與PTC的基因轉錄失調、p53蛋白誘導死亡結構域的下游通路和MET信號轉導通路。

圖2 生物過程GO富集通路熱圖：顏色越深富集越顯著，P值越小

2.3 關鍵基因的篩選應用STRING數據庫導出共表達基因編碼蛋白質間的相互作用網絡圖(圖3)，應用Cytoscape軟件、Cytohubba插件，采用MCC算法，篩選出10個關鍵基因，按節點排名分為三組：第一組為MET；第二集團為DPP4、ENTPD1、MUC1、SPINT1；第三集團為BHLHE40、CLDN1、GALNT7、TGFA、TYMS(圖4)。

圖3 共表達基因編碼蛋白質互作圖

圖4 關鍵基因排序及蛋白質互作圖：紅色為第一組基因；磚紅色為第二集團基因；黃色為第三集團基因

2.4 Oncomine數據庫分析關鍵基因的表達應用Oncomine數據庫依次對10個關鍵基因進行多數據集整合檢索，庫中所有包含PTC的5個數據集均被納入研究，得出每個基因的中位秩及P值，這10個基因在PTC中的表達均顯著高于正常甲狀腺組織(P<0.01)。

2.5 基于TCGA數據庫驗證PTC中關鍵基因的表達應用Ualcan網站挖掘10個關鍵基因在TCGA數據庫中PTC中的差異表達情況，相比于正常組織，PTC中MET、DPP4、ENTPD1、MUC1、SPINT1、BHLHE40、CLDN1、GALNT7、TGFA、TYMS基因的mRAN表達顯著增加(P<0.01)。除了ENTPD1、TGFA、TYMS基因，其他關鍵基因mRNA的表達與PTC淋巴結轉移相關(P<0.01)。

2.6 HPA數據庫分析PTC中關鍵基因編碼蛋白的表達從HPA網站(www.proteinatlas.org)瀏覽10個關鍵基因編碼蛋白在正常甲狀腺組織和PTC癌細胞中的表達圖像，MET、DPP4、ENTPD1、MUC1、SPINT1蛋白在PTC中呈陽性，表達于PTC細胞質/膜；BHLHE40、CLDN1、GALNT7、TGFA、TYMS蛋白在PTC中呈陽性，其中BHLHE40在癌細胞核中弱至中等強度表達，余蛋白表達于細胞質/膜；而在正常甲狀腺組織中均未見陽性染色。

3 討論

近30年來全球甲狀腺癌的發病率增加了3倍[2]，在長期隨訪中發現超過25%的PTC患者發生復發[3]，RET、RAS和BRAF三者是PTC最常見的驅動基因[4]，cBioPortal數據顯示高達83%的PTC檢測出MAPK通路的異常活化，它們彼此間幾乎不發生重疊突變[5]。近70%PTC發生BRAF、RAS或RET基因突變：36%～83%PTC發生BRAF突變[5]；13%～25% PTC發生RET基因重排[6]，近年其發生率逐步降低，低至5%[7]；約13%PTC發生RAS基因突變[5]，近幾年其突變發生率增長較快，尤其在濾泡型PTC中的突變率達25%[8]。進入二十一世紀以來腫瘤基因和蛋白質分子數據的爆發性增長和相應臨床病理資料采集的規模化、集約化，隨著生物信息技術在解析醫學大數據的交融和滲透，為進一步拓展RET/RAS/BRAF基因驅動PTC的上下游相關分子，本組以包含18例正常甲狀腺組織、18例BRAF+、8例RET+、1例RAS+PTC組織的基因表達譜數據集芯片為研究對象[9]，篩選出共表達基因34個：KCNJ2、GDF15、MET、SOX4、COL1A1、DPP4、MUC1、RYR1、CITED1、TENM1、RUNX1、TYMS、GALNT7、PTP4A3、SCEL、NELL2、HMGA2、PDZK1IP1、RXRG、TUSC3、MDK、SPINT1、BHLHE40、TGFA、PHLDA2、ENTPD1、NRCAM、TBC1D2、RAB27A、LRP4、CLDN1、TESC、SCG5、PSD3。基于Metascape大數據全面分析和解釋組學研究結果[10]，其中GO分析工具標準化展現了這些歸集基因參與皮膚的發育、骨骼系統的發育、對類固醇激素的反應、整合素介導細胞黏附的調節、對酸性化學物質的反應、發育性生長、細胞周期阻滯的正向調控、體液水平的調節、進入宿主細胞、解剖結構的成熟、蛋白質的糖基化、激素分泌的調節、衰老、細胞對非生物刺激的反應共14項生物過程調控功能，還涉及肝素與鈣離子的結合功能兩項分子功能調控。而KEGG、Canonical、Reactome通路分析顯示它們分別與轉錄失調、p53蛋白調控免疫球蛋白死亡結構域、MET信號轉導事件有關。STRING v11.0是目前在線分析蛋白質間相互作用的最新版本[11]，在STRING互作網絡中共發現2個關鍵基因模塊，分別包含8個、2個基因，進一步應用Cytoscape軟件、Cytohubba插件，把這10個關鍵基因按節點排名分為三組：第一組為MET；第二集團為DPP4、ENTPD1、MUC1、SPINT1；第三集團為BHLHE40、CLDN1、GALNT7、TGFA、TYMS，從Oncomine數據庫可視化的直觀結果可以看出10個關鍵基因在PTC中的表達均顯著高于正常甲狀腺組織。

MET編碼產物通過PI3K-AKT-mTOR和RAS-RAF-MEK-ERK途徑參與細胞增殖、凋亡、遷移等過程[12-13]，作為本組篩出的首選關鍵基因，在其他的研究也證實MET基因mRNA及編碼蛋白在PTC中的表達均明顯高于良性病變組[14]。DPP4蛋白又稱CD26，是一種細胞表面的絲氨酸蛋白酶，其在PTC組織中高表達可促進細胞增殖、遷移和浸潤，PTC患者血液中可溶性CD26水平顯著高于甲狀腺良性病變者及健康對照者[15]。MUC1基因編碼的黏蛋白主要表達于正常上皮細胞近官腔或腺腔面，而在PTC細胞會出現異常表達[16]，MUC1蛋白高表達與淋巴結癌轉移、甲狀腺外侵行為的發生密切相關[17]。CLDN1基因編碼產物是一種細胞間緊密連接跨膜蛋白，TGFA基因編碼的生長因子，作為表皮生長因子受體的配體，它們均可參與PTC的遷移與侵襲性生長[18-19]。據報道61%～68%PTC高表達SPINT1，其實SPINT1能調節細胞去適應活性基質金屬酶的溶解作用，而且還能誘導缺氧條件下的細胞遷移、侵襲[20]。這充分說明了本組挖掘的MET、DPP4、MUC1、CLDN1、TGFA、SPINT1基因在PTC中的表達及其機制得到了相應研究報道的佐證。除此之外，本文發現ENTPD1、BHLHE40、GALNT7、TYMS與PTC的發生、發展有關，目前尚無類似的研究報道。從Ualcan網站挖掘這10個基因在TCGA數據庫PTC中的差異表達情況，也得出同樣的結果；在HPA數據庫檢索上述10個基因編碼蛋白在PTC中均呈高表達。此外本組發現除了ENTPD1、TGFA、TYMS基因，其他7個基因mRNA表達與PTC淋巴結轉移相關。

本實驗應用GEO、Metascape、Cytoscape、STRING、Oncomine、Ualcan及HPA公共數據庫綜合分析RET/RAS/BRAF基因驅動PTC的10個共表達關鍵基因，結果表明它們相對應的mRNA及蛋白在PTC組織中的表達均顯著高于正常甲狀腺組織，且是PTC淋巴結轉移相關的關鍵基因。經檢索Pubmed、萬方數據知識平臺等，目前在PTC的研究中尚未見ENTPD1、BHLHE40、GALNT7和TYMS分子的報道，這將為參與RET/RAS/BRAF基因驅動PTC的分子家庭注入新成員，它們有望成為難治性PTC患者的潛在分子靶點。