金文東，陰慧娟，李迎新

據(jù)美國癌癥協(xié)會(huì)統(tǒng)計(jì)，癌癥的總死亡率在20世紀(jì)經(jīng)歷了增長后，從1991年至2018年下降了31%，這主要得益于腫瘤預(yù)防、早期發(fā)現(xiàn)和治療方法的改進(jìn)[1]。為進(jìn)一步了解癌癥關(guān)鍵轉(zhuǎn)變過程中腫瘤動(dòng)態(tài)生物系統(tǒng)中細(xì)胞間的相互作用，實(shí)現(xiàn)對癌癥的精確治療，美國國家癌癥研究所（National Cancer Institute，NCI）在癌癥登月計(jì)劃中提出建立包括腫瘤發(fā)生、局部擴(kuò)張、轉(zhuǎn)移和治療抗性等關(guān)鍵轉(zhuǎn)變的三維圖譜——人類腫瘤圖譜網(wǎng)絡(luò)（Human Tumor Atlas Network，HTAN）[2]。這項(xiàng)工作主要通過分子分析、空間分子分析、組織學(xué)分析和解剖學(xué)分析生成數(shù)據(jù)，不僅包括腫瘤樣本的收集，同時(shí)涵蓋細(xì)胞空間數(shù)據(jù)分析及可視化工具的開發(fā)。

免疫組織化學(xué)（immunohistochemistry，IHC）作為組織病理學(xué)分析的金標(biāo)準(zhǔn)，可以實(shí)現(xiàn)多種腫瘤標(biāo)志物的標(biāo)記。全切片圖像（whole slide images，WSIs）是通過圖像掃描技術(shù)獲得的全組織切片二維數(shù)字圖像，通過空間、時(shí)間維度的擴(kuò)展，可以實(shí)現(xiàn)單細(xì)胞圖譜的三維和四維繪制。目前，IHC分析圖像分割的金標(biāo)準(zhǔn)是顏色反卷積（colour deconvolution，CD）方法，常在ImageJ中實(shí)現(xiàn)[3]。針對WSIs這種大內(nèi)存圖像，Bankhead P等[4]開發(fā)了開源數(shù)字病理系統(tǒng)QuPath?；贑D的分割原理，ImageJ和QuPath均只能實(shí)現(xiàn)3種染色以下的明場圖像分析，而腫瘤生物學(xué)效應(yīng)常常具有多個(gè)腫瘤標(biāo)志物；同時(shí)，目前的可視化分析方法缺乏對腫瘤標(biāo)志物空間位置及組織與細(xì)胞間空間分布關(guān)系的分析。

Lab色彩空間是基于人類視覺原理，將顏色分成明度L*、紅綠a*和藍(lán)黃b*通道的一種色彩模式，比紅綠藍(lán)（red green blue，RGB）和青品紅黃黑（cyan magenta yellow black，CMYK）的色域更廣。K均值聚類是一種對目標(biāo)數(shù)據(jù)集劃分的無監(jiān)督學(xué)習(xí)方法，結(jié)合Lab色域可以模擬人眼對顏色分類?？臻g分布分析[5]是地理信息系統(tǒng)（geographic information system，GIS）中對實(shí)體或事件的空間定量研究方法，通過對空間數(shù)據(jù)和空間模型的聯(lián)合分析發(fā)掘空間對象的潛在信息，包括空間位置、分布、形態(tài)、形成、演變和空間對象間的空間關(guān)系等。據(jù)此，筆者提出對Lab色域的病理圖像進(jìn)行2次K均值聚類實(shí)現(xiàn)WSIs腫瘤生物標(biāo)志物的分割和閾值提取，結(jié)合空間分布分析方法實(shí)現(xiàn)量化和可視化分析。以小鼠乳腺癌光動(dòng)力治療（photodynamic therapy，PDT）為例，闡述其在腫瘤生物學(xué)效應(yīng)分析中的應(yīng)用，基于獲取的閾值對序列蘇木精-伊紅（hematoxylin-eosin，HE）、TdT介導(dǎo)的dUTP缺口末端標(biāo)記（TdT-mediated dUTP nick end labeling，TUNEL）和血小板內(nèi)皮細(xì)胞黏附因子（platelet endothelial cell adhesion molecule-1，PECAM-1/CD31）染色WSIs數(shù)據(jù)集進(jìn)行圖像分割，實(shí)現(xiàn)對腫瘤PDT后壞死、凋亡和血管等生物學(xué)效應(yīng)的空間分布分析。

1 材料與方法

1.1 實(shí)驗(yàn)材料

所用小鼠乳腺癌PDT3d后的HE、CD31和TUNEL染色WSIs數(shù)據(jù)集來自科學(xué)數(shù)據(jù)銀行（Science Data Bank，ScienceDB）[6]數(shù)據(jù)庫，圖像分辨率為0.5μm/pixel；由中國醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)工程研究所激光醫(yī)學(xué)實(shí)驗(yàn)室制作和采集。圖像閾值提取、分割和空間分布分析均在MATLAB R2020a中實(shí)現(xiàn)。

1.2 方法

1.2.1 閾值提取和圖像分割

實(shí)驗(yàn)算法（Kmeans-Lab）以Lab色彩空間K均值聚類實(shí)現(xiàn)圖像分割。閾值提取流程如圖1a所示，隨機(jī)選取3幅陽性WSIs，每幅截取3個(gè)1 000×1 000像素的陽性圖像拼接成一個(gè)3×3的陽性圖像（訓(xùn)練圖像），用K均值聚類進(jìn)行圖像分割后提取模板用于整張WSIs的圖像分割。HE、CD31和TUNEL染色切片均由3種顏色組分構(gòu)成：HE明場背景、紅染細(xì)胞質(zhì)和藍(lán)染細(xì)胞核；CD31明場背景、棕染細(xì)胞膜和藍(lán)染細(xì)胞核；TUNEL明場背景、棕染細(xì)胞核和藍(lán)染細(xì)胞核。K均值聚類均設(shè)置k=3，基于a*和b*兩個(gè)通道聚類后散點(diǎn)圖的邊界繪制a與b的關(guān)系曲線，作為WSIs的顏色閾值。L*通道上將聚類圖像分為背景、強(qiáng)陽性、弱陽性和陰性4類，2次K均值聚類獲取亮度閾值。

1.2.2 空間分布分析

如圖1b，基于顏色和亮度閾值制作掩膜分割WSIs，中值濾波（3×3像素）去除“椒鹽”噪聲并進(jìn)行空間分布分析。將被標(biāo)記的分子作為點(diǎn)對象，以其面積占腫瘤面積的百分比進(jìn)行量化分析，定義腫瘤生物效應(yīng)的空間分布密度。核密度估計(jì)[7]（kernel density estimation，KDE）是點(diǎn)模式最常用的一階效應(yīng)分析方法之一，其中二維高斯核密度函數(shù)具有旋轉(zhuǎn)對稱性，而且權(quán)重隨著中心點(diǎn)位置單調(diào)低減。考慮到細(xì)胞間（腫瘤細(xì)胞大小約為10μm）的相互作用，以21×21像素大小的高斯濾波器對目標(biāo)區(qū)域進(jìn)行卷積，并生成熱力圖，實(shí)現(xiàn)可視化分析。

圖1 WSIs圖像處理流程圖Fig.1 Image processing flow chart of WSIs

1.2.3 圖像分割性能驗(yàn)證

與訓(xùn)練圖像相同，從陽性WSIs中截取共9個(gè)1 000×1 000像素的陽性圖像，拼接成一個(gè)3×3的陽性圖像（測試圖像），用于圖像分割性能的驗(yàn)證，對比分析Kmeans-Lab與下列常用圖像分割方法的性能參數(shù)，包括歸一化互信息（normal mutual information，NMI）、Kappa系數(shù)、平均交互比（mean intersection over union，mIoU）、平均精準(zhǔn)率（mean precision，mPr）、平均召回率（mean recall，mRe）和平均準(zhǔn)確度（mean accuracy，mA）[8，9]。

1.2.3.1 基于色域閾值的圖像分割 基于色域閾值的圖像分割方法是常用的彩色圖像分割方法之一，目前比較常用的色域有RGB、色調(diào)飽和度明度（hue saturation value，HSV）、明度紅藍(lán)色度（luminance chroma-red chroma-blue，YCrCb）和Lab等。MATLAB的Color Thresholder應(yīng)用，通過手動(dòng)選取感興趣區(qū)域（region of interest，ROI）可以實(shí)現(xiàn)不同色域閾值的半自動(dòng)圖像分割。其中RGB和HSV色域無亮度通道，蘇木精著色區(qū)域與其他染料分割效果不佳，故僅比對MATLAB下Lab和YCrCb色域的閾值分割方法，分別命名為MATLAB-Lab和MATLAB-YCrCb。

1.2.3.2 基于機(jī)器學(xué)習(xí)和統(tǒng)計(jì)模型的圖像分割 Fiji中的IHC Toolbox插件[10]是目前較為常用的IHC分析工具，其原理是通過對訓(xùn)練圖像集中的ROI的機(jī)器學(xué)習(xí)，建立基于YCrCb色域的統(tǒng)計(jì)模型用于圖像分割，這里將其命名為Fiji-YCrCb。

1.2.3.3 基于CD的圖像分割 CD是目前IHC圖像顏色分割中應(yīng)用最為廣泛的算法，其原理是基于朗伯比爾定律，針對不同染料RGB分量光的特異性吸收，將明場圖像RGB通道轉(zhuǎn)換為最多3個(gè)代表單個(gè)染料光密度的通道。QuPath軟件以訓(xùn)練圖像中的不同染料著色區(qū)域計(jì)算各自的著色向量，用于測試圖像的CD，得到表征蘇木精、曙紅或二氨基聯(lián)苯胺（3，3’-diaminobenzidine，DAB）通道的圖像。Otsu算法[11]是圖像處理中常用的一種灰度圖像目標(biāo)和背景的自適應(yīng)分割方法，與CD聯(lián)合運(yùn)用實(shí)現(xiàn)不同染料著色區(qū)域的分割。這里將通過QuPath軟件中的CD結(jié)合Otsu實(shí)現(xiàn)圖像分割的方法命名為QuPath-CD。

1.3 統(tǒng)計(jì)學(xué)方法

采用GraphPad Prism軟件進(jìn)行Pearson相關(guān)性分析和線性回歸分析。Pearson相關(guān)系數(shù)為負(fù)數(shù)時(shí)，表明變量存在負(fù)相關(guān)關(guān)系；正數(shù)時(shí)存在正相關(guān)關(guān)系；絕對值越接近1相關(guān)性越強(qiáng)，絕對值大于0.5時(shí)具有強(qiáng)相關(guān)性。P<0.05為差異有統(tǒng)計(jì)學(xué)意義，P<0.01為差異有顯著統(tǒng)計(jì)學(xué)意義，P<0.001為差異有極其顯著的統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 TUNEL、HE和CD31染色圖像分割

腫瘤生物學(xué)效應(yīng)中凋亡、壞死和血管是療效評(píng)估中最為常見的幾個(gè)指標(biāo)，分別對應(yīng)著TUNEL、HE和CD31染色圖像中的棕染細(xì)胞核、紅色細(xì)胞質(zhì)及白色空泡、棕染細(xì)胞膜部分。故以TUNEL和CD31染色圖像提取的棕色目標(biāo)區(qū)域表征凋亡和血管的空間分布，HE染色圖像提取的紅色和白色區(qū)域表征壞死的空間分布。

2.1.1 顏色閾值提取

如圖2（a），兩種染料著色的TUNEL、HE和CD31染色圖像在a*、b*顏色通道上K均值聚類進(jìn)行圖像分割。TUNEL和CD31（k=3）均被分割成3個(gè)聚類圖像，包括明場背景（聚類1）、藍(lán)染細(xì)胞核（聚類2）和棕染細(xì)核/膜（聚類3）；HE在k=3時(shí)，部分紅染細(xì)胞膜被劃分到藍(lán)染細(xì)胞核，故選擇k=4，將圖像分割成明場背景（聚類1）、藍(lán)染細(xì)胞核（聚類2）及紅色深染和淺染的細(xì)核膜（聚類3+4）。其在a*、b*顏色通道上的散點(diǎn)圖分別被標(biāo)記為藍(lán)、綠、紅和黃色。聚類3和4為目標(biāo)區(qū)域，在散點(diǎn)圖中占據(jù)著紅色和黃色區(qū)域。TUNEL和CD31中，基于a*、b*顏色通道的K均值聚類將明場背景（白色）和部分深染細(xì)胞核（接近黑色）歸為一類，在聚類1中同時(shí)含有部分目標(biāo)區(qū)域；HE中白色空泡被認(rèn)為是壞死。故根據(jù)各種的散點(diǎn)圖，排除聚類2，提取出各自聚類1、3和4圖像的分布曲線作為顏色閾值（表1），由兩條閾值關(guān)系曲線組成。

表1 TUNEL、HE和CD31染色圖像分割閾值（聚類1+3+4）Tab.1 Segmentation threshold of TUNEL,HE and CD31 stained images(cluster 1+3+4)

2.1.2 亮度閾值提取

如圖2（b），基于L*通道對TUNEL和CD31染色圖像中的目標(biāo)區(qū)域（顏色聚類1+3）再次K均值聚類（k=4，CD31的k=5），濾除目標(biāo)區(qū)域中的背景和淺色陰性信號(hào)（亮度聚類3+4+5），得到深染的陽性信號(hào)（亮度聚類1+2）。對HE染色圖像（顏色聚類1+3+4）再次K均值聚類（k=4），濾除目標(biāo)區(qū)域中的深染細(xì)胞核信號(hào)（亮度聚類3+4），得到淺染的陽性信號(hào)（亮度聚類1+2）。得到的亮度閾值如表1所示，其中L=0為分割后黑色背景，故去除零點(diǎn)。

2.1.3 圖像分割性能驗(yàn)證

如圖3，F(xiàn)iji-YCrCb對3種染色圖像均有明顯欠分割或過分割，而Kmeans-Lab在TUNEL和CD31染色圖像中含有的淺染噪聲明顯少于其他4種方法。將5種方法分割結(jié)果分別作為真實(shí)分割圖像（1～5分別代表MATLAB-YCrCb、MATLAB-Lab、Fiji-YCr-Cb、QuPath-CD和Kmeans-Lab方法的真實(shí)圖像分割結(jié)果），對應(yīng)5種方法的預(yù)測分割圖像結(jié)果（A～E分別代表MATLAB-YCrCb、MATLAB-Lab、Fiji-YCrCb、QuPath-CD和Kmeans-Lab方法的預(yù)測圖像分割結(jié)果）的性能參數(shù)如圖4所示?？偟膩碚f，NMI、Kappa系數(shù)和mIoU的結(jié)果基本一致，MATLAB-YCrCb與MATLAB-Lab的結(jié)果相似度最高，NMI為0.55～0.76，Kappa系數(shù)為0.79～0.91，mIoU為0.82～0.91；MATLAB-Lab與Kmeans-Lab結(jié)果相似度次之，NMI為0.42～0.61，Kappa系數(shù)為0.68～0.84，mIoU為0.74～0.82；Fiji-YCrCb與其他方法的結(jié)果相似度均很低；QuPath-CD與除Fiji-YCrCb以外其他3種方法的NMI、Kappa系數(shù)和mIoU均無明顯差異。將基于實(shí)際染料著色的QuPath-CD分割結(jié)果作為真實(shí)值，則對HE染色圖像分割的mPr均為0.73，mRe為0.96～0.97，mA均為0.94；對TUNEL染色圖像分割的mPr為0.84～0.90，mRe為0.90～0.93，mA為0.94～0.95；對CD31染色圖像分割的mPr為0.80～0.91，mRe為0.81～0.89，mA為0.95～0.96。如圖5，以蒙太奇方式分析筆者算法與QuPath-CD分割結(jié)果的差異，白色區(qū)域?yàn)檎骊栃?，目?biāo)區(qū)域顏色特征明顯；紫色區(qū)域?yàn)榧訇幮?，HE染色圖像主要表現(xiàn)為細(xì)胞核與細(xì)胞質(zhì)的交界區(qū)域及過染的細(xì)胞質(zhì)，TUNEL染色圖像主要表現(xiàn)為淺染的細(xì)胞核及淺染的背景，CD31染色圖像表現(xiàn)為淺染的背景和細(xì)胞核；綠色區(qū)域?yàn)榧訇栃?，主要表現(xiàn)為“椒鹽”噪聲和污漬。

圖3 TUNEL、HE和CD31染色圖像不同算法分割結(jié)果比較Fig.3 Comparison of segmentation results of different algorithms for TUNEL,HE and CD31 stained images

圖4 TUNEL、HE和CD31染色圖像不同算法的性能參數(shù)熱力圖Fig.4 Thermodynamic diagram comparison of performance parameters of different algorithms for TUNEL,HE and CD31 stained images

圖5 TUNEL、HE和CD31染色圖像Kmeans-Lab與QuPath-YCrCb方法的分割結(jié)果比較Fig.5 Comparisonofimagesegmentationresultsbetween Kmeans-Laband QuPath-YCrCbmethodsforof TUNEL,HEand CD31stained images

2.2 腫瘤WSIs凋亡、壞死和血管的空間分布分析

如圖6，以上述提取的a*，b*通道上的顏色閾值曲線和L*通道上的亮度閾值對腫瘤光動(dòng)力數(shù)據(jù)集WSIs進(jìn)行凋亡、壞死和血管圖像分割和空間分布分析，1幅526.1 M內(nèi)存的圖像分割平均耗時(shí)約為28 s。通過筆者的圖像分割和空間分布分析算法，PDT小鼠腫瘤凋亡、壞死和血管的橫截分布和皮下組織深度的分布密度（面積百分比）如圖6（a）和6（b）所示。如圖6（c），對凋亡、壞死、血管密度和光通量經(jīng)線性插值后進(jìn)行Pearson相關(guān)性分析，發(fā)現(xiàn)腫瘤壞死密度與光通量線性相關(guān)程度最高，Pearson相關(guān)系數(shù)為0.88；進(jìn)行線性回歸分析，可得到關(guān)系式為：

圖6 PDT小鼠腫瘤生物學(xué)效應(yīng)空間分布分析Fig.6 Image analysis of spatial distribution of tumor biological effects in PDT mice

其中：Y為壞死密度；X為光通量。

3 討論

腫瘤生物學(xué)效應(yīng)空間數(shù)據(jù)分析和可視化是腫瘤圖譜中的重要組成部分，不同于傳統(tǒng)IHC分析，整體腫瘤生物標(biāo)志物的分布區(qū)域和分布密度將是其研究重點(diǎn)。筆者針對腫瘤圖譜中生物學(xué)效應(yīng)的評(píng)估，提出基于陽性區(qū)域圖像K均值聚類獲取的目標(biāo)區(qū)域的顏色和亮度閾值，實(shí)現(xiàn)WSIs的快速分割。采用a*和b*通道的關(guān)系曲線和L*通道閾值作為顏色閾值和亮度閾值，與傳統(tǒng)單通道閾值相比分割精度更高，更符合人視覺習(xí)慣。以K均值聚類方法對色彩進(jìn)行分割提取閾值，與常用的上述其他IHC圖像分割方法相比，不僅提高了閾值提取的穩(wěn)定性和重復(fù)性，同時(shí)解決了IHC基于染料著色CD分割的圖像顏色不超過3種的局限性。在計(jì)算速度上，筆者算法經(jīng)過1次色域矩陣變換和顏色、亮度閾值分割得到掩膜的時(shí)間頻度為T（15 n），而基于CD和灰度閾值分割得到掩膜的時(shí)間頻度為T（14 n），兩者時(shí)間復(fù)雜度均為O（n）；對大內(nèi)存WSIs的分割速率約為54.5 s/G，快于基于神經(jīng)網(wǎng)絡(luò)的圖像分割算法[12]（313.6 s/G）。在圖像分割性能上，如圖3~5所示，與其他分割算法相比，Kmeans-Lab對淺染背景噪聲的分割優(yōu)于其他圖像分割方法。與基于真實(shí)染料著色的QuPath-CD相比，HE、TUNEL和CD31染色圖像中著色明顯的區(qū)域分割結(jié)果基本一致；在DAB著色的TUNEL和CD31染色圖片中對背景著色的分割效果優(yōu)于QuPath-CD，但是對灰色污染和“椒鹽”噪聲的分割效果不佳，故對制作的切片樣本有一定的要求，圖像分割時(shí)需通過中值濾波去除“椒鹽”噪聲。

以PDT小鼠乳腺癌為例，實(shí)驗(yàn)運(yùn)用Kmeans-Lab方法實(shí)現(xiàn)WSIs圖像分割，結(jié)合GIS中的分布密度和KDE熱力圖對壞死、凋亡和血管量化和可視化。如圖6，Kmeans-Lab方法結(jié)合空間分布分析方法用于序列HE、TUNEL和CD31染色的WSIs顯示了PDT后壞死、凋亡和血管在腫瘤不同皮下組織深度橫截面中的分布位置，對分布密度進(jìn)行統(tǒng)計(jì)學(xué)分析判定和推算了各參數(shù)間的相關(guān)性和關(guān)系式，發(fā)現(xiàn)PDT引起的壞死與光通量線性相關(guān)，驗(yàn)證了筆者算法在腫瘤生物學(xué)效應(yīng)空間數(shù)據(jù)分析中的有效性。