基于單細胞轉錄組測序分析骨關節炎軟骨細胞分化的分子機制*

徐文飛 梅其杰 明春玉 段戡 袁長深 郭錦榮 胡琪 曾超 梁曉輝

骨關節炎是一種常見的慢性退行性疾病，常伴有關節腫脹和軟骨損傷等變化[1]；因其具有高發病率和致殘率，將會成為全球主要致殘疾病之一[2-4]。然而目前骨關節炎的發病機制仍存在爭議，闡明發病機制成為迫切需要解決的問題。而軟骨細胞對于骨關節炎具有至關重要的作用，其中軟骨細胞肥大或衰老擔任著重要角色[5]，故軟骨老化或過度使用將導致軟骨細胞穩態喪失，進而加重骨關節炎病情[6]。因此，深入研究骨關節炎與軟骨細胞之間關系，從分子水平闡述骨關節炎軟骨細胞發病機制成為新熱點和機遇。

單細胞轉錄組測序（single-cell RNA sequencing,scRNA-seq）是通過對組織中細胞進行獨立分析，來揭示細胞亞群之間相互作用的關系，以及細胞內基因之間相互調控的過程，較傳統的RNA 測序技術分析更加詳細，對疾病研究具有重要作用[7-9]。當前單細胞轉錄組測序已廣泛運用于多種疾?。òü顷P節炎），尚缺乏單細胞測序聯合生物信息學（基因芯片）關于骨關節炎軟骨細胞的分析，因此本研究在其他研究的基礎上將兩者相結合，共同探討骨關節炎軟骨細胞的作用機制，以期更好地治療骨關節炎。

1 材料與方法

1.1 基因芯片數據來源

通過GEO 數據庫（https://www.ncbi.nlm.nih.gov/geo/）下載骨關節炎軟骨細胞基因芯片數據集GSE104782，其所采用平臺為GPL20301，該芯片數據包括10 個骨關節炎樣本，1 464個軟骨單細胞測序信息。

1.2 質控與數據過濾

運用R 語言對GSE104782 進行質控與數據過濾，將重復基因取均值來進行數據質控，并轉換為“Seurat”對象，進一步通過“Percentage Feature Set”函數計算線粒體基因百分比，進而將數據過濾，從而獲得基因特征及測序深度，并對數據標準化，提取細胞間變異系數較大的基因。

1.3 主成分分析

主成分分析（principal component analysis, PCA）是一種最常用的線性降維方法，能確定可用維度并篩選相關基因。運用R語言對數據進行標準化預處理，進而使用PCA分析，可獲得每個PCA成分相關基因、P值及均勻分布情況。

1.4 TSNE聚類分析與Marker基因

TSNE 聚類分析是可視化高維數據工具，可在二維或三維空間里展示聚類結果。采用R語言對單細胞測序所得數據分組，優化標準模塊，進行可視化聚類分析；同時，通過尋找數據差異表達特征獲得Marker基因，并對各個成分中所含Marker基因進行分析。

1.5 細胞注型及分化軌跡分析

基因芯片數據集（GSE104782）所包含細胞類型均為軟骨細胞，采用R 語言中的SingleR、Monocle 包進行細胞注釋，進一步驗證該數據集細胞類型，并運用偽時間分析法來闡述軟骨細胞分化軌跡，同時查閱相關文獻深入分析骨關節炎軟骨細胞的具體分化途徑。

1.6 GO及KEGG分析

利用R 語言工具對軟骨細胞的Marker 基因進行GO 富集及KEGG信號通路分析，獲得軟骨細胞相關功能及信號通路，其中顯著性基因富集臨界值設定為P＜0.05。

1.7 蛋白互助網絡分析

選取細胞成分9 中所含Marker 基因進行蛋白互助網絡分析，將Marker 基因導入STRING（Version 11.0）在線分析軟件中（https://string-db.org/cgi/input.pl），獲得蛋白質-蛋白質相互作用（protein-protein interaction,PPI），并將其結果導入Cytoscape 3.7.2 軟件，運用插件cytoHubba 中Degree 算法獲得前10 位關鍵靶基因并構建PPI 網絡。

2 結果

2.1 質控與數據過濾結果

運用R 語言對基因芯片數據集GSE104782 進行分析，發現基因測序深度與基因特征具有顯著聯系，隨著測序不斷深入，基因特征也更加豐富（見圖1）；該數據集獲得基因17 380 種，篩選出COL1A1、IBSP、PRG4、TAC1、SPP1、HSPA6、COL10A1、CCL2等1 500 種重要基因（見圖2）；同時分析每種樣本基因特征及研究深度，發現基因特征主要2 500～5 000范圍內，而研究深度在0 ～ 500 000之間（見圖3）。

圖1 測序基本特征與深度關聯

圖2 測序基因表達

圖3 測序的基本特征與深度

2.2 PCA分析結果

通過PCA分析相關數據發現，軟骨細胞由20個主成分組成，按P＜0.05 進行統計學分析，結果顯示20 個主成分具有統計學差異，均可納入研究（見圖4）。以PC1 與PC2繪制主成分分析圖形，發現軟骨細胞主成分之間具有一定的關聯（見圖5）。同時，進一步分析PC1、PC2、PC3 及PC4 的相關基因，繪制主成分熱圖，展示樣本中主成分所含基因分布及聚類情況（見圖6、圖7），結果顯示PC1 中COL1A2基因、PC2中OGN基因、PC3中HLA-DRA基因及PC4中SOD2基因在樣品組織中顯著高表達。

圖4 主成分分析結果

圖5 PC1與PC2的關聯

圖6 前4位主成分基因分布曲線圖

圖7 前4位主成分基因分布熱圖

2.3 TSNE聚類分析結果

本研究進一步對單細胞測序所得的數據進行TSNE聚類分析，共獲得9 類不同顏色的細胞聚類，分別代表軟骨細胞中不同成分（見圖8）；同時以校正P=0.05 為標準進行差異分析，獲得1 886 個Marker 基因并以熱圖及氣泡圖的方式顯示有些基因在不同成分中的表達情況（見圖9、圖10），同時圍繞關鍵基因LECT1、TGFBI分析表達及分布情況，發現LECT1（logFC=2.33，P=0.045）、TGFBI（logFC=1.66，P=0.044）均在細胞成分9 種高表達（見圖11）。

圖8 軟骨細胞成分分布

圖10 Marker基因分布氣泡圖

圖11 關鍵基因的表達與分布

2.4 細胞注釋及細胞軌跡分析結果

運用R語言工具中SingleR包對獲得的9種細胞成分進行注釋，發現均為軟骨細胞，與研究內容相符；通過偽時間分析法進一步闡述軟骨細胞的分化軌跡（見圖12），發現共有6種分化途徑，其中顯示Cluster0細胞聚類單一，為最早分化細胞，查閱相關文獻發現Cluster0 可能為軟骨祖細胞，經關鍵調控因子分化為增殖軟骨細胞、肥大前軟骨細胞和肥大軟骨細胞等不同軟骨細胞成分[10]。

圖12 軟骨細胞的分化軌跡

2.5 GO和KEGG通路富集分析

GO 是一組大量的基因注釋術語，可用于識別高通量基因組或轉錄組數據的特征生物學屬性。采用R 語言工具對Marker 基因進行GO 功能富集分析，以校正P＜0.05 進行篩選，主要涉及軟骨發育、結締組織發育等生物過程（biological process, BP）；涉及含膠原蛋白的細胞外基質、內質網內腔等細胞組分（cellular component, CC）；涉及肝素結合、糖胺聚糖結合等分子功能（molecular function,MF）（見圖13）。KEGG 通路分析：主要與蛋白質消化吸收、ECM受體相互作用、人乳頭瘤病毒感染、粘著斑、癌癥中的蛋白多糖、補體與凝血級聯反應、PI3K-Akt信號通路、礦物吸收等信號通路有關（見圖14）。

圖13 GO富集分析

圖14 KEGG信號通路分析

2.6 差異基因的PPI網絡分析

STRING 數據庫是通過已知和預測的蛋白質相互作用數據所組成，本研究選取細胞成分9 種所含Marker 基因進行蛋白互助網絡分析，運用STRING 數據庫獲得PPI 網絡圖，并且將其所得結果導入Cytoscape 3.7.2 軟件，運用插件cytoHubba 中Degree 算法發現前5 位關鍵靶基因分別為COL1A1、COL2A1、VEGFA、COL1A2、DCN（見圖15）。

圖15 蛋白質互助網絡分析

3 討論

骨關節炎是一種由多種因素引起的慢性退行性疾病[11]，主要病理表現為軟骨細胞凋亡和細胞外基質丟失，故維持軟骨正常內環境穩定是治療骨關節炎的關鍵所在[12]。軟骨是一種無血管、淋巴和非神經支配的組織，而軟骨細胞作為軟骨唯一的細胞類型，具有內在修復潛力[13]；然而如何保護及修復軟骨細胞仍然十分艱難，故本研究采用單細胞轉錄組進一步探究軟骨細胞基因表達及分化途徑，獲得軟骨保護及修復方法，為更好治療骨關節炎提供新思路。

本研究通過單細胞測序及生物信息學分析發現17 380種軟骨細胞基因，其中包含LECT1、TGFBI等1 886種關鍵基因。其中LECT1在骨關節炎軟骨細胞中呈高表達狀態，其作為軟骨蛋白-1編碼可抑制軟骨細胞生長，進而影響骨關節炎病情[14]；而TGFBI不僅調節合成代謝和分解代謝標志物的表達，而且還降低炎癥因子[15]。同時進一步研究發現軟骨細胞共包含多種軟骨細胞類型及分化途徑，其中軟骨祖細胞為最早分化細胞，可經關鍵調控因子分化為增殖軟骨細胞、肥大前軟骨細胞和肥大軟骨細胞等不同成分，因此本研究推測干預軟骨祖細胞可能為修復及保護軟骨細胞提供新的方法。

采用GO、KEGG及蛋白互助網絡分析軟骨細胞具體機制，GO 富集分析發現軟骨細胞基因涉及軟骨發育、結締組織發育等生物過程（BP）；涉及含膠原蛋白的細胞外基質、內質網內腔等細胞組分（CC）；涉及肝素結合、糖胺聚糖結合等分子功能（MF）。其中，軟骨細胞及結締組織發育可起著修復軟骨作用，且關節軟骨退變與細胞外基質降解密切相關[16]；分泌型肝素結合肽通過刺激細胞外基質合成、減少降解基質金屬蛋白酶和誘導其抑制劑來刺激軟骨細胞聚集和增殖，增強軟骨形成[17]；糖胺聚糖鏈連接的蛋白聚糖存在于軟骨細胞表面和細胞周圍軟骨基質中，可參與細胞-細胞和細胞-基質的相互作用[18-19]。因此，本研究認為從生物過程、細胞組分、分子功能角度保護軟骨細胞可能會具有顯著的療效。

KEGG通路分析主要與蛋白質消化吸收、ECM受體相互作用、人乳頭瘤病毒感染、粘著斑、癌癥中的蛋白多糖、補體與凝血級聯反應、PI3K-Akt信號通路、礦物吸收等信號通路有關；其中發現黃芪甲苷可通過上調ECM受體相互作用來促進軟骨細胞的快速增殖，對關節軟骨退變發揮了保護作用[20]；而補體與凝血級聯及金黃色葡萄球菌感染主要參與天然免疫反應，可能提高軟骨細胞的修復作用[21-22]；同時，在骨關節炎軟骨細胞中PI3K-AKT 通路處于下調狀態，通過激活PI3K/AKT 可減少軟骨細胞損傷，進而促進骨關節炎軟骨細胞再生[23]。

通過Cytoscape 3.7.2 軟件獲得5 個關鍵靶基因，包含COL1A1、COL2A1、VEGFA、COL1A2、DCN。在軟骨受到損傷后，可促進膠原蛋白（COL1A1、COL2A1、COL1A2）下調，使得胞外囊泡-殼寡糖偶聯物逆轉軟骨損傷而起到保護作用[24]；血管內皮生長因子（vascular endo‐thelial growth factor, VEGF）是血小板衍生生長因子家族的重要成員之一，主要分布于關節的成骨細胞、軟骨細胞及滑膜內成纖維細胞等組織[25]，其主要通過促進血管內皮增殖、促進腫瘤淋巴管生長、上調抗凋亡蛋白Bcl-1表達促進骨關節炎的形成[26]；而核心蛋白聚糖（decorin, DCN）可延緩軟骨蛋白聚糖碎片的丟失和軟骨表面的纖維化，從而對改善軟骨退變起到保護作用[27]。

綜上所述，本研究運用單細胞測序聯合生物信息學方法研究骨關節炎軟骨細胞分化的具體機制，結果顯示軟骨細胞分為多種類型且具有多種分化途徑，可通過調節關鍵基因（LECT1、TGFBI）來干預軟骨細胞再分化，進而為保護軟骨細胞起著重要作用。然而本研究仍存在不足之處，首先未對結果采用試驗進行驗證；同時納入單細胞基因芯片樣本較少，可能存在一定偏倚，需后續擴大樣本進一步研究，因此期待后續能更深入研究軟骨細胞相關機制，為治療骨關節炎提供堅實基礎。