0 引言

       胶质瘤是最常见的原发性中枢神经系统恶性肿瘤，约占恶性脑肿瘤的80%，发病率和死亡率最高^{[1, 2]}。根据世界卫生组织标准，胶质瘤可分为较低级别（2/3级）和高级别（4级）^[3]。其中，4级胶质母细胞瘤（glioblastoma, GBM）是恶性程度最高的胶质瘤，中位生存时间仅为11~15个月，易复发，预后差^{[4, 5, 6]}。分子病理学研究发现异柠檬酸脱氢酶1（isocitrate dehydrogenase, IDH1）突变和染色体1p/19q共缺失等分子标记物是胶质瘤的关键预后因素^[7]，然而分子标记物的准确检测需要足够的手术标本、专业的技术人员和昂贵的设备材料，无法满足GBM预后生存期术前无创精准预测的临床需求^{[8, 9, 10]}。

       近年来，影像组学在临床诊断和预后预测方面取得了显著成果。通过将医学图像转换为可挖掘的高通量特征数据，提供了对肿瘤更全面的量化信息，基于这些特征数据构建机器学习预测模型可有效地辅助临床决策^{[11, 12]}。现有研究初步证明基于术前MRI和正电子发射断层扫描（positron emission tomography, PET）图像提取的影像组学特征与胶质瘤患者的总生存期有显著相关性^{[13, 14]}，影像组学能够反映肿瘤的恶性生物学特征^[15]，从而提供了一种潜在的GBM风险分层方法。然而，如何选择和提取有意义的预后相关影像特征，并将这些特征与患者的临床病理信息结合，实现GBM风险分层仍然是一个挑战性问题。

       本文以术前多序列MRI影像组学方法开发一种GBM总生存期预测模型，为临床提供更为精准的风险分层和预后评估工具，提高GBM患者的个体化诊疗水平，提升患者生存获益。

1 材料与方法

1.1 研究对象

       本研究遵守《赫尔辛基宣言》，经山西医科大学第一医院伦理委员会批准，免除受试者知情同意，批准文号：2021-K-K073。回顾性分析癌症基因组图谱（The Cancer Genome Atlas, TCGA）/肿瘤免疫图谱数据库（The Cancer Immunome Atlas, TCIA）数据库中三个数据集[TCGA/GBM，加利福尼亚大学旧金山术前弥漫性胶质瘤数据集（UCSF/PDGM），宾夕法尼亚大学GBM数据集（UPENN/GBM）]中GBM患者的影像、临床、病理和基因数据。纳入标准：（1）根据世界卫生组织组织学分类标准病理证实的4级GBM^[3]；（2）术前无肿瘤相关治疗史；（3）术前对比增强后T1加权（post-contrast enhanced T1-weighted, T1CE）和T2加权液体衰减反转恢复（T2-weighted fluid attenuation inversion recovery, T2 FLAIR）序列图像清晰完整；（4）总生存期（overall survival, OS）、临床、病理信息完整；（5）确诊O⁶-甲基鸟嘌呤-DNA甲基转移酶（O⁶-methylguanine-DNA methyltransferase, MGMT）状态；（6）随访时间超过2年或终点事件（满足OS分析要求），OS 定义为从手术病理确诊到患者死亡或最后一次随访的时间。排除标准：（1）MR图像质量差（明显的运动伪影）；（2）既往有脑部手术史；（3）OS、临床、病理信息不完整；（4）MGMT状态未确诊。最终将UCSF/PDGM中249例患者和TCGA/GBM中60例患者纳入研究，按7∶3比例随机划分为训练集和测试集（表1）。

点击查看表格

表1  基线患者特征
Tab. 1  Baseline patient characteristics

1.2 MRI数据采集与图像预处理

       每位患者的T1CE和T2 FLAIR两种序列的MR图像从TCIA数据库中获取，临床、病理和IDH等基因数据从TCGA数据库中获取。首先将图像信号强度归一化为（0~255），将所有MRI信号强度值转化为标准化的强度范围，避免异质性偏倚。然后利用开源软件3D_Slicer（5.2.2版，https://www.slicer.org/）对所有MR图像进行颅骨和头皮剥离。考虑不同机型和扫描协议带来的图像分辨率差异，采用基于Pyradiomics程序包（https://github.com/Radiomics/pyradiomics）的开源软件FAE（0.5.5版，http://github.com/salan668/FAE）将每层MR图像重新采样到3 mm×3 mm的均匀像素间距。

1.3 MRI肿瘤、水肿和坏死感兴趣区的分割

       本研究采用BraTS挑战赛获胜的分割算法对肿瘤、水肿和坏死区域进行感兴趣区（region of interest, ROI）自动分割。然后，由一名高年资放射科医生对自动分割结果进行手动校正。ROI包括瘤周水肿区域、增强肿瘤区域和囊变坏死核心区域。ROI分割和校正主要基于MRI不同序列的信号特征，具体如下：（1）瘤周水肿区域，指肿瘤周围的水肿组织区域，在T2 FLAIR序列上呈高信号；（2）增强肿瘤区域，在T1增强序列上呈高信号的肿瘤实质部分，包括实体肿瘤和微血管增生区域；（3）囊变坏死核心区域，指由于肿瘤生长过快导致内部供血不足而形成的坏死组织，在T1增强序列上呈低信号，在T2 FLAIR序列上呈混合信号。图1为一位GBM患者的T1增强和T2 FLAIR原始图像，ROI分割图像和校正后ROI分割图像，分别以红色、黄色和绿色标出坏死、肿瘤和水肿的分割结果。

点击查看大图

图1  男，54岁，GBM患者。从左到右分别为T1CE和T2 FLAIR原始图像、ROI分割图像和校正后ROI分割图像，坏死、增强肿瘤和水肿的注释分别以红色、黄色和绿色表示。GBM：胶质母细胞瘤；T1CE：对比增强后T1加权；T2 FLAIR：T2加权液体衰减反转恢复；ROI：感兴趣区。
Fig. 1  Male, a 54-year-old patient with GBM. From left to right are the T1CE and T2 FLAIR original images, the ROI segmentation images, and the corrected ROI segmentation images. Annotations for necrosis, enhancing tumor, and edema are represented in red, yellow, and green, respectively. GBM: glioblastoma; T1CE: post-contrast enhanced T1-weighted; T2 FLAIR: T2-weighted fluid attenuation inversion recovery; ROI: region of interest.

1.4 特征提取

       对每位患者的T1增强和T2 FLAIR两个序列的三种ROI区域利用基于Pyradiomics程序包（https://github.com/Radiomics/pyradiomics）的开源软件FAE（0.5.5版，http://github.com/salan668/FAE）提取影像组学特征。提取的特征可分为4组：（1）形状特征组（n=14），主要量化ROI的形状和大小，如最大3D直径、体积、表面积等；（2）一阶统计特征组（n=18），基于图像灰度直方图统计分布量化肿瘤信号的非均质性，例如平均值、标准偏差、偏度、峰度等；（3）纹理特征组，描述图像内部的纹理分布关系量化肿瘤微观异质性，包括GLCM（灰度共生矩阵，n=22）、GLRLM（灰度游程长度矩阵，n=16）、GLSZM（灰度大小区域矩阵，n=16）、GLDM（灰度依赖矩阵，n=14）和NGTDM（邻近灰度差异矩阵，n=5）；（4）滤波特征组（n=1583），对原始图像应用LoG、Wavelet、Square、SquareRoot、Logarithm、Exponential、Gradient、LocalBinaryPattern（3D）共计8种滤波器进行滤波处理，然后基于滤波后图像ROI提取一阶统计特征和纹理特征，每一个序列的每一个ROI提取1 688个影像组学特征，每位患者2个序列3个区域共计提取10 128个影像组学特征。

1.5 影像组学特征筛选与组学标签的构建和验证

       为了选择OS相关的预后特征，首先进行相关性分析初筛特征，选择相关性系数>0.7的特征，即与目标变量有较强关联的特征，利用初筛特征进行单因素Cox回归分析，选择P<0.05的特征，对筛选后特征继续进行主成分分析（principal component analysis, PCA）降维，计算每个主成分解释的方差百分比以及累计方差百分比，找到覆盖95%比例变异性的主成分，然后进行LASSO-Cox回归分析，最终筛选出与OS显著相关的主成分特征，并依据所选主成分特征的各自系数（β）加权的线性组合计算风险评分（Risk-score）作为影像组学标签，公式如下：

       然后利用R中“Surv_point”函数使用最显著差异选择最佳Cut-off值，将训练集和测试集统一划分为高、低风险组。为了评估Risk-score作为预后相关影像组学标签的有效性，通过Kaplan-Meier（KM）生存分析和Log-rank检验来比较高风险组和低风险组患者间的生存差异。

1.6 GBM总生存期个体化预测列线图构建

       首先对临床、病理、基因信息进行单因素Cox回归，得到与OS显著相关的风险因素。将Risk-score与临床、病理、基因信息相结合，采用多因素Cox回归构建临床-影像组学联合预测模型及列线图，同时构建临床预测模型，采用一致性指数（C-index）评估预后预测模型的判别能力并进行模型效能比较，采用校正曲线来评估列线图预测生存概率和实际生存概率之间的一致性。本研究影像组学分析流程见图2。

点击查看大图

图2  影像组学分析流程图。T1CE：对比增强后T1加权；FLAIR：液体衰减反转恢复；ROI：感兴趣区；PCA：主成分分析；LASSO：最小绝对收缩和选择算法；KM：Kaplan-Meier；MGMT：O⁶-甲基鸟嘌呤-DNA甲基转移酶。
Fig. 2  Radiomics workflow diagram. T1CE: post-contrast enhanced T1-weighted; FLAIR: fluid attenuation inversion recovery; ROI: region of interest; PCA: principal component analysis; LASSO: least absolute shrinkage and selection operator; KM: Kaplan-Meier; MGMT: O⁶-methylguanine-DNA methyltransferase.

1.7 统计学分析

       本研究中所有统计分析均在R 4.2.2版本（https://www.R-project.org/）进行。采用“Survival”包进行单变量和多变量Cox回归分析和KM生存分析。分别使用“glmnet”和“rms”包构建LASSO-Cox回归模型和绘制列线图。所有统计学检验均为双侧检验，以P<0.05为差异有统计学意义。

2 结果

2.1 临床风险因素和临床模型

       训练集和测试集患者的年龄分别为（61.38±12.44）岁和（59.75±13.81）岁；男性分别占58.33%和58.06%。关于基因状态，两组GBM中MGMT的描述性统计结果以甲基化MGMT居多。年龄（P=0.307）、性别（P>0.999）、MGMT（P=0.564）差异均无统计学意义。采用单因素Cox回归确定年龄，MGMT是影响GBM总生存期的显著风险因素。利用多因素Cox回归构建临床预测模型，采用C-index评估模型预测效能（训练集：C-index=0.659，测试集：C-index=0.653）。

2.2 影像组学特征筛选和标签建立

       首先通过相关性分析，将10 128个特征初筛到8 509个特征，再利用单因素Cox回归筛选到673个特征，之后进行PCA降维，找到覆盖95%比例变异性的41个主成分特征（图3），并对选择的41个主成分特征进行LASSO-Cox回归（图4），最终筛选出与OS显著相关的16个主成分特征，通过所选主成分的各自系数（β）加权线性组合计算风险评分：Risk-score=PC1×0.0163+PC2×（-0.0022）+PC4×（-0.0728）+PC5×（-0.0308）+PC7×0.0826+PC8×0.0874+PC11×（-0.0123）+PC12×（-0.0430）+PC14×0.0356+PC23×（-0.0218）+PC24 ×（-0.0245）+PC27×0.0164+PC31×（-0.0161）+PC32×（-0.0361）+PC35×0.0122+PC39×（-0.0597）。

       利用R中“Surv_point”函数使用最大统计学差异选择最佳Cut-off值，在训练集和测试集中，OS均可以通过Risk-score得到显著风险分层（图5）。依据Risk-score建立影像组学预测模型，采用C-index评估模型预测性能（训练集：C-index=0.744，测试集：C-index=0.710）。

点击查看大图

图3  经过主成分分析降维后主成分数量与覆盖变异性比例的关系，选择覆盖95%以上变异性的41个主成分，随着选择的主成分数量大于41，对变异性的影响不显著。
图4  最小绝对收缩和选择算法回归模型的调优参数（lambda）筛选。相对于log(lambda)生成部分似然偏差，最小的部分似然偏差对应于最优的特征数量。虚线表示根据1个标准误差标准和最小标准得到的最佳lambda值。当剩余16个主成分时，部分似然偏差最小。
图5  基于Risk-score对患者划分高、低风险组，训练集（5A）：P<0.0001，测试集（5B）：P<0.0001，二者均<0.05，差异具有统计学意义。
Fig. 3  The relationship between the number of principal components after principal component analysis dimensionality reduction and the proportion of variance covered. A selection of 41 principal components, covering over 95% of the variance. As depicted in the figure, the impact on variance becomes non-significant beyond selecting more than 41 principal components.
Fig. 4  The tuning parameter (lambda) selection for the least absolute shrinkage and selection operator regression model. The minimum partial likelihood deviance relative to log(lambda) generation corresponds to the optimal number of features. The dashed line represents the best lambda value obtained based on one standard error criterion and the minimum standard. The minimum partial likelihood deviance is achieved when 16 principal components remain.
Fig. 5  Based on the Risk-score, the patients were classified into high-risk and low-risk groups. The training set (5A): P<0.0001 and the test set (5B): P<0.0001, both with P-values less than 0.05, The differences are statistically significant.

2.3 临床-影像组学联合预测模型及列线图的构建

       利用单因素Cox回归和多因素Cox回归，最终选择风险评分、年龄、MGMT状态构建临床-影像组学联合预测模型和列线图预测GBM患者的6月、12月和24月的生存概率（表2），发现临床-影像组学联合预测模型性能（训练集：C-index=0.768，测试集：C-index=0.724）显著高于临床预测模型（训练集：C-index=0.659，测试集：C-index=0.653)。列线图在训练集和测试集中均显示出良好的校正效果（图6）。此外，决策曲线分析显示，在合理阈值概率的大部分范围内，临床-影像组学联合预测模型总体净收益高于影像组学预测模型和临床预测模型（图7）。

点击查看大图

图6  临床-影像组学联合预测模型和列线图。6A：基于训练集结合Risk-score、患者年龄和O⁶-甲基鸟嘌呤-DNA甲基转移酶（MGMT）来开发的列线图；6B：训练集6个月、12个月、24个月生存期校准曲线；6C：测试集6个月、12个月、24个月生存期校准曲线，灰色的对角线表示理想标准，越接近灰色对角线表示评估性能越好。
图7  决策曲线图。曲线代表净利润，y轴表示净收益。
Fig. 6  Clinical-radiomics combined prediction model and nomogram. 6A: Develop a Nomogram based on the training set, integrating Risk-score, patient age, and O⁶-methylguanine-DNA methyltransferase (MGMT) status; 6B: Calibration curves for the 6-month, 12-month, and 24-month survival periods in the training set; 6C: Calibration curves for the 6-month, 12-month, and 24-month survival periods in the test set, the gray diagonal line represents the ideal standard, and the closer the curves are to the gray diagonal line, the better the performance evaluation.
Fig. 7  Decision curve plot. The curve represents the net profit, with the y-axis representing the net income.

点击查看表格

表2  训练集Cox回归分析
Tab. 2  Training set Cox regression analysis

3 讨论

       本研究基于之前的研究成果^{[16, 17, 18]}，基于高通量影像组学特征构建了基于T1CE和T2 FLAIR联合序列的GBM总生存期预测模型，并与临床病理信息结合，模型效能得到了极大的改善。本研究在GBM层面，使用常规MRI影像组学方法，采用PCA降维建立Risk-score进而构建GBM总生存期预测模型为国内外被首次提出。该模型的构建将为GBM患者的个体化诊疗水平和总生存期预测准确性产生积极影响，切实提升患者生存获益。

3.1 MRI序列和ROI分割区域的选择

       在研究序列选择方面，有研究表明^[19]，多序列联合模型效能普遍高于单序列模型效能，单序列模型受限于其单一视角的数据，如果该视角的数据受到噪声干扰，可能会导致模型性能下降。而多序列模型通过融合多个视角的信息，能够更好地抵抗噪声干扰，同时，多序列联合模型通过整合不同序列的特征，能够建立起一个更具鲁棒性的模型，这对于解决GBM异质性问题非常关键。TAN等^[20]联合T1CE、T2 FLAIR序列，针对肿瘤区域和瘤周水肿区域研究表明，影像组学特征可以独立预测患者的OS，其中临床-影像组学联合模型性能（训练集：C-index=0.764，测试集：C-index=0.758）显著高于影像组学模型（训练集：C-index=0.707，测试集：C-index=0.711）。在ROI分割区域选择方面，YANG等^[21]对GBM瘤周水肿亚区提取影像组学特征建立预测模型，并结合临床风险因素构建联合模型，同样联合模型展现出最好的预测能力（训练集：C-index=0.770，测试集：C-index=0.785）。与本研究相比，该研究联合T1CE、T2 FLAIR序列，将MR图像更加细致地分割为增强肿瘤区域、水肿区域和坏死区域，提取了更加全面的影像组学特征，最终建立的影像组学模型具备良好的预测性能（训练集：C-index=0.744，测试集：C-index=0.710），结合临床风险因素所构建的联合模型预测效能在所建立的模型中表现最优（训练集：C-index=0.768，测试集：C-index=0.724）。

3.2 特征筛选方法

       在特征筛选过程中，高维数据通常包含大量不相关的、冗余的和有噪声的特征，这可能导致维度的破坏和模型的过拟合^[22]，自动编码器，稀疏表示和局部保持投影相结合，Kruskal-Wallis检验、Variance Threshold等方法在降低特征维数方面产生了不错的效果^{[22, 23, 24, 25]}。本研究针对每个患者分别提取了10 128个特征，特征维数过高，为此本研究采用了PCA，以消除特征间的多重共线性，提高模型的解释能力和预测精度。PCA通过寻找数据中的最大方差方向，并在这些方向上创建新的、线性无关的变量（主成分），从而有效地简化了数据集的复杂性。这些新的主要组成成分保留了原始数据的大部分信息，但具有更低的维度，使得数据可以在较低的计算成本下更易于处理和解释。PCA的应用使我们能够在减少计算复杂性和保证模型性能之间找到一个理想的平衡点，显示出强大的降维能力和数据解释能力。经过PCA降维之后构建的模型效能得到显著提高，同时过拟合问题得到了一定的解决。

3.3 风险分层依据及模型构建方法

       本研究根据风险评分进行分层，在之前的一些研究中，有许多依据风险评分中位值或中值进行分层的例子^{[26, 27, 28]}，本研究采用“R”包中“Surv_point”函数通过比较每个可能的切点下两组之间的生存差异来确定最佳切点。对分层后的训练集和测试集通过KM曲线和Log-rank检验，P均<0.05，差异具有统计学意义，并且低风险组的生存概率显著高于高风险组。在模型构建方面，一些研究尝试使用随机森林和其他机器学习分类器来预测患者生存期^{[29, 30]}，生存分析被作为一个二分类任务，将生存分析作为二分类任务的主要局限性在于忽视了时间因素的作用。二分类任务基于特定时间点进行预测，但在实际的临床生存分析中，患者的生存时间是一个连续的变量，其分布可能受到年龄，性别，疾病状态等因素的影响，忽视这些因素可能会导致模型的预测能力下降。本研究利用Cox回归进行OS模型构建，将风险变量作为连续变量，开发了一种不依赖任何阈值的模型。

3.4 本研究的局限性

       本研究存在一些局限性。首先，数据来自不同的数据中心，MR仪器的型号、参数存在一定差异，这可能对结果造成一定的误差。其次，由于对每位患者提取了包括一阶、纹理、形状与大小以及8种滤波器方式的海量特征，在特征选择过程中进行了PCA降维，其局限性在于经过PCA降维后的主成分是若干原始特征的线性组合，这让对其生物学解释变得困难。同时，对于多中心数据，需要依据OS、临床、病理信息完整的入组标准，导致一些信息不能被选用，如KPS评分、IDH状态等。有研究表明，将影像组学联合更多的临床、病理、基因信息对模型效能的提高有显著的影响^{[31, 32]}。在特征提取过程中，使用传统方法提取特征的效率低于深度学习方法^{[33, 34, 35]}。在下一项目的研究中，团队将在保持数据量充足的同时收集丰富的临床、病理、遗传信息，并进行学习方法的迁移和组合，力求解决本研究存在的限制。

4 结论

       综上所述，本研究通过基于多序列和多区域的MRI影像组学分析，联合临床病理信息开发GBM总生存期预测模型具有较好的预测性能，通过这一模型，临床医生能够更有效精准地进行GBM患者的风险分层评估，从而对患者进行个体化的治疗方案制订。