基于綜合生物信息學和機器學習算法構(gòu)建衰老相關(guān)分泌表型的骨關(guān)節(jié)炎預測模型

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

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

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

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

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

1 材料與方法

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

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

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

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

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

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

1.4 CIBERSORT免疫浸潤

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

1.5 OA預測模型構(gòu)建

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

1.6 微RNA（microRNA，miRNA）-轉(zhuǎn)錄因子（transcription factor，TF）-mRNA調(diào)控網(wǎng)絡預測

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

1.7 動物實驗

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

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

取正常組（n=6）和OA組（n=6）大鼠各50 mg膝關(guān)節(jié)軟骨組織，通過TRIzol法提取膝關(guān)節(jié)軟骨組織總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 結(jié)果

2.1 SASP相關(guān)OA基因

首先，對GSE48556數(shù)據(jù)集進行預處理，共得到19 613個OA相關(guān)基因。其次，將OA基因與SASP基因相結(jié)合，分離出125個與SASP相關(guān)的OA基因。

2.2 候選基因篩選結(jié)果

LASSO和SVM-RFE均使用10折交叉驗證。LASSO篩選出31個基因作為OA候選預測標志物（圖1）。SVM-RFE在特征基因數(shù)量為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 特征基因數(shù)量為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），結(jié)果顯示，正常樣本與OA樣本間漿細胞浸潤水平存在顯著差異（P=0.001 3）。

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

2.4 OA預測模型構(gòu)建

使用CGs構(gòu)建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的發(fā)生率呈正相關(guān)，CXCL8低表達（P=0.00 169）與OA的發(fā)生率呈正相關(guān)（圖5）。其中，TNFRSF1A為預測模型中最優(yōu)基因。預測模型的AUC值為0.891，說明模型具有良好的預測能力。

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

2.5 miRNA-TF-mRNA調(diào)控網(wǎng)絡構(gòu)建

利用miRTarBase、Starbase和TargetScan數(shù)據(jù)庫進行CGs-miRNA預測，結(jié)果顯示共有29個交集miRNA，其次在Enrichr數(shù)據(jù)庫中進行CGs-TF預測，結(jié)果顯示共有35個TF。通過mRNA-miRNA和mRNA-TF預測后，得到71個miRNA-TF-mRNA調(diào)控關(guān)系。通過Cytoscape建立調(diào)控網(wǎng)絡（圖6），包括29個miRNA，35個TF，7個mRNA。

圖6 miRNA-TF-mRNA調(diào)控網(wǎng)絡圖Fig.6 Diagram of miRNA-TF-mRNA regulatory network

2.6 正常組與OA組大鼠膝關(guān)節(jié)軟骨組織中TNFRSF1A的表達

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

3 討論

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

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

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

本研究存在一定的局限性。由于數(shù)據(jù)庫資源和樣本量的限制，預測模型的準確性可能存在偏差。此外，本研究預測相關(guān)的miRNA-TF-mRNA調(diào)控網(wǎng)絡未來需要深入的實驗驗證。

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