基于綜合生物信息學和機器學習算法構建衰老相關分泌表型的骨關節炎預測模型

劉孝生，魏東升，2，何信用，方策

（1.遼寧中醫藥大學研究生學院，沈陽 110847；2.遼寧中醫藥大學中醫臟象理論及應用教育部重點實驗室，沈陽 110847；3.撫順市中醫院骨傷一科，遼寧 撫順 113008）

骨關節炎（osteoarthritis，OA）是一種以關節軟骨損傷、骨贅形成為表現的退行性疾病，可導致關節疼痛、畸形和功能障礙。OA的致病因素繁多，增齡相關因素被認為是導致OA的重要原因。目前，臨床上仍通過體格檢查與影像學相結合的方式來診斷OA，但OA早期階段癥狀隱匿，X線檢查往往難以發現。臨床上缺乏針對OA早期的有效篩查手段，使得OA無法得到有效的早期預防和治療［1］。因此，構建OA預測模型，對OA進行早期預測和干預尤為重要。

衰老是一種細胞損傷導致的不可逆的細胞周期永久停滯形式。衰老相關分泌表型（senescenceassociated secretory phenotype，SASP）由編碼細胞因子、趨化因子和金屬蛋白酶的基因共同組成，介導與衰老相關的生物學效應［2］。軟骨細胞維持著關節軟骨的修復、再生、代謝平衡和結構完整性［3］。許多研究已經證實了軟骨細胞衰老參與了OA的進展［4］，并且衰老的軟骨細胞數量會隨著年齡的增長而逐漸增多［5］。當過量的軟骨細胞發生衰老后，軟骨細胞細胞外基質分解與合成的代謝平衡被打破，關節軟骨遭到破壞，加速了OA的發生［6］。

機器學習算法將大量的歷史數據導入系統，通過特定的算法，識別歷史數據中的模塊特征和趨勢，通過計算機在一定精度范圍內對結局指標進行預測。近年來，隨著計算機科學與醫學領域的交叉融合，越來越多的算法被用于臨床疾病的診斷、預測和預后［7］。因此，本研究選擇機器學習算法篩選OA的預測基因。然而，單獨的算法可能在一定程度上增加數據的過擬合性，而多種機器學習算法聯合使用能在一定程度上避免這種情況的發生。因此，本研究使用3種機器學習算法來避免數據的過擬合性，提高預測的精度。

1 材料與方法

1.1 基因表達綜合（gene expression omnibus，GEO）數據庫獲取OA數據集

從GEO數據庫中獲取OA微陣列數據集GSE48556，數據集樣本來源皆為人外周血單核細胞。GSE48556中共有139個樣本，其中OA患者外周血單核細胞樣本（以下簡稱OA樣本）106個，正常人外周血單核細胞樣本（以下簡稱正常樣本）33個。從文獻［8］中收集125個SASP基因。

1.2 數據處理

在R語言4.1.3中對數據集進行背景校正、歸一化和log2轉換處理。隨后，篩選OA數據集中SASP相關基因及其表達值。

1.3 3種機器學習算法篩選候選預測標志物

使用最小絕對收縮和選擇算子（least absolute shrinkage and selection operator，LASSO）、支持向量機遞歸特征消除（support vector machine-recursive feature elimination，SVM-RFE）和隨機森林（random forest，RF）3種機器學習算法篩選OA候選預測標志物。LASSO是一種利用收縮方式的回歸分析算法，能提高其生成的統計模型的預測準確性、可解釋性和防止過擬合。SVM-RFE是一種廣泛用于微陣列數據分析的技術，可以識別具有價值的特征，從而進行結局指標的分組。RF是一種基于遞歸劃分構建二叉樹的算法，由一組決策樹組成，并考慮每個決策樹的隨機特征，可以有效地規避高維數據集小樣本的過擬合。最后，為了提高特征的準確性，將3種機器學習算法分別篩選出的候選預測標志物取交集，得到共同基因（common genes，CGs）。

1.4 CIBERSORT免疫浸潤

CIBERSORT可用于定量分析基因集中相關免疫細胞和功能。本研究采用CIBERSORT評估OA樣本與正常樣本的免疫浸潤水平。使用“pheatmap”包將OA樣本與正常樣本的免疫浸潤水平顯示在熱圖中。免疫浸潤分析中P＜0.05為差異有統計學意義。

1.5 OA預測模型構建

使用CGs構建OA預測模型，采用受試者操作特征（receiver operating characteristic curve，ROC）曲線的曲線下面積（area under curve，AUC）值評價模型的預測能力，并選取預測模型中最優基因（P＜0.001）進行動物實驗驗證。

1.6 微RNA（microRNA，miRNA）-轉錄因子（transcription factor，TF）-mRNA調控網絡預測

利用miRTarBase、Starbase和TargetScan數據庫預測CGs靶向miRNA。為了提高預測的準確性，只保留3個數據庫中共同預測的miRNA。使用Enrichr數據庫預測CGs的TF，以P＜0.05為閾值。在獲得miRNA-TF-mRNA的調控關系后，使用Cytoscape可視化miRNA-TF-mRNA的調控網絡。

1.7 動物實驗

將12只SD大鼠適應性喂養7 d后，隨機分為正常組和OA組，每組6只。OA組采用前交叉韌帶切斷法構建OA模型，2%戊巴比妥鈉0.02 mL/kg腹腔注射麻醉，麻醉后對大鼠右膝關節進行備皮、碘伏消毒。在膝關節內側做長約1 cm的切口，分離皮下筋膜層并暴露關節囊和髕韌帶。剪斷髕韌帶與肌肉連接后，將髕骨和髕韌帶向外側剝離，暴露關節間隙。剪斷前交叉韌帶，行前抽屜實驗，陽性表示已剪斷前交叉韌帶。術后每只大鼠腹腔注射青霉素（80 000 U/d），連續3 d，以預防感染。30 d后處死大鼠，腹主動脈采血，切斷腹主動脈，造成急性失血死亡。切開右膝關節，肉眼觀察關節軟骨表面退變情況，刮取關節面軟骨組織至凍凝管，液氮快速冷凍。本研究獲得遼寧中醫藥大學倫理委員會批準（21000042023054）。

1.8 實時定量PCR（real-time quantitative PCR，RTqPCR）

取正常組（n=6）和OA組（n=6）大鼠各50 mg膝關節軟骨組織，通過TRIzol法提取膝關節軟骨組織總RNA，使用SYBR法進行RT-qPCR，引物采用Primer5軟件設計。引物序列：TNFRSF1A，正向5’-CCTCCTCAGTGGGTTTCT-3’，反向5’-CGCCTTTC TATGCTTGTCC-3’；GAPDH，正向5’-TGCGACTTCA ACAGCAACTC-3’，反向5’-ATGTAGGCCATGAGGTC CAC-3’。以TNFRSF1A與GAPDH的相對比值作為其表達量，采用2-ΔΔCt法計算相對比值。為確保實驗準確性，每組每只大鼠作3個復孔。

2 結果

2.1 SASP相關OA基因

首先，對GSE48556數據集進行預處理，共得到19 613個OA相關基因。其次，將OA基因與SASP基因相結合，分離出125個與SASP相關的OA基因。

2.2 候選基因篩選結果

LASSO和SVM-RFE均使用10折交叉驗證。LASSO篩選出31個基因作為OA候選預測標志物（圖1）。SVM-RFE在特征基因數量為30個時得到最佳SVMRFE模型（圖2）。RF選取重要性排名前10位的基因作為候選預測標志物（圖3）。將3種機器學習算法分別篩選出的候選預測標志物取交集，共得到7個CGs。

圖1 LASSO模型中的特征選擇Fig.1 Feature selection in least absolute shrinkage and selection operator（LASSO）model

圖2 特征基因數量為30時得到最佳SVM-RFE模型Fig.2 The optimal support vector machines recursive feature elimination（SVM-RFE）model was obtained at 30 featured genes

圖3 RF模型圖和RF預測基因重要性排序圖Fig.3 Random forest（RF）model diagram and importance ranking diagram

2.3 CIBERSORT免疫浸潤分析

使用CIBERSORT分析正常樣本與OA樣本免疫浸潤水平差異（圖4），結果顯示，正常樣本與OA樣本間漿細胞浸潤水平存在顯著差異（P=0.001 3）。

圖4 正常樣本與OA樣本中免疫浸潤水平差異性廂圖Fig.4 Differences in immune infiltration levels between normal and osteoarthritis（OA）samples

2.4 OA預測模型構建

使用CGs構建OA預測模型，生成列線圖并得到AUC值。列線圖中，NAP1L4（P=0.005 479）、TNFRSF1A（P=0.000 875）、白細胞介素（interleukin，IL）-1β（P=0.012 166）、IL-32（P=0.150 475）、CD55（P=0.150 475）和腫瘤壞死因子（tumor necrosis factor，TNF）（P=0.616 024）高表達與OA的發生率呈正相關，CXCL8低表達（P=0.00 169）與OA的發生率呈正相關（圖5）。其中，TNFRSF1A為預測模型中最優基因。預測模型的AUC值為0.891，說明模型具有良好的預測能力。

圖5 OA預測模型列線圖Fig.5 Nomogram for osteoarthritis（OA）prediction model

2.5 miRNA-TF-mRNA調控網絡構建

利用miRTarBase、Starbase和TargetScan數據庫進行CGs-miRNA預測，結果顯示共有29個交集miRNA，其次在Enrichr數據庫中進行CGs-TF預測，結果顯示共有35個TF。通過mRNA-miRNA和mRNA-TF預測后，得到71個miRNA-TF-mRNA調控關系。通過Cytoscape建立調控網絡（圖6），包括29個miRNA，35個TF，7個mRNA。

圖6 miRNA-TF-mRNA調控網絡圖Fig.6 Diagram of miRNA-TF-mRNA regulatory network

2.6 正常組與OA組大鼠膝關節軟骨組織中TNFRSF1A的表達

通過RT-qPCR檢測大鼠膝關節軟骨組織中TNFRSF1A的表達，正常組和OA組TNFRSF1A表達水平分別為1.00±0.14和4.08±1.21，OA組TNFRSF1A表達水平明顯高于正常組（P＜0.0001），與OA預測模型分析結果一致，表明TNFRSF1A對OA具有潛在的預測價值。

3 討論

OA是一種高致殘率的疾病，早期預測對限制其致殘率至關重要。許多研究［9］發現，在OA早期，炎癥通路的激活會加速SASP的分泌，誘導軟骨細胞衰老。因此，SASP相關基因能否作為潛在的OA預測標志物，是本研究主要關注的問題。本研究使用3種機器學習算法和預測模型篩選出1個最優基因TNFRSF1A，并利用CIBERSORT探究正常樣本與OA樣本間免疫浸潤水平的差異。

隨著對OA研究的深入，許多研究開始關注免疫調控在預防和治療OA中的作用［10］。然而，對于OA相關的免疫調控仍存在許多未解之謎。因此，本研究通過CIBERSORT進行免疫浸潤分析，結果顯示，在正常樣本和OA樣本中漿細胞浸潤水平存在顯著差異。許多研究表明，關節軟骨細胞在受損或應激后，會釋放蛋白酶和膠原酶，這些降解酶在OA軟骨細胞中表達失調，其表達和活性的增加是OA發生、發展過程中軟骨降解的主要因素。關節軟骨在損傷后巨噬細胞被激活成M1型巨噬細胞并浸潤到損傷局部，分泌大量促炎性因子，如TNF-α、IL-1β、IL-6。其中IL-6作為炎癥的重要調節因子，不僅可以與其他炎性因子相互作用，加速軟骨降解與細胞外基質的破壞，更可以誘導B細胞分化為漿細胞［11］。漿細胞是一種合成與儲存抗體的細胞，可以產生抗瓜氨酸化的蛋白抗體，這種抗體會直接刺激破骨細胞，促進破骨細胞的生成，進一步加重OA［12］。

機器學習算法和列線圖結果顯示，TNFRSF1A表現出最優秀的預測精度?，F階段的研究［13］表明，炎癥通路激活是OA的發病基礎，隨著炎癥浸潤程度的加深，軟骨細胞發生不可逆的炎癥性衰老，最終發展為OA。TNFRSF1A作為誘導核因子κB（nucler factor-κB，NF-κB）通路活性的主要介質［14］，也是細胞因子TNF-α的關鍵細胞表面受體。TNFRSF1A通過結合TNF-α，誘導轉錄因子刺激NF-κB通路的激活。該通路的持續激活會導致體內炎癥水平升高，并參與軟骨退化、軟骨下硬化等OA的早期病理過程［15］。本研究通過動物實驗驗證了TNFRSF1A在OA大鼠膝關節軟骨組織中高表達。因此，TNFRSF1A極有可能成為潛在的OA預測標志物。

本研究存在一定的局限性。由于數據庫資源和樣本量的限制，預測模型的準確性可能存在偏差。此外，本研究預測相關的miRNA-TF-mRNA調控網絡未來需要深入的實驗驗證。

綜上所述，本研究基于機器學習算法構建了一個包含7個預測候選標志物的OA預測模型。機器學習算法、免疫浸潤和miRNA-TF-mRNA調控網絡，可能為研究OA提供新的方向。