基于機器學習探討不同年齡階段半月板損傷鐵死亡作用機制和免疫浸潤研究

黃柯琪，陳躍平，曾 浩，李加根，陳尚桐

（廣西中醫藥大學附屬瑞康醫院，廣西 南寧 530000）

近期研究表明半月板損傷與鐵死亡之間存在密切關聯。鐵死亡是一種鐵依賴性的細胞死亡方式，其發展受到脂質過氧化的積累影響。鐵死亡在多種疾病中扮演重要角色，包括癌癥、壞死性炎癥性疾病以及器官損傷、手術和生殖系統的器官變化等［6-8］。研究發現，半月板軟骨細胞通常在炎癥和鐵超負荷條件下發生鐵死亡。越來越多的證據表明，鐵死亡可能成為治療軟骨損傷的新靶點［9-11］。然而，尚未明確不同年齡段半月板損傷診斷和治療中與鐵死亡相關的潛在生物標志物。因此，有必要積極開展差異鑒別研究，以期未來能提供更加個性化和有針對性的診療策略。

高通量技術的進展已經使得基因芯片分析成為鑒定差異表達基因的重要方法。已經有研究使用基因微陣列分析結合各種生物信息學分析和機器學習來識別多種疾病發病機制的關鍵基因，這些關鍵基因有可能成為疾病治療和診斷的潛在生物標志物。在本研究中采用了免疫浸潤分析來更準確地識別兩個年齡組之間半月板損傷的差異。通過綜合生物信息學分析和LASSO 機器學習，并結合免疫浸潤分析，鑒定了與鐵死亡相關的基因，這些基因有望成為不同年齡段半月板損傷治療的新靶點。這一發現為精準化治療半月板損傷提供了新的方向，為未來的實驗研究和臨床治療提供了有價值的線索。

1 資料與方法

1.1 數據采集

從GEO 篩選微陣列數據集。搜索關鍵字“meniscusinjury”，搜 索 方 法 如 下：“meniscusinjury”［MeSH 術語］，過濾掉非人類相關的基因芯片，最后得到GSE191157 基因芯片。GSE191157 包含8例行關節鏡半月板修復手術切除的損壞半月板表達譜，無明顯軟骨變性、骨關節炎、類風濕性關節炎等代謝性疾病的患者，其中年輕組4 例患者（年齡＜50）和老年組4 例患者（年齡≥50）。

1.2 差異分析

GSE191157 使用Perl 軟件整理相關矩陣數據，使用 R 語言中“limma”包對表達矩陣進行標準化處理和差異基因分析，以logFC 的絕對值＞1 且P＜0.05 篩選差異基因，并使用ggplot2 繪制火山圖和熱圖。

1.3 鐵死亡相關基因及GO、KEGG 功能分析

從“FerrDb”數據庫中收集和整理鐵死亡相關基因，使用DEGs 和鐵死亡相關基因取交集，通過維恩圖可視化。然后整理DE-FRGs 相關基因在GSE191157 的 表 達 矩 陣，使 用R 包“clusterProfiler”進行GO 和KEGG 功能富集分析，將相關富集分析結果可視化。

1.4 WGCNA 分析

對標準化處理過后的基因表達矩陣使用R 語言 中“WGCNA” 包 和“limma”包 進 行WGCNA 分析，分別進行刪除離群樣本，樣本聚類構建拓撲重疊矩陣和無尺度網絡，動態模塊分析及合并，以基因與模塊相關性＞0.9 為過濾條件，找出臨床信息相關性最強的基因模塊，整理模塊中的基因為下一步篩選核心基因做準備。

1.5 交集取核心基因及功能富集分析

將GSE191157 矩陣分析的差異基因、鐵死亡相關基因和WGCNA 中與臨床最相關模塊的基因取交集作為不同年齡段半月板損傷的核心基因，用venn 圖可視化。

1.6 LASSO 鑒定核心基因

將得到的venn 圖的交集基因數據整合，數據分為訓練集和測試集。訓練集用于訓練LASSO 模型，而測試集用于評估模型的性能。使用交叉驗證選擇合適的lambda 值，根據選擇的lambda 值構建LASSO 模型。從訓練好的LASSO 模型中提取非零系數，這些非零系數對應于被選擇的核心基因。

心肌組織中ATP等能量物質容易在各種酶的作用下被破壞，本研究從心肌組織離體后凍存（液氮迅速冷凍，-80℃冰箱暫存）、樣本前處理（預冷處理液、冰水浴下均漿、4℃離心）和上機檢測（柱溫15℃）全過程、關鍵實驗環節都保持在低溫下進行，最大限度地保持心肌組織中能量物質的穩定性，為后續研究奠定物質基礎［2，6］。這也是本研究建立一種HPLC法同時測定大鼠心肌組織中能量代謝物質的含量的前提條件和關鍵點，其保證了檢測結果的準確性和可靠性。

1.7 GSEA 富集分析

此分析使用R 語言中的“fgsea”包執行。所使用的基因集為KEGG 公共數據庫， 基因根據其在樣本中的表達水平進行排序，計算一個累積型的得分。通過對基因表達矩陣進行隨機化來生成多個隨機化數據集，形成一個隨機分布。將計算出的得分與隨機分布進行比較，計算基因集合的富集分數，用于衡量基因集合在樣本中的富集程度。通過比較富集分數分布和真實數據的富集分數，設P＜0.05 提示基因集合在樣本中顯著富集。然后將相關的通路可視化。

1.8 SSGSEA

ssGSEA 算法基于29 個免疫基因集，包括不同免疫細胞類型、功能、途徑和檢查點相關的基因。通過R 包（GSVA，GSEABase 和limma）采用ssGSEA 算法來全面評估研究中包含的每個樣本的免疫學特征和每個免疫細胞的相對浸潤豐度。

2 結果

2.1 年輕組和老年組的對照差異分析

使用“limma”包分析年輕組和老年組的RNA-seq 數 據，設 計FDR 閾 值＜0.05 和|log2FC|使用＞1，總共通過篩查確定了2 517 個DEG，其中661個被下調和1 856 個是上調基因（圖1A）。然后從“FerrDb”數據庫中提取了鐵死亡相關基因，與上述差異基因取交集并用venn 圖可視化，結果可見左邊的計數（2 456 個基因）是指GSE191157 的DEG，中間的計數（61 個基因）是指鐵死亡相關的DEG，右邊的計數（504 個基因）指的是獨特的鐵死亡基因（圖1B）。然后將鐵死亡相關DEG 基因在年輕組和老年組的表達情況用熱圖可視化（圖1C）。

圖1 差異基因火山圖、熱圖與鐵死亡相關Venn 圖Fig 1 Differential gene volcano map， heat map， and Venn map related to iron death

2.2 鐵死亡相關基因及GO、KEGG 功能分析

不同年齡段半月板損傷鐵死亡相關基因的GO 富集分析結果顯示，共富集到 616 個條目，其中生物過程（biological process，BP）579 條，細胞組分（cell components，CC）14 條，分子功能（molecularfunction，MF）24 條（圖2A）。生物過程主要涉及細胞代謝調整、細胞生長和分裂、細胞死亡、信號傳導、基因表達調控，細胞組分主要涉及過氧化物酶體、微小體、頂端質膜、過氧化物酶體膜、微小體膜、細胞質應激顆粒，分子功能主要涉及抗氧化活性、蛋白質絲氨酸激酶活性、核受體活性、配體激活的轉錄因子活性、過氧化物酶活性、氧化還原酶活性。KEGG 富集分析結果顯示，共映射出37 條通路，主要與病毒蛋白與細胞因子和細胞因子受體的相互作用、細胞因子-細胞因子受體相互作用、趨化因子信號通路、Toll 樣受體信號傳導途徑、IL-17 信號通路和腫瘤壞死因子信號通路等有關（見圖2B）。

圖2 鐵死亡相關基因及GO、KEGG 功能分析Fig 2 Analysis of iron death related genes and GO， KEGG functions

2.3 共表達模塊的構建與關鍵性狀模塊基因的鑒定

WGCNA 是發現基因與臨床表型的好方法，通過結合DEGs 挖掘出與臨床疾病表型相關的關鍵基因。GSE191157 沒有異常值，數據集中未刪除任何樣本。對本數據集中位絕對偏差的前75%基因進行了WGCNA 篩選（圖3B）。當軟閾值功率等于6時，構建了β=6 軟閾值能力的無標度網絡，并根據表達譜將基因分類為不同的模塊。通過平均鏈接聚類（圖3A）。模塊的顏色越深，提示該模塊相關性越強。在分析了正相關系數后，選擇了相關系數最強的模塊。結果顯示青綠色模塊（r=0.93，P＜0.01）和黃色模塊（r=-0.9，P=0.002）具有最顯著的基因意義，其他沒有顯著的模塊。最后分別提取了青綠色模塊中的2 072 個基因和黃色模塊中的397 個基因。這2 469 個基因被認為是與臨床表型相關的關鍵基因（圖3C）。

圖3 通過共表達網絡構建及特異性模塊鑒定分析相關的潛在鐵死亡基因Fig 3 Construction of a co expression network and identification and analysis of potential iron death genes through specific modules

2.4 LASSO 機器學習篩查核心基因

WGCNA 結果得到的青綠色模塊和黃色模塊共有2 469 個基因，再去與鐵死亡相關的半月板的61 個差異基因（F-DEG）取交集顯示在維恩圖中，兩者的交集部分共得到46 個基因（圖3D）。提取46 個基因的表達譜，用于構建LASSO 模型。在 LASSO 模型中確定了產生最小分類誤差的最佳λ（λ=4）?；?于λ值4，篩 選 出F-DEG 的LASSO 系 數 譜（圖4A）。隨后鑒定出4 個非零系數的核心基因，包括BEX1、MMP13、PTPN18和SLC38A1確定為不同年齡段半月板損傷核心基因（圖4B）。

圖4 LASSO 機器學習和GESA 通路富集分析Fig 4 LASSO machine learning and GESA pathway enrichment analysis

2.5 GESA

GSEA 分析顯示，富集途徑主要與Toll 樣受體信號通路、FOXO 信號通路、mTOR 信號通路、ECM-受體相互作用和信號通路、免疫球蛋白G 生產者信號通路、草酰乙酸和乙酰輔酶A 的降解有關（圖4C）。差異途徑主要集中細胞代謝和能量供給、細胞外基質與細胞受體的相互作用、細胞黏附、遷移和信號傳遞、蛋白質在內質網的折疊、修復、加工和蛋白質合成能在細胞應激和修復中起作用。GSEA 結果表明，上述信號通路可能參與半月板損傷的進展，并提示了這些關鍵通路的預后價值，可能會在調節半月板損傷患者進展和不良預后中起重要作用。

2.6 免疫浸潤分析

SSGSEA 算法用于檢查老年組半月板與年輕組對照在免疫細胞浸潤差異方面的關聯。利用R語言對GSE191157 數據集進行的29 種免疫細胞相關性分析。創建了一個熱圖來描述免疫細胞之間的差異相關模式，如圖5A 所示，紅色越深代表二者負相關性越高，藍色越深代表二者正相關性越高，正相關性較高的組合有效應記憶CD8 T 細胞與靜CD56dim 自然殺傷細胞，記憶B 細胞與未成熟樹突狀細胞等，負相關性較高的組合有自然殺傷T 細胞和激活的CD4 T 細胞，中央記憶CD4 T 細胞和效應記憶CD8 T 細胞。此外，通過探索核心基因相關性分析證實，說明樞紐基因表達與免疫細胞浸潤之間的相關性，MMP13在效應記憶CD8 T 細胞、CD56dim 自然殺傷細胞、中央記憶CD4 T 細胞中呈現負相關，而SLC38A1、BEX1、PTPN18這三個核心基因則呈正相關（圖5C）。

圖5 免疫浸潤分析Fig 5 Immunoinfiltration analysis

差異熱圖顯示老年組半月板損傷的中激活的B細胞、激活的CD4 T 細胞和效應記憶CD8 T 細胞等免疫細胞大部分浸潤豐度高于年輕組（圖5B）。免疫差異細胞結果顯示，老年組中激活的B 細胞，激活的CD4 T 細胞，效應記憶CD8 T 細胞在老年組中顯著低于年輕組，未成熟樹突狀細胞、漿細胞樣樹突狀細胞顯著高于年輕組（圖5D）。

3 討論

半月板是人體膝關節內具有彈性和韌性彎曲纖維軟骨結構，主要功能包括增強膝關節的穩定性，減少關節受力，提供緩沖作用，同時在膝關節屈伸、旋轉和側向運動中發揮重要作用［12］。隨著年齡的增長，身體供血、關節內微環境、半月板的磨損程度都會出現不同程度的改變，因此可知不同年齡段的半月板損傷存在著差異［13］。面對兩人群半月板損傷的差異，需要采用不同針對性的治療策略。鐵死亡是由鐵依賴性脂質過氧化介導的一種新型調節性細胞死亡形式［14］。鐵死亡與多種人類疾病密切相關，不僅在許多人類惡性腫瘤中觀察到鐵死亡，同時還與骨關節炎、類風濕關節炎、骨質疏松疾病有密切關系［15］。既往研究發現鐵死亡參與軟骨細胞凋亡的進展［16］。低濃度的鐵離子可以促進成骨細胞前體細胞（MC3T3-E1）的生長，而高濃度的鐵離子會抑制它們的生長并增加ROS 水平［17］。鐵超負荷可以在一定程度上抑制成骨細胞的活性，從而影響它們的分化過程。同時，它還可以激活破骨細胞的分化并引起骨破壞［18］。但是具體的作用方式尚不清楚，靶向鐵死亡或許能提供針對不同人群半月板損傷的有效治療方法。文章在GEO 數據庫中檢索得到年輕組和老年組半月板損傷的表達基因，鐵死亡基因相關數據從FerrDb 數據庫的得到，經過生物信息學分析后，兩者取交集得到不同人群半月板損傷的差異基因，并通過機器學習篩選出鐵死亡相關的標志基因，然后進行免疫浸潤分析，對鐵死亡在不同人群半月板損傷中的作用機制進行初步探討，明確不同年齡段半月板損傷鐵死亡的生物標志物。

鐵死亡相關基因GO 分析顯示，生物過程主要涉及內質網應激誘導的內源性凋亡信號通路的調控、細胞代謝調整細胞死亡、信號傳導和基因表達調控。其中最顯著的分組過程是涉及內質網應激的調控，研究發現過度和延長的內質網應激可導致軟骨細胞凋亡，而輕度內質網應激可通過GRP78 途徑激活自噬并保護軟骨細胞免于凋亡［19］。還有研究發現內質網應激可以通過調節破骨細胞的分化和功能以及骨細胞的存活來調節骨穩態［20］。另外在骨髓間充質干細胞中內質網應激過度和自噬功能障礙可促進炎癥介導的骨質流失［21］。細胞組分主要涉及過氧化物酶體、微小體、頂端質膜、過氧化物酶體膜、微小體膜、細胞質應激顆粒。其中多個顯著的細胞組分與氧化物酶體相關，研究發現氧化物酶體與軟骨細胞的死亡有著密切關系。在對小鼠軟骨細胞的實驗中，發現過氧化物酶體增殖物激活的受體γ（PPAR-γ）表達增加可以緩解由高同型半胱氨酸引起的促炎介質和基質金屬蛋白酶13（MMP13）的表達增加，進而抑制軟骨細胞死亡［22］。另外一個小鼠實驗發現過氧化物酶體增殖物激活受 體γ 共 激 活 劑1α 和FoxO3A 通 過AMP 活 化 蛋 白激酶活化蛋白激酶介導軟骨保護［23］。KEGG 富集分析顯示，主要與P13K 信號通路、趨化因子信號通路、Toll 樣受體信號傳導途徑、IL-17 信號通路和腫瘤壞死因子信號通路等。研究發現miR-7/EGFR/MEGF9 軸通過激活PI3K/AKT/mTOR 信號通路引起骨關節炎軟骨降解［24］。在小鼠實驗中，關節內注射IL-17 中和抗體可減少關節軟骨變性并降低衰老標志物Cdkn1a的表達，主要與關節中的IL-17 的降 低 和IL-4 表 達 的 增 加 有 關［25，26］。上 述 解 釋 有 關于鐵死亡在不同人群半月損傷中的作用機制，本研究還通過WGCNA 共表達網絡和機器學習算法進一步篩選出BEX1、MMP13、PTPN18和SLC38A1四個關鍵基因，然后對其進行免疫浸潤相關性分析。

MMP13 是靶向軟骨降解的主要酶，通過其切割Ⅱ型膠原蛋白的特殊能力參與軟骨降解的主要基質金屬蛋白酶［27］。關節軟骨合成代謝和分解代謝之間的平衡受復雜的因素網絡調節，但主要由MMPs 及其金屬蛋白酶（TIMPs）的內源性組織抑制 劑 維 持［28］。有 研 究 發 現CRISPR-Cas9 介 導 的MMP13 位點的DNA 編輯可以降低MMP13 的表達，可用于有效地建立人關節軟骨細胞（健康和骨關節炎）的細胞群，顯示出軟骨保護作用［29］。 Wang等［30］證明關節內注射摻入果素的熱凝膠通過下調MMP13 的表達來抑制OA 兔的關節炎癥。通過上述研究得知MMP13 可能是治療半月板損傷的關鍵靶基因。而PTPN18、BEX1、SLC38A1 在半月板、軟骨等方面相關的研究文獻不多，但是與鐵死亡關系密切。PTPN18 敲低增加了細胞內ROS 水平，并下調了GPX4 和xCT 表達，PTPN18 的沉默可能通過靶向p-p38/GPX4/xCT 軸誘導KLE 子宮內膜癌細胞中的鐵死亡［31］。PTPN18、BEX1、SLC38A1 能否成為不同年齡段半月板損傷診斷和治療靶點還需要更深入的實驗研究。

本研究采用SSGSEA 算法對老年組和年輕組的免疫細胞浸潤類型進行評估，發現多種免疫細胞亞型與半月板損傷的重要生物學過程密切相關。老年組中激活的B 細胞，激活的CD4 T 細胞，效應記憶CD8 T 細胞在老年組中顯著低于年輕組，未成熟樹突狀細胞、漿細胞樣樹突狀細胞顯著高于年輕組。 此 外，通 過 對BEX1、MMP13、PTPN18 和SLC38A1 和免疫細胞之間的相關性分析發現，MMP13 在 效 應 記 憶CD8 T 細 胞、CD56dim 自 然 殺傷細胞、中央記憶CD4 T 細胞中呈現負相關，而SLC38A1、BEX1、PTPN18這三個核心基因則呈正相關。T 細胞作為骨關節炎滑膜中的主要免疫細胞之一，可在早期被檢測出來［32，33］。但滑膜內襯層主要含有 CD4+ 和 CD8+ T 細胞，內膜層主要含有 CD4+ T 細胞，滑膜T 細胞浸潤可能在骨關節炎發病機制中承擔重要角色，誘發癥狀性骨關節炎［34，35］。T 細胞可能還通過 CD40/ CD40L 信號軸活化單核細胞和巨噬細胞，從而通過免疫細胞間的相互作用，影響骨關節炎的發生、發展［34］，此外，免疫應答的調節與細胞衰老亦關系密切［36］。因此軟骨細胞的衰老程度在介導的免疫反應方面存在差異，與不同年齡段半月板損傷之間存在差異相符，這值得進一步研究。

綜上所述，本研究運用生物信息學方法和機器學 習 算 法 共 篩 選 出BEX1、MMP13、PTPN18和SLC38A1四個鐵死亡相關的關鍵基因，這些關鍵基因可能與不同年齡段半月板損傷的差異表達有著緊密關系，同時分析結果提示鐵死亡在不同年齡段半月板損傷免疫浸潤方面發揮重要作用，對基礎研究具有一定的參考意義。但此次研究仍存在一定的局限性：第一，芯片數據樣本量較少，所以結果可能出現一定的偏差，第二，雖然不同年齡段半月板損傷相關的鐵死亡基因已經篩出，但對于篩出的基因研究較少，其具體的作用機制仍然有待解釋，這仍要后續研究進一步深入挖掘。

作者貢獻度說明：

黃柯琪：文章設計、數據分析及手稿撰寫；曾浩：相關文獻搜集與整理，數據分析；李加根、陳尚桐：部分文獻搜集與文章格式檢查；陳躍平：提供整體思路和修改意見。

所有作者聲明不存在利益沖突關系。