基于基因芯片技術結合CIBERSORT反卷積算法研究骨關節炎的免疫機制及潛在中藥靶向預測①

楊雷 伍搏宇 熊輝齊 新宇 簡功輝 楊少鋒 聶穎 馬露

（湖南中醫藥大學，長沙 410208）

骨關節炎（osteoarthritis，OA）是一種以關節破壞、骨質增生和關節軟骨退行性病變為特征的慢性疾病。目前機制尚不清楚，多種因素均可導致OA［1-3］。既往認為OA是一種“磨損”疾病，而目前卻認為OA是一種免疫相關疾病，免疫復合物介導的補體激活、破骨細胞激活以及細胞因子平衡失調在OA病情進展中處于重要地位，備受關注［4-5］。早期OA患者影像學改變并不明顯，而滑膜炎卻是所有關節炎最早出現的共有特征［6］。滑膜組織作為OA免疫系統的損傷對象，是持續性關節退行性變的主要驅動因素，尤其免疫細胞浸潤，包括漿細胞、M2型巨噬細胞、活化樹突狀細胞等多種免疫細胞滑膜浸潤參與并介導了OA的自身免疫應答，誘導促炎細胞因子、趨化因子及蛋白酶分泌，擾亂免疫平衡進而加速軟骨和骨侵蝕，在OA發生發展、惡化過程中起重要作用［7-10］。因此，從免疫浸潤角度研究OA發病機制是針對OA早期治療、防止惡化的先決條件，可為OA精準靶向免疫治療提供新途徑。

生物信息學是分子生物學和信息技術的新興交叉學科，可揭示疾病分子機制。基因芯片作為一種新興技術，已被用于生物信息的高效、大規模獲取，可廣泛收集疾病表達譜數據，尤其對多種免疫細胞間的動態調節研究具有指導作用。本文采用生物信息學工具對通用基因芯片數據庫（gene expression omnibus，GEO）中OA患者和正常人群滑膜數據進行分析，通過CIBERSORT反卷積算法分析OA滑膜中的免疫浸潤情況，并根據差異的免疫浸潤基因預測相關免疫調節中藥，為OA免疫機制研究和免疫調節相關中藥篩選提供參考。

1 材料與方法

1.1 芯片數據采用GEO（www.ncbi.nlm.nih.gov/geo/）數據庫搜索OA滑膜基因芯片數據，篩選GSE55235數據，由Thomas H?upl提供（包含10個正常人群滑膜組織、10個OA患者滑膜組織），設正常滑膜為對照組，OA患者滑膜數據為實驗組，根據GPL96-57554（［HG-U133A］Affymetrix Human Genome U133A Array）平臺將探針轉換為基因名稱。

1.2 差異表達基因獲取通過R語言程序“limma”安裝包對基因探針進行兩樣本t檢驗，并根據篩選條件|log2FC|≥1和校正P<0.05進行差異基因分析，其中log2FC<0代表基因表達下調，log2FC>0代表基因表達上調。

1.3 核心基因獲取和注釋將兩組差異基因通過在線平臺STRING（www.string-db.org/cgi/input.pl）獲取蛋白互作（protein-protein interaction，PPI）網絡模型，并將結果導入Cytoscape軟件，采用Cytoscape插件工具cyto Hubba中的MCC方法識別重要基因，將重要基因進行可視化，并列出排名前10的作為核心基因。采用R語言cluster Profiler對差異基因進行基因本體論（genetic ontology，GO）分析與京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析，并探究其潛在生物學過程（biological process，BP）、細胞成分（cellular component，CC）、分子功能（molecular function，MF）及重要信號通路。

1.4 免疫細胞浸潤矩陣獲取采用R語言程序并鏈接CIBERSORT反卷積算法進行活化記憶性CD4 T細胞、漿細胞、活化NK細胞、調節性T細胞等22種免疫細胞轉錄特征的模擬計算［11］。為了結果準確，設定模擬次數為1 000次，采用Kruskal-Wallis秩和檢驗對P<0.05的數據進行后續分析。

1.5 免疫細胞相關性和差異分析為分析各類免疫細胞相關性，對CIBERSORT反卷積算法以P<0.05基礎上篩選的可信樣本數據中進行不同免疫細胞間Pearson相關系數計算，并采用秩和檢驗比較OA組與對照組差異。

1.6 OA免疫浸潤過程潛在中藥預測將兩組差異的核心靶基因和GO富集與免疫相關BP導入在線平臺（www.coremine.com/medical/），篩選具有潛在免疫調節作用的藥物，采用秩和檢驗，以P<0.01表示差異具有統計學意義，并對中藥進行功效、歸經統計。

1.7 主要觀察指標通過GEO數據庫獲取正常人群滑膜和OA患者滑膜組織差異基因及免疫細胞浸潤分布，并進行潛在中藥預測，主要觀察指標：①正常滑膜組織和OA滑膜組織免疫細胞浸潤的分布差異及相關性；②正常滑膜組織和OA滑膜組織差異的靶基因；③OA滑膜組織免疫浸潤相關BP和可能信號通路；④挖掘調節OA滑膜組織局部免疫微環境的潛在治療中藥。

2 結果

2.1 OA患者差異表達基因結果對來自GSE55235的微陣列進行篩選，共得到10個健康志愿者和10個OA患者滑膜組織芯片，正常滑膜組織和OA滑膜組織差異基因共716個，其中上調基因382個，下調基因334個（圖1）。

圖1 OA滑膜組織基因芯片差異表達基因熱圖和火山圖Fig.1 Heat map and volcano map of OA synovium-related differential expression genes

2.2 核心靶基因獲取結果將716個差異表達基因導入STRING繪制PPI，并將分析結果的TSV文件導入Cytoscape軟件進行可視化展示（圖2），經cyto Hubba計算列取前10位作為核心基因，分別為IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9、ITGB2、FOS、APOB、CXCL12（表1）。

表1 PPI模型中OA組和正常組差異核心靶基因Tab.1 Core differential target genes in PPI networks of OA and normal groups

圖2 OA滑膜組織差異基因PPI及核心蛋白基因模塊Fig.2 OA synovium tissue differential gene PPI networks map and Hub genes

2.3 差異基因GO和KEGG富集分析采用R軟件cluster Profiler對OA差異基因進行BP和KEGG分析，上調和下調的差異基因分別為382、334個，GO分析得到BP 374、879條，上調基因共富集到34條通路，下調基因共富集到39條通路，對比上調基因和下調基因GO和KEGG富集分析結果，發現上調基因與免疫密切相關，GO主要涉及白細胞正向調節、吞噬作用、T細胞呈遞激活等BP（圖3A）。KEGG富集的相關通路主要涉及吞噬小體、TNF信號通路、IL-17信號通路等相關信號通路（圖3B），而下調基因GO主要涉及對氧水平的反應、對氧水平下降的反應、節律過程等非免疫BP（圖3C），相關信號通路主要涉及毒性心肌炎相關通路、利什曼病相關信號通路等轉導機制（圖3D）。

圖3 OA差異基因GO及KEGG富集分析Fig.3 GO and KEGG enrichment analysis of differential genes related to OA

2.4 免疫浸潤細胞分布采用CIBERSORT反卷積法以P<0.05為篩選條件對芯片進行篩選，共得到20個可信樣本，熱圖左側10個為正常滑膜組織，右側10個為OA滑膜組，未活化的記憶CD4 T細胞和活化的肥大細胞含量在OA滑膜組織中不同程度減少，而M0型巨噬細胞、未活化的肥大細胞含量在OA滑膜中均出現不同程度增加（圖4A），圖4B進一步展示了22種免疫細胞在各樣本中的含量分布。

圖4 正常人群和OA患者免疫細胞浸潤分布Fig.4 Infiltration of immune-associated cells in normal and OA samples

2.5 免疫細胞相關性及差異分析未活化的CD4記憶T細胞和活化的肥大細胞（r=0.66）、未活化的樹突狀細胞和M1型巨噬細胞（r=0.54）、M1型巨噬細胞和活化的CD4記憶T細胞（r=0.49）呈正相關；而未活化的肥大細胞和活化的肥大細胞（r=-0.78）、記憶B細胞和活化的B細胞（r=-0.67）、未活化的CD4記憶T細胞和M2型巨噬細胞（r=-0.62）呈負相關（圖5A）。通過小提琴圖較為直觀地展現正常人群和OA患者滑膜組織免疫浸潤細胞差異，以P<0.05表示差異有統計學意義，漿細胞（P=0.019）、M0型巨噬細胞（P=0.038）和未活化的肥大細胞（P<0.001）在OA滑膜組織中含量較高，而未活化的CD4記憶T細胞（P<0.001）、活化的NK細胞（P=0.026）、活化的肥大細胞（P<0.001）在OA滑膜組織中均出現不同程度降低（圖5B）。

圖5 OA組織中免疫浸潤細胞相關性及差異小提琴圖Fig.5 Correlation analysis and violin plot of differences of immune cells in OA group

2.6 干預OA滑膜組織免疫浸潤機制的中藥挖掘通過Coremine Medical預測和挖掘潛在干預OA免疫浸潤的中藥，以P<0.01作為篩選條件，共得到123味中藥，與白細胞黏附、白細胞增殖、吞噬作用、抗原呈遞與表達、T細胞呈遞等BP及IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9等核心靶基因密切相關，將結果導入Cytoscape進行可視化分析，共得到節點137個，邊226條，青風藤、姜黃、雷公藤、藤黃等度值較高，提示與OA免疫浸潤機制密切相關，將這些中藥進行歸經，功效分析顯示OA免疫浸潤相關藥物均歸屬肝、腎、肺經，藥物功效以清熱類、補虛類、祛風濕類、活血化瘀類、解表類等為主，提示以上藥物可能作為OA免疫治療的潛在候選中藥（圖6）。

圖6 OA免疫浸潤機制潛在中藥功效及雷達圖Fig.6 Potential traditional Chinese medicine efficacy map and radar map of osteoarthritis immune infiltra?tion mechanism

3 討論

滑膜組織是OA免疫系統的主要入侵對象，本研究將正常人群和OA患者滑膜組織差異基因進行PPI分析，篩選出10個與OA免疫相關的核心靶基因IL-6、CXCL8、JUN、VEGFA、IL-1β、MMP9、ITGβ2、FOS、APOB、CXCL12，除ITGβ2和APOB基因未見報道以外，其他均發現參與OA滑膜免疫浸潤過程。VEGFA可介導早期OA滑膜中性粒細胞、淋巴細胞、肥大細胞浸潤誘導血管炎癥［12］；JUN和FOS在OA滑膜中高表達，參與炎癥細胞彌散浸潤［13］；MMP9參與CD4+T細胞調控巨噬細胞浸滑膜病程；IL在OA患者滑膜組織中主要由白細胞、巨噬細胞、滑膜纖維細胞等產生［14］，OA患者滑液、滑膜、軟骨下骨和軟骨中均不同程度增高［15］；此外，過多IL產生可抑制軟骨中Ⅱ型膠原生成，誘導基質金屬蛋白酶類產生，加速軟骨侵蝕［16］。本研究預測出的IL-6、IL-1β在OA進程中備受關注，已展開大量相關研究［17-18］。而趨化因子（CXCL）是一類具有化學趨化作用的細胞因子，主要由滑膜巨噬細胞產生，可對中性粒細胞、淋巴細胞、單核細胞等細胞產生趨化作用，且與血管炎發展密切相關，介導早期OA出現滑膜炎癥［19］。CXCL8可影響JAK-STAT、NF-κB和JNK-MAPK信號通路并增強其他促炎細胞因子表達促進細胞凋亡并抑制軟骨細胞增殖［20］；CXCL12可驅動樹突狀細胞從血液募集至發生炎癥反應的滑膜誘導OA慢性炎癥［21］。本研究篩選出的大部分核心靶基因與目前OA免疫相關研究結果相符，進一步證實研究結果的可靠性。對OA滑膜組織差異基因進行GO和KEGG富集分析，結果顯示上調基因與OA免疫機制密切相關，主要涉及多種免疫細胞，進一步說明OA是由多種炎癥細胞及因子共同參與的免疫過程，存在錯綜復雜的轉導過程，而KEGG富集預測出大部分信號通路均介導了OA免疫過程，這些通路在誘導OA滑膜炎癥因子滲出、免疫細胞浸潤、補體系統激活、體液免疫和細胞免疫過程中均發揮重要作用［22］。

為深入研究免疫細胞在OA滑膜組織的作用，采用CIBERSORT反卷積法對兩組免疫相關差異基因進行篩選，CIBERSORT反卷積算法是基于線性支持向量回歸（SVR）的可靠機器學習方法，廣泛用于評估22種免疫細胞相對含量及動態調節過程，在噪聲、未知混合物含量和密切相關的細胞類型方面優于其他方法。最終結果顯示漿細胞、M0型巨噬細胞和未活化的肥大細胞在OA滑膜組織中含量較高，而未活化的CD4記憶T細胞、活化的NK細胞、活化的肥大細胞在OA滑膜組織不同程度降低，以上研究結果部分得到相關基礎研究證實，如OA滑膜組織中肥大細胞浸潤程度較高，與OA患者結構損傷有關，其機制可能與肥大細胞釋放大量TNF-α及IL-6等細胞因子誘導炎癥反應，刺激滑膜增生和軟骨侵蝕［23］；記憶CD4+T細胞積累是關節局部炎癥反應過程中的常見現象，可能與慢性OA形成有關［24］；調節性T細胞在OA滑膜中含量豐富，其水平與炎癥因子水平相關［25］；而滑液中過量的CD19+CD24hiC7+B細胞可加速軟骨破壞［26］；M1型巨噬細胞產生炎癥因子與產生抗炎因子的M2型巨噬細胞比例失衡是促進OA進展的關鍵機制之一，調節M2型巨噬細胞有望對延緩OA起積極治療作用［27-28］。上述文獻證據結合課題組分析結果表明巨噬細胞、漿細胞、肥大細胞等在OA中發揮重要作用，應成為進一步研究重點。

OA屬于中醫學“痹癥”范疇，中醫對其認識較早，長時間臨床實踐顯示中藥具有良好有效性及安全性，且中藥具有多成分、多靶點、多途徑作用機制，從中藥中尋找針對OA免疫機制整體調節的藥物應作為研究重點［29-31］。OA中醫病機屬于年老體衰，肝腎虧虛，筋骨失養，外感風寒濕邪氣，痹阻脈絡，流注關節［32］，證屬“本虛標實”，以肝腎虧虛為主，風寒濕邪為標實。本次預測藥物的歸經和功效分析，大部分藥物歸經多屬于肺、肝腎經且藥物多為補虛類、祛風濕類、活血化瘀類、解表類等，進一步說明本研究藥物預測結果符合OA中醫病機，且本研究預測藥物以藤類藥物居多，中醫學受“取類比象”思維影響，認為藤類藥物形條達，纏繞蔓延，如同人體經絡縱橫交錯走行，具有通暢經絡、舒通筋脈作用，善治痹證久治不愈，甚則變形之頑痹，完全符合OA病癥特點。本研究從免疫浸潤角度將OA差異基因中免疫相關BP及作用于PPI核心靶基因的中藥進行反向預測，篩選出青風藤、姜黃、雷公藤、藤黃等中藥在免疫調節干預OA過程中處于核心地位，現代藥理學也證實以上藥物可通過調節免疫細胞發揮OA治療作用，如青風藤可抑制TNF-α、IL-β、IL-6等炎癥因子產生，減少炎癥反應對關節軟骨的損傷［33］。此外，青風藤還可抑制滑膜組織血管新生和血管炎癥，減少免疫細胞浸潤對滑膜的刺激；姜黃中有效成分姜黃素具有強大的抗氧化和抗炎作用，可抑制IL-1β誘導的NF-κB介導的炎癥和細胞凋亡，保護軟骨［34］；雷公藤除抗炎、鎮痛、抗血管新生外，還對滑膜有免疫調節功能，抑制T細胞增殖、誘導T細胞凋亡、恢復Th1/Th2細胞平衡［35-36］。以上文獻支持本研究結果的可信度，也為中藥干預骨關節炎免疫機制候選中藥提供參考。

目前美國風濕病學會（ACR）、歐洲抗風濕聯盟（EULAR）和中國2018版《骨關節炎診療指南》均推薦針對OA治療的藥物，包括非甾體類抗炎藥（NSAIDs）、阿片類藥、糖皮質激素、抗焦慮藥、軟骨保護藥等，促使消化道出血、腦、腎和心血管疾病、關節感染、藥物依賴性等風險性增加，長期使用飽受爭議，治療目的主要是減輕或消除疼痛，糾正畸形，改善或恢復關節功能，但卻未對免疫系統進行調節，遠期效果仍不理想［37-40］。對于晚期OA患者關節軟骨缺損，組織工程學方法的自體軟骨移植修復是目前研究熱點［41］。但關節周圍組織免疫微環境引發的自身免疫排斥反應也是組織工程材料關節植入需考慮的關鍵問題。中醫藥是中國人民幾千年來的智慧結晶，其獨特魅力受到廣泛關注［42-43］。因此，從中草藥中尋找干預OA免疫失衡的潛在中藥，對早期OA患者減緩病程具有重要意義，對于晚期OA患者關節組織工程材料的特性選擇也具有指導意義。

本研究依賴于GEO數據庫，是既往發布數據集的二次挖掘和分析，實驗結果可能與既往實驗結論不同，原因可能為樣本量過少導致數據分析偏倚。此外，CIBERSORT反卷積算法分析是基于有限的遺傳數據，這些數據可能會偏離細胞異型相互作用、疾病的誘發因素以及疾病表型的可塑性，導致結果不準確。雖然本研究結果得到了多數文獻研究支持，但研究結果的可靠性仍需進一步實驗驗證。