金文東，陰慧娟，李迎新

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

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

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

1 材料與方法

1.1 實驗材料

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

1.2 方法

1.2.1 閾值提取和圖像分割

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

1.2.2 空間分布分析

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

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

1.2.3 圖像分割性能驗證

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

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

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

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

1.3 統(tǒng)計學方法

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

2 結(jié)果

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

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

2.1.1 顏色閾值提取

如圖2（a），兩種染料著色的TUNEL、HE和CD31染色圖像在a*、b*顏色通道上K均值聚類進行圖像分割。TUNEL和CD31（k=3）均被分割成3個聚類圖像，包括明場背景（聚類1）、藍染細胞核（聚類2）和棕染細核/膜（聚類3）；HE在k=3時，部分紅染細胞膜被劃分到藍染細胞核，故選擇k=4，將圖像分割成明場背景（聚類1）、藍染細胞核（聚類2）及紅色深染和淺染的細核膜（聚類3+4）。其在a*、b*顏色通道上的散點圖分別被標記為藍、綠、紅和黃色。聚類3和4為目標區(qū)域，在散點圖中占據(jù)著紅色和黃色區(qū)域。TUNEL和CD31中，基于a*、b*顏色通道的K均值聚類將明場背景（白色）和部分深染細胞核（接近黑色）歸為一類，在聚類1中同時含有部分目標區(qū)域；HE中白色空泡被認為是壞死。故根據(jù)各種的散點圖，排除聚類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染色圖像中的目標區(qū)域（顏色聚類1+3）再次K均值聚類（k=4，CD31的k=5），濾除目標區(qū)域中的背景和淺色陰性信號（亮度聚類3+4+5），得到深染的陽性信號（亮度聚類1+2）。對HE染色圖像（顏色聚類1+3+4）再次K均值聚類（k=4），濾除目標區(qū)域中的深染細胞核信號（亮度聚類3+4），得到淺染的陽性信號（亮度聚類1+2）。得到的亮度閾值如表1所示，其中L=0為分割后黑色背景，故去除零點。

2.1.3 圖像分割性能驗證

如圖3，F(xiàn)iji-YCrCb對3種染色圖像均有明顯欠分割或過分割，而Kmeans-Lab在TUNEL和CD31染色圖像中含有的淺染噪聲明顯少于其他4種方法。將5種方法分割結(jié)果分別作為真實分割圖像（1～5分別代表MATLAB-YCrCb、MATLAB-Lab、Fiji-YCr-Cb、QuPath-CD和Kmeans-Lab方法的真實圖像分割結(jié)果），對應5種方法的預測分割圖像結(jié)果（A～E分別代表MATLAB-YCrCb、MATLAB-Lab、Fiji-YCrCb、QuPath-CD和Kmeans-Lab方法的預測圖像分割結(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均無明顯差異。將基于實際染料著色的QuPath-CD分割結(jié)果作為真實值，則對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ū)域為真陽性，目標區(qū)域顏色特征明顯；紫色區(qū)域為假陰性，HE染色圖像主要表現(xiàn)為細胞核與細胞質(zhì)的交界區(qū)域及過染的細胞質(zhì)，TUNEL染色圖像主要表現(xiàn)為淺染的細胞核及淺染的背景，CD31染色圖像表現(xiàn)為淺染的背景和細胞核；綠色區(qū)域為假陽性，主要表現(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*通道上的亮度閾值對腫瘤光動力數(shù)據(jù)集WSIs進行凋亡、壞死和血管圖像分割和空間分布分析，1幅526.1 M內(nèi)存的圖像分割平均耗時約為28 s。通過筆者的圖像分割和空間分布分析算法，PDT小鼠腫瘤凋亡、壞死和血管的橫截分布和皮下組織深度的分布密度（面積百分比）如圖6（a）和6（b）所示。如圖6（c），對凋亡、壞死、血管密度和光通量經(jīng)線性插值后進行Pearson相關(guān)性分析，發(fā)現(xiàn)腫瘤壞死密度與光通量線性相關(guān)程度最高，Pearson相關(guān)系數(shù)為0.88；進行線性回歸分析，可得到關(guān)系式為：

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

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

3 討論

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

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