蔡江瑜 盛旦丹 蔣佳 陳世益

復(fù)旦大學(xué)附屬華山醫(yī)院運(yùn)動(dòng)醫(yī)學(xué)科（上海 200040）

前交叉韌帶離斷聯(lián)合半月板切除誘發(fā)的骨關(guān)節(jié)炎早期模型相關(guān)基因芯片數(shù)據(jù)的生物信息學(xué)分析

蔡江瑜 盛旦丹 蔣佳 陳世益

復(fù)旦大學(xué)附屬華山醫(yī)院運(yùn)動(dòng)醫(yī)學(xué)科（上海 200040）

目的：探討前交叉韌帶（anterior cruciate ligament，ACL）離斷聯(lián)合半月板切除誘發(fā)的骨關(guān)節(jié)炎（osteoarthritis，OA）的基因表達(dá)和生物過程的改變，為OA分子機(jī)制的進(jìn)一步研究提供生物信息學(xué)依據(jù)。方法：從GEO數(shù)據(jù)庫下載Appleton等構(gòu)建的ACL離斷聯(lián)合半月板切除誘發(fā)的OA早期大鼠模型的基因芯片數(shù)據(jù)集，應(yīng)用生物信息學(xué)方法篩選差異表達(dá)基因（differentially expressed genes，DEGs），并進(jìn)行基因本體論（Gene ontology，GO）和日本京都基因和基因組百科全書代謝通路數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析。結(jié)果：共篩選出170個(gè)DEGs，其中包括97個(gè)上調(diào)基因和73個(gè)下調(diào)基因。上調(diào)基因主要富集在細(xì)胞外基質(zhì)，并與細(xì)胞外基質(zhì)交互通路等信號(hào)通路關(guān)系密切。而下調(diào)基因主要富集在肌肉收縮的生物學(xué)功能中，并與PPAR通路等信號(hào)通路關(guān)系密切。結(jié)論：細(xì)胞外基質(zhì)和肌肉收縮的改變對(duì)OA的發(fā)生發(fā)展起著關(guān)鍵作用。細(xì)胞外基質(zhì)受體交互通路和PPAR信號(hào)通路與OA密切相關(guān)，值得進(jìn)一步深入研究。

骨關(guān)節(jié)炎；差異表達(dá)基因；生物信息學(xué)

前交叉韌帶（anterior cruciate ligament，ACL）損傷是臨床常見的運(yùn)動(dòng)損傷之一。ACL損傷常會(huì)引起關(guān)節(jié)松弛、半月板損傷、股四頭肌肌力下降以及膝關(guān)節(jié)力學(xué)的變化，造成關(guān)節(jié)功能損害，最終導(dǎo)致骨關(guān)節(jié)炎（osteoarthritis，OA）的發(fā)生[1]。Lohmander等[2]指出，至少50%的ACL損傷患者最終會(huì)發(fā)生膝關(guān)節(jié)OA。另外，一篇系統(tǒng)綜述中強(qiáng)調(diào)，ACL損傷合并半月板損傷時(shí)OA的發(fā)生概率將進(jìn)一步增加[3]。OA的病理改變包括軟骨的進(jìn)行性減少和破壞，軟骨下骨增厚，骨贅形成，不同程度的滑膜炎癥反應(yīng)，韌帶和半月板組織退化，以及關(guān)節(jié)囊肥大，嚴(yán)重影響患者的日常生活和工作[4]。然而，目前關(guān)于ACL損傷引起OA發(fā)病和進(jìn)展的分子機(jī)制尚不清楚。生物信息學(xué)是近年來生命科學(xué)領(lǐng)域的新興學(xué)科，可為疾病的分子機(jī)制研究提供新的思路和研究前景。本文通過生物信息學(xué)相關(guān)方法對(duì)Appleton等構(gòu)建的ACL離斷聯(lián)合半月板切除誘發(fā)的骨關(guān)節(jié)炎早期大鼠模型的基因芯片數(shù)據(jù)集進(jìn)行重新分析，以探討ACL離斷聯(lián)合半月板切除誘發(fā)OA的基因表達(dá)和生物代謝過程的改變，為OA分子機(jī)制的進(jìn)一步研究提供生物信息學(xué)依據(jù)。

1 材料與方法

1.1 芯片數(shù)據(jù)

基因芯片數(shù)據(jù)集GSE8077來源于NCBI的GEO（Gene Expression Omnibus）數(shù)據(jù)庫，基于GPL1355 平臺(tái)（[Rat230_2]Affymetrix Rat Genome 230 2.0 Array）。根據(jù)其研究方案可知，實(shí)驗(yàn)對(duì)象為5只大鼠，通過右膝（實(shí)驗(yàn)側(cè)）前交叉韌帶離斷和半月板部分切除手術(shù)來構(gòu)建OA的早期模型，左膝（健側(cè)）不做處理作為自身對(duì)照。本研究通過生物信息學(xué)方法分析其術(shù)后4周獲取的兩側(cè)軟骨標(biāo)本的基因芯片數(shù)據(jù)。

1.2 數(shù)據(jù)預(yù)處理

利用R軟件Affymetrix工具包的RMA算法（robust multiarray average algorithm）將原始的CEL文件進(jìn)行質(zhì)量控制和標(biāo)準(zhǔn)化處理，并轉(zhuǎn)化為探針表達(dá)矩陣，然后根據(jù)GPL1355平臺(tái)文件將探針名轉(zhuǎn)化為基因名。

1.3 差異基因的篩選

通過R軟件limma工具包（Linear Models for Microarray data）篩選出兩組之間的差異表達(dá)基因（differentially expressed genes，DEGs）。 log2fold change（log2FC）用于評(píng)價(jià)基因的表達(dá)水平，若P＜0.05且|log2FC|＞1可認(rèn)為具有統(tǒng)計(jì)學(xué)差異。利用pheatmap工具包繪制熱圖，直觀地展示每個(gè)差異基因在每個(gè)樣本中的表達(dá)情況。

1.4 基因本體論和通路富集分析

基因本體論（Gene ontology，GO）是用來注釋基因及其產(chǎn)物的常用方法，有利于集中研究感興趣的方向，發(fā)現(xiàn)新的現(xiàn)象。日本京都基因和基因組百科全書代謝通路數(shù)據(jù)庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）是用來分析生物功能和代謝通路的數(shù)據(jù)庫。我們利用DAVID（Database for Annotation，Visualization and Integrated Discovery）在線工具（http：//david.abcc.ncifcrf.gov/）對(duì)DEGs進(jìn)行GO和KEGG通路富集分析[5]。若P＜0.05則認(rèn)為差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 數(shù)據(jù)處理和DEGs篩選

基因芯片數(shù)據(jù)經(jīng)標(biāo)準(zhǔn)化后以箱圖形式呈現(xiàn)（圖1）。圖1b黑線基本在同一水平，表明一致性較高。經(jīng)數(shù)據(jù)預(yù)處理后，R軟件共篩選出170個(gè)DEGs，其中包括97個(gè)上調(diào)基因和73個(gè)下調(diào)基因（圖2）。DEGs表達(dá)的熱圖如圖3所示，可見兩組樣本基因表達(dá)具有顯著差異。

圖1 標(biāo)準(zhǔn)化之前（a）和之后（b）的基因表達(dá)數(shù)據(jù)

圖2 樣本基因芯片的火山圖分析

2.2 GO和KEGG通路富集分析

GO可分為生物過程（biological process，BP），細(xì)胞組成（cellular component,CC）和分子功能（molecular function,MF）。分別選取三類中排列前三的分析結(jié)果如表1所示，上調(diào)基因主要包括COL4A1、TNC、TNN、LAMB1、MYLPF、ACTN3等，下調(diào)基因有FABP4、ADIPOQ、FABP5、PCK1等。KEGG通路富集分析結(jié)果如表2所示。結(jié)合兩者結(jié)果可以發(fā)現(xiàn)，上調(diào)基因主要富集在細(xì)胞外基質(zhì)，并與細(xì)胞外基質(zhì)交互通路等信號(hào)通路關(guān)系密切。而下調(diào)基因主要富集在肌肉收縮的生物學(xué)功能中，并與PPAR通路等信號(hào)通路關(guān)系密切。

圖3 DEGs熱圖分析

表1 DEGs的GO分析

表2 DEGs的KEGG通路富集分析

3 討論

近年來，關(guān)于OA樣本的基因組學(xué)、轉(zhuǎn)錄組學(xué)和蛋白組學(xué)研究逐漸成為分子機(jī)制研究的熱點(diǎn)[6-8]。本研究的基因芯片數(shù)據(jù)集來源于GEO數(shù)據(jù)庫，我們利用生物信息學(xué)方法比較實(shí)驗(yàn)側(cè)和健側(cè)樣本的差異基因表達(dá)，共篩選出170個(gè)DEGs，其中包括97個(gè)上調(diào)基因和73個(gè)下調(diào)基因。

上調(diào)基因 COL4A1、TNC、TNN、LAMB1、MYLPF、ACTN3主要富集在細(xì)胞外基質(zhì)。功能富集分析顯示，OA相關(guān)基因與細(xì)胞外基質(zhì)受體交互通路有關(guān)。之前的研究曾指出，此信號(hào)通路在OA的進(jìn)展中起著重要作用。Koelling等[9]曾報(bào)道在取自O(shè)A晚期軟骨的祖細(xì)胞中，與細(xì)胞外基質(zhì)受體交互通路相關(guān)的基因表達(dá)顯著升高。Fei等[10]通過生物信息學(xué)分析了OA患者的關(guān)節(jié)液樣本，發(fā)現(xiàn)在KEGG分析中，細(xì)胞外基質(zhì)受體交互通路是最為相關(guān)的一條通路。同時(shí)，作者也發(fā)現(xiàn)COL4A1富集于此通路中。軟骨細(xì)胞合成并分泌基質(zhì)成分,構(gòu)成了細(xì)胞生活的微環(huán)境。因此，細(xì)胞外基質(zhì)成分和結(jié)構(gòu)的改變可引起軟骨代謝平衡穩(wěn)態(tài)的破壞，從而導(dǎo)致OA的發(fā)生和進(jìn)展。COL4A1在OA的基因芯片研究中有見報(bào)道[11-13]，但是，對(duì)于其具體導(dǎo)致OA發(fā)生的分子機(jī)制有待進(jìn)一步研究。

下調(diào)基因FABP4、ADIPOQ、FABP5、PCK1主要富集在肌肉收縮的生物學(xué)功能中。普遍認(rèn)為，在膝關(guān)節(jié)中，肌肉組織起到重要的運(yùn)動(dòng)和感覺功能，如活動(dòng)關(guān)節(jié)、維持關(guān)節(jié)穩(wěn)定性、減震緩沖，以及作為本體感受器等[14]。其中股四頭肌與OA的發(fā)生和進(jìn)展密切相關(guān)。有學(xué)者提出ACL損傷造成股四頭肌肌力下降的原因是存在關(guān)節(jié)源性的肌肉抑制（arthrogenic muscle inhibi-tion，AMI），即ACL損傷后個(gè)體無法徹底自主活動(dòng)股四頭肌，因此容易使股四頭肌肌力下降及萎縮[15]。股四頭肌肌力下降會(huì)造成主動(dòng)肌和拮抗肌的協(xié)調(diào)作用紊亂，使膝關(guān)節(jié)負(fù)荷加重，可作為OA發(fā)生的獨(dú)立危險(xiǎn)因素[16]。而OA所引起的疼痛和膝關(guān)節(jié)功能紊亂又會(huì)加重股四頭肌的萎縮[17]。因此，肌肉收縮的改變與OA發(fā)生和發(fā)展互為因果。此外，KEGG通路富集分析顯示，PPAR信號(hào)通路在OA中起著一定作用。最新的研究發(fā)現(xiàn)，過氧化物酶體增生物激活受體（peroxisome proliferators-activated receptor，PPAR）可能與OA的發(fā)病機(jī)制密切相關(guān)[18，19]。PPAR 包括 PPARα、PPAR β/δ、PPARγ 3個(gè)亞型，3個(gè)亞型在OA的發(fā)病中均有一定作用[20-22]。特別是PPARγ可以抑制IL-1β誘導(dǎo)的蛋白多糖的降解過程[23]。Vasheghani等[24]發(fā)現(xiàn)軟骨特異性（cartilage-specific）PPARγ敲除小鼠能自發(fā)形成OA。因此，PPAR有望成為潛在的藥物治療靶點(diǎn)。

4 總結(jié)

綜上所述，細(xì)胞外基質(zhì)和肌肉收縮的改變對(duì)OA的發(fā)生發(fā)展起著關(guān)鍵作用。細(xì)胞外基質(zhì)受體交互通路和PPAR信號(hào)通路與OA密切相關(guān)，值得進(jìn)一步深入研究。

[1]Allen CR，Livesay GA，Wong EK，et al.Injury and reconstruction of the anterior cruciate ligament and knee osteoarthritis.Osteoarthritis Cartilage，1999，7（1）：110-121.

[2]Lohmander LS，Englund PM，Dahl LL，et al.The longterm consequence of anterior cruciate ligament and meniscus injuries：osteoarthritis.Am J Sports Med，2007，35（10）：1756-1769.

[3]Oiestad BE，Engebretsen L，Storheim K，et al.Knee osteoarthritis after anterior cruciate ligament injury：a systematic review.Am J Sports Med，2009，37（7）：1434-1443.

[4]Loeser RF，Goldring SR，Scanzello CR，et al.Osteoarthritis：a disease of the joint as an organ.Arthritis Rheum，2012，64（6）：1697-1707.

[5]Sherman BT，Huang Da W，Tan Q，et al.DAVID Knowledgebase：a gene-centered database integrating heterogeneous gene annotation resources to facilitate high-throughput gene functional analysis.BMC Bioinformatics，2007，8：426.

[6]Steinberg J and Zeggini E.Functional genomics in osteoarthritis：Past，present，and future.J Orthop Res，2016，34（7）：1105-1110.

[7]Wang Q，Rozelle AL，Lepus CM，et al.Identification of a central role for complement in osteoarthritis.Nat Med，2011，17（12）：1674-1679.

[8]Ikeda D，Ageta H，Tsuchida K，et al.iTRAQ-based proteomics reveals novel biomarkers of osteoarthritis.Biomarkers，2013，18（7）：565-572.

[9]Koelling S，Kruegel J，Irmer M，et al.Migratory chondrogenic progenitor cells from repair tissue during the later stages of human osteoarthritis.Cell Stem Cell，2009，4（4）：324-335.

[10]Fei Q，Lin J，Meng H，et al.Identification of upstream regulators for synovial expression signature genes in osteoarthritis.Joint Bone Spine，2016，83（5）：545-551.

[11]Olex AL，Turkett WH，F(xiàn)etrow JS，et al.Integration of gene expression data with network-based analysis to identifysignalingand metabolicpathwaysregulated during the development of osteoarthritis.Gene，2014，542（1）：38-45.

[12]Lamas JR，Rodriguez-Rodriguez L，Vigo AG，et al.Largescale gene expression in bone marrow mesenchymal stem cells：a putative role for COL10A1 in osteoarthritis.Ann Rheum Dis，2010，69（10）：1880-1885.

[13]He P，Zhang Z，Liao W，et al.Screening of gene signatures for rheumatoid arthritis and osteoarthritis based on bioinformatics analysis.Mol Med Rep，2016，14（2）：1587-1593.

[14]Hurley MV.The role of muscle weakness in the pathogenesisofosteoarthritis.Rheum DisClin North Am，1999，25（2）：283-298，vi.

[15]Palmieri-Smith RM and Thomas AC.A neuromuscular mechanism of posttraumatic osteoarthritis associated with ACL injury.Exerc Sport Sci Rev，2009，37（3）：147-153.

[16]Roos EM，Herzog W，Block JA，et al.Muscle weakness，afferent sensory dysfunction and exercise in knee osteoarthritis.Nat Rev Rheumatol，2011，7（1）：57-63.

[17]Bennell KL，Wrigley TV，Hunt MA，et al.Update on the role of muscle in the genesis and management of knee osteoarthritis.Rheum Dis Clin North Am，2013，39（1）：145-176.

[18]Qu Y，Zhou L，and Wang C.Mangiferin Inhibits IL-1beta-Induced Inflammatory Response by Activating PPAR-gamma in Human Osteoarthritis Chondrocytes.Inflammation，2017，40（1）：52-57.

[19]Sacitharan PK and Vincent TL.Cellular ageing mechanisms in osteoarthritis.Mamm Genome，2016，27（7-8）：421-429.

[20]Ratneswaran A，Leblanc EA，Walser E，et al.Peroxisome proliferator-activated receptor delta promotes the progression of posttraumatic osteoarthritis in a mouse model.Arthritis Rheumatol，2015，67（2）：454-464.

[21]Vasheghani F，Zhang Y，Li YH，et al.PPARgamma deficiency results in severe，accelerated osteoarthritis associated with aberrant mTOR signalling in the articular cartilage.Ann Rheum Dis，2015，74（3）：569-578.

[22]Clockaerts S，Bastiaansen-Jenniskens YM，F(xiàn)eijt C，et al.Peroxisome proliferator activated receptor alpha activation decreases inflammatory and destructive responses in osteoarthriticcartilage.OsteoarthritisCartilage，2011，19（7）：895-902.

[23]Francois M，Richette P，Tsagris L，et al.Activation of the peroxisome proliferator-activated receptor alpha pathway potentiatesinterleukin-1 receptorantagonistproduction in cytokine-treated chondrocytes.Arthritis Rheum，2006，54（4）：1233-1245.

[24]Vasheghani F，Monemdjou R，F(xiàn)ahmi H，et al.Adult cartilage-specific peroxisome proliferator-activated receptor gamma knockout mice exhibit the spontaneous osteoarthritis phenotype.Am J Pathol，2013，182（4）：1099-1106.

Bioinformatics Analysis of Gene Arrays in an Early Osteoarthritis Model Induced by Anterior Cruciate Ligament Transection and Partial Medial Meniscectomy

Cai Jiangyu，Sheng Dandan，Jiang Jia，Chen Shiyi
Department of Sports Medicine，Huashan Hospital，F(xiàn)udan University，Shanghai 200040，China

Chen Shiyi，Email：cshiyi@163.com

ObjectiveTo explore the changes in gene expression and biological process of the osteoarthritis（OA）induced by anterior cruciate ligament（ACL）transection and partial medial meniscectomy，so as to provide bioinformatic basis for further studying the molecular mechanism of OA.MethodsThe gene chip datasets of a rat model of early 0A induced by ACL transection and partial medial meniscectomy were downloaded from GEO databases（submitted by Appleton,et al.）.The differential expression genes（DEGs）were identified，and the Gene ontology（GO） as well as the Kyoto Encyclopedia of Genes and Genomes（KEGG）pathway enrichment analyses for DEGs were conducted using bioinformatic methods.ResultsA total of 170 DEGs including 97 up-regulated genes and 73 down-regulated genes were identified.The up-regulated genes were mainly enriched in the extracellular matrix（ECM）and were closely related to the ECM-receptor interaction，while the down-regulated genes were mainly enriched in the biological function of muscle contraction and were linked with the peroxisome proliferators-activated receptor（PPAR） signaling pathway.ConclusionThe changes of ECM and muscle contraction play a key role in the occurrence and development of OA.The ECM-receptor interaction and PPAR signaling pathway are strongly associated with OA and worthy of further study.

osteoarthritis，differentially expressed genes，bioinformatics

2017.03.21

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2015AA033703)；國家自然科學(xué)基金（81572108，81370052）；國家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2016YFC1100300）

第1作者：蔡江瑜，caijiangyu1@126.com；

陳世益，Email：cshiyi@163.com