0 前言

       骨质疏松症是一种与年龄有关的以骨量减少、骨质量受损及骨强度减低，导致骨脆性增加，易发生骨折为特征的全身性疾病。研究表明，2016年中国60岁以上的老年人骨质疏松症患病率为36%。据预测，至2050年，我国骨质疏松骨折患病人数将达599万，相应的医疗支出高达1745亿元^[1]，造成严重的社会医疗负担，威胁人民健康。《中国老年骨质疏松症诊疗指南（2018）》^[2]明确指出了双能X线吸收法（dual X-ray absorptiometry, DXA）是国际和国内公认测量骨质密度的“金标准”，对临床诊断骨质疏松症具有重大意义，但因扫描仪普及性不够、辐射暴露等原因使其远不能满足临床需求^{[3, 4, 5]}，该指南也指出了MRI具有足够的空间分辨率和高对比度分辨率、无创伤、无电离辐射等优点，在诊断腰椎疾病方面具有其独特的优势，并且在临床使用广泛^{[6, 7, 8]}。目前已有MR波谱成像、水-脂分离成像、高分辨率成像、动态对比增强成像、扩散加权成像，超短回波时间成像用于评估骨质疏松的研究^{[9, 10, 11]}，并且随着未来MR设备和技术的进一步完善与优化，有望提供骨质疏松诊断或筛查的新策略，但这些方法都需要在常规MRI检查基础上额外增加扫描序列，延长扫描时间，不适用于因疼痛就诊、依从性较低的患者。

       影像组学是由荷兰学者Lambin于2012年提出的一种新型的定量分析技术，它使用计算机算法从医学影像图像中提取高通量的定量、特征性信息并分析建模^[12]，且应用先进的数据模型算法来检测视觉无法识别的放射组学特征，有助于早期无创地明确病灶的性质、评估疾病的疗效、判断预后，以及实施精准医疗、个体化医疗等^{[13, 14, 15]}。

       国外亦有基于计算机断层扫描（computed tomography, CT）影像组学分析检测骨质疏松症的研究^{[16, 17]}，其主要临床应用局限性表现在辐射暴露性，T2WI序列作为腰椎MRI检查最基本的序列，易操作且图像获取度高，并且临床以腰背部疼痛为主诉就诊行MR检查后最终诊断为骨质疏松症的患者日益增多，但基于腰椎T2WI序列影像组学诊断骨质疏松症的研究鲜有报道。因此，本研究旨在探讨基于T2WI序列的腰椎MRI放射组学分析在检测骨质疏松症中的潜力。

1 材料与方法

1.1 研究对象

       本研究遵守《赫尔辛基宣言》，经孝感市中心医院伦理委员会批准，免除受试者知情同意，批准文号：XGLY2022-11-07。回顾性分析了2022年12月至2023年3月期间因下腰部疼痛于我院行MRI检查的患者病例共300例，排除9例，最终291例入组。纳入标准：（1）14天内同时进行了腰椎MRI和DXA检查；（2）患者年龄区间为45-83岁。排除标准：（1）图像质量差影响成像评价；（2）有恶性肿瘤或血液疾病、糖尿病、甲状腺相关疾病等病史。依据DXA结果将患者分为骨质疏松组、非骨质疏松组（图1）。记录所有患者的临床资料，包括年龄、性别、身体质量指数（body mass index, BMI）。

点击查看大图

图1  研究人群选择流程图。
Fig. 1  Flow chart of study population selection.

1.2 影像检查方法

       采用SIEMENS 3.0 T超导型MR扫描仪（MAGNETOM Verio, Germany），用脊柱相控阵体线圈行腰1至腰5椎体T2WI序列图像采集，扫描参数：TR 2400 ms，TE 100 ms，层厚3 mm，层间距0.3 mm，矩阵384×268，视野320 mm×320 mm。

1.3 影像组学过程

1.3.1 图像获取

       所有数据来源于科室的图像存档和通信系统（picture archiving and communication system, PACS），以匿名DICOM格式进行下载。

1.3.2 感兴趣区的勾画、分割及特征提取

       将腰椎矢状位T2WI图像导入开源软件3D Slicer（V5.0.3,https://www.slicer.org），由一名具有3年工作经验的影像科主治医师逐层手动勾画腰1至腰5椎体感兴趣区域（region of interest, ROI），含皮质骨与松质骨，注意避开椎间盘、终板，共15层。利用3D Slicer软件中的“PyRadiomics”提取常见的7个特征组：一阶特征（First order features）、形状特征（shape features）、灰度依赖矩阵（gray level dependence matrix, GLDM）、灰度共生矩阵（gray-level co-occurrence matrix, GLCM）、灰度运行长度矩阵（gray-level run length matrix, GLRLM）、灰度大小区域矩阵（gray level size zone matrix, GLSZM）和邻域灰度色差矩阵（neighborhood gray-tone difference matrix, NGTDM）。由另一名5年工作经验的影像科主治医师随机抽取大约1/3的图像独立勾画ROI并提取相应特征，采用组内相关系数（intra-class correlation coefficient, ICC）比较两位医师所勾画的特征数据，当ICC≥0.75时，则认为该影像组学特征具有一定的可靠性。

1.3.3 特征筛选

       我们对所有ICC≥0.75的放射学特征进行了Mann-Whitney U统计检验和特征筛选。特征筛选基于两步：（1）引入并使用特征选择工具类（Feature Selector）进行初步特征选择；（2）对于重复性高的特征，还使用Spearman等级相关系数计算特征参数间的冗余性，选取0.9为冗余阈值，并使用最小绝对收缩和选择算子（least absolute shrinkage and selection operator, LASSO）算法根据最优调节权重λ进一步减少冗余。保留的非零系数特征用于回归模型拟合，并组合成放射组学特征。

1.3.4 机器学习算法选择

       本研究共使用了八种传统的机器学习算法：逻辑回归（logistics regression, LR）、支持向量机（support vector machine, SVM）、K邻近算法（K-nearestbor neighbor, KNN）、随机森林（RandomForest）、极度随机树（extremely randomized trees, ExtraTrees）、极度梯度提升树（eXtreme gradient boosting, XGBoost）、梯度提升决策树（light gradient boosting machine, LightGBM）、多层感知机（multi-layer perception, MLP）。采用受试者工作特征（receiver operating characteristic, ROC）曲线下面积（area under the curve, AUC）、准确度、敏感度及特异度对每个机器学习方法预测性能进行评估，并进行5折交叉验证筛选最佳机器学习方法应用于后续模型构建。

1.3.5 诊断模型的构建

       （1）临床诊断模型：对骨质疏松症与非骨质疏松症组间患者年龄、性别、BMI进行比较，选取具有统计学意义（P＜0.05）的临床指标建立临床诊断模型；（2）影像组学模型：应用特征筛选保留的影像组学特征建立影像组学模型；（3）联合模型：将筛选出的临床特征和影像组学特征联合，建立联合模型。

1.3.6 模型的评估

       通过各个模型的特异度、敏感度、准确度及AUC等参数分析模型的预测价值，并对效能最佳的模型使用校正曲线、Hosmer-Lemeshow检验验证模型的拟合情况，再利用DeLong检测对各模型之间的AUC进行统计学差异性比较，最后采用决策曲线分析（decision curve analysis, DCA）评估各个模型的临床净收益（图2）。

点击查看大图

图2  影像组学流程图。
Fig. 2  Flowchart of radiomics.

1.4 统计学分析

       采用IBM SPSS 20.0（version 20.0, USA）进行统计分析。采用Kolmogorov-Smirnov法对计量资料进行正态性检验，符合正态分布的计量资料采用独立样本t检验，不符合正态分布的计量资料采用Mann-Whitney U检验。计数资料以例数或构成比表示，两组间比较采用χ²检验。采用R软件（version 4.2.1,https://www.r-project.org）筛选纹理特征，首先采用ICC选择观察者间ICC≥0.75的特征，然后采用Spearman相关性检验剔除特征间r＞0.9的特征，最后采用LASSO进行特征降维。绘制八种机器学习方法的ROC曲线并获得AUC，选取最佳机器学习方法联合临床特征、影像组学特征建立模型，对各类模型绘制ROC曲线评估诊断效能，计算AUC、准确率、敏感度及特异度。利用Hosmer-Lemeshow检验、校正曲线验证模型的拟合情况。使用DCA评价模型的临床获益度。使用DeLong检验比较AUC的诊断差异。以P＜0.05为差异具有统计学意义。

2 结果

2.1 一般资料

       本研究共纳入291例受试者病例，其中男97例（33.3%），女194例（66.7%），年龄为43-85（62±10）岁。年龄、BMI在训练组及测试组的骨质疏松与非骨质疏松组间具有统计学意义（P＜0.05）。随后，采用8∶2的比例将所有病例随机分配到训练组（n=233）和测试组（n=58）（表1）。

点击查看表格

表1  入组研究患者的基线特征
Tab. 1  The baseline characteristics of enrolled patients

2.2 影像组学特征选择结果

       训练组共233例患者1165个椎体，骨质疏松和非骨质疏松患者分别为91例455个椎体及142例710个椎体；测试组共58例患者290个椎体，骨质疏松和非骨质疏松患者分别为23例115个椎体及35例175个椎体。本研究最终选取23个非零系数的最佳放射组学特征，分别为一阶特征5个、GLCM特征5个、GLRLM特征5个、GLSZM特征3个和NGTDM特征5个（图3）。

点击查看大图

图3  影像组学特征筛选。横坐标表示权重，纵坐标表示特征。
Fig. 3  Feature selection for radiomics. The abscissa represents the coefficients, and the ordinate represents the feature.

2.3 分类诊断模型在测试集的效能

       基于测试集数据及入选特征，各分类诊断模型的AUC从高到低排列如下：LR（0.913）、MLP（0.884）、SVM（0.872）、LightGBM（0.795）、XGBoost（0.787）、ExtraTrees（0.773）、KNN（0.772）、RandomForest（0.769）。各模型的分类诊断模型ROC曲线详见表2、图4。

点击查看大图

图4  各分类诊断模型鉴别骨质疏松与非骨质疏松的受试者工作特征（ROC）曲线。LR：逻辑回归；AUC：ROC曲线下面积；SVM：支持向量机；KNN：K邻近算法；RandomForest：随机森林；ExtraTrees：极度随机树；XGBoost：极度梯度提升树；LightGBM：梯度提升决策树；MLP：多层感知机。
Fig. 4  The receiver operating characteristic (ROC) curves for each classification diagnostic model to discriminate osteoporosis from non-osteoporosis. LR: logistics regression; AUC: area under the curve; SVM: support vector machine; KNN: K-nearestbor neighbor; RandomForest: random forest; ExtraTrees: extremely randomized trees; XGBoost: eXtreme gradient boosting; LightGBM: light gradient boosting machine; MLP: multi-layer perception.

点击查看表格

表2  各分类诊断模型鉴别骨质疏松与非骨质疏松的效能
Tab. 2  The effectiveness of each classification diagnostic model in differentiating osteoporosis from non-osteoporosis

2.4 诊断模型构建及应用效能

       基于LR的分类模型在鉴别骨质疏松与非骨质疏松方面表现最为优异。（1）临床诊断模型：临床模型由年龄、BMI基于LR分类模型构建而成；（2）影像组学诊断模型：提取的23个非零系数的最佳放射组学特征基于LR分类模型构建而成；（3）联合模型：该模型由临床特征及影像组学特征基于LR分类模型联合构建而成。临床模型预测训练集和测试集判定是否有骨质疏松的AUC值分别为0.791（95% CI: 0.733-0.849）及0.805（95% CI: 0.676-0.935）；影像组学模型预测训练集和测试集判定是否有骨质疏松的AUC值分别为0.879（95% CI: 0.833-0.925）及0.913（95% CI: 0.841-0.985）；基于临床模型及放射组学模型构建的联合模型，其预测训练集和测试集判定是否有骨质疏松的AUC值分别为0.893（95% CI: 0.853-0.934）及0.904（95% CI: 0.825-0.984）（表3、图5）。

       绘制临床模型、影像组学模型及联合模型的校正曲线，Hosmer-Lemeshow检验（P=0.250, 0.753, 0.575）验证了模型预测值与观察值之间存在一致性，具有较好的拟合情况。DCA结果显示，影像组学模型、联合模型预测骨质疏松的临床价值均优于临床模型（图6）。

点击查看大图

图5  不同模型区分非骨质疏松与骨质疏松的受试者工作特征曲线（ROC）。5A：训练组中的模型比较；5B：测试组中的模型比较。
图6  模型的校正曲线（6A）、决策曲线（6B）。Clinical Signature：基于临床特征的模型；Rad Signature：基于T2WI图像的模型；Nomogram：基于临床联合T2WI图像的联合模型；AUC：ROC曲线下面积；CI：置信区间。
Fig. 5  The receiver operating characteristic (ROC) curve of different models to distinguish between non-osteoporosis and osteoporosis. 5A: Comparison of models in the training group; 5B Comparison of models in the test group.
Fig.6  Calibration plot (6A) and decision curve (6B) diagram for model. Clinical Signature: a model based on clinical features; Rad Signature: a model based on T2WI images; Nomogram: a model based on combined model; AUC: area under the curve; CI: confidence interval.

点击查看表格

表3  各模型预测骨质疏松的效能
Tab. 3  Effectiveness of each model in predicting osteoporosis

3 讨论

       本研究通过对291例骨质疏松及非骨质疏松患者病例进行回顾性分析，并基于临床特征、腰椎T2WI序列影像组学特征建立相应模型，创新性地构建多种不同机器学习算法，结果显示基于LR的腰椎T2WI的影像组学模型对骨质疏松症有好的预测能力，它可以在常规诊断的基础上提供额外的骨质疏松与否的信息，这种方法可以潜在地减少医疗费用和辐射暴露，为患者提供巨大的临床获益。

3.1 影像组学研究方法的优势及相关研究比较

       影像组学是一种新兴的领域，能够无创地提取数字医学图像中肉眼无法观察到的高维数据，可以反映特定区域的异质性，已被广泛用于无创性评估肿瘤的临床诊断和预后^{[18, 19]}。目前，已有研究者将影像组学用于预测骨质疏松症，并取得了不错的效果，其中，JIANG等^[20]利用基于CT腰椎图像提取了包括一阶特征5个、GLCM特征1个、GLSZM特征4个和NGTDM特征2个在内的12个组学特征用于腰椎术前检测骨质疏松症，并取得了很好的预测价值（AUC=0.96）；ZHANG等^[21]利用CT腰椎图像提取包含形状特征（Shape）2个、一阶特征4个、GLCM特征1个、GLSZM特征4个和GLRLM特征3个在内的14个组学特征用于预测急性与慢性骨质疏松性椎体骨折，并取得了良好的预测效果（AUC=0.93）。在本研究中，一阶特征、纹理特征得以入选可能与在骨质疏松椎体中骨量丢失引发骨小梁变细、数量减少、间隙增大、黄骨髓增加导致骨质异质性有关，这与先前研究^{[20, 22, 23]}结论一致。而形状特征直接反映骨骼完整性、病变形态和像素情况^{[24, 25, 26]}，本次特征筛选并未入选，可能与形状特征自身不需要数据矩阵变换特征有关，也可能与骨质疏松症患者椎体的三维大小、空间几何特性等一阶特征差异无统计学意义有关^{[27, 28]}。并且，本研究中，年龄和BMI差异具有统计学意义（P＜0.001），这可能与年龄增大引起体内激素水平变化导致破骨细胞活跃、BMI值较高的人群骨骼负重更多而促进成骨形成等因素有关，最终两者作为临床特征参与临床模型建模，这与部分研究^{[17, 29]}选取的临床特征建模一致，但性别差异无统计学意义（P＞0.05），因而未能入组，笔者认为可能与本研究人群均来源于因腰部疼痛就诊的有症状人群有关。

3.2 影像组学特征参数意义

       基于T2WI图像的影像组学诊断模型对海量T2WI影像潜在的、有价值的数据信息进行更深层次挖掘、预测和分析，有助于对骨质疏松临床结果的预测。本研究得到的影像组学特征为一阶特征、GLCM、GLRLM、GLSZM、NGTDM 5类。其中，一阶特征可以反映所有体素的直方图特性，选取的偏度（skewness）、中值（median）、均方根（root mean squared, RMS）、第十百分位（10th percentile）4个指标，说明此4个直方图参数在鉴别骨质疏松与否方面具有重要贡献。GLCM特征不仅反映灰度的分布特征，也反映具有同样灰度或者接近灰度的像素之间的位置分布特性，IMC1参数量化纹理的复杂性；集群突出是衡量偏度和不对称性的指标，值越高意味着不对称性越大；相关性参数衡量灰度共生矩阵元素在行或列方向或角度上的相似程度；IDN是度量图像的同质性的一种方法。GLRLM特征可量化图像中的像素值的分布，短行程强调（short run emphasis, SRE）较大表示越短的运行长度以及更精细的纹理；短行程低灰度强调（short run low gray level emphasis, SRLGLE）测量较低灰度值的较短运行长度的联合分布。GLSZM特征可以量化图像中连续像素值的区域，归一化灰度不均匀性（gray level non-uniformity normalized, GLNN）的值越低、强度值的相似性越高；区域大小不均匀性（size-zone non-uniformity, SZN）的值越低、图像中的均匀性越大；灰度方差（gray level variance, GLV）测量区域的灰度强度方差。NGTDM特征量化了一个灰度值与其相邻距离内的平均灰度值之间的差异，对比度（contrast）度量空间强化变化；繁忙度（busyness）度量从像素到其相邻像素的变化；复杂度（complexity）越大、图像越不均匀；强度（strength）度量图像中体素。本研究中的23个特征参数，在不同程度上反映了骨质疏松患者与非骨质疏松患者在T2WI图像灰度值分布、纹理特征及空间异质性方面的差异，笔者认为可能是骨质疏松患者的骨密度和骨质量减低导致，骨密度由骨内的矿物质决定，占骨强度的70%，骨质量包括骨的微观组织结构及分子水平变化，占骨强度的30%，这些都会导致椎体成分含量、空间排布及骨小梁的微空间几何结构形态、数量等变化，从而产生纹理特征的差异性。

3.3 基于T2WI序列MRI影像组学诊断模型鉴别腰椎骨质疏松与非骨质疏松的临床价值

       骨质疏松患者骨丢失及早期没有明显临床症状，往往因骨质疏松性骨折（脆性骨折）这样的严重临床后果就医^[30]，既往有研究表明^{[31, 32, 33]}，MRI主要是通过测量脂肪含量，或通过获取骨灌注、微循环以及骨质量等信息评估骨质疏松程度，这些方法都需要在常规序列基础上增加额外扫描序列或加以后处理，增加了扫描时间、影响了日常工作进度，因此更加需要探索出一种简便可靠、易于使用、可量化及个性化地鉴别骨质疏松与非骨质疏松的影像学方法。常规T2WI序列作为腰椎MR检查最基本的序列，信号强度取决于细胞丰富度、组织含水量、纤维含量及水肿等特征^{[34, 35, 36]}，但并不能很好地应用于骨质疏松的判断。本研究中基于LR建立的T2WI腰椎影像组学在训练组和测试组中均取得了较好的效能，改善了影像医师判断的主观性和肉眼难以识别的微小信号强度变化，相较于普通影像学资料对预测骨质疏松的结果更具有说服力。

3.4 本研究的局限性

       （1）样本量相对较少，对于机器学习来说，病例数越多，每次模型运行的结果越稳定、过拟合越小；（2）缺乏外部验证组，这可能会导致过拟合；（3）本研究的序列存在一定局限，只纳入了T2WI序列，未涉及到T1WI及脂肪抑制序列；（4）临床特征纳入不全，对于绝经与否、体表面积、抽烟喝酒等临床特征未进行收集，导致临床模型预测效能欠佳。因此，未来的研究中，应该使用更大的样本量和多中心数据来验证这些发现，再基于多序列的MRI联合临床模型，进一步探讨临床模型、影像组学模型及联合模型的预测效能。

4 结论

       综上所述，本研究的结果证明了基于腰椎T2WI的影像组学分析在判断非骨质疏松和骨质疏松患者方面的潜力，为临床决策提供依据。