釀酒酵母單細(xì)胞形態(tài)參數(shù)精準(zhǔn)提取算法的研究?

劉 可王穎瀛耿楊燁肖 秦朱 真?

(1.東南大學(xué)電子科學(xué)與工程學(xué)院，江蘇 南京210096；2.東南大學(xué)微電子學(xué)院，江蘇南京210096)

細(xì)胞衰老是一個(gè)復(fù)雜的生理過程，伴隨著細(xì)胞生理結(jié)構(gòu)和功能的逐漸喪失。釀酒酵母(Saccharomyces cerevisiae)作為常用的模式細(xì)胞，具有繁殖快、壽命短、培養(yǎng)成本低等優(yōu)點(diǎn)。此外，釀酒酵母細(xì)胞與哺乳動(dòng)物細(xì)胞、人類具有相似的保守生化機(jī)制，其存活率曲線與人類壽命類似[1]，常被用于生命體衰老機(jī)制的研究[2]。釀酒酵母的復(fù)制壽命(Replicative lifespan，RLS)是酵母細(xì)胞衰老研究的主要指標(biāo)，定義為一個(gè)酵母細(xì)胞從新生到死亡過程中產(chǎn)生子代細(xì)胞的個(gè)數(shù)[3]。研究表明，在酵母細(xì)胞的生長、增殖、衰老過程中，細(xì)胞的截面面積、體積等形態(tài)參數(shù)變化用于表征細(xì)胞生長速率，并且與細(xì)胞周期調(diào)控和細(xì)胞復(fù)制壽命緊密相關(guān)[4]，但相關(guān)機(jī)理尚未被完全揭示。因此，需要進(jìn)一步研究細(xì)胞表型特征參數(shù)及其在細(xì)胞完整生命周期內(nèi)的動(dòng)態(tài)變化。微流控芯片利用微結(jié)構(gòu)捕獲并固定酵母單細(xì)胞，并通過流體剪切力實(shí)現(xiàn)子細(xì)胞的高效去除[2]，是酵母細(xì)胞復(fù)制衰老研究的重要平臺(tái)。在酵母細(xì)胞實(shí)驗(yàn)中，實(shí)驗(yàn)者可以借助微流控芯片捕獲、固定、培養(yǎng)酵母母細(xì)胞、并去除子代細(xì)胞、對細(xì)胞進(jìn)行時(shí)序顯微成像。釀酒酵母細(xì)胞平均復(fù)制壽命為20代~30代[5]，每代持續(xù)時(shí)間大約為90 min，若以10 min為時(shí)序成像時(shí)間間隔，單個(gè)酵母細(xì)胞完整生命周期的時(shí)序監(jiān)測將產(chǎn)生數(shù)百幀顯微圖像。對于具有多個(gè)單細(xì)胞捕獲結(jié)構(gòu)及成像位點(diǎn)的高通量陣列式微流控細(xì)胞培養(yǎng)芯片，單次實(shí)驗(yàn)需分析上萬幀圖像，因此急需提出一種精確高效、自動(dòng)化的細(xì)胞圖像處理算法，以提取酵母細(xì)胞的形態(tài)參數(shù)。

目前已經(jīng)有幾種成熟的算法和軟件用于細(xì)胞的圖像分割和目標(biāo)跟蹤[6]。主流算法包括主動(dòng)輪廓算法、分水嶺算法等傳統(tǒng)圖像分割算法的組合及卷積神經(jīng)網(wǎng)絡(luò)等[7－8]。細(xì)胞圖像處理的相關(guān)軟件包括ImageJ、CellProfiler、CellStar、CellX等。ImageJ使用最普遍，該軟件是基于Java開發(fā)的圖像處理軟件，具有開放式結(jié)構(gòu)和可擴(kuò)展性，支持一些用戶編寫的插件或宏如BudJ[9]，其優(yōu)點(diǎn)是步驟靈活，但是難以實(shí)現(xiàn)圖像的批量處理。CellProfiler是用于識別和量化細(xì)胞表型參數(shù)的相對完善的圖像分析軟件，用戶無需使用編程語言就能通過交互界面自行設(shè)計(jì)不同的算法流程，然而其更適用于高分辨率、視野中存在多個(gè)細(xì)胞的圖像，且不適合處理延時(shí)圖像[10]。Cristian等人提出了適用于明場酵母細(xì)胞圖像的處理軟件CellStar，并建立了酵母細(xì)胞圖像處理軟件的評估平臺(tái)，綜合比較了各種算法或軟件的適用范圍和效果，且平臺(tái)仍在持續(xù)更新[6]。大部分已有的算法和軟件要求圖像具有均一、明亮的背景，且多數(shù)針對細(xì)胞成片生長的情況進(jìn)行細(xì)胞分割?；跈C(jī)器學(xué)習(xí)的算法相較于傳統(tǒng)方法雖然具有更強(qiáng)的通用性，但其復(fù)雜度更高，運(yùn)行和調(diào)試成本也更高[11]。

然而，微流控芯片上細(xì)胞捕獲微結(jié)構(gòu)的幾何輪廓造成了圖像的背景不均一，且母細(xì)胞可能與微結(jié)構(gòu)邊緣產(chǎn)生交疊，影響細(xì)胞輪廓的判別[12－13]，已有的圖像處理軟件不適用。由于微流控芯片結(jié)構(gòu)和成像條件的差異，不同微流控芯片上采集的實(shí)驗(yàn)圖像通常需要設(shè)計(jì)專用的算法進(jìn)行分析處理。本文中的高通量的酵母單細(xì)胞捕獲－培養(yǎng)－解剖微流控芯片上設(shè)有陣列式的捕獲微結(jié)構(gòu)[14]，因此采用通過檢測形狀規(guī)則的細(xì)胞捕獲微結(jié)構(gòu)確定被捕獲細(xì)胞的位置作為圖像分割的思路。我們針對此芯片設(shè)計(jì)了專用的圖像處理算法，該算法結(jié)合Hough變換[15]、分水嶺算法、形態(tài)學(xué)膨脹[16]等步驟，對微流控芯片上被捕獲的酵母單細(xì)胞進(jìn)行母細(xì)胞區(qū)域分割和截面面積計(jì)算，并進(jìn)一步探究了酵母細(xì)胞衰老過程中截面面積與其復(fù)制壽命的關(guān)聯(lián)性。

1 微流控芯片和圖像采集

本工作實(shí)驗(yàn)中采用的微流控芯片設(shè)有1 100個(gè)“瓶頸式”微結(jié)構(gòu)，呈交錯(cuò)陣列式排布如圖1(a)，每行22個(gè)微結(jié)構(gòu)，共50行。如圖1(b)所示，兩側(cè)對稱的微柱形成向上游開口的“瓶身”，整個(gè)捕獲單元寬15 μm，深7 μm，高8 μm[14]，中間形成3 μm寬的狹窄孔口即“瓶頸”。陣列中的捕獲單元橫向間距為34 μm，縱向間距為30 μm，相鄰行橫向間距為17 μm。利用流體動(dòng)力原理實(shí)現(xiàn)酵母單細(xì)胞在捕獲單元處的穩(wěn)定捕獲。如圖1(c)所示，帶芽細(xì)胞經(jīng)過旋轉(zhuǎn)使小芽穩(wěn)定固定在下游瓶頸處，成熟后的子細(xì)胞能夠被流體力剪切去除。實(shí)驗(yàn)中，每隔10 min對細(xì)胞進(jìn)行一次明場顯微成像，記錄酵母細(xì)胞衰老過程中的形態(tài)變化，整個(gè)實(shí)驗(yàn)過程持續(xù)50 h左右。

圖1 微流控芯片

2 酵母細(xì)胞面積提取算法

細(xì)胞圖像分割是圖像處理領(lǐng)域中的一大難點(diǎn)。由于細(xì)胞形態(tài)各異，且細(xì)微結(jié)構(gòu)復(fù)雜，容易出現(xiàn)過度分割和分割不準(zhǔn)確等問題。本文中的微流控芯片上各個(gè)捕獲單元分別對應(yīng)單個(gè)被捕獲細(xì)胞，為實(shí)現(xiàn)酵母復(fù)制壽命檢測，被捕獲母細(xì)胞產(chǎn)生的子細(xì)胞成熟后即被流體剪切去除，因此圖像主要由捕獲單元兩側(cè)的微柱和被捕獲細(xì)胞構(gòu)成。捕獲單元的形狀規(guī)則，易于通過傳統(tǒng)的Hough變換方法檢測其邊緣直線，準(zhǔn)確率較高且運(yùn)行速度快。然后利用被捕獲細(xì)胞內(nèi)部區(qū)域灰度的相似性，及細(xì)胞邊緣與內(nèi)部區(qū)域灰度的不連續(xù)性，以捕獲位點(diǎn)為種子點(diǎn)，分割出細(xì)胞區(qū)域。結(jié)合酵母細(xì)胞近似橢球形的形態(tài)特點(diǎn)，利用區(qū)域圓形度進(jìn)行預(yù)判，最后進(jìn)行子－母細(xì)胞區(qū)域的分割。該圖像處理算法結(jié)合了微流控芯片的設(shè)計(jì)，且經(jīng)過實(shí)驗(yàn)驗(yàn)證，效果較好，改進(jìn)后可用于相似的微流控芯片上細(xì)胞圖像的形態(tài)參數(shù)提取。

面積提取算法的流程如圖2所示，細(xì)胞陣列時(shí)序圖像經(jīng)過取反和二值化，進(jìn)行Hough變換，進(jìn)行旋轉(zhuǎn)和平移的仿射變換完成配準(zhǔn)。配準(zhǔn)后的陣列圖像分割成為單個(gè)捕獲單元圖像，再次進(jìn)行Hough變換確定被捕獲細(xì)胞的中心位置，通過魔棒函數(shù)分割出細(xì)胞區(qū)域。然后結(jié)合區(qū)域圓形度利用分水嶺算法繼續(xù)分割圖像得到母細(xì)胞，最后進(jìn)行形態(tài)學(xué)膨脹，得到母細(xì)胞ROI(Region of interest)掩模。整個(gè)過程主要分為時(shí)序圖像的配準(zhǔn)、中心－魔棒法分割細(xì)胞區(qū)域、圓形度－分水嶺分割母細(xì)胞區(qū)域三個(gè)步驟，下面分別闡述其詳細(xì)步驟和具體原理。

圖2 母細(xì)胞截面面積提取的算法流程圖

2.1 時(shí)序圖像的配準(zhǔn)

由于載物臺(tái)上芯片放置造成顯微圖像角度存在偏移，為方便后續(xù)直線檢測，需要對全體時(shí)序圖像進(jìn)行水平旋轉(zhuǎn)，即角度配準(zhǔn)。首先，對圖像取反，并用Ostu算法即最大類間方差閾值選擇法二值化[17]。接著，對二值圖像進(jìn)行Hough變換[15]，檢測微柱邊緣直線，如圖3(a)所示，根據(jù)直線的傾斜角確定逆時(shí)針旋轉(zhuǎn)角θ＝π/2－β，對圖像進(jìn)行旋轉(zhuǎn)仿射變換，完成角度配準(zhǔn)。旋轉(zhuǎn)的坐標(biāo)變換公式如下:

式中:(v，w)是原圖像的像素坐標(biāo)，(x，y)是旋轉(zhuǎn)變換后圖像的像素坐標(biāo)。

此外，由于時(shí)序顯微成像持續(xù)時(shí)間長，實(shí)驗(yàn)過程中存在一定水平漂移(在圖3(b)所示的xy平面內(nèi))，需要以目標(biāo)捕獲單元為約束點(diǎn)對圖像進(jìn)行位移配準(zhǔn)。選取最鄰近原點(diǎn)處的捕獲單元作為位移配準(zhǔn)的基準(zhǔn)，通過Hough變換檢測其下邊緣直線，將此直線中點(diǎn)作為不同幀圖像的統(tǒng)一約束點(diǎn)。計(jì)算參考幀與待配準(zhǔn)幀約束點(diǎn)的坐標(biāo)差，作為待配準(zhǔn)幀的平移坐標(biāo)矢量，進(jìn)行平移變換，實(shí)現(xiàn)位移配準(zhǔn)。如圖3(b)所示，參考幀約束點(diǎn)坐標(biāo)記為(x，y)，待配準(zhǔn)幀約束點(diǎn)坐標(biāo)記為(x′，y′)，則待配準(zhǔn)幀的平移矢量為A＝(x－x′，y－y′)。為便于后續(xù)跟蹤處理，將目標(biāo)位置以矩陣形式記錄并編號，將捕獲單元的圖像逐個(gè)分割后對應(yīng)保存。

圖3 時(shí)序圖像配準(zhǔn)示意圖

2.2 中心－魔棒法分割細(xì)胞區(qū)域

基于前述通過檢測捕獲單元確定被捕獲細(xì)胞位置的思路，我們提出了中心－魔棒法用以細(xì)胞分割。由于原始圖像分辨率較低，為削弱鋸齒效應(yīng)并減小區(qū)域分割誤差，首先通過雙三次插值[18]將圖像尺寸擴(kuò)大兩倍；然后對圖像進(jìn)行取反和二值化，并做Hough變換檢測微柱的左右邊緣直線；再通過兩側(cè)直線的位置確定母細(xì)胞的捕獲位點(diǎn)。以捕獲位點(diǎn)作為種子點(diǎn)，應(yīng)用魔棒函數(shù)[19]提取細(xì)胞區(qū)域。類似于Photoshop中的魔棒工具，魔棒函數(shù)用于提取包含指定種子點(diǎn)的連通區(qū)域，且該區(qū)域內(nèi)所有像素點(diǎn)灰度值與種子點(diǎn)處灰度值的差值不大于指定的容差。令U表示整幅圖像像素占據(jù)的空間區(qū)域，Rc表示魔棒函數(shù)提取的空間區(qū)域，(x0，y0)表示種子點(diǎn)坐標(biāo)，Δ表示灰度值容差，滿足:①Ri?U，i＝1，2，…，n；②Ri是一個(gè)連通域；③?(x，y)∈Ri，｜f(x，y)－f(x0，y0)｜≤Δ；

由于成像條件和細(xì)胞形態(tài)差異，細(xì)胞內(nèi)部和細(xì)胞邊界處的像素差值不固定，因此容差需具有自適應(yīng)性。具體地，我們逐漸增大容差值，計(jì)算每次迭代時(shí)魔棒提取區(qū)域的面積增加量，大于細(xì)胞平均面積則判斷已經(jīng)越過細(xì)胞邊界，取上一次迭代時(shí)的容差，這樣使得提取的細(xì)胞區(qū)域達(dá)到最大。將背景像素記為0，細(xì)胞區(qū)域記為1，此二值圖像作為細(xì)胞分割的掩模。

圖4 中心－魔棒法分割細(xì)胞區(qū)域示意圖

2.3 圓形度－分水嶺算法分割母細(xì)胞區(qū)域

由于原始圖像中母細(xì)胞和微柱或子細(xì)胞存在交疊，提取的細(xì)胞區(qū)域除了包含母細(xì)胞，還可能包含子細(xì)胞和微柱。為避免過度分割，需對所提取區(qū)域進(jìn)行圓形度預(yù)判。由于區(qū)域圓形度取值范圍為(0，1)，為了增加判定區(qū)間，我們?nèi)A形度的倒數(shù)記為C，其取值范圍為(1，＋∞)。因此，規(guī)定C取值越接近1，則區(qū)域越近似圓形，計(jì)算公式如下:

式中:C為圓形度倒數(shù)，p為區(qū)域周長，S為區(qū)域面積。

考慮樣本中母細(xì)胞拖拽變形和與子細(xì)胞交疊的情況，確定圓形度倒數(shù)C的經(jīng)驗(yàn)閾值為1.7。小于該閾值則表示區(qū)域近似圓形，無需分割；大于閾值則使用分水嶺算法分割。將掩模圖像取反，使目標(biāo)區(qū)域值為0，背景區(qū)域值為1，如圖5(b)所示，做出其歐氏距離變換圖[20]。再利用分水嶺算法分割[21]，如圖5(c)所示，得到分割完成的標(biāo)記矩陣；如圖5(d)所示，只保留包含種子點(diǎn)的子區(qū)域作為母細(xì)胞。然后對掩模進(jìn)行孔洞填充，補(bǔ)全細(xì)胞內(nèi)部區(qū)域。最后對提取的母細(xì)胞進(jìn)行形態(tài)學(xué)膨脹，補(bǔ)償輪廓邊緣。

圖5 圓形度－分水嶺算法分割母細(xì)胞區(qū)域示意圖

3 結(jié)果與討論

3.1 單個(gè)細(xì)胞圖像處理結(jié)果

算法對實(shí)驗(yàn)圖像適用性良好，按照預(yù)期依次完成了配準(zhǔn)、分割和面積的提取。如圖6所示為典型的細(xì)胞圖像處理過程:圖6(a)為實(shí)驗(yàn)獲取的原始灰度圖像，首先檢測到微柱兩側(cè)直線并確定中心點(diǎn)位置，對種子點(diǎn)應(yīng)用魔棒函數(shù)分割細(xì)胞，接著計(jì)算區(qū)域圓形度的倒數(shù)C(C為1.96大于經(jīng)驗(yàn)閾值1.7)，然后利用分水嶺算法分割母細(xì)胞，最后進(jìn)行形態(tài)學(xué)膨脹補(bǔ)償輪廓，將母細(xì)胞區(qū)域疊加在原圖上并標(biāo)記像素面積。

圖6 單個(gè)細(xì)胞圖像處理及母細(xì)胞截面面積計(jì)算結(jié)果

3.2 多個(gè)細(xì)胞的面積變化趨勢

計(jì)算一個(gè)世代內(nèi)母細(xì)胞面積的平均值，如圖7所示是其隨世代增長的變化曲線?？梢钥吹剑m然母細(xì)胞的截面面積局部存在下降和波動(dòng)，但其整體呈上升趨勢。如圖7(a)展示了兩個(gè)細(xì)胞生命周期內(nèi)的面積－世代曲線及其標(biāo)準(zhǔn)差，其中面積增長較快的一個(gè)母細(xì)胞復(fù)制壽命較短，而另一個(gè)面積增長慢的母細(xì)胞復(fù)制壽命較長。如圖7(b)所示，匯總十個(gè)細(xì)胞的面積變化曲線(A1~A10)，其中酵母細(xì)胞的復(fù)制壽命最短15代，最長31代。如圖所示，單個(gè)酵母細(xì)胞截面面積約在第10代之后隨世代增加呈明顯增大的趨勢，且在最后數(shù)代增長速率明顯加快。研究結(jié)果表明，細(xì)胞的尺寸與酵母衰老具有一定相關(guān)性，是限制酵母壽命的潛在因素之一，與前期研究結(jié)論一致[22]。

圖7 酵母細(xì)胞的截面面積－世代曲線

4 總結(jié)

用于酵母細(xì)胞圖像分析的已有算法或軟件特異性較強(qiáng)，主要適用于背景均一、細(xì)胞成片生長或熒光蛋白圖像等情況。對于本文實(shí)驗(yàn)中用于酵母單細(xì)胞捕獲－培養(yǎng)－解剖微流控芯片上的細(xì)胞時(shí)序圖像，存在效率低、適用性差等問題。本文針對“瓶頸式”微結(jié)構(gòu)陣列的微流控芯片結(jié)構(gòu)，設(shè)計(jì)了完整的分析處理算法，依次實(shí)現(xiàn)了時(shí)序圖像的角度及位移配準(zhǔn)、單元區(qū)域分割，母細(xì)胞區(qū)域分割以及截面面積計(jì)算，并探究了酵母細(xì)胞復(fù)制衰老過程中的截面面積的動(dòng)態(tài)變化及其與復(fù)制壽命的關(guān)聯(lián)性，并由已有的生物學(xué)結(jié)論進(jìn)行了佐證。

本文提出的酵母母細(xì)胞截面面積提取算法推動(dòng)了基于微流控芯片的高通量酵母衰老壽命研究中時(shí)序圖像的自動(dòng)化分析處理，為后續(xù)計(jì)算細(xì)胞多種形態(tài)參數(shù)及細(xì)胞生長速率、并分析這些參數(shù)與酵母復(fù)制壽命的影響機(jī)制奠定了基礎(chǔ)。進(jìn)一步地，本文工作為利用傳統(tǒng)圖像處理方法分析單細(xì)胞圖像提供了新的思路，對與相似微流控單細(xì)胞圖像的自動(dòng)化處理也具有重要意義。