基于多模態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 患者的精確診斷、化療使用決策及生存狀況預測提供臨床輔助價值。

作者利益沖突聲明：全體作者均聲明無利益沖突。

作者貢獻聲明：吳磊設計本研究的方案，對稿件重要內容進行了修改，獲得了吳階平臨床科研專項基金的資助；焦凱劍起草和撰寫稿件，獲取、分析或解釋本研究的數據；楊波、陳文、方鈺輝、陳亞琳獲取、分析或解釋本研究的數據，對稿件重要內容進行了修改，其中楊波獲得了湖北省自然科學基金的資助；全體作者都同意發表最后的修改稿，同意對本研究的所有方面負責，確保本研究的準確性和誠信。