多芯片聯(lián)合分析揭示肌腱粘連的分子機(jī)制

佟春曉，雷則鳴

（1.中國(guó)醫(yī)科大學(xué)附屬第一醫(yī)院產(chǎn)科，沈陽 110001；2.沈陽醫(yī)學(xué)院附屬中心醫(yī)院手外五科，沈陽 110024；3.中國(guó)醫(yī)科大學(xué)附屬盛京醫(yī)院骨外科，沈陽 110004）

肌腱粘連是常見且棘手的并發(fā)癥，可導(dǎo)致肢體活動(dòng)嚴(yán)重受限，甚至功能喪失。目前肌腱粘連的發(fā)生機(jī)制仍未完全闡明，主要學(xué)說認(rèn)為是由肌腱內(nèi)源性和外源性愈合不平衡引起［1］。內(nèi)源性愈合主要是腱細(xì)胞經(jīng)增殖及遷移致肌腱的愈合。外源性愈合是肌腱損傷斷面形成肉芽組織，腱外膜的成纖維細(xì)胞增殖侵襲，由周圍組織向肌腱斷端生長(zhǎng)，腱周纖維化，膠原蛋白和纖維蛋白沉積，形成肌腱與周圍組織的瘢痕組織。外源性愈合占據(jù)主導(dǎo)時(shí)會(huì)導(dǎo)致肌腱粘連，減弱外源性愈合及炎癥參與，促進(jìn)內(nèi)源性愈合，是當(dāng)前解決肌腱粘連問題和提高肌腱術(shù)后或損傷后腱愈合預(yù)后的方法。鑒于目前對(duì)天然生態(tài)位和動(dòng)態(tài)肌腱愈合過程的了解有限，需要進(jìn)一步研究肌腱損傷、愈合、移除過程中的肌腱生物學(xué)變化。本研究利用生物信息學(xué)技術(shù)對(duì)基因表達(dá)數(shù)據(jù)庫（Gene Expression Omnibus，GEO）中含有肌腱病和轉(zhuǎn)化生長(zhǎng)因子β（transforming growth factor-β，TGF-β）處理的成纖維細(xì)胞的表達(dá)譜芯片進(jìn)行聯(lián)合分析，挖掘肌腱粘連中的關(guān)鍵基因、重要通路、調(diào)控網(wǎng)絡(luò)及生物代謝過程，為探索肌腱粘連的分子機(jī)制提供依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)來源

基因芯片數(shù)據(jù)集GSE26051和GSE1724均來源于GEO數(shù)據(jù)庫。GSE26051芯片基于GPL570平臺(tái)（Affymetrix Human Genome U133 Plus 2.0 Array），GSE1724芯片基于GPL8300平臺(tái)（Affymetrix Human Genome U95 Version 2 Array）。GSE26051研究對(duì)象為23例慢性肌腱病患者肌腱病變組織和肌腱非病變組織。在GSE1724數(shù)據(jù)集中選取研究對(duì)象為特發(fā)性肺纖維化患者的肺成纖維細(xì)胞非處理組和TGF-β處理組，各3例，比較基因表達(dá)譜。聯(lián)合分析2個(gè)數(shù)據(jù)集的基因芯片數(shù)據(jù)。

1.2 差異表達(dá)基因（different expression genes，DEGs）的篩選

采用R語言“l(fā)imma”包篩選組間DEGs。對(duì)樣本質(zhì)檢表達(dá)量的差異倍數(shù)（fold change，F(xiàn)C）取對(duì)數(shù)，用于評(píng)價(jià)基因的表達(dá)水平，若P＜0.05且|logFC|＞1，認(rèn)為有統(tǒng)計(jì)學(xué)意義。采用Pheatmap軟件包繪制熱圖，分析DEGs在樣本中的表達(dá)。

1.3 韋恩圖繪制

利用韋恩圖繪制工具（https：//bioinformatics.psb.ugent.be/webtools/Venn/）制作韋恩圖，將GSE26051和GSE1724數(shù)據(jù)集篩選出的DEGs取交集。

1.4 DEGs的京都基因和基因組數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集和基因本體（gene ontology，GO）功能富集分析

為了探索DEGs的主要功能和途徑，通過R語言“Bioconductor”包對(duì)DEGs進(jìn)行GO富集和KEGG通路富集，選出P＜0.05的結(jié)果；最后通過“ggplot2”包對(duì)富集結(jié)果分別進(jìn)行可視化。

1.5 預(yù)測(cè)DEGs可能結(jié)合的微RNA（microRNA，mi-RNA）

利用R語言結(jié)合miRDB數(shù)據(jù)庫（https：//mirdb.org）、ENCORI數(shù)據(jù)庫（https：//starbase.sysu.edu.cn）、miRTarBase數(shù)據(jù)庫（https：//mirtarbase.mbc.nctu.edu.tw）、miRWalk數(shù)據(jù)庫（https：//mirwalk.umm.uni-heidelberg.de）、RNA22數(shù)據(jù)庫（https：//cm.jefferson.edu/rna22）、RNAInter數(shù)據(jù)庫（https：//www.rnainter.org）、Target-Miner數(shù)據(jù)庫（https：//www.isical.ac.in/～bioinfo_miu/targetminer20.htm）、TargetScan數(shù)據(jù)庫（https：//www.targetscan.org），聯(lián)合預(yù)測(cè)DEGs可能結(jié)合的miRNA，結(jié)果用R語言“UpSetR”展示。

1.6 蛋白質(zhì)-蛋白質(zhì)相互作用（protein-protein interaction，PPI）網(wǎng)絡(luò)分析

將DEGs導(dǎo)入STRING數(shù)據(jù)庫（https：//string-db.org/）進(jìn)行PPI分析，并設(shè)定篩選閾值為0.4，獲得共有差異蛋白之間的互作圖。

2 結(jié)果

2.1 DEGs的篩選

對(duì)GSE26051芯片數(shù)據(jù)集（肌腱病患者病灶與非病灶，各23 例）及GSE1724芯片數(shù)據(jù)集（特發(fā)性肺纖維化患者的肺成纖維細(xì)胞TGF-β處理組與成纖維細(xì)胞對(duì)照組各3例）進(jìn)行分析，根據(jù)P＜0.05且|logFC|＞1篩選組間DEGs。按照|logFC|從大到小排列，GSE26051和GSE1724篩選出的前十位DEGs分別為TRIM54、BE156417、LRP1B、JSRP1、BF514585、ACADL、BF680284、ESRRG、DOK7、BG541187以及ESM1、IL11、NGF、HBEGF、SOCS1、PCYT2、MEOX1、IER3、TXNIP、NEDD9。

2.2 火山圖和聚類圖繪制分析差異表達(dá)譜

根據(jù)芯片篩選出的DEGs繪制火山圖和熱圖。火山圖（圖1）反映總體基因的表達(dá)情況，紅色和綠色區(qū)域代表P＜0.05且|logFC|＞1的基因。根據(jù)基因表達(dá)水平，將檢測(cè)的基因和標(biāo)本進(jìn)行聚類分析，采用Pheatmap軟件包繪制聚類分析熱圖，紅色代表基因表達(dá)量相對(duì)高，綠色代表基因表達(dá)量相對(duì)低。由聚類分析熱圖發(fā)現(xiàn)2組樣本中基因的表達(dá)存在差異。

圖1 GSE1724和GSE26051芯片中DEGs的火山圖Fig.1 Volcano map of differentially expressed genes in the GSE1724 and GSE26051 gene microarray

2.3 韋恩圖

利用韋恩圖繪制工具繪制韋恩圖（圖2），將GSE26051和GSE1724數(shù)據(jù)集分別篩選出的45個(gè)和526個(gè)DEGs取交集，共篩選出2個(gè)基因（EPB41、ETV6）。

圖2 GSE1724與GSE26051數(shù)據(jù)集DEGs的韋恩圖Fig.2 Mapping of the differential expression of gene between GSE1724 and GSE26051 data sets

2.4 DEGs的miRNA預(yù)測(cè)

利用miRDB數(shù)據(jù)庫、ENCORI數(shù)據(jù)庫、miRTarBase數(shù)據(jù)庫、miRWalk數(shù)據(jù)庫、RNA22數(shù)據(jù)庫、RNAInter數(shù)據(jù)庫、TargetMiner數(shù)據(jù)庫、TargetScan數(shù)據(jù)庫預(yù)測(cè)篩選出的差異基因（EPB41、ETV6）可能結(jié)合的miRNA，結(jié)果用R語言“UpSetR”包可視化。結(jié)果顯示，EPB41、ETV6分別最可能與hsa-miR-30e-5p、hsamiR-129-5p結(jié)合，EPB41、ETV6聯(lián)合在以上數(shù)據(jù)庫篩選出最可能與hsa-miR-181b-5p、hsa-miR-181c-5p、hsa-miR-181d-5p結(jié)合。

2.5 GO功能富集和KEGG通路富集分析

由于GSE26051篩選出的差異基因數(shù)目較少，不對(duì)其富集。通過R軟件對(duì)GSE1724篩選出的526個(gè)DEGs進(jìn)行富集分析，以了解其生物學(xué)功能，結(jié)果顯示，有13個(gè)GO項(xiàng)目和12個(gè)KEGG項(xiàng)目，對(duì)富集結(jié)果進(jìn)行可視化處理。其中與肌腱粘連密切相關(guān)的GO分析主要富集在轉(zhuǎn)錄激活與抑制、G蛋白耦聯(lián)受體、生長(zhǎng)因子活性、受體配體活性等（圖3A）；與肌腱粘連高度相關(guān)的 KEGG通路有PI3K/AKT、MAPK、鈣通道等通路（圖3B）。

圖3 GSE1724篩選出的差異基因功能富集Fig.3 Functional and signaling pathway analysis of GSE1724 DEGs

2.6 PPI網(wǎng)絡(luò)

將EPB41、ETV6導(dǎo)入STRING數(shù)據(jù)庫，在獲得的PPI網(wǎng)絡(luò)中，DLG1、CASK、GYPC、SF3A1、SF3A3、SPTAN1是已被驗(yàn)證過的與EPB41關(guān)系密切且相互作用的蛋白（圖4A）；BRCA1、SUMO1、MDC1、TP-53BP1、SIN3A、NCOR1、ETV7、SOCS1是已被驗(yàn)證過的與ETV6關(guān)系密切且相互作用的蛋白（圖4B）。

圖4 差異基因的PPI網(wǎng)絡(luò)Fig.4 The protein-protein interaction map of the differential gene

3 討論

如何改善及預(yù)防術(shù)后粘連是外科醫(yī)生關(guān)注的重要問題。急性或慢性肌腱損傷后，如治療效果不理想，愈合時(shí)間長(zhǎng)，最終會(huì)導(dǎo)致粘連。肌腱粘連被認(rèn)為是繼發(fā)于炎癥過程的纖維化組織反應(yīng)，持續(xù)的炎癥會(huì)加重纖維化［2］。目前研究較多的TGF-β被認(rèn)為與肌腱粘連的纖維化最相關(guān)［3］。肌腱損傷后，巨噬細(xì)胞被募集到愈合部位，釋放TGF-β，與成纖維細(xì)胞表面TGF-β受體結(jié)合，促進(jìn)成纖維細(xì)胞增殖及分化，膠原合成增加，最終導(dǎo)致粘連形成［4-5］。而阻斷TGF-β途徑被認(rèn)為是避免肌腱粘連形成且不干預(yù)肌腱愈合強(qiáng)度的有效方法［6-7］。除TGF-β外，其他通路也有被研究，如核轉(zhuǎn)錄因子-κB（nuclear factor-kappa B，NF-κB）通路及絲裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）通路也被認(rèn)為與肌腱粘連的形成相關(guān)。NF-κB途徑主要參與肌腱修復(fù)過程中炎癥的初始階段［8］。MAPK途徑則參與外源性成纖維細(xì)胞增殖、炎癥反應(yīng)及調(diào)節(jié)膠原合成，抑制該通路或蛋白磷酸化可阻止肌腱粘連的發(fā)生［9］。

通過納米顆粒、水凝膠、電紡纖維膜、羊膜等［1］物理屏障搭載藥物治療的方式是肌腱粘連主要的治療策略。從分子水平上尋找有效的藥物靶點(diǎn)目前仍是研究難題。本研究檢索人類肌腱粘連疾病的相關(guān)數(shù)據(jù)集無結(jié)果，考慮肌腱粘連為肌腱損傷未愈導(dǎo)致的疾病，故選用人類慢性肌腱損傷病灶組織芯片數(shù)據(jù)集GSE26051。而TGF-β與腱外膜中的成纖維細(xì)胞功能的關(guān)系與肌腱粘連的發(fā)生發(fā)展高度相關(guān)，且肌腱粘連［1］與特發(fā)性肺纖維化［10］的標(biāo)志都是豐富的膠原沉積和成纖維細(xì)胞增殖，故搜索到數(shù)據(jù)集GSE1724，試圖從TGF-β及成纖維細(xì)胞的角度探討慢性肌腱損傷導(dǎo)致的肌腱粘連疾病背后的分子機(jī)制。通過生物信息學(xué)分析，篩選出2個(gè)差異基因，即EPB41和ETV6。ETV6屬于ETS轉(zhuǎn)錄因子家族，作為轉(zhuǎn)錄抑制劑，在造血及癌癥的惡性轉(zhuǎn)化等方面發(fā)揮功能［11］。EPB41通過結(jié)構(gòu)域與各種跨膜蛋白結(jié)合，介導(dǎo)肌動(dòng)蛋白-膜蛋白的連接，在細(xì)胞黏附、細(xì)胞遷移、細(xì)胞膜架構(gòu)等生理過程中起重要作用［12］。研究［13］表明，在纖維連接蛋白上培養(yǎng)的EPB41-/-角質(zhì)形成細(xì)胞不能形成肌動(dòng)蛋白應(yīng)力纖維，不產(chǎn)生局灶性粘連。本研究從TGF-β處理的成纖維細(xì)胞及肌腱病灶數(shù)據(jù)集中篩選出的EPB41是否可能通過影響成纖維細(xì)胞的增殖、遷移、侵襲水平而改變肌腱粘連的結(jié)局，有待進(jìn)一步研究。

miRNA在肌腱粘連的分子機(jī)制中起關(guān)鍵作用，或可搭載載體治療肌腱粘連。如骨髓巨噬細(xì)胞可分泌直接靶向Smad7的外體miR-21-5p，從而激活肌腱細(xì)胞中的纖維生成，參與腱周纖維化［14］。有研究［15］通過將環(huán)氧合酶（COX-1和COX-2）工程化的miRNA質(zhì)粒嵌入透明質(zhì)酸水凝膠中，在屈肌肌腱愈合的早期階段，通過炎癥反應(yīng)下調(diào)肌腱組織中環(huán)氧合酶的表達(dá)，減少肌腱粘連。本研究通過預(yù)測(cè)DEGs可能結(jié)合的miRNA，擬構(gòu)建較完整的miRNA-mRNA網(wǎng)絡(luò)。為了預(yù)測(cè)DEGs可能結(jié)合的miRNA，分別將EPB41和ETV6的6個(gè)miRNA數(shù)據(jù)庫的預(yù)測(cè)結(jié)果取交集，結(jié)果聚焦于hsa-miR-181家族。研究［16］表明，miR-181b可介導(dǎo)表皮生長(zhǎng)因子受體依賴的血管細(xì)胞黏附分子-1表達(dá)，介導(dǎo)單核細(xì)胞通過整合素α4β1與多形性膠質(zhì)母細(xì)胞瘤之間的黏附。miR-181家族參與包括MYC、TGF-β、Wnt、STAT3、NOTCH、NF-κB、PI3K、MAPK等多條信號(hào)通路，可調(diào)節(jié)細(xì)胞增殖、凋亡、代謝、上皮-間質(zhì)轉(zhuǎn)化、遷移、侵襲、細(xì)胞周期、自噬等細(xì)胞功能［17］。miR-181b通過影響TGF-β通路在細(xì)胞外基質(zhì)重建功能方面發(fā)揮作用［18］。miR-181家族是否能靶作用于EPB41，影響TGF-β通路的表達(dá)，從而影響肌腱愈合階段的細(xì)胞外基質(zhì)重建，則有待后續(xù)深入研究。在功能富集方面，GO分析主要富集涉及轉(zhuǎn)錄激活與抑制、生長(zhǎng)因子活性等，KEGG通路富集涉及PI3K/AKT、MAPK、鈣通道、神經(jīng)激活配體受體等通路，這些是否與本研究的miRNA-mRNA網(wǎng)絡(luò)、PPI網(wǎng)絡(luò)相關(guān)，還需進(jìn)一步挖掘。

綜上所述，本研究通過分析基因芯片數(shù)據(jù)集，篩選肌腱粘連病變中的重要基因和通路，構(gòu)建潛在的miRNA-mRNA、PPI網(wǎng)絡(luò)，為闡明肌腱粘連的分子機(jī)制提供了科學(xué)依據(jù)。