基于多序列MRI 影像組學的膠質(zhì)母細胞瘤風險分層預測研究

牛文舉，徐懷文，高宇翔，王效春2,，譚艷2,，張輝2,，楊國強2,*

0 引言

膠質(zhì)瘤是最常見的原發(fā)性中樞神經(jīng)系統(tǒng)惡性腫瘤，約占惡性腦腫瘤的80%，發(fā)病率和死亡率最高[1-2]。根據(jù)世界衛(wèi)生組織標準，膠質(zhì)瘤可分為較低級別（2/3 級）和 高 級 別（4 級）[3]。其 中，4 級 膠 質(zhì) 母 細 胞瘤（glioblastoma, GBM）是惡性程度最高的膠質(zhì)瘤，中位生存時間僅為11～15個月，易復發(fā)，預后差[4-6]。分子病理學研究發(fā)現(xiàn)異檸檬酸脫氫酶1（isocitrate dehydrogenase, IDH1）突變和染色體1p/19q共缺失等分子標記物是膠質(zhì)瘤的關鍵預后因素[7]，然而分子標記物的準確檢測需要足夠的手術標本、專業(yè)的技術人員和昂貴的設備材料，無法滿足GBM 預后生存期術前無創(chuàng)精準預測的臨床需求[8-10]。

近年來，影像組學在臨床診斷和預后預測方面取得了顯著成果。通過將醫(yī)學圖像轉換為可挖掘的高通量特征數(shù)據(jù)，提供了對腫瘤更全面的量化信息，基于這些特征數(shù)據(jù)構建機器學習預測模型可有效地輔助臨床決策[11-12]?，F(xiàn)有研究初步證明基于術前MRI和正電子發(fā)射斷層掃描（positron emission tomography, PET）圖像提取的影像組學特征與膠質(zhì)瘤患者的總生存期有顯著相關性[13-14]，影像組學能夠反映腫瘤的惡性生物學特征[15]，從而提供了一種潛在的GBM 風險分層方法。然而，如何選擇和提取有意義的預后相關影像特征，并將這些特征與患者的臨床病理信息結合，實現(xiàn)GBM風險分層仍然是一個挑戰(zhàn)性問題。

本文以術前多序列MRI影像組學方法開發(fā)一種GBM 總生存期預測模型，為臨床提供更為精準的風險分層和預后評估工具，提高GBM 患者的個體化診療水平，提升患者生存獲益。

1 材料與方法

1.1 研究對象

本研究遵守《赫爾辛基宣言》，經(jīng)山西醫(yī)科大學第一醫(yī)院倫理委員會批準，免除受試者知情同意，批準文號：2021-K-K073?；仡櫺苑治霭┌Y基因組圖譜（The Cancer Genome Atlas, TCGA）/腫瘤免疫圖譜數(shù)據(jù)庫（The Cancer Immunome Atlas, TCIA）數(shù)據(jù)庫中三個數(shù)據(jù)集[TCGA/GBM，加利福尼亞大學舊金山術前彌漫性膠質(zhì)瘤數(shù)據(jù)集（UCSF/PDGM），賓夕法尼亞大學GBM 數(shù)據(jù)集（UPENN/GBM）]中GBM 患者的影像、臨床、病理和基因數(shù)據(jù)。納入標準：（1）根據(jù)世界衛(wèi)生組織組織學分類標準病理證實的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）確診O6-甲基鳥嘌呤-DNA 甲基轉移酶（O6-methylguanine-DNA methyltransferase, MGMT）狀態(tài)；（6）隨訪時間超過2 年或終點事件（滿足OS 分析要求），OS 定義為從手術病理確診到患者死亡或最后一次隨訪的時間。排除標準：（1）MR 圖像質(zhì)量差（明顯的運動偽影）；（2）既往有腦部手術史；（3）OS、臨床、病理信息不完整；（4）MGMT 狀態(tài)未確診。最終將UCSF/PDGM 中249 例 患 者 和TCGA/GBM 中60例患者納入研究，按7∶3比例隨機劃分為訓練集和測試集（表1）。

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

1.2 MRI數(shù)據(jù)采集與圖像預處理

每位患者的T1CE 和T2 FLAIR 兩種序列的MR圖像從TCIA 數(shù)據(jù)庫中獲取，臨床、病理和IDH 等基因數(shù)據(jù)從TCGA 數(shù)據(jù)庫中獲取。首先將圖像信號強度歸一化為（0～255），將所有MRI 信號強度值轉化為標準化的強度范圍，避免異質(zhì)性偏倚。然后利用開源軟件3D_Slicer（5.2.2 版，https://www.slicer.org/）對所有MR圖像進行顱骨和頭皮剝離?？紤]不同機型和掃描協(xié)議帶來的圖像分辨率差異，采用基于Pyradiomics程序包（https://github.com/Radiomics/pyradiomics）的開源軟件FAE（0.5.5 版，http://github.com/salan668/FAE）將每層MR 圖像重新采樣到3 mm×3 mm 的均勻像素間距。

1.3 MRI腫瘤、水腫和壞死感興趣區(qū)的分割

本研究采用BraTS 挑戰(zhàn)賽獲勝的分割算法對腫瘤、水腫和壞死區(qū)域進行感興趣區(qū)（region of interest,ROI）自動分割。然后，由一名高年資放射科醫(yī)生對自動分割結果進行手動校正。ROI 包括瘤周水腫區(qū)域、增強腫瘤區(qū)域和囊變壞死核心區(qū)域。ROI分割和校正主要基于MRI不同序列的信號特征，具體如下：（1）瘤周水腫區(qū)域，指腫瘤周圍的水腫組織區(qū)域，在T2 FLAIR 序列上呈高信號；（2）增強腫瘤區(qū)域，在T1 增強序列上呈高信號的腫瘤實質(zhì)部分，包括實體腫瘤和微血管增生區(qū)域；（3）囊變壞死核心區(qū)域，指由于腫瘤生長過快導致內(nèi)部供血不足而形成的壞死組織，在T1 增強序列上呈低信號，在T2 FLAIR 序列上呈混合信號。圖1 為一位GBM 患者的T1 增強和T2 FLAIR 原始圖像，ROI 分割圖像和校正后ROI 分割圖像，分別以紅色、黃色和綠色標出壞死、腫瘤和水腫的分割結果。

圖1 男，54 歲，GBM 患者。從左到右分別為T1CE 和T2 FLAIR 原始圖像、ROI 分割圖像和校正后ROI 分割圖像，壞死、增強腫瘤和水腫的注釋分別以紅色、黃色和綠色表示。GBM：膠質(zhì)母細胞瘤；T1CE：對比增強后T1加權；T2 FLAIR：T2加權液體衰減反轉恢復；ROI：感興趣區(qū)。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 區(qū)域利用基于Pyradiomics 程序包（https://github.com/Radiomics/pyradiomics）的開源軟件FAE（0.5.5 版，http://github.com/salan668/FAE）提取影像組學特征。提取的特征可分為4 組：（1）形狀特征組（n=14），主要量化ROI 的形狀和大小，如最大3D 直徑、體積、表面積等；（2）一階統(tǒng)計特征組（n=18），基于圖像灰度直方圖統(tǒng)計分布量化腫瘤信號的非均質(zhì)性，例如平均值、標準偏差、偏度、峰度等；（3）紋理特征組，描述圖像內(nèi)部的紋理分布關系量化腫瘤微觀異質(zhì)性，包括GLCM（灰度共生矩陣，n=22）、GLRLM（灰度游程長度矩陣，n=16）、GLSZM（灰度大小區(qū)域矩陣，n=16）、GLDM（灰度依賴矩陣，n=14）和NGTDM（鄰近灰度差異矩陣，n=5）；（4）濾波特征組（n=1583），對原始圖像應用LoG、Wavelet、Square、SquareRoot、Logarithm、Exponential、Gradient、LocalBinaryPattern（3D）共 計8 種濾波器進行濾波處理，然后基于濾波后圖像ROI提取一階統(tǒng)計特征和紋理特征，每一個序列的每一個ROI 提取1 688 個影像組學特征，每位患者2 個序列3個區(qū)域共計提取10 128個影像組學特征。

1.5 影像組學特征篩選與組學標簽的構建和驗證

為了選擇OS 相關的預后特征，首先進行相關性分析初篩特征，選擇相關性系數(shù)＞0.7 的特征，即與目標變量有較強關聯(lián)的特征，利用初篩特征進行單因素Cox 回歸分析，選擇P＜0.05 的特征，對篩選后特征繼續(xù)進行主成分分析（principal component analysis,PCA）降維，計算每個主成分解釋的方差百分比以及累計方差百分比，找到覆蓋95%比例變異性的主成分，然后進行LASSO-Cox 回歸分析，最終篩選出與OS 顯著相關的主成分特征，并依據(jù)所選主成分特征的各自系數(shù)（β）加權的線性組合計算風險評分（Risk-score）作為影像組學標簽，公式如下：

然后利用R中“Surv_point”函數(shù)使用最顯著差異選擇最佳Cut-off 值，將訓練集和測試集統(tǒng)一劃分為高、低風險組。為了評估Risk-score 作為預后相關影像組學標簽的有效性，通過Kaplan-Meier（KM）生存分析和Log-rank 檢驗來比較高風險組和低風險組患者間的生存差異。

1.6 GBM總生存期個體化預測列線圖構建

首先對臨床、病理、基因信息進行單因素Cox 回歸，得到與OS 顯著相關的風險因素。將Risk-score與臨床、病理、基因信息相結合，采用多因素Cox 回歸構建臨床-影像組學聯(lián)合預測模型及列線圖，同時構建臨床預測模型，采用一致性指數(shù)（C-index）評估預后預測模型的判別能力并進行模型效能比較，采用校正曲線來評估列線圖預測生存概率和實際生存概率之間的一致性。本研究影像組學分析流程見圖2。

圖2 影像組學分析流程圖。T1CE：對比增強后T1 加權；FLAIR：液體衰減反轉恢復；ROI：感興趣區(qū)；PCA：主成分分析；LASSO：最小絕對收縮和選擇算法；KM：Kaplan-Meier；MGMT：O6-甲基鳥嘌呤-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: O6-methylguanine-DNA methyltransferase.

1.7 統(tǒng)計學分析

本研究中所有統(tǒng)計分析均在R 4.2.2版本（https://www.R-project.org/）進行。采用“Survival”包進行單變量和多變量Cox 回歸分析和KM 生存分析。分別使用“glmnet”和“rms”包構建LASSO-Cox 回歸模型和繪制列線圖。所有統(tǒng)計學檢驗均為雙側檢驗，以P＜0.05為差異有統(tǒng)計學意義。

2 結果

2.1 臨床風險因素和臨床模型

訓練集和測試集患者的年齡分別為（61.38±12.44）歲和（59.75±13.81）歲；男性分別占58.33%和58.06%。關于基因狀態(tài)，兩組GBM 中MGMT的描述性統(tǒng)計結果以甲基化MGMT 居多。年齡（P=0.307）、性別（P＞0.999）、MGMT（P=0.564）差異均無統(tǒng)計學意義。采用單因素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 個主成分特征，通過所選主成分的各自系數(shù)（β）加權線性組合計算風險評分：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）。

圖3 經(jīng)過主成分分析降維后主成分數(shù)量與覆蓋變異性比例的關系，選擇覆蓋95%以上變異性的41 個主成分，隨著選擇的主成分數(shù)量大于41，對變異性的影響不顯著。 圖4 最小絕對收縮和選擇算法回歸模型的調(diào)優(yōu)參數(shù)（lambda）篩選。相對于log(lambda)生成部分似然偏差，最小的部分似然偏差對應于最優(yōu)的特征數(shù)量。虛線表示根據(jù)1個標準誤差標準和最小標準得到的最佳lambda值。當剩余16個主成分時，部分似然偏差最小。 圖5 基于Risk-score對患者劃分高、低風險組，訓練集（5A）：P＜0.0001，測試集（5B）：P＜0.0001，二者均＜0.05，差異具有統(tǒng)計學意義。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.

利用R中“Surv_point”函數(shù)使用最大統(tǒng)計學差異選擇最佳Cut-off 值，在訓練集和測試集中，OS 均可以通過Risk-score 得到顯著風險分層（圖5）。依據(jù)Risk-score建立影像組學預測模型，采用C-index評估模型預測性能（訓練集：C-index=0.744，測試集：C-index=0.710）。

2.3 臨床-影像組學聯(lián)合預測模型及列線圖的構建

利用單因素Cox回歸和多因素Cox回歸，最終選擇風險評分、年齡、MGMT 狀態(tài)構建臨床-影像組學聯(lián)合預測模型和列線圖預測GBM 患者的6 月、12 月和24 月的生存概率（表2），發(fā)現(xiàn)臨床-影像組學聯(lián)合預測模型性能（訓練集：C-index=0.768，測試集：C-index=0.724）顯著高于臨床預測模型（訓練集：C-index=0.659，測試集：C-index=0.653)。列線圖在訓練集和測試集中均顯示出良好的校正效果（圖6）。此外，決策曲線分析顯示，在合理閾值概率的大部分范圍內(nèi)，臨床-影像組學聯(lián)合預測模型總體凈收益高于影像組學預測模型和臨床預測模型（圖7）。

圖6 臨床-影像組學聯(lián)合預測模型和列線圖。6A：基于訓練集結合Risk-score、患者年齡和O6-甲基鳥嘌呤-DNA 甲基轉移酶（MGMT）來開發(fā)的列線圖；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 O6-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-18]，基于高通量影像組學特征構建了基于T1CE 和T2 FLAIR 聯(lián)合序列的GBM 總生存期預測模型，并與臨床病理信息結合，模型效能得到了極大的改善。本研究在GBM 層面，使用常規(guī)MRI影像組學方法，采用PCA降維建立Risk-score 進而構建GBM 總生存期預測模型為國內(nèi)外被首次提出。該模型的構建將為GBM 患者的個體化診療水平和總生存期預測準確性產(chǎn)生積極影響，切實提升患者生存獲益。

3.1 MRI序列和ROI分割區(qū)域的選擇

在研究序列選擇方面，有研究表明[19]，多序列聯(lián)合模型效能普遍高于單序列模型效能，單序列模型受限于其單一視角的數(shù)據(jù)，如果該視角的數(shù)據(jù)受到噪聲干擾，可能會導致模型性能下降。而多序列模型通過融合多個視角的信息，能夠更好地抵抗噪聲干擾，同時，多序列聯(lián)合模型通過整合不同序列的特征，能夠建立起一個更具魯棒性的模型，這對于解決GBM 異質(zhì)性問題非常關鍵。TAN 等[20]聯(lián)合T1CE、T2 FLAIR 序列，針對腫瘤區(qū)域和瘤周水腫區(qū)域研究表明，影像組學特征可以獨立預測患者的OS，其中臨床-影像組學聯(lián)合模型性能（訓練集：C-index=0.764，測試集：C-index=0.758）顯著高于影像組學模型（訓練集：C-index=0.707，測試集：C-index=0.711）。在ROI 分割區(qū)域選擇方面，YANG 等[21]對GBM 瘤周水腫亞區(qū)提取影像組學特征建立預測模型，并結合臨床風險因素構建聯(lián)合模型，同樣聯(lián)合模型展現(xiàn)出最好的預測能力（訓練集：C-index=0.770，測試集：C-index=0.785）。與本研究相比，該研究聯(lián)合T1CE、T2 FLAIR 序列，將MR 圖像更加細致地分割為增強腫瘤區(qū)域、水腫區(qū)域和壞死區(qū)域，提取了更加全面的影像組學特征，最終建立的影像組學模型具備良好的預測性能（訓練集：C-index=0.744，測試集：C-index=0.710），結合臨床風險因素所構建的聯(lián)合模型預測效能在所建立的模型中表現(xiàn)最優(yōu)（訓練集：C-index=0.768，測試集：C-index=0.724）。

3.2 特征篩選方法

在特征篩選過程中，高維數(shù)據(jù)通常包含大量不相關的、冗余的和有噪聲的特征，這可能導致維度的破壞和模型的過擬合[22]，自動編碼器，稀疏表示和局部保持投影相結合，Kruskal-Wallis 檢驗、Variance Threshold等方法在降低特征維數(shù)方面產(chǎn)生了不錯的效果[22-25]。本研究針對每個患者分別提取了10 128個特征，特征維數(shù)過高，為此本研究采用了PCA，以消除特征間的多重共線性，提高模型的解釋能力和預測精度。PCA 通過尋找數(shù)據(jù)中的最大方差方向，并在這些方向上創(chuàng)建新的、線性無關的變量（主成分），從而有效地簡化了數(shù)據(jù)集的復雜性。這些新的主要組成成分保留了原始數(shù)據(jù)的大部分信息，但具有更低的維度，使得數(shù)據(jù)可以在較低的計算成本下更易于處理和解釋。PCA 的應用使我們能夠在減少計算復雜性和保證模型性能之間找到一個理想的平衡點，顯示出強大的降維能力和數(shù)據(jù)解釋能力。經(jīng)過PCA 降維之后構建的模型效能得到顯著提高，同時過擬合問題得到了一定的解決。

3.3 風險分層依據(jù)及模型構建方法

本研究根據(jù)風險評分進行分層，在之前的一些研究中，有許多依據(jù)風險評分中位值或中值進行分層的例子[26-28]，本研究采用“R”包中“Surv_point”函數(shù)通過比較每個可能的切點下兩組之間的生存差異來確定最佳切點。對分層后的訓練集和測試集通過KM 曲線和Log-rank 檢驗，P均＜0.05，差異具有統(tǒng)計學意義，并且低風險組的生存概率顯著高于高風險組。在模型構建方面，一些研究嘗試使用隨機森林和其他機器學習分類器來預測患者生存期[29-30]，生存分析被作為一個二分類任務，將生存分析作為二分類任務的主要局限性在于忽視了時間因素的作用。二分類任務基于特定時間點進行預測，但在實際的臨床生存分析中，患者的生存時間是一個連續(xù)的變量，其分布可能受到年齡，性別，疾病狀態(tài)等因素的影響，忽視這些因素可能會導致模型的預測能力下降。本研究利用Cox 回歸進行OS 模型構建，將風險變量作為連續(xù)變量，開發(fā)了一種不依賴任何閾值的模型。

3.4 本研究的局限性

本研究存在一些局限性。首先，數(shù)據(jù)來自不同的數(shù)據(jù)中心，MR 儀器的型號、參數(shù)存在一定差異，這可能對結果造成一定的誤差。其次，由于對每位患者提取了包括一階、紋理、形狀與大小以及8 種濾波器方式的海量特征，在特征選擇過程中進行了PCA降維，其局限性在于經(jīng)過PCA 降維后的主成分是若干原始特征的線性組合，這讓對其生物學解釋變得困難。同時，對于多中心數(shù)據(jù)，需要依據(jù)OS、臨床、病理信息完整的入組標準，導致一些信息不能被選用，如KPS 評分、IDH 狀態(tài)等。有研究表明，將影像組學聯(lián)合更多的臨床、病理、基因信息對模型效能的提高有顯著的影響[31-32]。在特征提取過程中，使用傳統(tǒng)方法提取特征的效率低于深度學習方法[33-35]。在下一項目的研究中，團隊將在保持數(shù)據(jù)量充足的同時收集豐富的臨床、病理、遺傳信息，并進行學習方法的遷移和組合，力求解決本研究存在的限制。

4 結論

綜上所述，本研究通過基于多序列和多區(qū)域的MRI 影像組學分析，聯(lián)合臨床病理信息開發(fā)GBM 總生存期預測模型具有較好的預測性能，通過這一模型，臨床醫(yī)生能夠更有效精準地進行GBM 患者的風險分層評估，從而對患者進行個體化的治療方案制訂。

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

作者貢獻聲明：楊國強設計了本研究的方案，對稿件的重要內(nèi)容進行了修改；牛文舉起草和撰寫稿件，獲取、分析或解釋本研究的數(shù)據(jù)；徐懷文、高宇翔、王效春、譚艷、張輝獲取、分析或解釋本研究的數(shù)據(jù)，對稿件重要內(nèi)容進行了修改，王效春、張輝獲得了國家自然科學基金項目的資助，楊國強獲得了山西省基礎研究計劃面上項目的資助；全體作者都同意發(fā)表最后的修改稿，同意對本研究的所有方面負責，確保本研究的準確性和誠信。