蜂窩狀數字PCR 微陣列熒光圖像的信息提取

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

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

1 引 言

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

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

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

2 實驗材料

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

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

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

圖1 數字PCR 芯片結構Fig. 1 Structure of digital PCR chip

3 圖像處理

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

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

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

3.1 三通道圖像配準

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

3.2 樣點尋址

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

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

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

（3）增強樣點孔內外的對比度、二值化；

（4）計算統計樣點的位置信息。

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

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

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

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

算法主要步驟如下：

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

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

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

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

（5）將圖像從HSV 空間轉化到RGB 空間，得到I ′RGB。

二維γ函數公式為：

為了增強樣點與周圍環境，即芯片上孔內外的對比度，本文提出一種基于頂帽變換與底帽變換增強對比度的算法。頂帽變換與底帽變換都是基于數學形態學的圖像處理方法［28］。

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

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

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

3.3 樣點分割

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

3.4 信息提取

計算每個六邊形區域內所有像素點的均值，作為樣點的熒光強度值。統計得到所有樣點的熒光信號值后，本文采用基于K 均值的分類方法對數據分類，得到芯片的陰陽性點數目信息。

4 結果與討論

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

圖3 樣點定位結果Fig. 3 Result of spot addressing

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

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

4.1 三通道圖像配準

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

4.2 對比度增強的作用

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

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

圖6 展示了整體對比度增強前后，圖片的灰度直方圖分布。原圖的灰度值集中在0～100，圖片整體亮度非常暗，而且幾乎分辨不出芯片區域的邊界；處理后圖片的灰度值均勻地分布在0～255，邊界非常明顯，為后續ROI 的提取打下了很好的基礎。

4.3 光照不均勻對結果的影響

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

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

4.4 去除無效樣點的必要性

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

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

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

4.5 定位時間及精度

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

為了驗證算法的定位效率，本文在一塊完整的數字PCR 芯片上隨機截取4 組局部圖片，每組圖片分別有100 張，4 組圖片分別含有約10 個，100 個，1 000 個，10 000 個樣點。對每張圖片，分別用Voronoi 定位方法與本文提出的AHMAA定位算法進行計算，得出每張圖片運算的時間，并對每組數據計算其均值與標準差，畫出圖9 的誤差條形圖。當計算約有10 個樣點的圖片時，Voronoi 耗 時 約0. 297 s，AHMAA 耗 時 約0. 022 s；當計算約有100 個樣點的圖片時，Voronoi 耗時約0. 968 s，AHMAA 耗時約0. 028 s；當計算約有1 000 個樣點的圖片時，Voronoi 耗時約17. 099 s，AHM AA 耗 時 約0. 100 s；當 計 算 約 有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

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

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

表2 拷貝數Tab. 2 Number of copies

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

5 結 論

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