游 婧，傅 瑜，單 杰，張 杰，林曉虎，江裕程，卞成玥，羅頤辰

成釉細(xì)胞瘤(ameloblastoma)是多發(fā)于頜骨內(nèi)的較為常見的牙源性上皮性腫瘤，主要臨床表現(xiàn)為頜骨進(jìn)行性、無痛性膨大，由于其具有局部侵襲性，能夠?qū)е禄颊哳M面部畸形[1-2]。目前，對(duì)于成釉細(xì)胞瘤的主要治療方式為外科手術(shù)，但術(shù)后復(fù)發(fā)率高并存在惡變風(fēng)險(xiǎn)，影響患者的生存質(zhì)量[3-4]。因此，深入研究成釉細(xì)胞瘤的發(fā)生發(fā)展機(jī)制，探索有效的治療手段具有重要意義。

lncRNA是一類長(zhǎng)度大于200個(gè)核苷酸的非編碼RNA，不具有或很少具有蛋白編碼功能，近年來研究顯示lncRNA在細(xì)胞周期、細(xì)胞分化、表觀遺傳等生物學(xué)過程中發(fā)揮重要調(diào)控作用[5-7]。lncRNA及其調(diào)控基因的異常表達(dá)可能與多種疾病的發(fā)生發(fā)展有關(guān)[8]。Davanian等[9]的研究首次報(bào)道成釉細(xì)胞瘤與正常組織中l(wèi)ncRNA的表達(dá)差異以及l(fā)ncRNA的表達(dá)與成釉細(xì)胞瘤瘤體大小的關(guān)系，但關(guān)于lncRNA及相關(guān)mRNA在成釉細(xì)胞瘤發(fā)生發(fā)展中的作用尚未深入分析。本研究通過對(duì)公共數(shù)據(jù)庫基因表達(dá)數(shù)據(jù)庫(gene expression omnibus, GEO)中的兩組成釉細(xì)胞瘤患者組織的芯片數(shù)據(jù)集GES38494及GES132472進(jìn)行合并分析，篩選差異表達(dá)的基因及l(fā)ncRNA并構(gòu)建lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)，對(duì)異常表達(dá)的lncRNA及mRNA進(jìn)行生物信息學(xué)分析，尋找與成釉細(xì)胞瘤發(fā)生發(fā)展相關(guān)的靶基因，為成釉細(xì)胞瘤的分子機(jī)制研究及臨床治療提供理論基礎(chǔ)。

1 資料與方法

1.1 研究對(duì)象

本研究所選成釉細(xì)胞瘤患者組織芯片表達(dá)譜數(shù)據(jù)集GES38494及GES132472下載自美國(guó)國(guó)立生物信息中心(National Center for Biotechnology Information, NCBI)公共數(shù)據(jù)庫基因表達(dá)數(shù)據(jù)庫(Gene Expression Omnibus, GEO)，芯片分析均取材于手術(shù)中的相關(guān)組織，芯片數(shù)據(jù)GES38494由Heikinheimo等研究者上傳[10]，采用Affymetrix Human Genome U133 Plus 2.0 Array芯片平臺(tái)，芯片數(shù)據(jù)包括15個(gè)成釉細(xì)胞瘤組織樣本及4個(gè)正常組織樣本；GES132472由Kondo等研究者上傳[11]，采用Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381芯片平臺(tái)，芯片數(shù)據(jù)包括8個(gè)成釉細(xì)胞瘤組織樣本及8個(gè)正常組織樣本。

1.2 分析方法

1.2.1 數(shù)據(jù)預(yù)處理 使用R-studio軟件對(duì)芯片數(shù)據(jù)進(jìn)行分析。應(yīng)用R語言“affyPLM”及“affy”程序包進(jìn)行背景矯正，RMA算法對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化預(yù)處理，合并兩數(shù)據(jù)集并應(yīng)用“sva” 程序包的ComBat算法對(duì)合并的數(shù)據(jù)集進(jìn)行處理以矯正因?qū)嶒?yàn)儀器、人員、試劑等因素造成的批次效應(yīng)。

1.2.2 篩選差異表達(dá)基因及l(fā)ncRNA 應(yīng)用R語言“l(fā)imma”程序包對(duì)經(jīng)過預(yù)處理的數(shù)據(jù)集進(jìn)行分析，以改變倍數(shù)在2倍以上(|Fold Change| ≥2)且矯正后P值(adj.P. Value)< 0.05為篩選條件，得到差異表達(dá)基因。應(yīng)用R語言“ggplot2”程序包繪制火山圖，選取排名前100的差異表達(dá)基因應(yīng)用“pheatmap”程序包繪制聚類分析圖。在GENCODE數(shù)據(jù)庫(https://www.gencodegenes.org)下載lncRNA注釋包，應(yīng)用R語言“AnnotationHub”程序包進(jìn)行ID轉(zhuǎn)換，獲得差異表達(dá)基因中的lncRNA，即差異表達(dá)lncRNA。

1.2.3 lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)構(gòu)建 對(duì)差異表達(dá)lncRNA與mRNA的表達(dá)量進(jìn)行皮爾森相關(guān)分析，選取皮爾森相關(guān)系數(shù)(Pearson correlation coefficient，PCC)的絕對(duì)值(|PCC|)>0.85且P<0.05的lncRNA-mRNA對(duì)，采用Cytosacape軟件繪制lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)圖。

1.2.4 GO及KEGG通路富集分析 應(yīng)用R語言“clusterProfiler”程序包對(duì)lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中的mRNA進(jìn)行GO及KEGG通路富集分析，以P<0.05為標(biāo)準(zhǔn)進(jìn)行篩選。

1.2.5 GESA信號(hào)通路富集分析 以應(yīng)用R語言“clusterProfiler”程序包對(duì)差異表達(dá)基因進(jìn)行GSEA分析，以標(biāo)準(zhǔn)化P<0.05、FDR<0.25為標(biāo)準(zhǔn)篩選陽性基因集并使用“ridgeplot”及“gseaplot2”程序包進(jìn)行可視化分析。

1.2.6 PPI網(wǎng)絡(luò)分析 為評(píng)估lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中的mRNA之間的交互關(guān)系，采用STRING(https://string-db.org/cgi)數(shù)據(jù)庫以confidence≥0.4 為閾值條件構(gòu)建蛋白質(zhì)相互作用(Protein-protein interaction, PPI)網(wǎng)絡(luò)，利用Cytosacape軟件進(jìn)行可視化分析，應(yīng)用cytoHubba插件獲得核心基因。

2 結(jié) 果

2.1 差異表達(dá)基因的篩選結(jié)果

通過對(duì)合并數(shù)據(jù)信息進(jìn)行處理，共有17 631個(gè)基因被納入后續(xù)分析。以|Fold Change|≥2且adj.P. Value<0.05為標(biāo)準(zhǔn)，共篩選出801個(gè)差異表達(dá)的基因，其中416個(gè)在成釉細(xì)胞瘤組織中表達(dá)上調(diào)，385個(gè)表達(dá)下調(diào)(圖1A)。取排名前100的差異基因繪制熱圖(圖1B)。

圖1 差異表達(dá)基因的火山圖(A)及聚類分析圖(B)

2.2 lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)分析

通過對(duì)差異表達(dá)基因進(jìn)行l(wèi)ncRNA重注釋發(fā)現(xiàn)其中包含6個(gè)有意義的lncRNA，分別為SOX2-AS1、HCG22、LOC100506098、MAGI2-AS3、MEG3、MIR100HG，篩選與差異表達(dá)lncRNA有關(guān)聯(lián)的mRNA，發(fā)現(xiàn)283個(gè)mRNA與6個(gè)lncRNA具有高度共表達(dá)關(guān)系。構(gòu)建lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)進(jìn)行分析，方框節(jié)點(diǎn)代表lncRNA，圓形節(jié)點(diǎn)代表mRNA，6個(gè)lncRNA與283個(gè)mRNA共組成了305個(gè)相關(guān)鏈接，其中160個(gè)呈正相關(guān)，145個(gè)呈負(fù)相關(guān)。MAGI2-AS3居于共表達(dá)網(wǎng)絡(luò)的核心位置，與236個(gè)mRNA具有共表達(dá)關(guān)系，且關(guān)聯(lián)度高的mRNA多呈正相關(guān)。23個(gè)mRNA分別與2個(gè)lncRNA相關(guān)聯(lián)(圖2)。

方框代表lncRNA，圓形代表mRNA，直線兩端連接的mRNA與lncRNA存在共表達(dá)關(guān)系。實(shí)線表示正相關(guān)，虛線表示負(fù)相關(guān)。紅色代表在成釉細(xì)胞瘤中高表達(dá)，綠色代表低表達(dá)。方框的大小代表lncRNA在成釉細(xì)胞瘤中的表達(dá)程度，圓形的大小代表共表達(dá)關(guān)聯(lián)度

2.3 GO及KEGG通路富集分析

GO富集分析結(jié)果顯示lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中的mRNA參與的生物過程主要富集于上皮發(fā)育(GO:0008544)、皮膚發(fā)育(GO:0043588)、表皮細(xì)胞分化(GO:0009913)、角化細(xì)胞分化(GO:0030216)、角化作用(GO:0031424)等，細(xì)胞組成主要定位于含膠原的細(xì)胞外基質(zhì)(GO:0062023)、角化細(xì)胞套膜(GO:0001533)、膜的錨定成分(GO:0031225)、頂端結(jié)合復(fù)合體(GO:0043296)、細(xì)胞橋粒(GO:0030057)等，顯著富集的分子功能包括肽酶調(diào)節(jié)因子活性(GO:0061134)、肽鏈內(nèi)切酶抑制因子活性(GO:0004866)、肽酶抑制因子活性(GO:0030414)、肽鏈內(nèi)切酶調(diào)節(jié)因子活性(GO:0061135)等(圖3A)。

KEGG通路富集分析發(fā)現(xiàn)lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中的mRNA主要參與干細(xì)胞多能性調(diào)節(jié)通路(hsa04550)、胰腺分泌(hsa04972)、視黃醇代謝(hsa00830)、急性骨髓性白血病(hsa05221)、糖酵解和糖異生(hsa00010)等通路(圖3B)。

圖3 Go分析(A)和KEGG通路富集分析(B)

2.4 GSEA富集分析

采用GSEA對(duì)差異表達(dá)基因進(jìn)行KEGG通路富集分析，發(fā)現(xiàn)共有13個(gè)通路與成釉細(xì)胞瘤相關(guān)，選取富集分?jǐn)?shù)排名前十位的通路分析(圖4A)，其中代謝通路(Metabolic pathways)與成釉細(xì)胞瘤成正相關(guān)(P=0.001)，其他通路則為負(fù)相關(guān)。按照富集分?jǐn)?shù)進(jìn)行排序，軸突導(dǎo)向(Axon guidance)、細(xì)胞外基質(zhì)-受體相互作用(ECM-receptor interaction)以及黏著斑(Focal adhesion)相關(guān)通路(P=0.003)為排名前三位的通路(圖4B)。

A：KEGG通路富集；B：上圖：正相關(guān)通路；下圖：以富集分?jǐn)?shù)排名前三位的負(fù)相關(guān)通路

2.5 蛋白質(zhì)-蛋白質(zhì)相互作用(protein-protein interaction, PPI)分析

基于String數(shù)據(jù)庫中的信息進(jìn)行PPI分析，包含280個(gè)節(jié)點(diǎn)，280條連接線，平均每個(gè)節(jié)點(diǎn)的度為2。利用Cytoscape進(jìn)行可視化分析(圖4)，對(duì)關(guān)鍵基因進(jìn)行篩選排序，篩選出節(jié)點(diǎn)連接度(degree)排名前10位的核心基因，其中連接最廣泛的是外膜蛋白(involucrin，IVL)(表1)。

圖5 PPI分析圖

表1 連接度Top 10核心基因

3 討 論

lncRNA是在真核生物中發(fā)現(xiàn)的一類長(zhǎng)度大于200個(gè)核苷酸的RNA，多不具備編碼蛋白質(zhì)功能[12]，廣泛參與染色質(zhì)的修飾、轉(zhuǎn)錄、表觀遺傳調(diào)控以及免疫監(jiān)控等生物學(xué)過程，其轉(zhuǎn)錄與功能的失調(diào)可能是許多疾病發(fā)生發(fā)展的重要因素[13]。lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)的建立是研究lncRNA功能和調(diào)控機(jī)制的重要方法[14]。近年來越來越多的研究者將對(duì)成釉細(xì)胞瘤的研究聚焦于lncRNA對(duì)其生物學(xué)特征的調(diào)控[15-16]，但是關(guān)于lncRNA及相關(guān)mRNA在成釉細(xì)胞瘤發(fā)生發(fā)展中的作用尚未見深入分析。在本研究中，我們通過生物信息學(xué)方法分析研究了lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)的功能，為后續(xù)實(shí)驗(yàn)研究提供基礎(chǔ)。

研究表明lncRNA的調(diào)控機(jī)制主要包括以下四個(gè)方面：通過堿基互補(bǔ)配對(duì)的方式調(diào)控mRNA；吸附miRNA間接抑制miRNA對(duì)靶基因的調(diào)控作用；與轉(zhuǎn)錄因子結(jié)合形成復(fù)合體調(diào)控基因轉(zhuǎn)錄活性；改變DNA/RNA甲基化狀態(tài)、染色體結(jié)構(gòu)和修飾狀態(tài)控制相關(guān)基因表達(dá)[12-13，17]。本研究通過對(duì)差異表達(dá)基因進(jìn)行注釋篩選出6個(gè)差異表達(dá)的lncRNA，其中SOX2-AS1在成釉細(xì)胞瘤中表達(dá)上調(diào)最明顯。根據(jù)皮爾斯相關(guān)分析(|PCC|> 0.85)在差異基因表達(dá)譜中篩選出283個(gè)與差異表達(dá)lncRNA高度相關(guān)的mRNA，構(gòu)建lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)。在此網(wǎng)絡(luò)中我們發(fā)現(xiàn)一個(gè)lncRNA可能調(diào)控多個(gè)mRNA，而一個(gè)mRNA也可能與多個(gè)lncRNA相關(guān)，lncRNA及與其關(guān)聯(lián)度高的mRNA在成釉細(xì)胞瘤中的表達(dá)趨勢(shì)大多一致，推測(cè)lncRNA可能通過調(diào)控與其關(guān)聯(lián)的mRNA參與成釉細(xì)胞瘤的發(fā)生發(fā)展。然而，共表達(dá)網(wǎng)絡(luò)中l(wèi)ncRNA與mRNA的具體作用關(guān)系及相關(guān)機(jī)制仍有待進(jìn)一步研究。

目前研究認(rèn)為，成釉細(xì)胞瘤是一種上皮來源的牙源性腫瘤[18-19]。本研究中GO分析發(fā)現(xiàn)lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中的差異mRNA主要影響上皮發(fā)育、皮膚發(fā)育、表皮細(xì)胞分化等GO生物學(xué)功能，上皮發(fā)育等生物學(xué)功能與成釉細(xì)胞瘤的發(fā)生發(fā)展密切相關(guān)。KEGG通路富集分析顯示共表達(dá)網(wǎng)絡(luò)中的差異mRNA在干細(xì)胞多能性調(diào)節(jié)通路、胰腺分泌、視黃醇代謝等通路富集程度較高，而這些通路在成釉細(xì)胞瘤的發(fā)生機(jī)制中研究較少，為發(fā)病機(jī)制研究提供了新的方向?；蚣患治?gene set enrichment analysis，GSEA)是一種基于基因集的富集分析方法，其相對(duì)于GO/KEGG富集分析的優(yōu)勢(shì)在于能夠避免遺漏部分差異表達(dá)不顯著但卻具有重要生物學(xué)意義的基因[20-21]。本研究中，我們通過GSEA分析共表達(dá)網(wǎng)絡(luò)中的差異mRNA的KEGG通路富集情況，發(fā)現(xiàn)了13條與KEGG通路富集分析結(jié)果不同的通路信息，其中12條通路與成釉細(xì)胞瘤呈負(fù)相關(guān)，只有代謝通路(hsa01100)的基因集在成釉細(xì)胞瘤中表達(dá)上調(diào)。GSEA結(jié)果為成釉細(xì)胞瘤中異常表達(dá)基因富集通路分析進(jìn)行補(bǔ)充以及完善。

另外，本研究對(duì)lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)中差異表達(dá)的基因進(jìn)行蛋白質(zhì)互作網(wǎng)絡(luò)分析，并篩選出互作網(wǎng)絡(luò)中的核心基因，發(fā)現(xiàn)核心基因外膜蛋白(involucrin，IVL)連接度最高。IVL是角質(zhì)形成細(xì)胞終分化的早期指標(biāo)，在表皮角化包膜形成過程中發(fā)揮支架作用。有研究顯示，IVL與角化過度的手濕疹[22]、角化病[23]、接觸性皮炎[24]有關(guān)。該結(jié)果進(jìn)一步為成釉細(xì)胞瘤的治療提供了研究方向。

綜上所述，本研究利用生物信息學(xué)分析工具，篩選出成釉細(xì)胞瘤患者異常表達(dá)的基因及l(fā)ncRNA，構(gòu)建lncRNA-mRNA共表達(dá)網(wǎng)絡(luò)，對(duì)網(wǎng)絡(luò)中的mRNA參與的生物過程及通路進(jìn)行分析，并根據(jù)網(wǎng)絡(luò)中的差異表達(dá)基因譜分析相關(guān)的蛋白相互作用關(guān)系，列出連接度最高的10個(gè)核心基因，為成釉細(xì)胞瘤發(fā)生發(fā)展的機(jī)制研究及臨床治療靶基因的選擇提供理論基礎(chǔ)。