基于多模态MRI 影像组学的胶质母细胞瘤生境亚区预测MGMT启动子甲基化表达

焦凯剑，杨波，陈文，方钰辉，陈亚琳，吴磊

0 前言

胶质母细胞瘤（glioblastoma, GBM）是最具侵袭性和致死性的原发脑肿瘤[1-2]。目前，GBM 主流治疗方法仍为手术治疗及术后化疗，化疗药物通常使用烷化剂替莫唑胺（temozolomide, TMZ）[3]。O6-甲基鸟嘌呤-DNA甲基转移酶（O6-methylguanine-DNA methyltransferase,MGMT）基因启动子的甲基化状态可以评价神经GBM患者对TMZ 的敏感性，并预测TMZ 的疗效[4-5]。这成为选择特异性个体治疗的重要依据。对于对TMZ 无效的患者，可以避免进行化疗，从而减少过度治疗和化疗带来的毒副作用。目前，评判MGMT 基因启动子甲基化状态仍依赖于有创的免疫组织化学方法，需基于手术标本[6]。

然而，近年来，影像组学（Radiomics）的兴起为术前非侵入性评估GBM MGMT启动子的甲基化状态提供了新的方法[7]。已有多项研究探讨利用纹理特征来预测MGMT 甲基化状态的有效性[8-10]。相关研究表明，相较于肿瘤区域，水肿区域也包含肿瘤微环境的信息[11-12]。先前的关于GBM 的工作中，没有系统性地理解肿瘤生境栖息地不同亚区（增强区域、坏死区域和水肿区域）的放射组学特征与肿瘤病理生理信息的相关性。多数研究仍是以整体肿瘤为感兴趣区进行研究，对于GBM在空间上不同区域的图像价值未进行完善的研究。弥散技术相较于结构MRI除了能够提供肿瘤位置形态之外，还能提供与肿瘤基因相关的更全面的信息。此外先前研究多使用来自单中心的研究队列，无法进一步确保模型的泛化能力和鲁棒性。

因此，在本项回顾性多中心研究中，我们基于自动分割将GBM 分为三个生境亚区（坏死、增强和瘤周水肿），对比了不同肿瘤生境亚区下多模态MRI 图像的放射组学特征，以评估其对预测GBM 中MGMT 启动子甲基化状态的能力和有效性，旨在构建一个具有较高鲁棒性和准确性的多生境区域多模态MRI 放射组学模型，术前无创预测GBM 患者的MGMT 启动子甲基化状态。

1 材料与方法

1.1 一般资料

本研究回顾性分析的数据来自湖北医药学院附属太和医院（Taihe Hospital, the Affiliated Hospital of Hubei University of Medicine, HBMUTH）、宾法尼西亚大学（University of Pennsylvania,UPENN）以及加州大学弗朗西斯科分校（University of California San Francisco, UCSF）[13-14]。HBMUTH数据收集通过湖北医药学院附属太和医院伦理委员会批准，免除受试者知情同意，批准文号：2023KS11，收集日期为2014 年10 月至2023 年7 月。UPENN 和UCSF 数据为公共数据库资源，二者的完整数据已分别获得UPENN 和UCSF 的伦理委员会批准且均可在癌症影像档案（The Cancer Imaging Archive, TCIA）获得，免除受试者知情同意。此外，UPENN 数据收集日期为2006 年至2018 年，UCSF 数据收集日期为2015 年至2021 年。本研究的纳入标准：（1）符合2016、2021 WHO 中枢神经系统肿瘤分级为GBM 的患者；（2）术前MRI 图像包括对比增强T1 加权成像（contrast enhanced T1-weighted imaging,T1WI-CE）序列、T2 液体衰减反转恢复（T2 fluid attenuation inversion recovery, T2-FLAIR）序列和弥散张量成像（diffusion tensor imaging, DTI）的各向异性分数（fractional anisotropy, FA）参数图，并且序列完整、图像清晰；（3）完整的MGMT分子信息及临床信息。排除标准：排除低级别胶质瘤。最终，600 例患者病例符合纳入标准，其中HBMUTH、UPENN和SPPH数据集分别有13、214和373例。

1.2 MRI扫描

HBMUTH 数据采用48 通道阵列线圈的3.0 T Scanner （Discovery 750W, SIGNA Architect; GE Healthcare）进行MRI 扫描。采集协议具体采用T1WI-CE 序列（TR/TE 2000-2600 ms/20-22 ms，FOV 240 mm×240 mm，层厚5.0 mm，矩阵512×512）、T2-FLAIR 序列（TR/TE 9000-9900 ms/130-150 ms，FOV 240 mm×240 mm，层厚5.0 mm，矩阵512×512）和DTI序列（TR/TE 9500 ms/98 ms，FOV 240 mm×240 mm，层厚5.0 mm，矩阵256×256），由GE 后处理工作站对DTI图像进行后处理得到FA图。注射0.1 mmol/kg的钆喷替酸葡甲胺（Gd-DPTA）后获得T1WI-CE 图像。UPENN 术 前MRI 图 像 使 用1.5 T 和3.0 T 扫 描 仪（Siemens）获得，所有患者术前均扫描T1WI-CE 序列（矩阵192×256×192，分辨率0.98×0.98×1.00；TR/TE 1760 ms/3.1 ms），T2-FLAIR 序列（矩阵192×256×60，分辨率0.94×0.94×3.00，TR/TE 9420 ms/141 ms）和DTI 序列（矩阵128×128×40，分辨率1.72×1.72×3.00，TR/TE 3100-8000 ms/104 ms）。注射0.3 mL/kg 钆贝酸二葡甲胺后获得T1WI-CE 图像。由Siemens 后处理工作站对DTI 图像进行后处理得到FA 图。UCSF 术前MRI 图像均使用3.0 T 扫描仪（Discovery 750, GE Healthcare）和专用的8 通道头部线圈（Invivo，Gainesville），扫描T2-FLAIR序列（TR/TE 5700 ms/115 ms，层厚1.2 mm，矩阵256×256）、T1WI-CE 序列（TR/TE 6 ms/2.3 ms，层厚1.0 mm，矩阵256×256）、DTI 序列（TR/TE 8400 ms/73 ms，层厚2 mm，矩阵128×128）；在研究期间，使用两种钆基对比剂：0.1 mL/kg 钆布醇（gadobutrol,Gadovist, Bayer）和0.2 mL/kg 钆 特 酸 葡 胺（gadoterate, Dotarem, Guerbet）。DTI 数据经涡流校正并使用来自FSL版本6.0.2（FMRIIB,https://fsl.fmrib.ox.ac.uk/fsl/fslwiki）的涡流和DTIFIT 模块处理，产生FA图像。

1.3 肿瘤感兴趣区勾画

在预处理过程中，包括去除颅骨、配准、重采样以及图像强度归一化处理等步骤。针对经过预处理后的图像，我们使用先前多模态脑肿瘤分割挑战赛（Multimodal Brain Tumor Segmentation Challenges,BraTS）中的先进分割算法集成模型进行自动分割[15-17]。然后，由一位具有15年神经放射科工作经验的医生利用ITK-SNAP 软件（http://www.itksnap.org）进行手动校正。图像的分割结果涵盖了三个主要肿瘤区域：增强肿瘤、中心非增强和/或坏死肿瘤以及周围FLAIR异常（由非增强肿瘤和相关水肿区域组成），如图1所示。

图1 肿瘤生境亚区分割图。第一行为79 岁男性胶质母细胞瘤O6-甲基鸟嘌呤-DNA 甲基转移酶（MGMT）启动子甲基化患者，第二行为52 岁男性胶质母细胞瘤非甲基化患者。上下两行从左到右依次为对比增强T1WI序列（T1WI-CE）、液体衰减反转恢复（FLAIR）序列和各向异性分数（FA）参数图。红色亚区为坏死区域，黄色亚区为增强区域，绿色亚区为水肿区域。Fig.1 Diagram of tumor habitat subsegmentation.The first line is a 79-year-old male glioblastoma O6-methylguanine-DNA methyltransferase(MGMT) promoter methylation patient, and the second line is a 52-year-old male glioblastoma non-methylated patient.The upper and lower rows, from left to right, are the contrast-enhanced T1-weighted imaging (T1WI-CE)sequence, fluid attenuation inversion recovery sequence (FLAIR), and fractional anisotropy (FA) map.The red subregion is necrotic region, the yellow subregion is the enhanced region, and the green subregion is the edema region.

1.4 MRI特征提取

为了从MRI 图像中高通量地提取肿瘤部分的影像组学特征，我们将使用Python 3.10 来调用PyRadiomic 包。该包可以从T1WI-CE、FLAIR 和FA 三个序列的3 个亚区域中提取影像组学特征，总共有9 种区域进行特征提取。这些特征可以分为三组：几何形状、强度和纹理。几何特征描述了肿瘤的三维形状特征。强度特征描述了肿瘤内体素强度的一阶统计（Firstorder）分布。纹理特征描述了强度的模式或二阶和高阶空间分布。我们将使用多种方法来提取纹理特征，包括灰度共现矩阵（gray level co-occurrence matrix, GLCM）、灰度运行长度矩阵（gray level run length matrix, GLRLM）、灰度大小区域矩阵（gray level size zone matrix, GLSZM）、灰度依赖矩阵（gray level dependence matrix,GLDM）和邻域灰度差异矩阵（neighborhood gray tone difference matrix, NGTDM）。对于每个患者的MRI单序列中的每个区域，我们将提取2153个特征：14个形状特征、18个一阶特征、75个纹理特征以及基于多种滤波变换所获得的一阶特征和纹理特征。这些特征可以对肿瘤的不同维度进行量化分析，以提取其特征和特性。更详细的特征介绍请参考PyRadiomic官网（https://pyradiomics.readthedocs.io/en/latest/features.html#）。

1.5 影像组学特征筛选和模型构建

本研究将HBMUTH 数据、UPENN 数据和UCSF 数据整合，并以7∶3比例随机分为训练集和测试集。在进行特征筛选之前，对训练集特征进行了归一化处理，将不同特征缩放到相同数量级。接着，对训练集进行了特征筛选，具体分为以下四步：（1）采用Pearson相关性分析来衡量特征与甲基化状态的相关性，并排除相关性绝对值小于0.05的特征。（2）采用最小冗余最大相关算法（minimum redundancy maximum relevance, MRMR）将复杂的变量系统分解成一组最相关的变量，以减少冗余。（3）通过对变量之间的相关性进行分析，选择了TOP 100 特征，并确保这些特征之间没有重叠。（4）采用Boruta 算法来选择“confirm”特征。Boruta是一种用于选择相关特征的包装器算法，它通过比较原始特征的重要性和随机特征的重要性来自上而下地搜索相关特征。在每次迭代中检查当前真实特征是否比最好的阴影特征具有更高的重要性，并逐步剔除那些被视为不重要的特征，以最大程度地减少误差[18]。基于这些步骤，我们选择了XGBoost来构建预测模型。XGBoost是一种强大的梯度提升算法，以其出色的预测性能、准确的分类能力以及对不平衡数据的有效管理而备受青睐。同时，该算法还应用了L1 和L2 正则化，以避免过拟合[19]。目前，它在许多分类竞赛中得到广泛应用。根据各区域的影像学特征，我们使用XGBoost开发组学模型，以预测GBM中MGMT基因启动子甲基化的有效性。以受试者工作特征（receiver operating characteristic,ROC）曲线及其曲线下面积（area under the curve,AUC）、准确度、敏感度和特异度等指标评价模型诊断效能，DeLong检验对比模型差异性。

1.6 统计分析

进行卡方检验以确定两组间性别的差异性，使用独立样本t检验评估年龄分布的差异，P＜0.05 认为差异具有统计学意义。本研究特征筛选部分采用R软件4.1.3版本（www.R-project.org），模型构建采用Python 中的scikit-learn 包（scikit-learn 版本0.23，Python 版本3.10），统计分析采用 SPSS 软件27.0版本（https://www.ibm.com/analytics/spssstatistics-software）实现。

2 结果

2.1 临床资料

MGMT 基因启动子的两种分子亚型在训练集和测试集中的临床资料统计描述如表1 所示。在训练集及测试集中，两种分子亚型在年龄和性别方面P均＞0.05，表示在训练集和测试集上两种分子亚型之间的差异没有统计学意义。

表1 训练集和测试集临床资料的比较Tab.1 Comparison of clinical data in training set and validation set

2.2 影像组学结果

2.2.1 特征筛选结果

通过Pearson 相关性分析、MRMR 和Boruta 进行降维处理后，我们筛选出了每组的最佳特征。在增强区域FA+FLAIR+T1WI-CE 序列得到了8 个特征。在坏死区域FA+FLAIR+T1WI-CE 序列得到了7 个特征。在瘤周水肿区域中FA+FLAIR+T1WI-CE 序列得到了10 个特征。最终，多模态多区域的特征集上得到10 个特征，其中2 个属于增强区域，3 个属于水肿区域，5个属于坏死区域，具体结果见表2和图2。

表2 Brouta降维处理后影像组学特征及其重要性平均值Tab.2 Important features and importance mean values in each group after Brouta dimensionality reduction

图2 Brouta特征筛选图。Fig.2 Features selection map of Brouta.

2.2.2 模型诊断效能评估

GBM 三个肿瘤生境亚区（坏死，水肿，增强）的影像组学模型对MGMT基因启动子甲基化状态的预测中均表现出较优的性能（图3）。在每个区域中，多模态（T1WI-CE_FA_FLAIR）组合构建影像组学模型是优于单一序列模型及结构序列组合（T1WI-CE_FLAIR）。所构建的具有3 个模态和3 个区域的多模态多区域联合模型在训练集上的AUC为0.874（0.841-0.906），在测试集上的AUC 为0.899 （0.853-0.944）（表3）。经DeLong 检验，测试集上多区域多模态组合模型与增强和水肿的联合模型存在显著差异（Z=2.65,P＜0.05），与其他模型间不存在显著差异（P＞0.05）。

图3 训练集和测试集多模态多区域影像组学模型受试者工作特征曲线。Fig.3 Receiver operating characteristic curves of the multi-region multi-modal radiomics model for the training and validation sets.

3 讨论

在本研究中，我们利用多模态预处理MRI 图像（T1WI-CE、FLAIR和FA）从肿瘤生境的不同亚区（增强肿瘤、瘤周水肿和非增强肿瘤区域以及肿瘤坏死区域）获得的放射组学特征，建立一个用于预测GBM的MGMT甲基化状态的多模态多区域影像组学模型。结果显示模型具有良好的预测性能（训练集和测试集的AUC分别为0.874 和0.899）。本研究在GBM 层面使用常规MRI和弥散MRI参数FA结合的影像组学方法，预测MGMT甲基化状态为国内外首次被提出。MGMT甲基化状态的准确预测为GBM患者分子分型的精准诊断、TMZ的临床治疗决策及生存期预测提供了重要的临床辅助价值。

3.1 生境成像对于GBM的评估和应用

生境一般指物种或物种群体赖以生存的自然生态环境。肿瘤本身作为一个复杂的类似于自然生态系统的生态部落，其内部的异质性亚区类似于自然生态系统的栖息地，每个亚区遵循适者生存原则，在各种压力作用下进行生长增殖和治疗抵抗，与肿瘤的进展和预后息息相关。先前通过反复穿刺活检以证实肿瘤内部异质性，但小样本组织并不能全面揭示，且受限于患者的配合程度，操作者的熟练程度等条件。生境成像根据放射组学理念，基于病理学和生物学等差异，利用定量影像学标志可以全面非侵入性地表征肿瘤环境，可视化和量化肿瘤内部异质性。对于GBM，放射组学的能力已经扩展到捕获肿瘤栖息地中的疾病异质性。先前常用无监督算法识别肿瘤生境亚区，我们根据Brats脑肿瘤分割比赛的排名top分割算法，将肿瘤自动划分为三个具有可能导致不同生物学行为区域：水肿区域、增强区域和坏死区域[17,20-22]。水肿区域通常位于肿瘤周边，反映肿瘤正常组织对肿瘤的生物学反应。增强区域反映肿瘤的血供情况，其形状和分布可以提供关于肿瘤边缘和边界的信息。坏死区域常由于血供不足导致缺血和坏死，对肿瘤的生长和扩散产生影响。

我们基于多模态MRI 图像每个亚区可获取6000 多个高维特征。因此，特征筛选是构建影像组学模型前的关键步骤。本研究中所采取的三步法从不同方面进行特征筛选。基于随机森林的Boruta算法筛选出所有与因变量具有相关性的特征集合，帮助我们更全面理解因变量的影响因素，得到更优的影像组学特征[23]。经过特征筛选我们得到10 个特征以描述肿瘤异质性：2个特征来自FA序列增强区域；3个水肿区域的特征分别来自于本研究所纳入的三个序列；5 个坏死区域的特征来自于FA 和FLAIR 序列。所获取的特征以滤波处理后的图像强度和纹理特征为主，通过滤波变换挖掘图像更深层次与肿瘤生物学变化相关信息。VERMA等[21]从肿瘤的增强、坏死、水肿生境中提取影像组学特征，特征筛选后所保留的纹理特征和强度特征与肿瘤细胞浸润及血管增生、坏死增强交界的缺氧区域密切相关。此外，来自于DTI参数图FA的特征占大多数提示功能序列可能包含更多的与GBM MGMT 甲基化等相关的分子标记物的状态变化的特征信息，在之后要增加这一部分的研究。

3.2 影像组学模型对比分析

结合多序列多模态MRI 和生境成像技术可以帮助我们更好地理解GBM 的病理生理变化。在本研究中从多模态MRI 序列（T1WI-CE、FLAIR、FA）提取肿瘤三个生境亚区的影像组学特征，建立用于预测MGMT甲基化状态的影像组学模型。结果显示在训练集中的AUC 为0.874，在测试集中的AUC 为0.899，多模态多区域联合的影像组学模型更具有鲁棒性。先前几项进行影像组学的研究多是基于结构序列进行胶质瘤MGMT 甲基化状态预测的研究。SHA 等[24]基于498 名异柠檬酸脱氢酶突变型胶质瘤患者的T1CE 和FLAIR图像预测MGMT启动子甲基化状态，在测试集和独立验证集中AUC 分别为0.916 和0.866。此外，两位学者基于T1CE 的组学模型在预测MGMT 启动子甲基化状态方面的AUC 为0.9 以上[25-26]。XI 等[27]研究证明T1、T2 以及T1CE 序列放射组学特征作为预测GBM中MGMT 启动子甲基化潜在影像学标记，训练集准确度为0.86，测试集准确度为0.80。部分研究在多参数多序列组合时选择纳入功能序列相关图像（如ADC）对MGMT 甲基化的表达进行预测，也表现出优于结构序列的性能[28-30]。这些研究多是以整个肿瘤作为高通量信息的获得区域，这样就无法与潜在的肿瘤微环境相联系起来。不同模态的MRI 序列所拥有的信息是不相同的，每个生境都包含了独特的微环境信息，结合二者在GBM患者的肿瘤病理生理变化评估和帮助新辅助治疗等方面具有广阔前景。

3.3 本研究的优越性与局限性

多模态MR 生境成像研究相较于以往传统结构MRI 不仅表征了肿瘤解剖结构，还进一步评估了患者的病理生理情况和潜在的肿瘤微环境相联系。Brats分割算法相较于手动分割能提供更为准确的生境亚区，更准确地描述胶质瘤复杂的空间异质性。因此在此基础上所构建的影像组学模型能够获取到更多的信息层面以及更全面地了解GBM 特征，有助于改进GBM的治疗和诊断策略。此外，特征选择方面，特征组合带来的非线性关系更适合于Boruta算法以获取鲁棒特征。XGBoost在许多机器学习竞赛和实际应用中表现出更高的鲁棒性，同时减少了过拟合的风险。

当然，本研究依然属于回顾性研究，样本量虽然整体较多并来自于多个中心，但是当地的国内的患者样本量与公开数据库中的样本量处于极不平衡的状态，下一步还需多与国内的其他研究中心或医院合作，开展多中心研究，进一步增加病例数。此外，MRI中其他具有临床意义的功能序列如灌注成像、弥散成像相关的其他参数及更高级的弥散技术也将被纳入研究。除了机器学习外，深度学习也越来越火热，在胶质瘤分割分类方面也可进行这方面的研究[31-34]。

4 结论

综上所述，本研究通过构建多模态生境亚区影像组学模型术前无创预测GBM 的MGMT 启动子甲基化状态，该模型能够带来鲁棒性更高的预测性能，为GBM 患者的精确诊断、化疗使用决策及生存状况预测提供临床辅助价值。

作者利益冲突声明：全体作者均声明无利益冲突。

作者贡献声明：吴磊设计本研究的方案，对稿件重要内容进行了修改，获得了吴阶平临床科研专项基金的资助；焦凯剑起草和撰写稿件，获取、分析或解释本研究的数据；杨波、陈文、方钰辉、陈亚琳获取、分析或解释本研究的数据，对稿件重要内容进行了修改，其中杨波获得了湖北省自然科学基金的资助；全体作者都同意发表最后的修改稿，同意对本研究的所有方面负责，确保本研究的准确性和诚信。