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

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

作者單位：

1. 華東師范大學物理與材料科學學院，上海市磁共振重點實驗室，上海200062

2. 上海康達卡勒幅醫(yī)療科技有限公司，上海 200062

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

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

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

1 影像組學流水線實現(xiàn)

1.1 影像組學研究流程

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

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

影像組學分析中，可以使用很多不同類型的特征，包括信號強度、幾何形狀、紋理等不同類別的特征。信號強度特征是基于圖像中像素點計算的統(tǒng)計值。幾何形狀特征是指腫瘤的形狀以及腫瘤的大小。紋理特征又稱為高階特征，是對圖像中多個像素灰度關系進行統(tǒng)計，Kumar等[2]研究提到紋理特征會比信號強度特征有更好的潛在價值。本研究使用的紋理特征包括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]特征等。此外，可能影響模型的臨床或個人信息也可以加入特征空間。

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

1.2 影像組學流水線的總體設計

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

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

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

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

圖2 流水線要求的數(shù)據(jù)集的存儲方式。圖中Sample1、Sample2、…等表示不同的樣本文件夾；data1.nii、data2.nii、…表示不同模態(tài)的圖像；ROI.nii表示分割圖像；Label.xlsx表示當前樣本的標記數(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；如果不指定端口，那么就使用默認輸入“Input”或輸出“Output”。

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

1.3 Python支持

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

Python是一種面向對象的解釋型語言，它具有簡單易學、運行速度快、免費開源、良好移植性、面向對象特性、擴展性和豐富的第三方庫等諸多優(yōu)點。它可以方便地調用C/C++、Java等語言編寫的模塊，所以也有“膠水語言”稱號。Python在科學計算方面有強大的第三方庫，例如：Numpy、scikit-learn等，在大數(shù)據(jù)與人工智能火熱的背景下，Python得到了迅速的發(fā)展，在IEEE Spectrum2017年發(fā)布的編程語言排行榜上高居榜首[20]，絕大多數(shù)的機器學習平臺也都支持Python編程。

圖3 影像組學流水線在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 實現(xiàn)的影像組學流程的YAP腳本代碼 圖5 Python服務中提供的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ā)效率。由于在同一進程內無法使用多個Python解釋器，更不用說同時使用不同版本的Python解釋器，所以在YAP的系統(tǒng)服務中，以IPython接口的形式，將調用Python腳本進行數(shù)據(jù)處理的能力提供給處理器開發(fā)人員。需要指出的是，目前筆者在IPython中只提供了專門用于YAP流水線數(shù)據(jù)處理的高級接口，并不能調用Python的所有功能，相應代碼如圖5所示。

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

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

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

2 應用影像組學流水線對膠質瘤分級

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

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

圖6 Python引擎處理數(shù)據(jù)流程。PythonProcessor：調用Python服務的處理模塊，PythonProcessor中的extract表示腳本中被調用的函數(shù)；Python Interpreter：Python解析器；radiomics.py：被調用的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 腫瘤分級的結果。A：不同的特征數(shù)量得到的模型準確率，選擇特征數(shù)為12時，可以獲得最佳的平均準確率94.5%。B：ROC曲線，其中驗證集的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也提供了腫瘤的分級標記，即高級別膠質瘤(higher grade glioma，HGG)和低級別膠質瘤(lower grade glioma，LGG)標記。數(shù)據(jù)集中包括72例LGG、163例HGG，共235個樣本。按照流水線的解析要求，把數(shù)據(jù)預先整理如圖2所示的形式，然后將路徑設置給SampleIterator處理器的Path屬性。

2.2 實驗方法

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

2.3 結果和討論

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

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

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

3 討論

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