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

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

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

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

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

1 材料與方法

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

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

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

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

1.3 3種機(jī)器學(xué)習(xí)算法篩選候選預(yù)測標(biāo)志物

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

1.4 CIBERSORT免疫浸潤

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

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

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

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

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

1.7 動物實(shí)驗(yàn)

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

1.8 實(shí)時定量PCR（real-time quantitative PCR，RTqPCR）

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

2 結(jié)果

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

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

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

LASSO和SVM-RFE均使用10折交叉驗(yàn)證。LASSO篩選出31個基因作為OA候選預(yù)測標(biāo)志物（圖1）。SVM-RFE在特征基因數(shù)量為30個時得到最佳SVMRFE模型（圖2）。RF選取重要性排名前10位的基因作為候選預(yù)測標(biāo)志物（圖3）。將3種機(jī)器學(xué)習(xí)算法分別篩選出的候選預(yù)測標(biāo)志物取交集，共得到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預(yù)測基因重要性排序圖Fig.3 Random forest（RF）model diagram and importance ranking diagram

2.3 CIBERSORT免疫浸潤分析

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

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

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

使用CGs構(gòu)建OA預(yù)測模型，生成列線圖并得到AUC值。列線圖中，NAP1L4（P=0.005 479）、TNFRSF1A（P=0.000 875）、白細(xì)胞介素（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）高表達(dá)與OA的發(fā)生率呈正相關(guān)，CXCL8低表達(dá)（P=0.00 169）與OA的發(fā)生率呈正相關(guān)（圖5）。其中，TNFRSF1A為預(yù)測模型中最優(yōu)基因。預(yù)測模型的AUC值為0.891，說明模型具有良好的預(yù)測能力。

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

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

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

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

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

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

3 討論

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

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

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

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

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