兒童腎上腺皮質癌表達譜時序分析用于篩選相關靶點的研究

韓媛媛，李笑各，王書煥，劉戈力

腎上腺皮質癌（adrenocortical tumor，ACT）是一種罕見的內分泌腫瘤，發病率約5/10萬～20/10萬，ACT往往是致命的惡性腫瘤，5年生存率只有20%～35%［1］。ACT在兒童癌癥中的占比（1.3%）遠高于在成年人（0.02%～0.20%）［2］。ACT可以分為4個分級（stage）：stage 1和stage 2定義為腎上腺皮質瘤，stage 3會浸潤到周圍組織或局部區域淋巴結，stage 4則會表現出轉移性。在過去的10年中，通過DNA微陣列分析基因表達譜為許多實體瘤的分類和預后預測提供了一種新技術［3-4］。既往的基因表達譜研究發現，基因L13RA2、HTR2B、CCNB2、RARRES2和SLC16A9可用于區分ACT的良惡性［5］，且IGF-Ⅱ基因集的高表達和steroidogenesis基因集的低表達可作為惡性ACT的預測指標［6］。兒童ACT的MicroRNA表達譜研究則揭示了miR-99a和miR-100表達下調使得胰島素樣生長因子-Ⅰ受體（IFG-1R）、雷帕霉素靶蛋白（MTOR）和mTOR復合體調控相關蛋白（RPTOR）的過表達［7］。另外，可以根據有絲分裂活性和細胞周期調節基因的表達水平幫助確定ACT的亞型［8］。對于ACT生存期相關基因的研究表明，DLG7和PINK1聯合表達水平可預測無病生存情況，而BUB1B和PINK1聯合表達水平則可用來預測總體存活率［9］。絲裂原活化蛋白激酶相關基因和生長因子受體信號通路相關基因在兒童ACT中的表達下調，表明其可以一定程度上反映治療狀況［10］。雖然這些研究已經推動了這一領域的發展，然而目前有關ACT的致病機制仍不明確，尤其對兒童ACT生存期相關研究更少。因此，本研究選取了兒童ACT相關的基因芯片數據，對各個分級的兒童ACT進行研究，以期對兒童ACT獲得更多了解。

1 材料與方法

1.1 材料 從Gene Expression Omnibus（https://www.ncbi.nlm.nih.gov/geo/）選取了兒童ACT相關RNA芯片數據集GSE75415，該數據集基于Affymetrix Human Genome U133A Array（GPL96）芯片平臺進行信使RNA（mRNA）表達譜的測定。共包含31個樣本：5個腎上腺皮質腺瘤組織樣本、18個ACT組織樣本、1個未確定樣本及7個正常腎上腺皮質組織樣本。本研究只保留了18個ACT組織樣本（實驗組）以及7個正常腎上腺皮質組織樣本（對照組）。18個ACT樣本處于stage 1的樣本4個、stage 2的樣本6個、stage 3的樣本5個、stage 4的樣本3個。

1.2 方法

1.2.1 差異表達基因篩選 使用R語言Affy包中的RMA函數對各樣本表達量進行去背景噪聲，標準化，得到標準化表達值矩陣。引入limma包進行差異表達基因的篩選，以實驗組相對于對照組變化倍數（Fold Change，FC）的對數值log2|Fold Change|＞1以及校正值P＜0.05（Benjamini&Hochberg方法）作為選取差異表達基因的閾值。分別將18個ACT樣本的各分期的表達情況與對照組的7個正常腎上腺組織樣本的表達情況進行輪流比較，篩選各個分級相對于對照組的差異表達基因（DEGs），并取得各分級間DEGs的交疊部分作為重合差異基因（OLDEGs）。

1.2.2 功能富集分析 使用R語言的clusterProfile包對各分級的DEGs進行KEGG通路分析以及GO功能富集分析。以假陽性率FDR＜0.05作為選取富集功能項的閾值。此外，對各分級間的OLDEGs用DAVID（Database for Annotation,Visualization and Integrated Discovery,http://david.abcc.ncifcrf.gov）進行GO生物過程（Biological Process）富集分析，使用cytoscape中的enrichment Map插件進行OLDEGs的GO富集可視化結果以觀察各個生物過程之間的相互聯系。

1.2.3 中心基因篩選將各癌癥分級中的DEGs的log2（FC）值進行可視化，選取癌癥4個階段中均有同樣差異表現的DEGs作為中心基因，這些中心基因在各個癌癥階段中皆具有表達共性，可以作為ACT的特異性基因。基因表達水平高低通過log2（FC）的符號來判斷，當log2（FC）為正時，表示在ACT中高表達，log2（FC）為負時，表示在ACT中低表達。

1.2.4 生存分析使用R語言的cgdsr包下載cbioportal（http://www.cbioportal.org/）中腎上腺癌的TCGA數據集，其中包括79個ACT樣本的二代測序數據。引入Survival包對OLDEGs進行生存分析，檢驗方法為Log-rankχ2檢驗，保留表達量與生存時間存在顯著相關性的生存結果，P＜0.05為差異有統計學意義。

2 結果

2.1 差異表達基因篩選 相對于對照組的7個正常樣本，stage1～4樣本組分別存在248、334、315和561個DEGs，如圖1所示?？梢姼鞣旨墭颖窘M與對照組均聚類成獨立的2組，存在明顯差異表達。

4個分級樣本組存在73個OLDEGs，這73個OLDEGs在各分級樣本組表達量熱圖見圖2A；各個分級之間的差異表達基因重合情況見圖2B?？梢姼鞣旨壷g有明顯的表達差異，也有部分基因在各分級中表達情況一致。

2.2 功能富集分析 各分級間DEG的KEGG通路富集見圖3A，stage 1～4樣本組分別有2、7、6和7條通路富集。各分級間DEG的GO富集結果見圖3B，可見Stage 1～4樣本組分別在5、9、7和8個過程富集，其中molecular_function，protein binding和binding這3個過程在各個分級中均表現明顯富集。圖3C展示了各個分級間DEG的KEGG富集結果重合情況，其中各個分級均包含cell cycle這一通路。圖3D展示了各個分級間DEG的GO富集結果重合情況。

對73個OLDEGs的GO富集結果通路進行串話分析見圖4，所有的通路大概分為4個主要的功能模塊：免疫相關功能、細胞周期相關功能、磷酸化相關功能以及凝血相關功能。另外，對73個OLDEGs進行KEGG富集分析，得到46個富集通路，其中多條通路與代謝、感染和信號通路轉導有關，見表1。

2.3 中心基因篩選 對各個分級中的DEGs進行分析，篩選得到12個中心基因，其中HSPA13、GARS、STXBP1、AKIRIN1、TUBB3這5個基因在各分級中均上調，ADH1B、DCN、RASSF2、PDGFRA、PLAT、C3、FOS這7個基因在各分級中均下調，見表2。

2.4 生存分析 對TCGA數據集中79個ACT樣本的二代測序數據進行73個OLDEGs生存分析發現，9個基因與ACT的生存相關，它們的生存曲線見圖5。其中基因XPO1、RACGAP1、PDGFD、NR4A2和MXRA5的高表達會導致ACT患者的生存時間縮短；而VPS51、TMED3、NDFIP1和CDKN1C高表達的ACT患者生存時間延長（P＜0.05）。

Fig.3 Functional enrichment analysis圖3 功能富集

3 討論

數據集GSE75415曾經被用于兒童ACT的表達譜分析，以方便對兒童ACT診斷和預后的了解，但相關研究中缺乏相關DEGs的生存期分析［11］。不同于既往研究，本研究旨在通過對正常樣本和4個分級ACT的表達譜進行分析與比較，探索ACT的發病機制，并通過對相關基因的生存期分析，篩選與ACT治療相關的靶點。

Fig.4 Crosstalk analysis of GO terms of 73 OLDEGs圖4 73個OLDEGs顯著富集通路之間的通路串化分析

Fig.5 Kaplan-Meier curve of ACT survival associated genes圖5 ACT生存相關基因的生存分析

3.1 DEG篩選及功能富集分析 與對照組正常樣本相比，在stage1～4的ACT樣本中，分別有248、334、315和561個基因發生了差異表達。對stage1～4樣本組的DEG進行KEGG功能富集分析發現，各個分級均包含cell cycle、molecular_function、protein binding和binding通路。4個分級樣本組間存在73個OLDEGs，12個中心基因中HSPA13、GARS、STXBP1、AKIRIN1、TUBB3在4個分級中均上調，基因ADH1B、DCN、RASSF2、PDGFRA、PLAT、C3、FOS在4個分級中均下調。對73個OLDEGs進行富集功能分析發現，富集通路大多分布在免疫相關功能，細胞周期相關功能，磷酸化相關功能模塊以及凝血相關功能，且多條KEGG富集通路與代謝，感染和信號通路轉導有關。

Tab.1 Significantly enriched KEGG pathways of 73 OLDEGs表1 73個OLDEGs的KEGG通路富集結果

續表

Tab.2 The twelve core genes表2 12個中心基因

基因GARS、PLAT和C3與腫瘤免疫功能密切相關。GARS編碼的蛋白是人類自身免疫性疾病、多發性肌炎或皮肌炎的自身抗體的目標。C3基因編碼的補體C3在補體系統的激活中起著決定性作用。PLAT和C3均參與Complement and coagulation cascades途徑。Anurag等［12］研究表明，Complement and coagulation cascades途徑在ACT中下調，這可能與PLAT和C3的表達下調有關。因此，筆者推測GARS、PLAT和C3通過影響免疫相關功能影響ACT的發生與發展。

基因TUBB3、DCN、RASSF2和FOS與腫瘤細胞周期調控密不可分。TUBB3在腦腫瘤（85%～100%）、肺癌（35%～80%）、胰腺癌（50%）、腎細胞癌（15%～80%）和惡性黑色素瘤（77%）中均有高表達［13］。研究表明，TUBB3的過表達與侵襲性膀胱癌的遺傳不穩定性增加有關，而p53改變和細胞增殖迅速有關［14］。DCN低表達可促進腎臟癌細胞的增殖和轉移［15］。RASSF2在賁門癌、乳腺癌、食管鱗狀細胞癌、上皮性卵巢癌，尤因肉瘤和宮頸癌中是一種新的抑癌基因［16-21］。研究表明，c-fos的表達缺失與卵巢癌的腫瘤進展有關［22］。因此，筆者推測，ACT細胞周期相關功能的富集與TUBB3、DCN、RASSF2和FOS的差異表達有密切關系。另外，DCN蛋白還可通過與多種細胞表面受體結合，發揮抑制血管生成和腫瘤發生的作用。

基因PDGFRA下調可影響腫瘤細胞的磷酸化過程。癌細胞內若有多條信號傳導途徑發生變化，則多與蛋白質磷酸化有關。研究顯示，酪氨酸磷酸化是控制細胞增殖、凋亡、遷移和侵襲的幾個關鍵的致癌信號通路的一個精致開關［23］。PDGFRA編碼細胞表面酪氨酸激酶受體，在器官發育、傷口愈合和腫瘤進展中起著重要作用。Pastor等［24］對6例無血緣關系的ACT患者進行研究發現，PDGFRA在患者中呈低表達。另外，HSPA13編碼熱休克蛋白70（HSP70）超級蛋白，可能通過影響蛋白分泌、清除變性或不正常折疊蛋白等過程，從而參與ACT的發生和發展［12］。ADH1B編碼的醇脫氫酶則可能通過影響各種底物代謝通路和藥物代謝參與ACT的發?。?4］。

綜上所述，Binding相關功能在各個分級中富集可能與STXBP1和AKIRIN1過度表達相關。NCBI數據庫關于STXBP1和AKIRIN1的介紹中提到，STXBP1編碼一種合成結合蛋白，與各種類型的binding密切相關，而AKIRIN1的表達與protein binding密切相關。

3.2 生存分析 本研究生存分析顯示，XPO1、RACGAP1、PDGFD、NR4A2和MXRA5的高表達，以及VPS51、TMED3、NDFIP1和CDKN1C的低表達與ACT患者的生存期密切相關。

XPO1編碼的細胞周期調控蛋白可通過控制cyclin B、MPAK和MAPKAP激酶2的表達，參與部分細胞的生長過程。RACGAP1的高表達提示乳腺癌、胃癌和直腸癌的預后較差［25-27］，本研究亦顯示RACGAP1高表達的ACT患者預后不良。PDGFD基因參與Gap junction通路，相比正常細胞和良性細胞，ACT細胞間的gap junction數量顯著降低［28］。在大約60%的ACT患者均呈現出腎上腺類固醇激素過量現象，這可能與NR4A2編碼的激素-維生素A表達增加有關，且該受體參與Aldosterone synthesis and secretion通路［29］。另外有研究顯示，在腎小管細胞中MXRA5編碼蛋白表現出免疫活性［30］。而MXRA5編碼的蛋白包含12個免疫球蛋白樣的c2型域，因此筆者推測該基因可能通過影響免疫功能來發揮作用。

表達下調基因主要通過影響物質運輸與蛋白分泌影響預后。VPS51編碼的蛋白是高爾基體相關逆行蛋白復合物（GARP）的組成部分。TMED3編碼一種穿膜蛋白，參與從內質網到高爾基體的囊泡運輸。DNFIP1編碼的蛋白屬于一組具有3個跨膜結構域的進化保守蛋白，這種蛋白質被認為是高爾基體膜蛋白家族的一部分。目前這3個基因的作用機制還不清楚，根據3個基因的功能推測，它們應該與物質的運輸和蛋白分泌有關系。另外研究顯示，TMED3基因的敲除增加了直腸癌的轉移生長［31］，本研究亦顯示，TMED3的低表達導致ACT的預后不良。CDKN1C編碼的蛋白對于幾種cyclin/CDK復合物有很強的抑制作用，并抑制細胞增殖，在成人ACT中也檢測到該基因的表達有所下調［32］。

綜上所述，通過ACT和正常組織樣本的表達譜分析，本研究篩選出了與ACT發生和發展相關的靶點，并通過對重合差異基因的生存分析，確定了與ACT預后有密切關系的基因，這將有助于進一步了解ACT的致病機制、治療和預后。

Fig.1 Heatmap of expression values of the DEGs圖1 差異表達基因的表達量熱圖

Fig.2 Analysis of the 73 OLDEGs圖2 73個OLDEGs的表達量分析