釀酒酵母單細胞形態參數精準提取算法的研究?

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

(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 總結

用于酵母細胞圖像分析的已有算法或軟件特異性較強，主要適用于背景均一、細胞成片生長或熒光蛋白圖像等情況。對于本文實驗中用于酵母單細胞捕獲－培養－解剖微流控芯片上的細胞時序圖像，存在效率低、適用性差等問題。本文針對“瓶頸式”微結構陣列的微流控芯片結構，設計了完整的分析處理算法，依次實現了時序圖像的角度及位移配準、單元區域分割，母細胞區域分割以及截面面積計算，并探究了酵母細胞復制衰老過程中的截面面積的動態變化及其與復制壽命的關聯性，并由已有的生物學結論進行了佐證。

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