基于数字图像处理的初级纤毛定量分析

李 金 杨鹏飞 王菊芳 何进鹏 马 伟 周 珩∗

1(兰州大学 核科学与技术学院，兰州 730000)

2(中国科学院近代物理研究所 甘肃省空间辐射生物学重点实验室 兰州 730000)

3(中国科学院大学 核科学与技术学院 北京 100049)

引言

纤毛是一种来源于细胞基体并突出于细胞表面的特殊细胞器，有3 种类型:初级纤毛(primary cilia)、运动纤毛(motile cilia)和节点纤毛(nodal cilia)[1]。 初级纤毛作为细胞的“天线”，通过调控Hedgehog[2]、Wnt[3]、TGF－β[4]、Notch[5]等多种信号通路感应物理、化学和生物信号的刺激，并传导至细胞内部。 初级纤毛的结构异常会导致多种疾病[6-7]，甚至影响肿瘤的发生和发展[7]。 通过化学药物抑制初级纤毛生长能够有效控制肿瘤细胞的恶性增殖[8]。

目前，多采用Photoshop 软件中的“标度尺”对初级纤毛进行测量，但存在准确度低、重复性差、分析效率低等问题，因此，建立一种基于数字图像处理技术的高效分析方法非常必要。

数字图像处理技术可实现包括噪声去除[9]、分割[10]、配准[11]、融合[12]等图像处理。 噪声去除是指消除、弱化或压制图像中的随机变化干扰。 这些干扰可能源于成像过程中的随机效应与测量误差。分割是将图像细分为它的组成要素，细分水平取决于要解决的问题，即感兴趣的物体被隔离出来时，就应该停止分割。 常见的分割方法有基于边缘分割、基于区域分割、基于特定理论方法和基于深度学习的分割算法等。

初级纤毛的数字图像处理的难点在于噪声去除与边界分割。 首先对相关医学图像降噪分割方法进行了调研。 陈炳权等[13]提出了一种基于离散小波变换和修正中值滤波的医学图像耦合去噪算法，用于CT 与超声波图像降噪。 Sarungbam Bonny等[14]提出了一种非局部均值滤波法用于超声图像去除散斑噪声，该方法能够保持目标边缘细节和结构信息。 Erwin 等[15]介绍了一种改进后的Prewitt算子用于视网膜血管分割，取得了较好的分割效果。 另外，还有一些对边缘细节保持和分割较好的算法如:中值滤波[13]、均值滤波[14]、Prewitt[15]、Sobel[16]、Roberts[17]、高斯－拉普拉斯(Log)[18]等。本研究将上述算法用于对初级纤毛形态特征的定量分析，对相应结果进行了对比分析，然后，选取了较优算法进行了二次改进，构建了基于较优算法的初级纤毛形态特征的自动分析方法。 该方法能自动对每个初级纤毛进行编号，自动计算其周长和面积等形态参数。 从而达到了客观、准确、快速、多元化的处理目的，为初级纤毛图像数字特征分析提供了有效的技术方案。

1 材料与方法

1.1 图像采集

首先进行图像采集，其质量将直接影响后续操作，如分割算法的难易以及数据结果和分析效率等。 本研究采用电离辐射后的M059K 细胞(人脑神经胶质瘤细胞)的初级纤毛图像作为案例，如图1所示。 X 射线由PXI Precision X-RAY225 系统产生，电压为225 kV，电流为13.3 mA，照射剂量为4 Gy。 图像分辨率为2 048×1 536像素，由激光扫描共聚焦显微镜(FLUOVIEW FV3000， Olympus)拍摄，图片保存为jpg 格式。 像素长度与真实长度间的比例尺为11.209 7 pixel/μm。

图1 X 射线照射M059K 细胞后初级纤毛免疫荧光图像Fig.1 Immunofluorescence image of primary cilia of M059K cells irradiated with X-ray

制片操作过程:首先利用4%多聚甲醛对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)染色细胞核后封片、采集制作数字图像。

1.2 图像分离与分割

图像分析设备为Intel core i5，3 GHz CPU，8 GB内存计算机，分析软件为Matlab 2017。 初级纤毛原始数字图像如图1 所示，采用RGB 彩色模式，每个像素点由R、G、B 等3 个字节构成。 经Arl13b(Proteintech，＃ 66739)标记染色的初级纤毛呈现红色；由γ-Tubulin(Sigma，＃SAB4503045)标记的细胞骨架呈现绿色；细胞核呈现蓝色。 根据染色特点，可以利用不同的色彩信息将细胞核与细胞骨架从图像中分离出来。 纤毛与背景中杂质的亮度不同，可依据灰度信息差异将纤毛从背景杂质中分割出来。 以此引入最大类间方差法进行分割。

1.3 图像色彩分离

RGB 图像可以看作由3 个灰度图像形成的“堆栈”，图像是一个M×N×3 大小的彩色像素数组，当发送到彩色监视器的红、绿、蓝输入端时，就会在屏幕上产生彩色图像。 通过对纤毛图像色彩分析，只需将三维RGB 图像提取出二维的R 分量，即可将细胞核与细胞骨架分离出来。

1.4 不同滤波方法对比

滤波是为了突出图像的空间信息，压抑或者去除其他无关的信息。 滤波方法包括空间滤波和频域滤波，空间滤波又包括线性空间滤波和非线性空间滤波。 线性空间滤波基于计算乘积和(即线性操作)，非线性空间滤波基于涉及领域像素内的非线性操作，例如通过领域内的像素最大值代替每个中心点的响应值。 本研究中的滤波结果如图2 所示，初级纤毛图像(a)分别经均值滤波(b)、Log 滤波(c)、Prewitt 滤波(d)、Sobel 滤波(e)、中值滤波(f)以及改进型滤波(g)处理。

图2 不同滤波方法处理后的纤毛图像。 (a)原始图像；(b)均值滤波；(c)Log 滤波；(d)Prewitt 滤波；(e)Sobel 滤波；(f)中值滤波；(g)改进型滤波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.

线性空间滤波也称“空间卷积”。 卷积本身是数学领域内的一种积分变化方法，在数字图像处理技术中，滤波方法如式(1)，滤波原理如图3 所示。

式中:x，y是像素在图片中的位置坐标；k，l是卷积核中的位置；m[x，y] 是滤波结果，即图3 中的结果40；f[k，l] 是滤波函数也称权重系数，即图示中的卷积核；I[x ＋k，y ＋l] 是与w[k，l] 相对应的图片像素值。

图3 图像滤波方法流程Fig.3 The filtering method flow

综合不同滤波结果得出结论:均值与中值滤波能够过滤部分细胞骨架，但也会导致初级纤毛边界模糊；Sobel 滤波算子能够在一个方向上增加纤毛的边界对比度；Prewitt 滤波算子能够增强整个纤毛的对比度，但当纤毛呈竖直排列时反而会减弱对比度；Log 滤波对比度没有Prewitt 滤波算子明显。 针对以上问题，将Prewitt 处理后的图片与原图相加，或者将Log 滤波后的图像取反与原图相加。 经过改进后的滤波方法使得初级纤毛图像的边界更加清晰，整体亮度提升明显。

1.5 OTSU′s 阈值处理方法

OTSU′s 最大类间方差法[19-20]是在最小二乘法原理为基础上进行的推导演变。 基本实现原理依据阈值k将图像分为C1 与C2 两个部分，最优阈值k的选择标准是C1 与C2 两部分之间的最大类间方差。 阈值k选取类似于最小二乘法，k的范围是C1中0 到L－1 内的整数。 OTSU′s 阈值处理方法具有良好的自适应分割性能，具体实现方法如下:

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

式中，n是图像中像素的总数，ni是灰度级为i的数目，L是图像中所有可能的灰度级数。 假设阈值k已经选定，C1 是灰度值从0 到k的像素，C2 是灰度值从k＋1 到L－1 的像素，此处的k均为正整数。

阈值k由Otsu′s 公式选定，即

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

如果设置k＝0，那么拥有为k赋值的任何像素C1 集合的概率为0。

Otsu′s 算法的优点是:方差(k) 较大，完全分割一幅图像的阈值将会更加接近。 另外，因为k是范围0 到L－1 内的整数，所以更加易于寻找(k)的最大值。 逐步通过k 的可能个L 值，计算每一步的(k) ，然后选择) 最大值为k值，即最佳阈值。 如果最大值不唯一，则选取的最优阈值即所有k值的平均值。

1.6 连通域标记与初级纤毛定量分析

一个坐标(x，y)的像素p具有两个垂直和两个水平的相邻像素，它们的坐标分别为(x＋1，y)，(x－1，y)，(x，y＋1)和(x，y－1)。 这个p的相邻像素集合记为N4(p)。p的4 个对角线相邻像素坐标分别为(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)。

现有另一像素q，若q∈N4(p)，则像素p和q称为4 邻接。 同样，若q∈N8(p)，则p和q称为8邻接。 若在前景像素p和q之间存在一条完全由前景像素组成的4 连接路径，则这两个前景像素称为4 连接。 若他们直接存在一条8 连接路径，则称为8连接。

由于在分辨率为2 048像素×1 536像素的图像中初级纤毛面积大于60 像素，利用最大类间方差法获得的二值图像，应进行形态学开运算，以便获取面积大于60 像素的区域。 本研究采用8 连接处理，并将编号标记在每个连通域的质心，结果如图4 所示。 初级纤毛像素周长等于图4 中被标记的连通域周长，像素面积等于连通域所含有像素点的个数。纤毛实际周长＝纤毛像素面积/比例尺；纤毛实际面积＝纤毛像素面积/比例尺2。

图4 连通域标记Fig.4 The connected domain markers

1.7 数据评价指标

根据初级纤毛定量分析结果，通过计算不同方法下的相对误差进行评价，计算公式如下:

式中，σA是本模型算法获取周长的标准差，σB是手动画线量取长度的标准差。 考虑两者应在同一水平下对比，纤毛周长CA应换算成长度A，换算公式如下:

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

2 实验结果

2.1 算法性能评价

在图4 中随机选取10 个初级纤毛，分别采用本研究模型和手动画线方法采集数据。 如表1 所示，自动栏显示本研究模型采集并计算出的10 个初级纤毛周长，手动栏显示直线“标度尺”手动分段测量方法采集的10 个对应初级纤毛的长度(单位:pixel)。 根据表1 数据分别算出本模型算法与手动画线量取的检测结果的相对误差ω，为了使衡量的结果更具有说服力，选择ε＝0.01，算出ω＝19.08% 。 从结果可以得出本模型相比手动画线方法算法离散程度更小，从一定程度可以说明测量结果更加准确。 从算法性能方面比较，手动画线方法获得上述一张图片的纤毛长度需要50～80 min，且操作过程复杂，精确度较低。 而应用本研究模型算法仅需2 s，即可获取较为准确的初级纤毛周长与面积参数。

表1 本模型与手动画线量取的初级纤毛长度Tab.1 The length of primary cilia measured automatically and manually

3 讨论和结论

本研究提出了一种基于数字图像处理技术的初级纤毛自动分析方法，该方法有助于大幅提高对初级纤毛的定量分析精确性。 与手动直线“标度尺”分段测量方法相比，本研究通过二次改进Prewitt、Log 算子以及结合图像色彩分离、连通域标记等技术，实现了对初级纤毛结构指标的全自动精确定量分析；由于算法性能的提升，分析时间较手动处理所需的50～80 min 减少至2 s 左右。 另外，本研究的分析离散程度也更低。 综上所述，本研究在技术上较大改进了初级纤毛检测的准确度、分析效率、重复性。 但由于本研究样本获取的多样性和样本量均不充分，因此仍需更充分的研究予以验证。 为此，相关后续研究将补充更多数据来源的图像，进一步更多的训练，来提高本研究算法的性能。

国内外对初级纤毛的研究多集中在纤毛介导的相关信号通路响应和外界刺激对纤毛增殖与形态功能改变的生物学效应影响。 其中，纤毛增殖与形态学研究因受限于传统测量方法，而影响其快速进展。 本研究对纤毛形态特征的自动分析，快速精确地测量纤毛的周长和面积参数，有助于促进初级纤毛的研究发展。