基于融合基因作用于橫紋肌肉瘤的miRNA-mRNA 調控網(wǎng)絡的生物信息學分析

趙志娟, 孟 蓮, 劉春霞,2

（1.石河子大學醫(yī)學院病理學系 石河子大學醫(yī)學院第一附屬醫(yī)院病理科，新疆 石河子 832002；2.廣州醫(yī)科大學附屬第二醫(yī)院病理科，廣東 廣州 510260）

橫紋肌肉瘤（rhabdomyosarcoma，RMS） 是兒童最常見的軟組織惡性腫瘤，主要包括腺泡狀RMS （alveolar RMS，ARMS）、 胚胎型RMS（embryonal RMS，ERMS） 和成人多形性RMS（polymorphic RMS，PRMS）。ARMS 和ERMS 分別占所有RMS 病例的20% 和60%，ARMS 經(jīng)常出現(xiàn)t （2； 13）（q35； q14） 或t （1； 13）（p36；q14） 染色體易位，分別產(chǎn)生配對盒基因3（paired box gene 3，PAX3） - 叉頭型轉錄因子（forkhead-type transcripition factor，F(xiàn)KHR） 或配對盒基因7 （paired box gene 7，PAX7）-FKHR融合基因［1］。在ARMS 中融合基因具有明確的診斷及預后價值，其代表一種具有明顯侵襲性的亞群，該亞群通常對常規(guī)治療無反應，且復發(fā)風險高［2］。

微小核糖核酸（microRNAs，miRNAs） 是小的核糖核酸分子，廣泛參與腫瘤細胞增殖、分化和凋亡的過程［3］。本課題組前期研究［4］顯示：miR-410-3p 高表達可以抑制RMS 細胞的增殖、侵襲和遷移，促進細胞凋亡，并影響上皮- 間質轉化（epithelial-mesenchymal transition，EMT） 進程，miRNA 與RMS 的發(fā)展和預后有密切關聯(lián)。miR-874 在RMS 中起抑癌作用，并通過調控鳥嘌呤核苷酸交換因子T （guanine nucleotide exchange factor T，GEFT） 抑制RMS 細胞的生長和轉移［5］。但miRNA 在RMS 中參與融合基因形成的相關研究較少。

本研究通過生物信息學方法尋找融合基因陽性和陰性RMS 組織中差異表達基因（differentially expressed genes，DEGs） 和差異表達miRNA，構建miRNA-mRNA 分子調控網(wǎng)絡，為RMS 的診斷、治療和預后提供新的思路。

1 資料與方法

1.1 資料來源從基因表達圖譜（Gene Expression Omnibus，GEO） 數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo/） 中下載GSE66533和GSE97553 數(shù)據(jù)集。其中GSE66533 為33 例融合基因陽性和25 例融合基因陰性RMS 組織的基因表達譜數(shù)據(jù)，GSE97553 為2 例融合基因陽性和7 例融合基因陰性RMS 異種移植模型組織的miRNA 微陣列數(shù)據(jù)。

1.2 篩選DEGs 和差異表達miRNA使用GEO2R 分析工具（https：//www.ncbi.nlm.nih.gov/geo/geo2r/） 分別對GSE66533 和GSE97553數(shù)據(jù)集進行分析，以P＜0.01 且|logFC|≥1 為標準篩選DEGs，以P＜0.05 且|logFC|≥3 為標準篩選差異表達miRNA。使用R 軟件ggplot2 和pheatmap程序包進行可視化分析，生成火山圖和熱圖。

1.3 miRNA 靶基因預測通過在線數(shù)據(jù)庫miRDB （http：//mirdb.org/）、 TargetScan（http：//www.targetscan.org/vert_71/） 和StarBase （http：//starbase.sysu.edu.cn/） 3 個平臺進行差異表達miRNA 靶基因預測，采用Bioinformatics Evolutionary Genomics （http：//bioinformatics.psb.ugent.be/webtools/Venn/） 將3 個平臺的預測結果與前面篩選的DEGs 取交集，得到目標基因。

1.4 基因本體（Gene Ontology，GO）功能和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析使用在線工具DAVID （https：//david.ncifcrf.gov/）對篩選的目標基因進行GO 功能和KEGG 通路富集分析，GO 功能分析包括生物學過程、細胞組成和分子功能，以P＜0.05 為有效。

1.5 蛋白-蛋白互作（protein-protein interaction，PPI）網(wǎng)絡構建和核心基因（hub 基因）篩選通過在線分析工具STRING （https：//string-db.org/） 預測并構建目標基因的PPI 網(wǎng)絡，將文件導入Cytoscape 3.6.1 軟件進行可視化，并使用插件cytoHubba 篩選Top10 hub 基因。

1.6 生存分析使用Kaplan-Meier Plotter 數(shù)據(jù)庫（http：//kmplot.com/analysis/） 評估Top10 hub基因與肉瘤患者總生存（overall survival，OS） 率的關系，繪制Kaplan-Meier 生存曲線。

1.7 統(tǒng)計學分析采用R Studio 軟件制圖并進行統(tǒng)計學分析。采用Kaplan-Meier 法（Log-rank 檢驗） 繪制生存曲線，進行生存分析。以P＜0.05 為差異有統(tǒng)計學意義。

2 結 果

2.1 融合基因陽性和陰性RMS 組織中DEGs 和差異表達miRNA對融合基因陽性和陰性RMS 組織的基因和miRNA 表達數(shù)據(jù)進行分析，以P＜0.01且|logFC|≥1 為標準 篩 出891 個DEGs ，其中512 個表達上調，379 個表達下調（圖1A 和B）。以P＜0.05 且|logFC|≥3 為標準篩出14 個差異表達miRNA （圖1C 和D），表達上調和表達下調各7 個（表1）。

表1 本研究篩選的差異表達miRNATab.1 Differentially expressed miRNA screened in this study

圖1 GSE66533 數(shù)據(jù)集中DEGs 和GSE97553 數(shù)據(jù)集中差異表達miRNAFig.1 DEGs in GSE66533 dataset and differentially expressed miRNA in GSE97553 dataset

2.2 差異表達miRNA 靶基因的預測采用TargetScan、miRDB 和StarBase 數(shù)據(jù)庫預測差異表達miRNA 的靶基因，并對3 個數(shù)據(jù)庫的結果進行重疊，共預測出1 654 個靶基因。與891 個DEGs 取交集得到115 個目標基因（圖2）。

圖2 DEGs 與預測的靶基因交集的韋恩圖Fig.2 Venn diagram of intersection of DEGs and predicted target genes

2.3 目標基因的GO 功能和KEGG 通路富集分析GO 功能富集分析顯示：目標基因富集于RNA聚合酶Ⅱ啟動子轉錄、細胞增殖的正向調節(jié)和多細胞生物的發(fā)展等生物學過程，富集于細胞質膜、細胞表面和轉錄因子復合體等細胞組成成分，富集于轉錄因子活性和DNA 序列特異性結合、蛋白結合和染色質結合等分子功能（圖3A～C）。KEGG 通路富集分析顯示：目標基因參與癌癥中的蛋白聚糖、黑色素生成和細胞外基質-受體相互作用等通路（圖3D）。

圖3 目標基因的GO 功能（A－C）和KEGG 通路（D）富集分析Fig.3 GO function(A-C) and KEGG pathway enrichment analysis(D) of target genes

2.4 PPI 網(wǎng)絡和hub 基因分析通過STRING 軟件繪制115 個目標基因的PPI 網(wǎng)絡，得到1 個包含74 個結點和81 條邊的PPI 網(wǎng)絡圖（圖4）。使用插件cytoHubba 篩選出Top10 hub 基因，分別為表皮生長因子受體4 （oncogene encoding human epidermal growth factor receptor 4，ERBB4）、受體型蛋白質酪氨酸磷酸酶D （protein tyrosine phosphate receptor type delta，PTPRD）、 胰島素受體底物1 （insulin receptor substrate 1，IRS1）、整合素金屬蛋白酶 10 （a disintegrin and metalloprotease 10，ADAM10）、 Yes 相關蛋白1（Yes associated protin1，YAP1）、轉錄因子AP-2A（transcription factor AP-2 alpha，TFAP2A）、細胞黏附分子1 （cell adhesion molecule 1，CADM1）、ELAV 類似RNA 結合蛋白2 （ELAV like RNA binding protein 2，ELAVL2）、 鋅指轉錄因子1（zinc finger protein 1，SNAI1） 和ERBB 受體反饋抑制劑1 （ERBB receptor feedback inhibitor 1，ERRFI1）（圖5A）。通過Cytoscape 軟件對hub 基因進行miRNA-mRNA 分子調控網(wǎng)絡的構建（圖5B）。

圖4 目標基因的PPI 網(wǎng)絡圖Fig.4 PPI network diagram of target genes

圖5 PPI 網(wǎng)絡的Top10 hub 基因（A）和對應的miRNA-mRNA 分子調控網(wǎng)絡（B）Fig.5 Top10 hub genes of PPI network（A）and corresponding miRNA-mRNA molecular regulation network（B）

2.5 Top10 hub 基因的生存分析因GEO 數(shù)據(jù)庫里無RMS 患者相關臨床信息，所以針對篩選出的Top10 hub 基因，在Kaplan-Meier Plotter 數(shù)據(jù)庫對其進行生存分析，該數(shù)據(jù)庫包括259 例肉瘤患者的數(shù)據(jù)，但是未對肉瘤進行具體分型，結果顯示：PTPRD （HR=1.79，P=0.005）、ADAM10 （HR=1.94，P=0.010）、ELAVL2 （HR=1.56，P=0.031） 和ERRFI1 （HR=2.05，P=0.005） 表達水平升高與RMS 患者的不良預后有關，而上述4 個基因在融合基因陽性的肉瘤患者腫瘤組織中表達上調（圖6）。其余各基因的表達差異對肉瘤患者的預后無明顯影響，差異無統(tǒng)計學意義（P＞0.05）。

圖6 與肉瘤患者預后相關的hub 基因Fig.6 Hub genes associated with prognosis of sarcoma patients

3 討 論

本研究采用生物信息學方法對融合基因陽性和陰性RMS 患者進行分析，共篩選出115 個目標基因。KEGG 富集分析顯示：目標基因主要富集在細胞外基質-受體相互作用、癌癥中蛋白聚糖等途徑中。細胞外基質是調節(jié)細胞發(fā)育、正常功能和體內(nèi)平衡的最重要的因素之一。研究［6］表明：細胞外基質和細胞外基質結合受體的異常組成和功能對癌癥發(fā)生發(fā)展具有重要意義。miR-215-5p 通過調節(jié)細胞外基質-受體相互作用和局灶性黏附等分子途徑，抑制結直腸癌細胞的轉移［7］。在腫瘤組織中，腫瘤細胞分泌的可溶性因子可被蛋白聚糖激活。由于其作為共受體的功能，蛋白聚糖可促進腫瘤細胞的增殖［6］。蛋白聚糖不僅在腫瘤基質中發(fā)揮作用，還可以對癌細胞的表型、特性和耐藥性的發(fā)展進行調節(jié)［8］。上述研究表明：細胞外基質-受體相互作用和癌癥中蛋白聚糖等途徑可能是融合基因在RMS中潛在的作用機制。

本研究利用cytoHubba 插件對115 個目標基因進行篩選，選出Top10 hub 基因，根據(jù)hub 基因查詢預測靶基因的miRNA，將hub 基因聯(lián)合miRNA得到miRNA-mRNA 分子調控網(wǎng)絡，結合生存分析篩 選 出miRNA-214-3p （PTPRD、 ERRFI1）、miR-375 （ELAVL2） 和miR-199a-3p （ADAM10、ELAVL2），可為研究RMS 的融合基因提供參考。

miRNA 是一類保守的非編碼核糖核酸，廣泛參與腫瘤的增殖、 侵襲和血管生成等過程。在骨肉瘤組織中miR-214-3p 的高表達能通過靶向泛醌-細胞色素c 還原酶復合物核心蛋白1 （ubiquinol cytochrome c reductase core protein 1，UQCRC1）促進骨肉瘤細胞增殖，為預防和治療骨肉瘤提供了潛在的治療靶點［9］。在食管鱗狀細胞癌中miR-214-3p為抑癌因子，其表達下調通過靶向生存素（survivin） 和CUG 結 合 蛋 白1 （CUG-binding protein 1，CUG-BP1） 促進食管癌細胞的化療耐藥性［10］。有研究［11］顯示：miR-214 在人RMS 細胞系中明顯下調，且miR-214 和神經(jīng)母細胞瘤RAS病毒致癌基因（neuroblastoma RAS viral oncogene，

N-ras） 調控還能抑制RMS 細胞生長和異種移植瘤的發(fā)生。miR-214-3p 在融合基因陽性的RMS 中下調，且PTPRD 和ERRFI1 為miR-214-3p 的靶基因。PTPRD 是蛋白質酪氨酸磷酸酶家族的成員之一，在膠質母細胞瘤和肺癌中已有通過缺失、突變或表觀遺傳甲基化使PTPRD 失活的報道［12］。PTPRD 在神經(jīng)母細胞瘤中通過Aurora 激酶A（Aurora kinase A，AURKA） 去磷酸化和去穩(wěn)定化以及下游編碼轉錄因子N-MYC （encoding N-MYC，MYCN） 蛋白的去穩(wěn)定化發(fā)揮抑瘤作用［13］。外源性表達PTPRD 在人膠質母細胞瘤、黑色素瘤和結腸癌中可以抑制腫瘤細胞生長和遷移，還可以通過誘導細胞凋亡降低細胞活性［14-15］。膠質瘤衍生的miR-148a-3p 通過抑制ERRFI1 激活表皮生長因子受體（epidermal growth factor receptor，EGFR） /絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK） 信號通路促進腫瘤血管生成［16］。長鏈非編碼RNA （long non- coding RNA，LncRNA） 中細胞周期激酶抑制因子4 基因座中反義非編碼RNA （antisense non-coding RNA in INK4 locus，ANRIL） 通過表觀遺傳抑制ERRFI1的表達促進膽管癌的惡性進展［17］。目前在RMS 中尚未發(fā)現(xiàn)miR-214-3p 與靶基因PTPRD 和ERRFI1關系的相關研究報道。

miR-375 在多種癌組織中表達下調，通過作用于重要的癌基因來調節(jié)癌癥相關的過程，如細胞增殖、凋亡、侵襲、遷移和自噬［18-22］。在骨肉瘤患者中miR-375 表達下調，其可以成為骨肉瘤診斷、預后和化療敏感性預測的標志物［23］。在膠質瘤中miR-375 表達下調且通過抑制結締組織生長因子-表皮生長因子受體（connective tissue growth factorepidermal growth factor receptor，CTGF-EGFR）信號通路影響神經(jīng)膠質瘤細胞的增殖、侵襲和遷移［24］。然而在某些類型癌組織中miR-375 表達上調，如前列腺癌，miR-375 的高表達影響前列腺癌的分期和淋巴結轉移［25］。本研究結果顯示：在融合基因陽性的RMS 中，miRNA-375 表達上調。ELAVL2 作為miRNA-375 的靶基因，在食管鱗狀細胞癌組織中過表達，且通過抑制細胞凋亡提高對紫杉醇的耐藥性［26］。在RMS 中上調的miRNA-375可結合ELAVL2 做進一步分析。

MiR-199a-3p 通過調節(jié)乙醇胺激酶1（ethanolamine kinase 1，ETNK1） 可以促進胃癌細胞的侵襲和遷移，影響胃癌患者的預后［27］。LncRNA 中DNA 損傷激活的非編碼RNA （noncoding RNA-activated by DNA damage，NORAD）通過與miR-199a-3p 內(nèi)源性競爭調節(jié)人骨肉瘤癌細胞的增殖和遷移［28］。本研究結果顯示：ADAM10為miR-199a-3p 的靶基因，但目前尚未見相關研究報道。研究［29］表明： 在三陰性乳腺癌患者中miR-365 通過調節(jié)ADAM10 抑制腫瘤細胞的侵襲。促紅細胞生成素產(chǎn)生肝細胞受體A8（erythropoietin-producing hepatocyte receptor A8，EphA8） 是一種致癌基因，在胃癌組織中通過調節(jié)ADAM10 導致胃癌患者預后不良［30］。ADAM10 可能在骨肉瘤的進展和成骨細胞骨肉瘤的病理發(fā)展過程中起重要作用［31］。

本研究通過生物信息學的方法對融合基因作用于RMS 的miRNA-mRNA 調控網(wǎng)絡進行分析，可為研究融合基因在RMS 發(fā)生發(fā)展中的作用提供新思路，為后續(xù)的實驗研究提供參考。