0 前言

       中轴型脊椎关节炎（axial spondyloarthritis, axSpA）是一类以腰背痛及骶髂关节和脊柱的慢性炎症为代表症状的疾病^[1]。axSpA初期症状隐匿且临床表现缺乏特异性，经常被延误诊断^[2]。随着axSpA病程的进展，其临床表现进行性加重，未及时治疗可导致不可逆的中轴骨畸形，致使活动受限^[3]。因此，早期诊断及早期治疗以改善axSpA的预后意义重大。

       影像技术在辅助诊断及鉴别诊断axSpA的过程中不可或缺。在《影像学技术在脊柱关节炎中应用的中国专家共识（2021年版）》^[4]中，骶髂关节MRI被认为是检测放射学阴性（即X线平片阴性）的疑诊患者和axSpA活动性病变的首选方法。骶髂关节MRI在axSpA的早期诊断中起到了积极的作用。然而，近年越来越多的研究却发现骨髓水肿缺乏特异性，可能存在于产后女性、入伍新兵乃至正常人群的骶髂关节MRI中^{[5, 6]}。这可能会导致假阳性表现，进而导致axSpA的过度诊断。另外，不同放射科医生对骶髂关节MRI图像的解读可能不同，这种观察者之间的变异性可能会影响诊断的准确性^[7]。2021年，JIBRI等^[8]对2019年更新的年国际脊柱关节炎评估协会（Assessment of SpondyloArthritis International Society, ASAS）axSpA骶髂关节阳性病变分类标准的阅读者间一致性进行了研究。该研究中，2名风湿领域影像学专家对94例活动性axSpA患者的骶髂关节MRI图像进行判读，发现对骶髂关节活动性炎症的判读呈中等一致性，骨侵蚀的一致性较差。采用新技术对骶髂关节MRI进行分析，可能会减少视觉解读带来的主观性，并允许对成像数据进行标准化评估，进而提高诊断准确性。

       影像组学是指从医学影像中对特定感兴趣区（region of interest, ROI）的定量特征进行高通量计算提取，将医学数字图像转换为深层数据进行定量分析，最终可用于疾病的辅助诊断和分类分级^[9]。利用该技术，可从MRI中提取众多细微的定量特征，建立提高临床诊断准确度的模型。影像组学在肿瘤学领域已有较丰富的研究^[10]，但由于包括axSpA在内的肌骨疾病影像表现多样，病变界限模糊，目前仅有少数研究使用影像组学方法分析骶髂关节MRI^{[11, 12]}，且这些研究多集中于对影像特征的区分，尚未建立有效的能够辅助临床诊断的模型。本研究探讨利用临床数据及骶髂关节MRI影像建立用于诊断axSpA的多模态融合影像组学模型，以期辅助提高临床医师诊断的效率。

1 材料与方法

1.1 研究对象

       本研究遵守《赫尔辛基宣言》，经北京大学第三医院伦理委员会批准，免除受试者知情同意，伦理批准号：M2022399，国家临床试验编号：MR-11-22-009236。

       回顾性分析于2019年4月至2021年9月间在北京大学第三医院接受骶髂关节MRI检查的患者临床及影像资料。纳入标准：（1）年满18岁；（2）因怀疑axSpA（患有下背痛）就诊；（3）无腰椎及骨盆手术史；（4）1年内无骶髂关节外伤史；（5）无骶髂关节感染史；（6）不合并其他结缔组织病史；（7）无恶性肿瘤病史。排除标准：（1）没有明确的临床诊断；（2）骶髂关节MRI图像质量差（如运动伪影重）、覆盖范围不完整或缺少所需序列；（3）临床资料不完整。结合医院电子病历系统中临床医师的记录，将患者分为确诊患者及健康对照组。所有axSpA确诊患者均符合2019年ASAS工作组的诊断标准。健康对照组与axSpA患者1∶1配对，为同期接受骶髂关节MRI检查，符合上述纳入及排除标准，因怀疑axSpA就诊而后经临床医师明确排除axSpA，且未患有其他疾病的患者，使两组年龄和性别一致。

1.2 临床数据

       从医院电子病历系统中收集患者的病历及临床数据。我们收集了所有患者和健康对照组的性别、年龄、红细胞沉降率（erythrocyte sedimentation rate, ESR）和C反应蛋白（C-reactive protein, CRP）结果。

1.3 影像检查方法

       所有的MRI图像从3台3.0 T磁共振仪中获取，2台型号为GE Discovery 750w（GE医疗，美国），1台型号为Discovery 750（GE医疗，美国），3台磁共振仪的扫描参数一致。纳入研究的扫描序列有斜冠状位（平行于骶骨长轴）快速自旋回波T1WI和斜冠状位快速自旋回波脂肪抑制（fat-supressed, FS）T2WI。T1WI序列扫描参数：重复时间450 ms，回波时间85 ms，层厚4 mm，层间隔0.5 mm，视野30 cm×30 cm，矩阵320×256，激发次数1；FS-T2WI序列扫描参数：重复时间3200 ms，回波时间85 ms，层厚4 mm，层间隔0.5 mm，视野30 cm×30 cm，矩阵320×256，激发次数4。

1.4 图像分割与特征提取

       将所有入组者的骶髂关节MRI斜冠状位T1WI及FS-T2WI导入ITK-SNAP软件（版本3.8.0，www.itksnap.org）^[13]。在高年资医师的指导下，由2名分别具有2年和5年经验的放射科医生，在6个连续骶髂关节中心层面上，沿着每侧骶髂关节间隙的10 mm范围勾画ROI，如图1所示。

       将每个层面的ROI合并为1个ROI之后，使用PyRadiomics平台（版本3.0.1）^[14]从每种序列勾画的ROI中提取了相应的组学特征。为了提取更多的特征，将多种滤波器应用于ROI图像，再从中提取特征。提取的特征遵循成像生物标志物标准化倡议定义^[15]，可分为几何形状特征、强度特征和纹理特征三组：几何特征捕捉了ROI的三维形状特征；强度特征代表ROI内体素强度的一阶统计分布；纹理特征描述强度的图案或二阶及高阶空间分布。

点击查看大图

图1  T1WI（1A）及FS-T2WI（1B）图像上勾画受试者工作特征（ROI）示意图。
Fig. 1  Regions of interest (ROI) segmented on T1WI (1A) and FS-T2WI (1B) images.

1.5 特征选择

       为了评估图像分割对组学特征稳定性的影响，利用组内相关系数（intra-class correlation coefficient, ICC）进行了同一评估者及2个不同评估者的重复性分析。同一评估者间隔2周对随机选择的30例患者图像进行重复两次分割，而2名评估者则对另外随机选择的30例患者图像进行独立分割。ICC被用于评估从这些ROI中提取的组学特征，其中ICC≥0.85被认为是不受分割影响的稳定特征。经使用ICC进行初步筛选后，使用Z分数方法对所有特征进行标准化，以确保特征的数量级一致。随后，使用t检验计算所有组学特征的P值。仅保留P＜0.05的影像组学特征用于进一步分析。

       为了识别高度相关的特征，使用皮尔逊相关系数对高度可重复特征进行进一步分析。为了避免冗余，当任意两种特征之间的相关系数超过0.9时只保留一个特征。为了在最大化特征具有代表性的同时使冗余最小化，采用贪心递归式消除策略过滤特征，每次迭代时删除当前集中冗余度最高的特征。

       最终，使用最小绝对收缩和选择算子（least absolute shrinkage and selection operator, LASSO）回归模型筛选构建影像组学标签的特征。根据正则化权重λ，LASSO将许多相关性低的特征系数设置为0。使用最小标准10倍交叉验证确定最佳λ值，选择具有最低平均标准误的λ。

1.6 列线图模型构建

       对临床特征，我们先后进行了单变量和多变量分析。选择多变量分析中差异具有统计学意义的特征用于构建列线图。对单变量分析中差异具有统计学意义的临床特征，使用逻辑回归（logistic regression, LR）、随机森林（RandomForest）、极端随机树（ExtraTrees）、极端梯度提升（XGBoost）、多变量逻辑回归（multivariate logistic regression）5种方式建立机器学习模型。建立这些模型时，采用了5倍交叉验证，并利用了网格搜索进行超参数调整。选择性能最好的模型为临床标签。

       对经LASSO筛选后的影像组学特征，利用与构建临床标签相同的机器学习模型，选择其中性能最好的作为影像组学标签，并用于构建列线图。本研究比较了T1WI和FS-T2WI两种序列的性能。为了评估这两种序列的融合是否优于单一序列，将两种序列的特征组合以获得融合特征集。该融合特征集特征选择和建立机器学习模型的过程与单一序列特征集相同。

       最终，整合选择的临床特征和影像组学标签以建立列线图模型。

1.7 统计学分析

       采用SPSS 24.0（IBM SPSS Inc., USA）软件进行统计学分析。采用Kolmogorov-Smirnov方法检验定量变量临床数据的正态性，符合正态分布的数据表示为平均值±标准差（x¯±s）并采用独立样本t检验，偏态分布的定量数据则用中位数（四分位数间距）[M (Q1, Q3)]表示，使用Mann-Whitney U检验，离散变量则使用卡方检验。P＜0.05为差异有统计学意义。利用准确性、受试者工作特征（receiver operating characteristic, ROC）曲线下面积（area under the curve, AUC）及其95%置信区间（confidence interval, CI）、敏感度、特异度、真阳性率、假阳性率以及F1分数在训练及测试集中评估各个机器学习模型两种标签及列线图的诊断性能。生成校准曲线来评估列线图模型的校准性能，利用决策曲线分析确定列线图模型的临床效用。

2 结果

2.1 临床特征

       本研究最终建立了一个纳入204例患者的研究队列，其中axSpA确诊患者102例，健康对照组102例。利用计算机随机数字生成程序进行简单随机化，将204例患者按照7∶3的比例分为训练集和测试集。具体统计分析结果见表1。

       对所有临床特征进行单变量分析，对单变量分析中P＜0.05的变量进行多变量分析，并获得每个特征的比值比和相应的P值，如表2所示。其中，选择多变量分析中的ESR（P＜0.05）来构建列线图。图2中显示了每个临床特征之间的相关性，表明CRP、ESR与是否确诊axSpA之间存在明显的线性关系。

       我们提取了单变量分析中P＜0.05的临床特征，包括CRP和ESR，并构建了多个机器学习模型。将测试集中AUC最高的LR模型（AUC＝0.884,　95% CI: 0.804-0.964）选为临床标签（表3、图3）。

点击查看大图

图2  C反应蛋白（CRP）、红细胞沉降率（ESR）与是否确诊中轴型脊柱关节炎的Spearman相关系数。
图3  训练集（3A）和测试集（3B）中基于临床特征的不同机器学习模型的受试者工作特征曲线。AUC：曲线下面积；CI：置信区间；LR：逻辑回归；RandomForest：随机森林；ExtraTrees：极端随机树；XGBoost：极端梯度提升；MLP：多变量逻辑回归。
Fig. 2  Spearman coefficients of C-reactive protein (CRP), erythrocyte sedimentation rate (ESR) and diagnosis of axial spondyloarthritis.
Fig. 3  The receiver operating characteristic (ROC) curves for different machine learning models based on clinical characteristics in the training set (3A) and testing set (3B). AUC: area under the curve; CI: confidence interval; LR: logistic regression; RandomForest: random forest ExtraTrees: Extremely randomized trees; XGBoost: eXtreme Gradient Boosting; MLP: multivariate logistic regression.

点击查看表格

表1  训练集和测试集中患者的基线特征
Tab. 1  Basic characteristics of the patients in the training set and testing set

点击查看表格

表2  临床特征的单变量和多变量分析
Tab. 2  Uni-variable and multi-variable analysis of clinical characteristics

点击查看表格

表3  基于临床特征的机器学习模型在训练集和测试集中的效能
Tab. 3  Performance of the machine learning models based on clinical characteristics in the training and testing set

2.2 组学特征

       从T1WI和FS-T2WI图像中提取了共计1834种影像组学特征。包括：360种一阶特征，14种形状特征，以及1460种纹理特征。为了增强特征集，我们合并了2种不同时间模态的特征，最终得到3368种融合特征。利用LASSO回归模型，从这些特征中分别筛选了36种融合特征、23种T1WI特征和17种FS-T2WI特征。所筛选的特征及λ系数如图4所示。

       对于所筛选的T1WI、FS-T2WI以及融合影像组学特征，利用了多种机器学习模型构建预测模型（表4）。其中利用T1W1和融合组学特征训练的机器学习模型中表现最佳的为LR，其测试集AUC分别为0.778（95% CI: 0.661-0.895）和0.905（95% CI: 0.829-0.981）；利用FS-T2WI特征训练的模型中表现最佳的则为XGBoost，其测试集AUC为0.907（95% CI: 0.832-0.981）。基于筛选的融合影像组学特征的LR模型诊断性能最佳，过度拟合低且AUC高（图5），这表明融合组学特征优于单序列特征。最终，将基于所筛选的融合影像组学特征的LR模型选为最终的影像组学标签（Rad_Sig）。

点击查看大图

图4  筛选的融合特征（4A），T1WI特征（4B）和脂肪抑制（FS）T2WI特征（4C）及其λ系数。
Fig. 4  Selected fusion features (4A), T1WI features (4B) and fat-suppressed (FS) T2WI features (4C), and λ coefficients.

点击查看大图

图5  训练集（5A）和测试集（5B）中基于T1加权、脂肪抑制T2加权以及融合影像组学特征的LR模型的受试者工作特征曲线。AUC：曲线下面积；CI：置信区间。
Fig. 5  The receiver operating characteristic (ROC) curves of LR model based on T1-weighted, fat-suppressed T2-weighted and fused radiomics features in the training set (5A) and testing (5B) set. AUC: area under the curve; CI: confidence interval.

点击查看表格

表4  基于所筛选的T1WI、脂肪抑制T2WI以及融合影像组学特征的不同机器学习模型的axSpA诊断性能
Tab. 4  The diagnostic performance in axSpA of different machine learning models based on selected T1WI, fat-suppressed T2WI and fused radiomics features

2.3 不同标签及列线图效能评价

       基于ESR和前述所得的Rad_Sig，我们构建了可用于辅助诊断的列线图模型，如图6所示。最终所得的列线图模型在训练集（AUC＝0.997, 95% CI: 0.992-1.000）和测试集（AUC＝0.944, 95% CI: 0.889-1.000）中表现出了优良的诊断性能。表5及图7展示了列线图模型同临床标签及影像组学标签的效能比较。图8所示校准曲线显示列线图模型的校准性能良好，决策曲线显示列线图模型具有更好的预测性能和临床应用价值。表6所示DeLong检验发现，训练集中临床标签与列线图及融合影像组学标签诊断axSpA的效能存在显著差异（P＜0.05）（表6）。

点击查看大图

图6  临床使用的列线图。ESR：红细胞沉降率；Rad_Sig：融合影像组学标签。
Fig. 6  Nomogram for clinical use. ESR: erythrocyte sedimentation rate; Rad_Sig: fusion radiomics label.

点击查看大图

图7  训练集（7A）及测试集（7B）中临床标签（Clinic_Sig）、融合影像组学标签（Rad_Sig）及列线图模型（nomogram）的受试者工作特征曲线。AUC：曲线下面积；CI：置信区间。
Fig. 7  The receiver operating characteristic (ROC) curves of clinical label (Clinic_Sig), fusion radiomics label (Rad_Sig) and nomogram model in training set (7A) and testing set (7B). AUC: area under the curve; CI: confidence interval.

点击查看大图

图8  测试集中临床标签（Clinic_Sig）、融合影像组学标签（Rad_Sig）及列线图模型（nomogram）的校准曲线（8A）和决策曲线（8B）。
Fig. 8  Calibration curve (8A) and decision curve (8B) of clinical label (Clinic_Sig), fusion radiomics label (Rad_Sig), and nomogram model in the testing set.

点击查看表格

表5  临床标签、融合影像组学标签及列线图模型诊断axSpA的指标
Tab. 5  Indicators for the diagnosis of axSpA using clinical label, fusion radiomics label and nomogram model

点击查看表格

表6  临床标签、融合影像组学标签及列线图模型之间的DeLong检验结果数值
Tab. 6  Results of DeLong test between clinical label, fusion radiomics label and nomogram model

3 讨论

       本研究中，我们结合广泛用于辅助诊断的临床特征及基于影像组学的融合特征标签，成功地建立了能用于协助诊断axSpA患者的列线图模型。该模型在训练集和测试集中的AUC分别为0.997（95% CI: 0.992-1.000）和0.944（95% CI: 0.889-1.000）。2019年，JONES等^[16]总结了近年来国内外将骶髂关节MRI用于诊断axSpA的研究发现，这些研究以骨髓水肿为主要标准，MRI的敏感度在0.35至0.91之间，特异度在0.42至0.97之间。本研究所得的影像组学标签及列线图模型的诊断效能表现优良，优于上述基于判读骶髂关节MRI的研究，也优于目前国内外其他类似利用影像组学技术分析骶髂关节MRI的研究。这表明影像组学技术及列线图模型优于人工判读影像，为更早的诊断axSpA带来了希望。长久以来，早期且准确诊断axSpA是一项艰巨的挑战，患者常常在首诊数年后方能确诊^[17]。本研究为进一步利用影像组学技术，建立能够协助临床医师判读骶髂关节MRI以及更早诊断axSpA的工具打下了基础。

3.1 临床特征

       ESR和CRP是少数在临床用于诊断axSpA以及预测其疾病进展的生物标志物^[18]。本研究对临床因素的多变量回归分析表明，ESR水平是诊断axSpA的独立预测因素，这与之前的研究一致^{[19, 20]}。尽管CRP和ESR在单变量分析中都具有显著性，CRP却并不是诊断的独立危险因素。非特异性腰痛患者和肥胖患者中的CRP水平也会升高，这可能导致其在本研究中对于axSpA诊断价值的下降^[21]。人类白细胞抗原B27（Human leukocyte antigen B27, HLA-B27）蛋白也是临床用于辅助诊断axSpA的生物标志物，其与axSpA的关联是迄今已知的任何人类疾病中最强的，大约85%的axSpA患者HLA-B27为阳性。但HLA-B27的特异性较低，阳性个体一生中患axSpA的概率仅为1%-2%^[22]。本研究受临床资料限制，未能完整收集入组患者的HLA-B27检查结果，因此未被纳入本研究中加以分析。此外，由于本研究入组的患者与健康对照组两组患者年龄和性别一致，因此本研究中患者的年龄与性别不是独立危险因素。

3.2 MRI影像特征

       2019年ASAS的最新定义认为，骶髂关节MRI显示的骨髓水肿是axSpA骶髂关节炎主要的诊断条件^[23]。另外，如关节囊炎、骨侵蚀部位炎症、附着点炎和关节间隙积液等骶髂关节的急性炎症改变，以及如骨侵蚀、脂肪沉积、骨侵蚀腔内脂肪化生、骨质硬化、骨性强直和骨芽等结构性变化也可辅助诊断^[24]。国内外指南推荐且临床常用的骶髂关节MRI扫描序列包括FS-T2WI、短时反转恢复（short tau inversion recovery, STIR）以及T1WI等^{[25, 26]}。FS-T2WI能够清晰显示关节面下骨髓水肿、囊变及关节积液等信号变化。STIR对炎症活动比较敏感，可检测axSpA的活动性，但其特异度欠佳^[27]。2019年，GREESE等^[28]的研究显示FS-T2WI检测骨炎的效能比STIR更好。T1WI主要显示关节软骨和软骨下区的形态结构，但不易分辨骨皮质、骨髓水肿及关节积液。本研究中建立的影像组学模型之所以能有良好的诊断性能，正是因为构建模型所用的影像组学特征是从上述MRI序列中包含上述急性炎症改变以及结构性变化的骶髂关节ROI中所提取的。

       本研究结果中，基于FS-T2WI建立的影像组学特征机器学习模型诊断axSpA的效能优于T1WI，而结合了两种序列的融合特征模型的效能又优于FS-T2WI和T1WI。这些结果印证了之前的相关研究。先前研究显示，T1WI主要显示的脂肪沉积和结构性病变对axSpA的诊断效能低于FS-T2WI显示的骨髓水肿^[29]。WEBER等^[30]研究发现，骶髂关节骨髓水肿诊断axSpA的敏感度为67%，而骨髓水肿联合骨质侵蚀，在保持良好特异度的同时，诊断的敏感度能提高到81%。

3.3 既往影像组学及列线图模型研究

       利用蕴藏在骶髂关节MRI中的凭人眼无法识别的大量的深层次数据，影像组学技术能挖掘骶髂关节MRI辅助诊断axSpA的潜能^[31]。到目前为止，国内外仅有少数研究团队探索了影像组学在axSpA诊断中的应用。TENÓRIO等^[32]使用47例axSpA患者的队列，将从骶髂关节MRI的频率衰减反转恢复和STIR序列中提取的影像组学特征与各种临床指标相关联，显示组学特征具有评估SpA骶髂关节病变的潜力，且构建的影像组学模型支持SpA临床评估^[13]。利用90例患者的队列，KEPP等^[33]揭示在区分骶髂关节炎和退行性变化方面，基于影像组学的纹理分析优于定性评估。然而，上述研究与本研究不同，本研究结合临床特征建立了能够辅助临床诊断的模型。基于638例患者（424例为axSpA，214例为非axSpA），YE等^[20]建立了一个结合影像组学特征与临床指标相的列线图模型。该模型在训练集和测试集中的AUC分别为0.90（95% CI: 0.86-0.93）和0.90（95% CI: 0.85-0.94）。YE等仅使用由影像组学标签中各个组学特征及其相应LASSO系数计算得出的Rad_Score构建列线图模型，而本研究在构建临床模型和组学模型时尝试了多种机器学习算法，由这些算法所得的最佳模型构建的列线图能显示出更好的性能。因此，本研究所得的列线图模型在训练集和测试集中的AUC均优于YE等的列线图模型，能够更进一步提高辅助诊断axSpA的效率。

3.4 局限性及展望

       本研究的局限性包括：（1）本研究为回顾性队列研究，尽管已选取性别和年龄匹配的对照组，但可能存在选择偏倚，需要前瞻性研究以进一步证实结果；（2）本研究样本量相对较少且来源于单一医疗机构，有待进行大样本、多个医疗中心的研究，以验证研究结果的可靠性；（3）本研究只采用了常规扫描序列提取的组学特征，从动态对比增强、扩散加权成像、Dixon方法及三维高分辨率MR等不常用于axSpA筛查的序列图像中^{[34, 35]}，可提取更多种类的影像组学特征，以构建效能更高的模型；（4）本研究使用手工勾画ROI既复杂又耗时，引入类似于KUCYBAŁA等^[36]使用的ROI自动分割技术，可提高ROI勾画的可靠性和再现性。

4 结论

       综上所述，骶髂关节MRI影像组学特征与临床参数结合建立的联合列线图模型能够有效区分axSpA患者和健康对照组，影像组学技术在辅助临床诊断axSpA方面具有重大潜能以及广泛的临床应用前景。