酿酒酵母单细胞形态参数精准提取算法的研究∗

刘 可王颖瀛耿杨烨肖 秦朱 真∗

(1.东南大学电子科学与工程学院，江苏 南京210096；2.东南大学微电子学院，江苏南京210096)

细胞衰老是一个复杂的生理过程，伴随着细胞生理结构和功能的逐渐丧失。酿酒酵母(Saccharomyces cerevisiae)作为常用的模式细胞，具有繁殖快、寿命短、培养成本低等优点。此外，酿酒酵母细胞与哺乳动物细胞、人类具有相似的保守生化机制，其存活率曲线与人类寿命类似[1]，常被用于生命体衰老机制的研究[2]。酿酒酵母的复制寿命(Replicative lifespan，RLS)是酵母细胞衰老研究的主要指标，定义为一个酵母细胞从新生到死亡过程中产生子代细胞的个数[3]。研究表明，在酵母细胞的生长、增殖、衰老过程中，细胞的截面面积、体积等形态参数变化用于表征细胞生长速率，并且与细胞周期调控和细胞复制寿命紧密相关[4]，但相关机理尚未被完全揭示。因此，需要进一步研究细胞表型特征参数及其在细胞完整生命周期内的动态变化。微流控芯片利用微结构捕获并固定酵母单细胞，并通过流体剪切力实现子细胞的高效去除[2]，是酵母细胞复制衰老研究的重要平台。在酵母细胞实验中，实验者可以借助微流控芯片捕获、固定、培养酵母母细胞、并去除子代细胞、对细胞进行时序显微成像。酿酒酵母细胞平均复制寿命为20代~30代[5]，每代持续时间大约为90 min，若以10 min为时序成像时间间隔，单个酵母细胞完整生命周期的时序监测将产生数百帧显微图像。对于具有多个单细胞捕获结构及成像位点的高通量阵列式微流控细胞培养芯片，单次实验需分析上万帧图像，因此急需提出一种精确高效、自动化的细胞图像处理算法，以提取酵母细胞的形态参数。

目前已经有几种成熟的算法和软件用于细胞的图像分割和目标跟踪[6]。主流算法包括主动轮廓算法、分水岭算法等传统图像分割算法的组合及卷积神经网络等[7－8]。细胞图像处理的相关软件包括ImageJ、CellProfiler、CellStar、CellX等。ImageJ使用最普遍，该软件是基于Java开发的图像处理软件，具有开放式结构和可扩展性，支持一些用户编写的插件或宏如BudJ[9]，其优点是步骤灵活，但是难以实现图像的批量处理。CellProfiler是用于识别和量化细胞表型参数的相对完善的图像分析软件，用户无需使用编程语言就能通过交互界面自行设计不同的算法流程，然而其更适用于高分辨率、视野中存在多个细胞的图像，且不适合处理延时图像[10]。Cristian等人提出了适用于明场酵母细胞图像的处理软件CellStar，并建立了酵母细胞图像处理软件的评估平台，综合比较了各种算法或软件的适用范围和效果，且平台仍在持续更新[6]。大部分已有的算法和软件要求图像具有均一、明亮的背景，且多数针对细胞成片生长的情况进行细胞分割。基于机器学习的算法相较于传统方法虽然具有更强的通用性，但其复杂度更高，运行和调试成本也更高[11]。

然而，微流控芯片上细胞捕获微结构的几何轮廓造成了图像的背景不均一，且母细胞可能与微结构边缘产生交叠，影响细胞轮廓的判别[12－13]，已有的图像处理软件不适用。由于微流控芯片结构和成像条件的差异，不同微流控芯片上采集的实验图像通常需要设计专用的算法进行分析处理。本文中的高通量的酵母单细胞捕获－培养－解剖微流控芯片上设有阵列式的捕获微结构[14]，因此采用通过检测形状规则的细胞捕获微结构确定被捕获细胞的位置作为图像分割的思路。我们针对此芯片设计了专用的图像处理算法，该算法结合Hough变换[15]、分水岭算法、形态学膨胀[16]等步骤，对微流控芯片上被捕获的酵母单细胞进行母细胞区域分割和截面面积计算，并进一步探究了酵母细胞衰老过程中截面面积与其复制寿命的关联性。

1 微流控芯片和图像采集

本工作实验中采用的微流控芯片设有1 100个“瓶颈式”微结构，呈交错阵列式排布如图1(a)，每行22个微结构，共50行。如图1(b)所示，两侧对称的微柱形成向上游开口的“瓶身”，整个捕获单元宽15 μm，深7 μm，高8 μm[14]，中间形成3 μm宽的狭窄孔口即“瓶颈”。阵列中的捕获单元横向间距为34 μm，纵向间距为30 μm，相邻行横向间距为17 μm。利用流体动力原理实现酵母单细胞在捕获单元处的稳定捕获。如图1(c)所示，带芽细胞经过旋转使小芽稳定固定在下游瓶颈处，成熟后的子细胞能够被流体力剪切去除。实验中，每隔10 min对细胞进行一次明场显微成像，记录酵母细胞衰老过程中的形态变化，整个实验过程持续50 h左右。

图1 微流控芯片

2 酵母细胞面积提取算法

细胞图像分割是图像处理领域中的一大难点。由于细胞形态各异，且细微结构复杂，容易出现过度分割和分割不准确等问题。本文中的微流控芯片上各个捕获单元分别对应单个被捕获细胞，为实现酵母复制寿命检测，被捕获母细胞产生的子细胞成熟后即被流体剪切去除，因此图像主要由捕获单元两侧的微柱和被捕获细胞构成。捕获单元的形状规则，易于通过传统的Hough变换方法检测其边缘直线，准确率较高且运行速度快。然后利用被捕获细胞内部区域灰度的相似性，及细胞边缘与内部区域灰度的不连续性，以捕获位点为种子点，分割出细胞区域。结合酵母细胞近似椭球形的形态特点，利用区域圆形度进行预判，最后进行子－母细胞区域的分割。该图像处理算法结合了微流控芯片的设计，且经过实验验证，效果较好，改进后可用于相似的微流控芯片上细胞图像的形态参数提取。

面积提取算法的流程如图2所示，细胞阵列时序图像经过取反和二值化，进行Hough变换，进行旋转和平移的仿射变换完成配准。配准后的阵列图像分割成为单个捕获单元图像，再次进行Hough变换确定被捕获细胞的中心位置，通过魔棒函数分割出细胞区域。然后结合区域圆形度利用分水岭算法继续分割图像得到母细胞，最后进行形态学膨胀，得到母细胞ROI(Region of interest)掩模。整个过程主要分为时序图像的配准、中心－魔棒法分割细胞区域、圆形度－分水岭分割母细胞区域三个步骤，下面分别阐述其详细步骤和具体原理。

图2 母细胞截面面积提取的算法流程图

2.1 时序图像的配准

由于载物台上芯片放置造成显微图像角度存在偏移，为方便后续直线检测，需要对全体时序图像进行水平旋转，即角度配准。首先，对图像取反，并用Ostu算法即最大类间方差阈值选择法二值化[17]。接着，对二值图像进行Hough变换[15]，检测微柱边缘直线，如图3(a)所示，根据直线的倾斜角确定逆时针旋转角θ＝π/2－β，对图像进行旋转仿射变换，完成角度配准。旋转的坐标变换公式如下:

式中:(v，w)是原图像的像素坐标，(x，y)是旋转变换后图像的像素坐标。

此外，由于时序显微成像持续时间长，实验过程中存在一定水平漂移(在图3(b)所示的xy平面内)，需要以目标捕获单元为约束点对图像进行位移配准。选取最邻近原点处的捕获单元作为位移配准的基准，通过Hough变换检测其下边缘直线，将此直线中点作为不同帧图像的统一约束点。计算参考帧与待配准帧约束点的坐标差，作为待配准帧的平移坐标矢量，进行平移变换，实现位移配准。如图3(b)所示，参考帧约束点坐标记为(x，y)，待配准帧约束点坐标记为(x′，y′)，则待配准帧的平移矢量为A＝(x－x′，y－y′)。为便于后续跟踪处理，将目标位置以矩阵形式记录并编号，将捕获单元的图像逐个分割后对应保存。

图3 时序图像配准示意图

2.2 中心－魔棒法分割细胞区域

基于前述通过检测捕获单元确定被捕获细胞位置的思路，我们提出了中心－魔棒法用以细胞分割。由于原始图像分辨率较低，为削弱锯齿效应并减小区域分割误差，首先通过双三次插值[18]将图像尺寸扩大两倍；然后对图像进行取反和二值化，并做Hough变换检测微柱的左右边缘直线；再通过两侧直线的位置确定母细胞的捕获位点。以捕获位点作为种子点，应用魔棒函数[19]提取细胞区域。类似于Photoshop中的魔棒工具，魔棒函数用于提取包含指定种子点的连通区域，且该区域内所有像素点灰度值与种子点处灰度值的差值不大于指定的容差。令U表示整幅图像像素占据的空间区域，Rc表示魔棒函数提取的空间区域，(x0，y0)表示种子点坐标，Δ表示灰度值容差，满足:①Ri⊆U，i＝1，2，…，n；②Ri是一个连通域；③∀(x，y)∈Ri，｜f(x，y)－f(x0，y0)｜≤Δ；

由于成像条件和细胞形态差异，细胞内部和细胞边界处的像素差值不固定，因此容差需具有自适应性。具体地，我们逐渐增大容差值，计算每次迭代时魔棒提取区域的面积增加量，大于细胞平均面积则判断已经越过细胞边界，取上一次迭代时的容差，这样使得提取的细胞区域达到最大。将背景像素记为0，细胞区域记为1，此二值图像作为细胞分割的掩模。

图4 中心－魔棒法分割细胞区域示意图

2.3 圆形度－分水岭算法分割母细胞区域

由于原始图像中母细胞和微柱或子细胞存在交叠，提取的细胞区域除了包含母细胞，还可能包含子细胞和微柱。为避免过度分割，需对所提取区域进行圆形度预判。由于区域圆形度取值范围为(0，1)，为了增加判定区间，我们取圆形度的倒数记为C，其取值范围为(1，＋∞)。因此，规定C取值越接近1，则区域越近似圆形，计算公式如下:

式中:C为圆形度倒数，p为区域周长，S为区域面积。

考虑样本中母细胞拖拽变形和与子细胞交叠的情况，确定圆形度倒数C的经验阈值为1.7。小于该阈值则表示区域近似圆形，无需分割；大于阈值则使用分水岭算法分割。将掩模图像取反，使目标区域值为0，背景区域值为1，如图5(b)所示，做出其欧氏距离变换图[20]。再利用分水岭算法分割[21]，如图5(c)所示，得到分割完成的标记矩阵；如图5(d)所示，只保留包含种子点的子区域作为母细胞。然后对掩模进行孔洞填充，补全细胞内部区域。最后对提取的母细胞进行形态学膨胀，补偿轮廓边缘。

图5 圆形度－分水岭算法分割母细胞区域示意图

3 结果与讨论

3.1 单个细胞图像处理结果

算法对实验图像适用性良好，按照预期依次完成了配准、分割和面积的提取。如图6所示为典型的细胞图像处理过程:图6(a)为实验获取的原始灰度图像，首先检测到微柱两侧直线并确定中心点位置，对种子点应用魔棒函数分割细胞，接着计算区域圆形度的倒数C(C为1.96大于经验阈值1.7)，然后利用分水岭算法分割母细胞，最后进行形态学膨胀补偿轮廓，将母细胞区域叠加在原图上并标记像素面积。

图6 单个细胞图像处理及母细胞截面面积计算结果

3.2 多个细胞的面积变化趋势

计算一个世代内母细胞面积的平均值，如图7所示是其随世代增长的变化曲线。可以看到，虽然母细胞的截面面积局部存在下降和波动，但其整体呈上升趋势。如图7(a)展示了两个细胞生命周期内的面积－世代曲线及其标准差，其中面积增长较快的一个母细胞复制寿命较短，而另一个面积增长慢的母细胞复制寿命较长。如图7(b)所示，汇总十个细胞的面积变化曲线(A1~A10)，其中酵母细胞的复制寿命最短15代，最长31代。如图所示，单个酵母细胞截面面积约在第10代之后随世代增加呈明显增大的趋势，且在最后数代增长速率明显加快。研究结果表明，细胞的尺寸与酵母衰老具有一定相关性，是限制酵母寿命的潜在因素之一，与前期研究结论一致[22]。

图7 酵母细胞的截面面积－世代曲线

4 总结

用于酵母细胞图像分析的已有算法或软件特异性较强，主要适用于背景均一、细胞成片生长或荧光蛋白图像等情况。对于本文实验中用于酵母单细胞捕获－培养－解剖微流控芯片上的细胞时序图像，存在效率低、适用性差等问题。本文针对“瓶颈式”微结构阵列的微流控芯片结构，设计了完整的分析处理算法，依次实现了时序图像的角度及位移配准、单元区域分割，母细胞区域分割以及截面面积计算，并探究了酵母细胞复制衰老过程中的截面面积的动态变化及其与复制寿命的关联性，并由已有的生物学结论进行了佐证。

本文提出的酵母母细胞截面面积提取算法推动了基于微流控芯片的高通量酵母衰老寿命研究中时序图像的自动化分析处理，为后续计算细胞多种形态参数及细胞生长速率、并分析这些参数与酵母复制寿命的影响机制奠定了基础。进一步地，本文工作为利用传统图像处理方法分析单细胞图像提供了新的思路，对与相似微流控单细胞图像的自动化处理也具有重要意义。