成釉細胞瘤患者lncRNA-mRNA共表達網絡分析及作用研究

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

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

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

1 資料與方法

1.1 研究對象

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

1.2 分析方法

1.2.1 數(shù)據預處理 使用R-studio軟件對芯片數(shù)據進行分析。應用R語言“affyPLM”及“affy”程序包進行背景矯正，RMA算法對數(shù)據進行標準化預處理，合并兩數(shù)據集并應用“sva” 程序包的ComBat算法對合并的數(shù)據集進行處理以矯正因實驗儀器、人員、試劑等因素造成的批次效應。

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

1.2.3 lncRNA-mRNA共表達網絡構建 對差異表達lncRNA與mRNA的表達量進行皮爾森相關分析，選取皮爾森相關系數(shù)(Pearson correlation coefficient，PCC)的絕對值(|PCC|)>0.85且P<0.05的lncRNA-mRNA對，采用Cytosacape軟件繪制lncRNA-mRNA共表達網絡圖。

1.2.4 GO及KEGG通路富集分析 應用R語言“clusterProfiler”程序包對lncRNA-mRNA共表達網絡中的mRNA進行GO及KEGG通路富集分析，以P<0.05為標準進行篩選。

1.2.5 GESA信號通路富集分析 以應用R語言“clusterProfiler”程序包對差異表達基因進行GSEA分析，以標準化P<0.05、FDR<0.25為標準篩選陽性基因集并使用“ridgeplot”及“gseaplot2”程序包進行可視化分析。

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

2 結 果

2.1 差異表達基因的篩選結果

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

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

2.2 lncRNA-mRNA共表達網絡分析

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

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

2.3 GO及KEGG通路富集分析

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

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

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

2.4 GSEA富集分析

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

A：KEGG通路富集；B：上圖：正相關通路；下圖：以富集分數(shù)排名前三位的負相關通路

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

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

圖5 PPI分析圖

表1 連接度Top 10核心基因

3 討 論

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

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

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

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

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