李樹力，李金澤，郭 振，3，朱文艷，3，周連群，3*，張芷齊，3*

（1. 上海大學(xué) 通信與信息工程學(xué)院，上海200444；2. 中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所 中國科學(xué)院生物醫(yī)學(xué)檢驗技術(shù)重點實驗室，江蘇 蘇州215163；3. 中國科學(xué)技術(shù)大學(xué)，安徽 合肥230026）

1 引 言

數(shù)字聚合酶鏈?zhǔn)椒磻?yīng)（digital Polymerase Chain Reaction，dPCR）技術(shù)在循環(huán)腫瘤DNA 檢測、稀有突變檢測和拷貝數(shù)變異分析等核酸定量分析中有著廣泛的應(yīng)用［1-2］。1999 年，Kinzler 和Vogelstein 等 提 出 了 數(shù) 字PCR 的 概 念［3］。 數(shù) 字PCR 實驗包含3 個環(huán)節(jié)，即樣本的分散、PCR 擴(kuò)增、熒光信號的采集與數(shù)據(jù)分析。 在數(shù)據(jù)分析時，利用泊松分布可以精確計算得到原始樣本的模板拷貝數(shù)［4-5］。影響最終結(jié)果的關(guān)鍵在于對熒光信號的采集與處理，即讀取各個反應(yīng)單元的熒光信號并判斷其陰性或陽性。

隨著光刻與MEMS 技術(shù)的發(fā)展［6］，芯片式數(shù)字PCR 得以廣泛應(yīng)用，解決了微升級樣品到納升級反應(yīng)單元的精準(zhǔn)化分問題［7］。 芯片式數(shù)字PCR 相比于液滴式數(shù)字PCR，圖像信息獲取更方便。數(shù)字PCR 芯片圖像的特征是由大量的形狀類似的樣點組成的陣列，它們之間的排列一般都是規(guī)則的，處理這種微陣列圖像一般有三個步驟：樣點定位、樣點分割和信號強(qiáng)度提?。?-10］。其中，樣點定位是微陣列圖像分析中最基本、最重要的步驟，準(zhǔn)確的樣點定位大大提高了后續(xù)分割與量化步驟的效率。Saeid 等［11］提出一種局部閾值網(wǎng)格化的方法。Nguyen 等［12］提出一種基于小波檢測、主動輪廓分割的方法。李鐵軍等［13］利用一種基于四階矩陣的對比度增強(qiáng)的方法，提高樣點定位的準(zhǔn)確性。Belean 等［14］提出一個基于偏微分方程的微陣列網(wǎng)格對齊的方法。牟穎等［15］使用最大類間方差法來提取信息。這些方法均基于投影法，容易實現(xiàn)，加以改進(jìn)也可以達(dá)到不錯的定位結(jié)果，但是它們只適用于矩形排列方式的微陣列芯片圖像。而近些年，出現(xiàn)了一些類似蜂窩狀的六邊形排布方式的數(shù)字PCR 芯片［16］。與矩形網(wǎng)格相比，蜂窩狀網(wǎng)格在芯片上的利用率更高，在相同的面積里，蜂窩狀的排布方式可以設(shè)計更多的樣點位。對于這種非正交排列的陣列圖像，投影法無法分隔所有的樣點。Genpix［17］和ScanAlyze［18］等軟件包使用固定模板的方法，預(yù)先設(shè)置了幾種不同排布方式的模板，分析數(shù)據(jù)前，讓用戶手動選擇，然后將模板對準(zhǔn)嵌套在采集到的芯片圖像上，利用設(shè)定好的坐標(biāo)定位每個樣點的位置。這種方法的局限性在于：圖片處理中需要人工參與，無法處理和預(yù)設(shè)模板排列方式不同的芯片。 Giannakeas 等［19］提出了一個生長同心六邊形算法，利用Voronoi 圖［20］對樣點定位。這種方法可以處理蜂窩狀排布方式的芯片圖像，但是由于是逐點定位，計算量過于龐大，不適用于處理數(shù)據(jù)量大的數(shù)字PCR 芯片圖像［21-22］。由相機(jī)拍攝采集到的芯片圖像除了存在光照分布不均勻的問題，其樣點的信號強(qiáng)度一般也非常低，這導(dǎo)致圖像整體灰度值集中在一個比較低的范圍，信號與背景的對比度低。如果不進(jìn)行特定的運算處理，直接采用全局固定閾值分割的方法，無法準(zhǔn)確地分割出每個有效樣點。有一些方法對圖像對比度進(jìn)行了改進(jìn)［23］，其中最簡單和有用的是自適應(yīng)直方圖均衡化方法［24］。然而，這些方法都試圖將原始圖像分割成小塊并調(diào)節(jié)增強(qiáng)程度，必然涉及到參數(shù)設(shè)置和閾值搜索，信號和亮度會被破壞。

針對目前任意排布方式的數(shù)字PCR 芯片，本文提出基于形態(tài)學(xué)的三通道蜂窩狀熒光圖像樣點尋址算法及數(shù)字PCR 圖像信息提取方法。該方法包含三通道圖像配準(zhǔn)、樣點的尋址定位（Automatic Hexagon Microarray Addressing Algorithm，AHMAA）、樣點分割、信息提取以及對樣點的篩選，能夠快速、準(zhǔn)確地獲取芯片多個通道的信息。

2 實驗材料

實驗采用Applied Biosystems?QuantStudio ?3D的數(shù)字PCR 芯片［25］。該芯片是一種基于硅基基片的表面疏水孔壁親水的毛細(xì)管微陣列芯片，使用硅基底板制造并在其上蝕刻有約20 000 個大小均一的微米級反應(yīng)孔。芯片上的微孔以蜂窩狀網(wǎng)格的方式排列，微孔直徑為60 μm，孔間隔15 μm，結(jié)構(gòu)如圖1 所示。使用數(shù)字PCR 檢測系統(tǒng)在明場下拍攝芯片整體圖像，芯片尺寸為10 mm×10 mm，其中局部放大部分是使用掃描電鏡拍攝的芯片圖像，可以清晰地觀察到六邊形結(jié)構(gòu)。

PCR 實驗使用到的試劑主要有質(zhì)粒DNA（生工生物工程股份有限公司），TaqManTMLiquid Biopsy dPCR Assay（Thermofisher），Quant-Studio ?3D Digital PCR Master Mix（Thermofisher）等。

芯片圖像由本實驗室自主研制的數(shù)字PCR檢測系統(tǒng)采集。本系統(tǒng)主要基于CMOS 成像原理，有多個熒光通道，可同時采集多個通道的熒光信號。本文使用的數(shù)據(jù)是5 塊不同的PCR 芯片 ，分 別 在FAM ，VIC ，ROX（校 正）3 個 通 道采集圖片，每張圖片的規(guī)格為8 位2 048×2 448 pixel。

圖1 數(shù)字PCR 芯片結(jié)構(gòu)Fig. 1 Structure of digital PCR chip

3 圖像處理

為了快速準(zhǔn)確地處理三通道蜂窩狀堆疊的微孔陣列芯片圖像，本文提出基于形態(tài)學(xué)的三通道蜂窩狀熒光圖像樣點尋址算法及數(shù)字PCR 圖像信息提取方法。圖2 展示了本文的整個數(shù)字PCR 芯片圖像處理流程。

該方法主要分為四個步驟，第一步：三通道的圖像配準(zhǔn)，針對不同通道的圖像進(jìn)行配準(zhǔn)融合，使樣點排布整齊；第二步：樣點尋址（AHMAA），通過增強(qiáng)圖像的對比度，選取有效樣點區(qū)域，基于改進(jìn)伽馬算法去除光照分布不均效應(yīng)，基于形態(tài)學(xué)算法對緊密排列樣點進(jìn)行識別，對微陣列芯片圖像進(jìn)行單元定位；第三步：樣點分割，以樣點坐標(biāo)為中心，選擇符合樣點結(jié)構(gòu)形狀的區(qū)域作為信息提取的區(qū)域；第四步：樣點的信息提取，提取每個樣點的生物分子熒光信息。

圖2 數(shù)字PCR 微陣列熒光圖像信息提取流程Fig. 2 Flow chart for information extraction of digital PCR microarray fluorescence image

3.1 三通道圖像配準(zhǔn)

實驗室研發(fā)的數(shù)字PCR 圖像采集系統(tǒng)可以同時采集芯片在3 個不同波段的熒光圖像，為了同時獲取多個通道的圖片信息，本文將不同通道的圖片進(jìn)行融合處理。由于圖像采集系統(tǒng)中存在光學(xué)誤差，同一個樣點在不同通道里的位置并不是完全重合的，有幾個像素的偏差，進(jìn)而導(dǎo)致后續(xù)獲取樣點數(shù)據(jù)的不準(zhǔn)確。為了解決這一問題，本文采用了基于互信息的多模態(tài)圖像配準(zhǔn)方法，將不同通道的樣點對齊?；バ畔⒘慷仁怯糜跍y量兩個變量之間的相關(guān)性的信息理論技術(shù)。算法使用來自兩個圖像的像素采樣的聯(lián)合概率分布來測量一組像素的值映射到另一圖像中的相似值的確定性，該信息是圖像相似程度的定量度量。高互信息意味著兩個分布之間的不確定性（熵）大大降低，表明圖像更好地被配準(zhǔn)［26］。

3.2 樣點尋址

圖像配準(zhǔn)完成后，得到一張包含多個通道信息的融合圖片，由于多個通道的圖像已經(jīng)配準(zhǔn)，所以對其中一個通道的圖像做尋址，得到樣點坐標(biāo)可以直接用于其他通道的圖像。為了準(zhǔn)確地定位樣點，針對芯片圖像存在的問題，本文對圖像進(jìn)行了以下幾步處理：

（1）提高圖像整體對比度，選擇ROI；

（2）去除光照不均影響；

（3）增強(qiáng)樣點孔內(nèi)外的對比度、二值化；

（4）計算統(tǒng)計樣點的位置信息。

圖像采集系統(tǒng)采集到的芯片圖像，除了芯片部分，還包含了芯片載具的部分，需要將芯片以外的部分都去掉。本文在提取邊緣之前，先對圖像直方圖進(jìn)行均值化處理，增強(qiáng)芯片區(qū)域與背景之間的對比度；然后，模糊處理，去除芯片區(qū)域內(nèi)細(xì)小邊緣的干擾；最后邊緣提取，得到準(zhǔn)確的外邊界。

針對上一步提取的芯片主體區(qū)域，本文使用基于形態(tài)學(xué)的方法，去除芯片內(nèi)部的連通區(qū)域；采用多次提取邊緣并減去的方法，將整個芯片主體區(qū)域的邊緣部分去掉。最終得到只含有有效樣點的區(qū)域，作為下一步處理的感興趣區(qū)域（Region of Interest，ROI）。

在采集芯片圖像的過程中，由于受到激發(fā)光源大小、鏡頭視野等外部因素的影響，芯片上的光照不均勻，主要表現(xiàn)為芯片中間區(qū)域的光照強(qiáng)度比芯片邊緣區(qū)域的光照強(qiáng)度更高［27］。這會導(dǎo)致樣點的分割結(jié)果不理想，陰陽性樣點的分類結(jié)果不準(zhǔn)確等問題。

本文使用改進(jìn)的基于二維γ 函數(shù)的光照校正算法來校正芯片的光照分布。該算法的主要思想是提取圖片的光照分量信息，通過校正使得光照過亮的區(qū)域整體亮度降低，光照過暗的地方亮度得到提高，提升整塊芯片的光照亮度均勻性。為了防止直接在RGB 3 個通道做校正導(dǎo)致圖像色彩失真，選擇在HSV 色彩空間進(jìn)行處理［26］，其中色調(diào)（H）、飽和度（S）、亮度（V）是相互獨立的。

算法主要步驟如下：

（1）輸入芯片圖像IRGB，將圖像從RGB 空間轉(zhuǎn)化到HSV 空間IHSV，單獨對亮度（V）處理；

（2）多尺度高斯卷積，賦予不同尺度不同的權(quán)值，得到芯片圖像的光照分量（gaus）；

（3）根據(jù)二維γ函數(shù)公式校正，得到校正后的亮度（V ′）；

（4）將校正后的亮度（V ′）與原圖的色調(diào)（H）、飽和度（S）重新合成，得到新的圖像I ′HSV；

（5）將圖像從HSV 空間轉(zhuǎn)化到RGB 空間，得到I ′RGB。

二維γ函數(shù)公式為：

為了增強(qiáng)樣點與周圍環(huán)境，即芯片上孔內(nèi)外的對比度，本文提出一種基于頂帽變換與底帽變換增強(qiáng)對比度的算法。頂帽變換與底帽變換都是基于數(shù)學(xué)形態(tài)學(xué)的圖像處理方法［28］。

灰度圖像f的頂帽變換定義為f減去其開運算結(jié)果，即：

在設(shè)置合適尺寸的構(gòu)造元素后，可以得到比較理想的孔內(nèi)與孔外的圖像信息。在得到頂帽與底帽變換之后的圖像后，將二者做差，可以大幅增強(qiáng)樣點孔內(nèi)與孔外的對比度。這種方法可以保證只有樣點被增強(qiáng)，而不受噪聲的干擾。

對上一步結(jié)果進(jìn)行二值化處理，去掉面積過小與過大的連通域，得到一幅包含所有有效樣點的二值化點陣。計算統(tǒng)計所有連通域的中心點的坐標(biāo)，將它作為樣點的位置信息，至此樣點的尋址定位步驟完成。

3.3 樣點分割

理論上，微孔陣列芯片的每個樣點孔都有固定的空間結(jié)構(gòu)和形狀，而使用AHMAA 算法得到的定位結(jié)果是每個樣點的中心點。當(dāng)我們得到了每個樣點的中心坐標(biāo)，則以定位點坐標(biāo)為中心，劃取一個樣點大小的正六邊形區(qū)域作為信息采集區(qū)域。這種直接分割的方法比先網(wǎng)格化定位，再對所有網(wǎng)格小碎片里的樣點分割的方法更快速、準(zhǔn)確。

3.4 信息提取

計算每個六邊形區(qū)域內(nèi)所有像素點的均值，作為樣點的熒光強(qiáng)度值。統(tǒng)計得到所有樣點的熒光信號值后，本文采用基于K 均值的分類方法對數(shù)據(jù)分類，得到芯片的陰陽性點數(shù)目信息。

4 結(jié)果與討論

圖3 是通過AHMAA 算法最終定位的結(jié)果。圖中視野內(nèi)所有樣點無論是陽性還是陰性，都被成功地識別并準(zhǔn)確地分割出來。

圖3 樣點定位結(jié)果Fig. 3 Result of spot addressing

圖4 展示了熒光通道1 與熒光通道2（FAMvs. VIC）的二維聚類散點圖。根據(jù)計算得出的閾值將圖像分為4 個象限，并將樣點分為4 類，分別代表FAM 陽性與VIC 陰性、FAM 陽性與VIC陽 性 、FAM 陰 性 與VIC 陰 性 、FAM 陰 性 與VIC陽性。

圖4 二維聚類散點Fig. 4 Two-dimensional clustering scatter plots

4.1 三通道圖像配準(zhǔn)

由于沒有一個準(zhǔn)確描述兩幅圖像對齊性的量，因此很難對配準(zhǔn)結(jié)果進(jìn)行定量分析。通常通過可視化的結(jié)果來定性地判斷配準(zhǔn)效果。圖5 顯示了三通道融合圖像配準(zhǔn)前后的效果（彩圖見期刊電子版）。在配準(zhǔn)前，由紅色和綠色代表的兩張不同通道的圖片有明顯的錯位現(xiàn)象，導(dǎo)致同一個樣點在不同通道圖片上的坐標(biāo)位置不同；配準(zhǔn)后，綠色的圓圈與紅色的圓重合了。

4.2 對比度增強(qiáng)的作用

為了提取芯片的主體區(qū)域，需要準(zhǔn)確地識別芯片的邊緣范圍，直接使用傳統(tǒng)的邊緣提取方法，會受到大量的樣點邊緣信息干擾。所以，需要先提高圖像的整體對比度。

圖5 三通道芯片的圖像配準(zhǔn)Fig. 5 Image registration of multi-channel chip

圖6 展示了整體對比度增強(qiáng)前后，圖片的灰度直方圖分布。原圖的灰度值集中在0～100，圖片整體亮度非常暗，而且?guī)缀醴直娌怀鲂酒瑓^(qū)域的邊界；處理后圖片的灰度值均勻地分布在0～255，邊界非常明顯，為后續(xù)ROI 的提取打下了很好的基礎(chǔ)。

4.3 光照不均勻?qū)Y(jié)果的影響

在實際實驗中，經(jīng)常遇到一些光照不均勻的圖片，如果不加以處理，則會對后續(xù)樣點的識別和定位有很大影響；除此之外，光照不均也會對獲取到的樣點信號值有影響，導(dǎo)致亮度高的區(qū)域信號值偏高，影響陰陽性分類結(jié)果。

圖7 顯示了去除光照不均影響前后的差別，圖7（a）是待處理的圖片，圖7（b）是沒有經(jīng)過處理的圖像分割結(jié)果，可以從放大的局部細(xì)節(jié)7（e），7（f）中看到，在原圖中過暗或過亮的區(qū)域，都無法很好地分割和識別樣點。圖7（c）是圖7（a）的光照分量圖，在做了去除光照不均的處理后，分割結(jié)果如圖7（d）所示，可以觀察到局部細(xì)節(jié)圖7（g）和圖7（h），分割結(jié)果得到了很好的改善。

4.4 去除無效樣點的必要性

圖6 灰度直方圖Fig. 6 Gray histograms

圖7 去除光照不均影響前后圖像的分割效果Fig. 7 Segmentation results before and after removing effects of uneven illumination

實驗中，芯片上經(jīng)常會出現(xiàn)如圖8 中的連通區(qū)塊。這導(dǎo)致大片的樣點呈陽性表征，這些點就是典型的偽陽性點，會對后續(xù)生物學(xué)統(tǒng)計結(jié)果產(chǎn)生巨大的誤差，所以在芯片圖像處理的過程中，需要去除掉這一部分的點。使用本文提出的算法可以有效地去除這些點。圖8 中顯示的是2 組圖片，其中，（a），（c）是采集到的原圖，（b），（d）是最終選擇的點的區(qū)域，可以清楚地看到，箭頭標(biāo)記的問題區(qū)域被識別并去除了。

4.5 定位時間及精度

圖8 去除無效的樣點Fig. 8 Removal of invalid spots

為了驗證算法的定位效率，本文在一塊完整的數(shù)字PCR 芯片上隨機(jī)截取4 組局部圖片，每組圖片分別有100 張，4 組圖片分別含有約10 個，100 個，1 000 個，10 000 個樣點。對每張圖片，分別用Voronoi 定位方法與本文提出的AHMAA定位算法進(jìn)行計算，得出每張圖片運算的時間，并對每組數(shù)據(jù)計算其均值與標(biāo)準(zhǔn)差，畫出圖9 的誤差條形圖。當(dāng)計算約有10 個樣點的圖片時，Voronoi 耗 時 約0. 297 s，AHMAA 耗 時 約0. 022 s；當(dāng)計算約有100 個樣點的圖片時，Voronoi 耗時約0. 968 s，AHMAA 耗時約0. 028 s；當(dāng)計算約有1 000 個樣點的圖片時，Voronoi 耗時約17. 099 s，AHM AA 耗 時 約0. 100 s；當(dāng) 計 算 約 有10 000 個樣點的圖片時，Voronoi 耗時約1 230. 163 s，AHMAA 耗時約0. 649 s。

圖10 樣點間距的一致性Fig. 10 Consistency of distance between spots

圖9 樣點定位耗時對比Fig. 9 Comparison of spots location time

為了探究算法定位的準(zhǔn)確性，本文計算了樣點之間的距離。 理論上，對于精加工的微孔芯片，每個相鄰孔之間的距離都是固定的，而實際上，通過光學(xué)系統(tǒng)拍攝得到圖像上樣點的間距也是不會有太大誤差的。所以，如果通過本文的算法計算得出的樣點位置坐標(biāo)足夠精確，由此計算得出的樣點之間的距離應(yīng)該也是一個相對固定的值。

如圖10 所示，在芯片中間區(qū)域隨機(jī)選取了1 000 片方形碎片，每塊碎片都包含約20 個樣點，使用AHMAA 算法取得每個樣點的位置信息，計算得出這20 個點與點之間的歐式距離。在這樣的小碎片中，點與點之間大約有24 種不同的距離。其中，相鄰點之間的距離最短，對角線兩端兩個點之間的距離最長。本文計算統(tǒng)計了所有的距離，然后將這1 000 塊碎片的相同距離的數(shù)據(jù)放在一起，計算每種距離的CV 值。最后得出結(jié)果，CV 值最大的是相鄰點之間的距離，為2. 69%，間隔一個點之間的距離CV 值為1. 43%，隨著距離的增大，CV 值逐漸減小。 圖10（b）是截取的前10 組距離的CV 值。 為了驗證最終采集數(shù)據(jù)的準(zhǔn)確性，本文對5 張數(shù)字PCR芯片樣品，分別使用標(biāo)準(zhǔn)儀器系統(tǒng)軟件和本實驗室自行研制的圖像采集系統(tǒng)及AHMAA 算法處理，得到有效填孔數(shù)與拷貝數(shù)，結(jié)果如表1～表2 所 示 。 這5 張 芯 片 中 ，編 號 為C0HGWS、C0HS51 的只含有FAM 染料，編號為C0HY06的只含有VIC 染料，編號為C0HUN8、C0I6MQ的是同時含有FAM 與VIC 染料的芯片。 與標(biāo)準(zhǔn)儀器的結(jié)果相比，通過AHMAA 算法處理得出的有效填孔數(shù)目結(jié)果的誤差率在0. 35%～3. 03%，平均相差 1. 208%，準(zhǔn)確率達(dá)到98. 792%；計算得到的拷貝數(shù)與之相比，誤差率在0. 57%～8. 31%，平均誤差率為3. 8%，準(zhǔn)確率達(dá)到96. 2%。

表2 拷貝數(shù)Tab. 2 Number of copies

表1 有效填孔數(shù)Tab. 1 Effective number of holes to be filled

5 結(jié) 論

為了解決蜂窩狀排列方式的數(shù)字PCR 芯片圖像對比度低、光照不均勻、樣點定位困難等問題，本文提出了基于形態(tài)學(xué)的圖像處理方法。該方法對一組包含約20 000 個樣點的數(shù)字PCR 芯片的三通道圖像，平均運算時間約為20 s，定位最小樣點間距離誤差CV 值小于2. 69%。對5 組不同拷貝數(shù)的數(shù)字PCR 芯片，應(yīng)用AHMAA 算法處理統(tǒng)計得出的有效樣點數(shù)目結(jié)果與標(biāo)準(zhǔn)儀器［1］相比，誤差率在0. 35%～3. 03% 之間，平均相差1. 208%；計算得到的拷貝數(shù)相比，誤差率在0. 57%～8. 31%，平均誤差率為3. 8%。 實驗結(jié)果表明，本文方法在樣點尋址的時候，使用了AHMAA 算法，可以準(zhǔn)確定位蜂窩狀排布方式的樣點，而其他方法使用基于投影法的算法，只能對矩形排布方式的樣點進(jìn)行尋址定位。與Voronoi 定位方法相比，處理相同數(shù)量樣點的圖片時，時間相差3 個數(shù)量級。 在計算的準(zhǔn)確度方面，與使用了預(yù)設(shè)模板的標(biāo)準(zhǔn)儀器相比，樣點識別數(shù)精度達(dá)到98. 79%，拷貝數(shù)計算準(zhǔn)確度達(dá)到96. 2%。該方法實現(xiàn)了對蜂窩狀排布方式圖像的快速準(zhǔn)確處理，能在盡可能短的時間內(nèi)為用戶提供準(zhǔn)確的樣點熒光強(qiáng)度數(shù)據(jù)，為后續(xù)拷貝數(shù)計算等生物統(tǒng)計學(xué)研究提供了精確的原始數(shù)據(jù)。