基于YAP流水線的影像組學(xué)研究流程的實(shí)現(xiàn)

李冬寶，宋陽，羅慶，謝海濱，2*，楊光，2*

作者單位：

1. 華東師范大學(xué)物理與材料科學(xué)學(xué)院，上海市磁共振重點(diǎn)實(shí)驗(yàn)室，上海200062

2. 上?？颠_(dá)卡勒幅醫(yī)療科技有限公司，上海 200062

醫(yī)療影像是臨床診斷中的重要工具，可以無創(chuàng)地提供患者信息來輔助臨床醫(yī)師進(jìn)行診斷。影像組學(xué)是[1-2]是一種常見的將影像和臨床診斷相關(guān)聯(lián)的方法，通過對大量的醫(yī)學(xué)圖像進(jìn)行分析得到隱含的圖像特征數(shù)據(jù)，并對特征數(shù)據(jù)與臨床診斷進(jìn)行關(guān)聯(lián)分析，從而建立從圖像到臨床診斷的模型。該模型可以用于分析醫(yī)學(xué)影像，輔助醫(yī)師進(jìn)行臨床診斷、術(shù)后干預(yù)等。相較于活體組織檢查，它不會對人體產(chǎn)生身體或心理上的傷害。磁共振成像作為臨床中常用的檢測方法，多角度、多參數(shù)、多序列成像，具有較高的軟組織對比度，不會造成放射性和侵入性的傷害，可以提供豐富的臨床特征，經(jīng)常用于影像組學(xué)研究[3-5]。

又一條流水線(yet another pipeline，YAP)是一個醫(yī)學(xué)圖像重建與后處理的開發(fā)框架[6]，采用基于接口的設(shè)計(jì)，可以通過插件對系統(tǒng)進(jìn)行擴(kuò)展。YAP具有以下特點(diǎn)：(1)處理模塊支持多個輸入輸出端口，流水線支持分支結(jié)構(gòu)；(2)提供了參數(shù)管理，可以向處理模塊傳遞伴隨數(shù)據(jù)的參數(shù)，也可以在模塊間傳遞參數(shù)；(3)支持靈活多樣的數(shù)據(jù)組織形式。此外，YAP還提供了基本的核心部件、基礎(chǔ)磁共振圖像處理插件庫、接口實(shí)現(xiàn)和使用的幫助類等。

由于影像組學(xué)很多開發(fā)工具是基于Python語言開發(fā)的，為了方便在YAP中搭建影像組學(xué)流水線，本研究為YAP增加了Python語言的支持，允許使用Python編寫流水線中的處理模塊，方便了算法的研究；同時利用基于Python的影像組學(xué)軟件包van Griethuysen等[7]在YAP上實(shí)現(xiàn)了影像組學(xué)處理流水線。本研究的最后一部分，使用影像組學(xué)流水線進(jìn)行多模態(tài)磁共振圖像腫瘤病理分級預(yù)測的方法。

1 影像組學(xué)流水線實(shí)現(xiàn)

1.1 影像組學(xué)研究流程

自Lambin等[1]和Kumar等[2]提出影像組學(xué)概念以來，這種方法越來越多被用來研究腫瘤病理分級，臨床診斷以及治療反饋。影像組學(xué)進(jìn)行分類的流程如圖1所示，其過程大致包括：獲取圖像、預(yù)處理、圖像分割、特征提取、特征選擇、分類以及結(jié)果分析。

其中，影像組學(xué)最常用圖像包括CT、PET以及MRI[8-9]。預(yù)處理過程常常包含圖像配準(zhǔn)、歸一化等過程，使得用于組學(xué)分析的圖像盡可能一致。圖像分割步驟可以采用自動分割或者手動分割，不同的分割算法會影響影像組學(xué)的結(jié)果[8]。用于研究的金標(biāo)準(zhǔn)則常常采用手工分割的方法，高質(zhì)量自動分割算法的匱乏可能是妨礙影像組學(xué)應(yīng)用的主要障礙之一。

影像組學(xué)分析中，可以使用很多不同類型的特征，包括信號強(qiáng)度、幾何形狀、紋理等不同類別的特征。信號強(qiáng)度特征是基于圖像中像素點(diǎn)計(jì)算的統(tǒng)計(jì)值。幾何形狀特征是指腫瘤的形狀以及腫瘤的大小。紋理特征又稱為高階特征，是對圖像中多個像素灰度關(guān)系進(jìn)行統(tǒng)計(jì)，Kumar等[2]研究提到紋理特征會比信號強(qiáng)度特征有更好的潛在價值。本研究使用的紋理特征包括Gray-Level Co-occurrence matrix (GLCM)[10]特征、Gray-Level Size Zone Matrix (GLSZM)[11]特征、Gray-Level Run Length Matrix (GLRLM)[12]特征、Gray Level Dependence Matrix (GLDM)[13]特征和Neighboring Gray Tone Difference Matrix (NGTDM)[14]特征等。此外，可能影響模型的臨床或個人信息也可以加入特征空間。

雖然使用更多的特征可能會導(dǎo)致更準(zhǔn)確的模型[15]，但過多的特征容易導(dǎo)致過擬合現(xiàn)象，因此需要通過特征選擇來去除冗余的或與分類結(jié)果無關(guān)的特征，減少分析所用的特征數(shù)量。特征選擇方法有很多，例如Pearson相關(guān)系數(shù)[16]、遺傳算法[17]以及遞歸特征消除方法[18]等。影像組學(xué)分析的最后一步是建立特征到臨床結(jié)論的預(yù)測模型，支持向量機(jī)(support vector machine，SVM)[19]就是最常用的分類器之一。

1.2 影像組學(xué)流水線的總體設(shè)計(jì)

由于在科研工作中，圖像分割常常由有經(jīng)驗(yàn)的影像科醫(yī)師手工完成，本研究實(shí)現(xiàn)的影像組學(xué)流水線假定已經(jīng)完成了圖像預(yù)處理與圖像分割工作。由于不同數(shù)據(jù)來源的原始數(shù)據(jù)組織形式不盡相同，為了使用本研究的流水線處理，應(yīng)使用手工或者自動的方式將圖像數(shù)據(jù)轉(zhuǎn)換成如圖2所示的文件目錄結(jié)構(gòu)中。

影像組學(xué)流水線的結(jié)構(gòu)如圖3所示。其中，SampleIterator表示樣本迭代器，它將數(shù)據(jù)文件夾中的每個樣本的文件夾一次送入ModalityIterator處理器處理。ModalityIterator表示模態(tài)迭代器，它將每個樣本中不同模態(tài)(序列)的圖像送入后續(xù)處理器處理。NiiReader處理器，實(shí)現(xiàn)對*.nii格式文件的讀取，它讀取的數(shù)據(jù)被送入FeatureExtractor進(jìn)行特征提取，提取的特征依次有ModalityCollector、SampleCollector收集到樣本特征矩陣。最后樣本特征矩陣匯總后，送入Classifier進(jìn)行特征選擇、分類等分析處理。實(shí)現(xiàn)上述流程的YAP腳本代碼如圖4所示。

其中第1行導(dǎo)入PythonRadiomics.dll模塊，使得腳本可以使用該模塊提供的處理器。影像組學(xué)流水線中所使用的所有處理器，都在此模塊中實(shí)現(xiàn)。第3～13行，定義了流水線中要使用的處理器對象。其語法類似C++中的對象定義，在構(gòu)造處理器對象的同時，可以指定處理器的屬性。

圖1 影像組學(xué)圖像分類處理流程Fig. 1 Flowchart of radiomics-based classification.

圖2 流水線要求的數(shù)據(jù)集的存儲方式。圖中Sample1、Sample2、…等表示不同的樣本文件夾；data1.nii、data2.nii、…表示不同模態(tài)的圖像；ROI.nii表示分割圖像；Label.xlsx表示當(dāng)前樣本的標(biāo)記數(shù)據(jù)Fig. 2 Required structure of folders containing dataset images. Sample1,Sample2, … means different sample folders. data1.nii, data2.nii, … means different modality images. ROI.nii means a segmentation image. Label.xlsx means current sample label data.

第15～22行，使用“-＞”操作符將處理器對象相互連接，在處理器對象Id后面可以指定具體的輸入或輸出端口Id；如果不指定端口，那么就使用默認(rèn)輸入“Input”或輸出“Output”。

第24行指定將整個流水線的輸入端口設(shè)置為sample處理器的Input端口。運(yùn)行時，直接向流水線中饋送數(shù)據(jù)對象即可啟動流水線處理。由于sample處理器可以從其屬性獲得待處理的文件夾路徑，所以用于啟動流水線工作的數(shù)據(jù)對象可以是個空數(shù)據(jù)對象。

1.3 Python支持

如前文所述，本研究影像組學(xué)流水線所用的所有處理器都在PythonRadiomics.dll中提供。在實(shí)現(xiàn)FeatureExtractor處理器時，筆者使用了PyRadiomics軟件包提供的功能，這是一個由第三方提供的用Python語言編寫的影像組學(xué)處理軟件包。

Python是一種面向?qū)ο蟮慕忉屝驼Z言，它具有簡單易學(xué)、運(yùn)行速度快、免費(fèi)開源、良好移植性、面向?qū)ο筇匦?、擴(kuò)展性和豐富的第三方庫等諸多優(yōu)點(diǎn)。它可以方便地調(diào)用C/C++、Java等語言編寫的模塊，所以也有“膠水語言”稱號。Python在科學(xué)計(jì)算方面有強(qiáng)大的第三方庫，例如：Numpy、scikit-learn等，在大數(shù)據(jù)與人工智能火熱的背景下，Python得到了迅速的發(fā)展，在IEEE Spectrum2017年發(fā)布的編程語言排行榜上高居榜首[20]，絕大多數(shù)的機(jī)器學(xué)習(xí)平臺也都支持Python編程。

圖3 影像組學(xué)流水線在YAP中的流程圖。圖中使用的模塊包括樣本迭代器(SampleIterator)、模態(tài)迭代器(ModalityIterator)、Nii文件讀取器(NiiReader)、特征提取器(FeatureExtractor)、模態(tài)收集器(ModalityCollector)、樣本收集器(SampleCollector)以及分類處理器(Classifier)Fig. 3 Flowchart of radiomics pipeline in YAP. Modules used in the flowchart includes SampleIterator, ModalityIterator, NiiReader,FeatureExtractor, ModalityCollector, SampleCollector and Classifier processor.

圖4 實(shí)現(xiàn)的影像組學(xué)流程的YAP腳本代碼 圖5 Python服務(wù)中提供的IPython接口函數(shù)簽名Fig. 4 Implementation of Radiomics flowchart’s YAP script code.Fig. 5 IPython Interface function signature in the Python service.

本研究在YAP中增加了對Python語言的支持，允許在編寫處理器時使用Python腳本，這樣就簡化了一些算法的研究工作，一些對于速度要求不高的處理器，也可以使用Python編寫，以提升開發(fā)效率。由于在同一進(jìn)程內(nèi)無法使用多個Python解釋器，更不用說同時使用不同版本的Python解釋器，所以在YAP的系統(tǒng)服務(wù)中，以IPython接口的形式，將調(diào)用Python腳本進(jìn)行數(shù)據(jù)處理的能力提供給處理器開發(fā)人員。需要指出的是，目前筆者在IPython中只提供了專門用于YAP流水線數(shù)據(jù)處理的高級接口，并不能調(diào)用Python的所有功能，相應(yīng)代碼如圖5所示。

在流水線中，調(diào)用Python處理器的流程如圖6所示。圖中虛線框內(nèi)表示Python引擎內(nèi)部運(yùn)行的大致過程，Python引擎將數(shù)據(jù)轉(zhuǎn)換為Python識別的數(shù)據(jù)類型，然后調(diào)用Python解析器(Python Interpreter)對腳本文件Radiomics.py(PythonProcessor的Module屬性制定，并通過IPython::Process()函數(shù)的參數(shù)傳遞)進(jìn)行解析，最后得到結(jié)果轉(zhuǎn)換得到C++形式數(shù)據(jù)類型，作為IPython::Process()函數(shù)的輸出，由PythonProcessor處理器傳往下一處理器。需要指出的是，Python引擎會在首次使用時啟動Python 解釋器。

在YAP中增加了流水線支持后，就可以使用Python語言編寫處理器，實(shí)現(xiàn)特定的處理步驟。Python的快速開發(fā)能力和強(qiáng)大的科學(xué)計(jì)算、可視化相關(guān)的第三方庫的支持，能給利用YAP進(jìn)行算法研究的人員帶來極大的便利。本研究工作中的FeatureExtractor利用了PyRadiomics中的腳本進(jìn)行特征提取；Classifier處理器，同樣利用Python腳本進(jìn)行影像組學(xué)的特征選擇和分類工作。

另一方面，我們的影像組學(xué)流水線中采用了醫(yī)學(xué)圖像常用的NIfTI格式，NiiReader處理器用于對該格式文件的支持，目前僅支持NIfTI-1和NIfTI-2版本的基本類型數(shù)據(jù)的解析，對＜*.img,*.hdr＞形式的文件對解析不支持。筆者用C++編寫了該處理器。不同語言編寫的處理器，可以在同一條流水線中應(yīng)用，體現(xiàn)了YAP的靈活性。

2 應(yīng)用影像組學(xué)流水線對膠質(zhì)瘤分級

2.1 使用的數(shù)據(jù)

為驗(yàn)證本研究影像組學(xué)流水線的可用性，筆者使用本研究的影像組學(xué)流水線對BRATS2017[21]公開數(shù)據(jù)集中的三維多模態(tài)膠質(zhì)瘤數(shù)據(jù)進(jìn)行了影像組學(xué)研究。該數(shù)據(jù)集中的成像序列包括FLAIR、T1CE、T1和T2。數(shù)據(jù)經(jīng)過標(biāo)準(zhǔn)化處理，提供了分割好的多種腫瘤感興趣區(qū)域，包括壞死或非增強(qiáng)部分、腫瘤周圍水腫部分和腫瘤增強(qiáng)區(qū)域。在本研究中，筆者將以標(biāo)記區(qū)域合并作為圖像感興趣區(qū)域。

圖6 Python引擎處理數(shù)據(jù)流程。PythonProcessor：調(diào)用Python服務(wù)的處理模塊，PythonProcessor中的extract表示腳本中被調(diào)用的函數(shù)；Python Interpreter：Python解析器；radiomics.py：被調(diào)用的Python腳本文件Fig. 6 Python engine processing data flow. PythonProcessor refers to a processor which calls Python service.Python Interpreter refers to Python script parser. Radiomics.py refers to Python script file called. ‘extract’ in the PythonProcessor refers to a called function in the ‘radiomics.py’ script.

圖7 BRATS2017 腫瘤分級的結(jié)果。A：不同的特征數(shù)量得到的模型準(zhǔn)確率，選擇特征數(shù)為12時，可以獲得最佳的平均準(zhǔn)確率94.5%。B：ROC曲線，其中驗(yàn)證集的AUC值為0.9650Fig. 7 Results for classification of BRATS2017 dataset. A: Model accuracy using different features number.The highest average accuracy of 94.5% was achieved when the number of selected features is 12. B: Receiver operator characteristic curve of the models. AUC for the validation data is 0.9650.

此外，BRATS2017也提供了腫瘤的分級標(biāo)記，即高級別膠質(zhì)瘤(higher grade glioma，HGG)和低級別膠質(zhì)瘤(lower grade glioma，LGG)標(biāo)記。數(shù)據(jù)集中包括72例LGG、163例HGG，共235個樣本。按照流水線的解析要求，把數(shù)據(jù)預(yù)先整理如圖2所示的形式，然后將路徑設(shè)置給SampleIterator處理器的Path屬性。

2.2 實(shí)驗(yàn)方法

本實(shí)驗(yàn)的特征提取利用了PyRadiomics工具包中的腳本，對圖像進(jìn)行預(yù)處理和特征提取。結(jié)合不同的模態(tài)，提取的可量化特征一共有383種。采用RFE (recursive feature elimination)[18]進(jìn)行特征篩選，采用線性SVM進(jìn)行模型建立，為了獲得較好的模型結(jié)果，對線性SVM分類器的分類懲罰參數(shù)C進(jìn)行基本范圍的遍歷調(diào)節(jié)，C值分別取0.001、0.005、0.01、0.05、0.1、0.25、0.5、1、2。10-Fold交叉驗(yàn)證(cross-validation，CV)進(jìn)行模型評估。評估方法包括準(zhǔn)確率、受試者操作特性(receiver operating characteristic，ROC)曲線和曲線下面積(area under curve，AUC)值。

2.3 結(jié)果和討論

使用RFE方法對特征排序后，觀察10次交叉驗(yàn)證實(shí)驗(yàn)的平均準(zhǔn)確率隨所選的特征數(shù)變化的情況，最終得到平均準(zhǔn)確率最高的情況出現(xiàn)在選擇特征數(shù)為12時，如圖7A所示。其中，得出當(dāng)C=0.25時，10次交叉驗(yàn)證的平均準(zhǔn)確率可以達(dá)到94.5%。除準(zhǔn)確率外，還使用ROC曲線的AUC值作為判定模型好壞。在圖7B中，表示訓(xùn)練和驗(yàn)證的ROC曲線以及對應(yīng)的AUC值為0.9650。

由實(shí)驗(yàn)過程和結(jié)果可見：通過YAP流水線，可以實(shí)現(xiàn)對影像組學(xué)處理過程的應(yīng)用。得到良好的分級結(jié)果，可以為臨床醫(yī)師提供很大幫助，也為該類型的疾病預(yù)測和治療分析提供幫助。

在影像組學(xué)流水線上，不同數(shù)據(jù)的研究，僅需要通過改變流水線腳本里Python處理器里的參數(shù)，例如圖6中的Module和Method等即可，而無需將項(xiàng)目進(jìn)行重新編譯。而且在流水線中，也可以很方便地添加一些其他處理模塊，例如利用YAP提供的基礎(chǔ)設(shè)施，可以方便地進(jìn)行特征的保存、可視化等操作，也可以將計(jì)算過程中發(fā)生的錯誤或關(guān)鍵數(shù)據(jù)記錄到日志。這些對算法開發(fā)人員和圖像研究人員都提供很大的便利。在影像組學(xué)流水線上，可以方便地開發(fā)并驗(yàn)證自己的算法，這為算法研究、醫(yī)學(xué)影像研究以及更大規(guī)模的應(yīng)用提供了便利。

3 討論

本研究在YAP框架中增加了對Python語言的支持，并在此基礎(chǔ)上實(shí)現(xiàn)了影像組學(xué)流水線。該流水線充分利用了Python語言中豐富的科學(xué)計(jì)算軟件包，也可以利用C++的速度和YAP流水線提供的基礎(chǔ)服務(wù)?；旌险Z言開發(fā)醫(yī)學(xué)影像后處理流水線，既可以利用Python的快速開發(fā)進(jìn)行算法研究，也可以在需要時，將Python處理器用C++處理器替代，以提高性能。我們利用影像組學(xué)流水線對BRATS2017公開數(shù)據(jù)集進(jìn)行了腦部膠質(zhì)瘤的影像組學(xué)的研究，得到較好的分類結(jié)果，這進(jìn)一步確認(rèn)使用YAP流水線形式進(jìn)行影像組學(xué)研究是可行的。