基于數(shù)字圖像處理的初級纖毛定量分析

李 金 楊鵬飛 王菊芳 何進(jìn)鵬 馬 偉 周 珩?

1(蘭州大學(xué) 核科學(xué)與技術(shù)學(xué)院，蘭州 730000)

2(中國科學(xué)院近代物理研究所 甘肅省空間輻射生物學(xué)重點實驗室 蘭州 730000)

3(中國科學(xué)院大學(xué) 核科學(xué)與技術(shù)學(xué)院 北京 100049)

引言

纖毛是一種來源于細(xì)胞基體并突出于細(xì)胞表面的特殊細(xì)胞器，有3 種類型:初級纖毛(primary cilia)、運(yùn)動纖毛(motile cilia)和節(jié)點纖毛(nodal cilia)[1]。 初級纖毛作為細(xì)胞的“天線”，通過調(diào)控Hedgehog[2]、Wnt[3]、TGF－β[4]、Notch[5]等多種信號通路感應(yīng)物理、化學(xué)和生物信號的刺激，并傳導(dǎo)至細(xì)胞內(nèi)部。 初級纖毛的結(jié)構(gòu)異常會導(dǎo)致多種疾病[6-7]，甚至影響腫瘤的發(fā)生和發(fā)展[7]。 通過化學(xué)藥物抑制初級纖毛生長能夠有效控制腫瘤細(xì)胞的惡性增殖[8]。

目前，多采用Photoshop 軟件中的“標(biāo)度尺”對初級纖毛進(jìn)行測量，但存在準(zhǔn)確度低、重復(fù)性差、分析效率低等問題，因此，建立一種基于數(shù)字圖像處理技術(shù)的高效分析方法非常必要。

數(shù)字圖像處理技術(shù)可實現(xiàn)包括噪聲去除[9]、分割[10]、配準(zhǔn)[11]、融合[12]等圖像處理。 噪聲去除是指消除、弱化或壓制圖像中的隨機(jī)變化干擾。 這些干擾可能源于成像過程中的隨機(jī)效應(yīng)與測量誤差。分割是將圖像細(xì)分為它的組成要素，細(xì)分水平取決于要解決的問題，即感興趣的物體被隔離出來時，就應(yīng)該停止分割。 常見的分割方法有基于邊緣分割、基于區(qū)域分割、基于特定理論方法和基于深度學(xué)習(xí)的分割算法等。

初級纖毛的數(shù)字圖像處理的難點在于噪聲去除與邊界分割。 首先對相關(guān)醫(yī)學(xué)圖像降噪分割方法進(jìn)行了調(diào)研。 陳炳權(quán)等[13]提出了一種基于離散小波變換和修正中值濾波的醫(yī)學(xué)圖像耦合去噪算法，用于CT 與超聲波圖像降噪。 Sarungbam Bonny等[14]提出了一種非局部均值濾波法用于超聲圖像去除散斑噪聲，該方法能夠保持目標(biāo)邊緣細(xì)節(jié)和結(jié)構(gòu)信息。 Erwin 等[15]介紹了一種改進(jìn)后的Prewitt算子用于視網(wǎng)膜血管分割，取得了較好的分割效果。 另外，還有一些對邊緣細(xì)節(jié)保持和分割較好的算法如:中值濾波[13]、均值濾波[14]、Prewitt[15]、Sobel[16]、Roberts[17]、高斯－拉普拉斯(Log)[18]等。本研究將上述算法用于對初級纖毛形態(tài)特征的定量分析，對相應(yīng)結(jié)果進(jìn)行了對比分析，然后，選取了較優(yōu)算法進(jìn)行了二次改進(jìn)，構(gòu)建了基于較優(yōu)算法的初級纖毛形態(tài)特征的自動分析方法。 該方法能自動對每個初級纖毛進(jìn)行編號，自動計算其周長和面積等形態(tài)參數(shù)。 從而達(dá)到了客觀、準(zhǔn)確、快速、多元化的處理目的，為初級纖毛圖像數(shù)字特征分析提供了有效的技術(shù)方案。

1 材料與方法

1.1 圖像采集

首先進(jìn)行圖像采集，其質(zhì)量將直接影響后續(xù)操作，如分割算法的難易以及數(shù)據(jù)結(jié)果和分析效率等。 本研究采用電離輻射后的M059K 細(xì)胞(人腦神經(jīng)膠質(zhì)瘤細(xì)胞)的初級纖毛圖像作為案例，如圖1所示。 X 射線由PXI Precision X-RAY225 系統(tǒng)產(chǎn)生，電壓為225 kV，電流為13.3 mA，照射劑量為4 Gy。 圖像分辨率為2 048×1 536像素，由激光掃描共聚焦顯微鏡(FLUOVIEW FV3000， Olympus)拍攝，圖片保存為jpg 格式。 像素長度與真實長度間的比例尺為11.209 7 pixel/μm。

圖1 X 射線照射M059K 細(xì)胞后初級纖毛免疫熒光圖像Fig.1 Immunofluorescence image of primary cilia of M059K cells irradiated with X-ray

制片操作過程:首先利用4%多聚甲醛對X 射線照射后的細(xì)胞于室溫固定10 min，再使用含0.5%Triton X－100 的磷酸鹽緩沖液(PBS)透膜5 min，然后在含5%山羊血清的PBS 室溫封閉2 h，之后利用Arl13b(1 ∶500)和γ－Tubulin(1 ∶500)一抗室溫孵育2 h，熒光二抗(1 ∶500)室溫避光孵育2 h，最后使用含0.5% Tween－20 的PBS 洗片，并通過DAPI(4′， 6－diamidino－2－phenylindole)染色細(xì)胞核后封片、采集制作數(shù)字圖像。

1.2 圖像分離與分割

圖像分析設(shè)備為Intel core i5，3 GHz CPU，8 GB內(nèi)存計算機(jī)，分析軟件為Matlab 2017。 初級纖毛原始數(shù)字圖像如圖1 所示，采用RGB 彩色模式，每個像素點由R、G、B 等3 個字節(jié)構(gòu)成。 經(jīng)Arl13b(Proteintech，＃ 66739)標(biāo)記染色的初級纖毛呈現(xiàn)紅色；由γ-Tubulin(Sigma，＃SAB4503045)標(biāo)記的細(xì)胞骨架呈現(xiàn)綠色；細(xì)胞核呈現(xiàn)藍(lán)色。 根據(jù)染色特點，可以利用不同的色彩信息將細(xì)胞核與細(xì)胞骨架從圖像中分離出來。 纖毛與背景中雜質(zhì)的亮度不同，可依據(jù)灰度信息差異將纖毛從背景雜質(zhì)中分割出來。 以此引入最大類間方差法進(jìn)行分割。

1.3 圖像色彩分離

RGB 圖像可以看作由3 個灰度圖像形成的“堆棧”，圖像是一個M×N×3 大小的彩色像素數(shù)組，當(dāng)發(fā)送到彩色監(jiān)視器的紅、綠、藍(lán)輸入端時，就會在屏幕上產(chǎn)生彩色圖像。 通過對纖毛圖像色彩分析，只需將三維RGB 圖像提取出二維的R 分量，即可將細(xì)胞核與細(xì)胞骨架分離出來。

1.4 不同濾波方法對比

濾波是為了突出圖像的空間信息，壓抑或者去除其他無關(guān)的信息。 濾波方法包括空間濾波和頻域濾波，空間濾波又包括線性空間濾波和非線性空間濾波。 線性空間濾波基于計算乘積和(即線性操作)，非線性空間濾波基于涉及領(lǐng)域像素內(nèi)的非線性操作，例如通過領(lǐng)域內(nèi)的像素最大值代替每個中心點的響應(yīng)值。 本研究中的濾波結(jié)果如圖2 所示，初級纖毛圖像(a)分別經(jīng)均值濾波(b)、Log 濾波(c)、Prewitt 濾波(d)、Sobel 濾波(e)、中值濾波(f)以及改進(jìn)型濾波(g)處理。

圖2 不同濾波方法處理后的纖毛圖像。 (a)原始圖像；(b)均值濾波；(c)Log 濾波；(d)Prewitt 濾波；(e)Sobel 濾波；(f)中值濾波；(g)改進(jìn)型濾波Fig.2 The primary ciliary images processed by different filtering methods.(a) Original image；(b)Mean filtering； (c)Log filtering； (d)Prewitt filtering； (e)Sobel filtering； (f)Median filtering；(g)Improved filtering.

線性空間濾波也稱“空間卷積”。 卷積本身是數(shù)學(xué)領(lǐng)域內(nèi)的一種積分變化方法，在數(shù)字圖像處理技術(shù)中，濾波方法如式(1)，濾波原理如圖3 所示。

式中:x，y是像素在圖片中的位置坐標(biāo)；k，l是卷積核中的位置；m[x，y] 是濾波結(jié)果，即圖3 中的結(jié)果40；f[k，l] 是濾波函數(shù)也稱權(quán)重系數(shù)，即圖示中的卷積核；I[x ＋k，y ＋l] 是與w[k，l] 相對應(yīng)的圖片像素值。

圖3 圖像濾波方法流程Fig.3 The filtering method flow

綜合不同濾波結(jié)果得出結(jié)論:均值與中值濾波能夠過濾部分細(xì)胞骨架，但也會導(dǎo)致初級纖毛邊界模糊；Sobel 濾波算子能夠在一個方向上增加纖毛的邊界對比度；Prewitt 濾波算子能夠增強(qiáng)整個纖毛的對比度，但當(dāng)纖毛呈豎直排列時反而會減弱對比度；Log 濾波對比度沒有Prewitt 濾波算子明顯。 針對以上問題，將Prewitt 處理后的圖片與原圖相加，或者將Log 濾波后的圖像取反與原圖相加。 經(jīng)過改進(jìn)后的濾波方法使得初級纖毛圖像的邊界更加清晰，整體亮度提升明顯。

1.5 OTSU′s 閾值處理方法

OTSU′s 最大類間方差法[19-20]是在最小二乘法原理為基礎(chǔ)上進(jìn)行的推導(dǎo)演變。 基本實現(xiàn)原理依據(jù)閾值k將圖像分為C1 與C2 兩個部分，最優(yōu)閾值k的選擇標(biāo)準(zhǔn)是C1 與C2 兩部分之間的最大類間方差。 閾值k選取類似于最小二乘法，k的范圍是C1中0 到L－1 內(nèi)的整數(shù)。 OTSU′s 閾值處理方法具有良好的自適應(yīng)分割性能，具體實現(xiàn)方法如下:

令一幅圖像的直方圖成分由下式所示:

式中，n是圖像中像素的總數(shù)，ni是灰度級為i的數(shù)目，L是圖像中所有可能的灰度級數(shù)。 假設(shè)閾值k已經(jīng)選定，C1 是灰度值從0 到k的像素，C2 是灰度值從k＋1 到L－1 的像素，此處的k均為正整數(shù)。

閾值k由Otsu′s 公式選定，即

式中，) 是方差，P1(k)，P2(k) 分別是像素小于等于閾值集合和像素大于閾值集合發(fā)生概率；m1k和m2k分別是像素小于等于閾值集合和像素大于閾值集合的平均灰度；mG全局均值。

如果設(shè)置k＝0，那么擁有為k賦值的任何像素C1 集合的概率為0。

Otsu′s 算法的優(yōu)點是:方差(k) 較大，完全分割一幅圖像的閾值將會更加接近。 另外，因為k是范圍0 到L－1 內(nèi)的整數(shù)，所以更加易于尋找(k)的最大值。 逐步通過k 的可能個L 值，計算每一步的(k) ，然后選擇) 最大值為k值，即最佳閾值。 如果最大值不唯一，則選取的最優(yōu)閾值即所有k值的平均值。

1.6 連通域標(biāo)記與初級纖毛定量分析

一個坐標(biāo)(x，y)的像素p具有兩個垂直和兩個水平的相鄰像素，它們的坐標(biāo)分別為(x＋1，y)，(x－1，y)，(x，y＋1)和(x，y－1)。 這個p的相鄰像素集合記為N4(p)。p的4 個對角線相鄰像素坐標(biāo)分別為(x＋1，y＋1)，(x＋1，y－1)，(x－1，y＋1)，(x－1，y－1)，這樣p的4 個對角像素被記為ND(p)。N4(p)和ND(p)的并集是p的8 個相鄰像素，記為N8(p)。

現(xiàn)有另一像素q，若q∈N4(p)，則像素p和q稱為4 鄰接。 同樣，若q∈N8(p)，則p和q稱為8鄰接。 若在前景像素p和q之間存在一條完全由前景像素組成的4 連接路徑，則這兩個前景像素稱為4 連接。 若他們直接存在一條8 連接路徑，則稱為8連接。

由于在分辨率為2 048像素×1 536像素的圖像中初級纖毛面積大于60 像素，利用最大類間方差法獲得的二值圖像，應(yīng)進(jìn)行形態(tài)學(xué)開運(yùn)算，以便獲取面積大于60 像素的區(qū)域。 本研究采用8 連接處理，并將編號標(biāo)記在每個連通域的質(zhì)心，結(jié)果如圖4 所示。 初級纖毛像素周長等于圖4 中被標(biāo)記的連通域周長，像素面積等于連通域所含有像素點的個數(shù)。纖毛實際周長＝纖毛像素面積/比例尺；纖毛實際面積＝纖毛像素面積/比例尺2。

圖4 連通域標(biāo)記Fig.4 The connected domain markers

1.7 數(shù)據(jù)評價指標(biāo)

根據(jù)初級纖毛定量分析結(jié)果，通過計算不同方法下的相對誤差進(jìn)行評價，計算公式如下:

式中，σA是本模型算法獲取周長的標(biāo)準(zhǔn)差，σB是手動畫線量取長度的標(biāo)準(zhǔn)差。 考慮兩者應(yīng)在同一水平下對比，纖毛周長CA應(yīng)換算成長度A，換算公式如下:

式中，ε取值范圍0.01～0.05。

2 實驗結(jié)果

2.1 算法性能評價

在圖4 中隨機(jī)選取10 個初級纖毛，分別采用本研究模型和手動畫線方法采集數(shù)據(jù)。 如表1 所示，自動欄顯示本研究模型采集并計算出的10 個初級纖毛周長，手動欄顯示直線“標(biāo)度尺”手動分段測量方法采集的10 個對應(yīng)初級纖毛的長度(單位:pixel)。 根據(jù)表1 數(shù)據(jù)分別算出本模型算法與手動畫線量取的檢測結(jié)果的相對誤差ω，為了使衡量的結(jié)果更具有說服力，選擇ε＝0.01，算出ω＝19.08% 。 從結(jié)果可以得出本模型相比手動畫線方法算法離散程度更小，從一定程度可以說明測量結(jié)果更加準(zhǔn)確。 從算法性能方面比較，手動畫線方法獲得上述一張圖片的纖毛長度需要50～80 min，且操作過程復(fù)雜，精確度較低。 而應(yīng)用本研究模型算法僅需2 s，即可獲取較為準(zhǔn)確的初級纖毛周長與面積參數(shù)。

表1 本模型與手動畫線量取的初級纖毛長度Tab.1 The length of primary cilia measured automatically and manually

3 討論和結(jié)論

本研究提出了一種基于數(shù)字圖像處理技術(shù)的初級纖毛自動分析方法，該方法有助于大幅提高對初級纖毛的定量分析精確性。 與手動直線“標(biāo)度尺”分段測量方法相比，本研究通過二次改進(jìn)Prewitt、Log 算子以及結(jié)合圖像色彩分離、連通域標(biāo)記等技術(shù)，實現(xiàn)了對初級纖毛結(jié)構(gòu)指標(biāo)的全自動精確定量分析；由于算法性能的提升，分析時間較手動處理所需的50～80 min 減少至2 s 左右。 另外，本研究的分析離散程度也更低。 綜上所述，本研究在技術(shù)上較大改進(jìn)了初級纖毛檢測的準(zhǔn)確度、分析效率、重復(fù)性。 但由于本研究樣本獲取的多樣性和樣本量均不充分，因此仍需更充分的研究予以驗證。 為此，相關(guān)后續(xù)研究將補(bǔ)充更多數(shù)據(jù)來源的圖像，進(jìn)一步更多的訓(xùn)練，來提高本研究算法的性能。

國內(nèi)外對初級纖毛的研究多集中在纖毛介導(dǎo)的相關(guān)信號通路響應(yīng)和外界刺激對纖毛增殖與形態(tài)功能改變的生物學(xué)效應(yīng)影響。 其中，纖毛增殖與形態(tài)學(xué)研究因受限于傳統(tǒng)測量方法，而影響其快速進(jìn)展。 本研究對纖毛形態(tài)特征的自動分析，快速精確地測量纖毛的周長和面積參數(shù)，有助于促進(jìn)初級纖毛的研究發(fā)展。