數字液滴PCR圖像的熒光信息提取

李姍姍，王子超，于成壯，魏春陽，李軍委

（1.河北工業大學 機械工程學院，天津300401；2.河北省智能傳感與人機融合重點實驗室，天津300130；3.電工裝備可靠性與智能化國家重點實驗室，天津300132；4.河北工業大學 生物物理研究所，天津300401）

1 引 言

數字聚合酶鏈式反應（Polymerase Chain Reaction，PCR）技術是應用于核酸檢測的主要手段之一。隨著微流控技術的發展，液滴技術以其體積小、精度高、反應腔之間完全隔離等優點，在核酸檢測中得到了廣泛的應用［1］。

dd PCR實驗包含4個環節，即待檢測樣本與液滴的制備、目標核酸的擴增、熒光信號的采集、實驗數據分析與處理。其中，熒光信號采集的精度會直接影響dd PCR實驗的結果［2］。目前，熒光信號采集方法有顯微圖像檢測法和基于流式細胞術［3］的熒光檢測法。相較于后者，顯微圖像檢測法的液滴破裂更少、精度更高、檢測速度更快。

ddPCR實驗中，液滴熒光顯微圖像由大量尺寸相似熒光強度不同的圓形亮點組成，影響實驗結果的關鍵在于對圖像中液滴的判別、分類與計數，即計算液滴的總數并讀取各個液滴的熒光信號，判斷信號是陰性還是陽性。圖像斑點檢測主要包括圖像去噪聲、增強目標斑點特征和斑點計數［4-5］。常規的圖像斑點計數通常使用圖像灰度閾值分割法將圖像二值化，再計算目標斑點對應區域的數量［6］，然而這種單一閾值的圖像二值化方法只適用于圖像背景均勻且目標點單一的情況；在多目標檢測中，這種方法的誤差會隨著檢測目標種類的增多而增大。因此，熵最小值法［7］和最大類間方差法等［8］多種閾值分割法被提出。這些方法改進了閾值分割法，通過數學運算自動確定用于圖像二值化的灰度閾值。這種方法僅依靠全局閾值劃分，不能區分目標的形狀和幾何特征，雖然適用于一部分實驗，但并不能滿足所有的dd PCR熒光檢測需求。

在ddPCR實驗中，液滴的數量多且排列緊密，圖像噪聲多［9-12］，需要采取對液滴特性具有針對性的檢測方法。Zhang等提出了多尺度方差穩定變換（MSVST）檢測器方法［13］。MSVST方法結合了方差穩定變換和各向同性非抽樣小波變換，在混合泊松-高斯噪聲濾除和熒光斑點檢測方面表現良好。然而，dd PCR實驗圖像中出現的亮斑噪聲、小氣泡以及圖像的光照不均勻效應都會導致液滴的錯誤判別。Basset等提出了針對圖像中不同光斑尺寸選擇最佳LoG尺度或多尺度的方法［14］，但各個實驗液滴圖像最優的尺度信息通常難以獲得，所以該方法并不適用于dd PCR的熒光檢測。

為解決ddPCR實驗中液滴難以計數的問題，本文提出了一種圖像處理方法，針對dd PCR圖像特性去除噪聲與增強對比度，并基于灰度遍歷法采集圖像，對液滴進行識別、分類和計數。人類gDNA的dd PCR核酸檢測實驗結果表明，本方法可準確獲取ddPCR實驗的圖像信息。

2 待檢測液滴的制備與圖像拍攝

液滴的制備采用自研制的多噴嘴階梯式微流控芯片，圖1為微流控芯片的結構與液滴生成示意圖。

圖1 數字液滴PCR芯片與液滴生成Fig.1 ddPCR chip and droplets formation

聚二甲基硅氧烷（Polydimethylsiloxane，PDMS）流道層的模板使用玻璃片通過SU-8光刻得到，陽模用二甲基二氯硅烷試劑硅烷化進行表面處理，將PDMS與固化劑按質量比10∶1混合，攪拌均勻澆注在微通道陽模上，80℃加熱1 h左右使PDMS固化，揭下帶有微通道的PDMS并在兩相流體的入口處打孔，完成流道層的制作。再使用等離子機將PDMS流道層與玻璃基底進行鍵合，便得到了芯片。

流道層頂部有兩個入口，通過氣泵和導管分別通入生物油和待檢測樣本溶液，生物油為連續相，待檢測樣本為分散相。在流體剪切力的作用下在流道右側生成大量均一液滴，直徑為90～110μm。芯片油相流道高度為200μm，樣本溶液流道高度為30μm。實驗檢測樣本為人類gDNA組樣本（杭州沃森生物）、熒光標記物羧基熒光素（杭州沃森生物），將其制成待檢測樣本溶液備用。在液滴制備完成后通過數字PCR儀（蘇州思納福Sniperdx DQ24）對液滴進行40次熱循環，使核酸分子大量擴增并與熒光探針相結合，并對液滴進行拍攝，獲取液滴的8位灰度圖像，每張圖片的尺寸為3 344×3 092 pixel。

3 圖像預處理

3.1 圖像去噪聲

在熒光顯微圖像拍攝過程中，顯微拍攝模塊與熒光激發光源會產生高斯噪聲，導致生成的液滴顯微圖像會出現很多獨立的像素塊和亮點。為了減小圖像噪聲的影響，這里使用高斯濾波法［15］去除圖像噪聲。

實驗中，熒光顯微圖像的主要噪聲來源是熒光激發光源光照的不均勻性［16-17］。在圖像采集過程中，光照環境或物體表面反光等原因造成的圖像整體光照不均勻，導致圖像局部灰度值差異較大，局部信息無法辨認。為了去除圖像光照不均勻效應，本文通過基于正態分布平均化的光照分量分離法去除圖像的光照不均勻效應。

二維圖像正態分布平均化方程為：

式中：σ是正態分布的標準偏差；x，y是圖像像素點坐標。首先，計算各個像素坐標點的正態分布權重，使用正態分布平均化對圖像像素進行變換，再用圖像中像素點的灰度值乘以該坐標點的正態分布權重，每個像素的值為相鄰像素值的加權平均。通過此方法濾去圖片中的高頻成分，即細節剩下的低頻成分為光照分量，如此得到光照分量的圖像，再用原圖與光照分量圖像作差便可有效去除圖像的光照不均勻效應。

3.2 圖像灰度對比度的增強

在經過高斯濾波與去除光照不均勻后，圖像對比度降低，液滴圖像的灰度值范圍小，不便于觀察，難以識別和分類。本文使用基于ddPCR圖像特點的圖像對比度增強方法，通過累積分布函數［18］對圖像進行灰度分布均衡化處理。式中：n為圖像中像素的總和，nk為灰度級為r k的像素個數，L為圖像當中的灰度級數。由于ddPCR圖像是暗場下拍攝得到的，所以圖像中灰度最大值為熒光液滴部分，灰度值0為暗場背景部分，因此將原圖像的灰度級從r k（k=0，1，2，…，L-1）映射到完整灰度級sk（k=0，1，2，…，255），能夠以最大限度提高圖像的整體灰度對比度，并且不會造成圖像信息的改變與丟失。

4 液滴計數

4.1 圖像二值化

圖像預處理完成后，需要通過灰度閾值分割法分離熒光液滴、空液滴和背景區域，并計算熒光液滴與全部液滴的數量。獲得準確的灰度閾值是液滴分類與計數的必要條件；本文通過灰度遍歷法以0～255的灰度點逐一作為圖像的二值化閾值，對圖像進行二值化處理，得到一組256幅的二值化圖像；由于實驗中的液滴數量多且排列緊密，在圖像二值化后，圖像中大部分液滴會發生粘連的情況，使用分水嶺算法［19］將得到的256幅圖像中粘連的液滴分割開，后續再基于每一幅圖像所提取的信息確定灰度閾值，并對液滴進行分類與計數。

4.2 液滴的識別條件

為了計算圖像中有效液滴的數量，需根據液滴的特性設置識別條件。因為液滴的幾何均勻性較好，尺寸和圓度差異較小，所以本文將圖像中達到尺寸范圍或圓度范圍的區域記為液滴，以提高液滴識別的準確度。

根據芯片中生成液滴的直徑范圍和圖像尺寸，通過顯微放大倍數換算得出單個液滴的尺寸為375～418 pixel。

液滴圓度的計算公式為：R=S1/S2，R為圓度，S1為圖形面積，S2為最大外接圓面積。實驗中的液滴圓度均一性較高且為標準圓形，但受到圖像分辨率的影響，液滴圖像的圓度略小于1。選取圖組的局部實驗圖像，以液滴尺寸375～418 pixel為約束，統計共1 650個液滴的圓度，統計圖如圖2所示。計算得到液滴圓度的期望值E（R）=0.95，設置液滴圓度為0.95～1。

圖2 液滴的圓度統計Fig.2 Roundness statistics of droplets

4.3 液滴的分類與計數

在dd PCR實驗中，由于各液滴核酸分子的含量不同，液滴所發出的熒光強度也不同，導致同類液滴的圖像灰度也有所差異，用于區別熒光液滴與空液滴的灰度閾值無法直接確定。因此，本文基于灰度遍歷法所采集的圖像與數據，通過微分分析法對液滴進行分類和計數。通過檢測累加器分別計算256幅二值化圖像中滿足液滴識別條件的封閉區域的數量，檢測結果即為識別到的液滴數量；定義液滴數-灰度閾值函數表達式為g(k)，k=0，1，2，…，255；通 過 拉 格 朗 日 插 值 法對灰度閾值0～255的計數結果散點圖進行擬合，得到擬合曲線并對擬合曲線做濾波處理，令得到的液滴計數結果擬合曲線的函數表達式為g1(k)，k∈(0，255]。由 于 圖 像 中 僅 存 在 暗 場 背景、空液滴和熒光液滴3部分，且同種類液滴的圖像灰度值較為接近，因此，在擬合曲線g1(k)中存在兩部分明顯平緩的區間，顯然熒光液滴與全部液滴的數量信息分別位于兩個區間的數據之中。為了確定這兩個數據，對擬合曲線g1(k)進行一階求導，設導數為z(k)，k∈(0，255]。在|z(k)|取極小值時，可認為該閾值下液滴數量最穩定，即最接近實際數量。

將點k=1，2，…，255代入導函數z(k)中，解出|z(k1)|≤|z(k2)|≤|z(k3)|≤…≤|z(k255)|。解出點k1，k2，即得出熒光液滴數量為min{g(k1)，g(k2)}，圖像中全部液滴數量為max{g(k1)，g(k2)}。

5 實驗與結果分析

5.1 實驗圖像

圖3為以人類gDNA組為樣本的dd PCR待處理實驗圖像，圖像中均勻分布著熒光液滴和空液滴，共包含約20 000個液滴。

圖3 液滴圖像Fig.3 Image of droplets

5.2 圖像預處理的作用與效果

圖像去除光照不均效果如圖4所示，如圖4（b）右側的局部放大圖所示，圖像右下方光照明顯較強，左上方光照較暗，導致不同區域的二值化效果差別較大。如圖4（d）局部放大圖所示，在消除了光照不均勻的影響后，圖像整體的二值化效果得到明顯改善。

圖4 消除光照不均勻前后的對比圖Fig.4 Contrast images before and after eliminating uneven illumination

經過高斯濾波與去除光照不均勻后，圖像的灰度值集中在0～120，不足灰度圖灰度范圍的1/2。圖5（a）展示了灰度對比度增強前后的圖像效果對比，圖5（b）為圖像對比度增強前后的圖像灰度分布率對比。在圖像對比度增強前圖像的灰度值集中在0～120，圖片整體亮度非常低，幾乎分辨不出液滴之間的界限。處理后圖片的灰度值分布在0～255，圖像灰度分布率的標準差由19.28增大到54.29，灰度變化趨勢明顯平緩，為后續液滴計數打下基礎。

圖5 圖像對比度增強效果Fig.5 Enhancement effect of image contrast

5.3 液滴計數結果

以灰度0～255逐一作為圖像的二值化閾值，對圖像組做二值化處理。通過液滴識別條件統計每幅圖所識別到的液滴數量。設液滴數-灰度閾值函數表達式為g(k)，k=0，1，2，…，255。液滴計數結果擬合曲線的函數表達式為g1(k)，k∈(0，255]。液滴的計數結果、計數結果擬合曲線及其求導函數z(k)如圖6所示。

圖6 液滴計數結果Fig.6 Counting results of droplets

解出點k1=38，k2=137，得出圖像中熒光液滴數量為g(k1)=2 782，全部液滴數量為g(k2)=21 923。

5.4 耗時分析

一組完整的ddPCR檢測實驗包括樣本配制、液滴生成、熱循環、圖像拍攝與圖像處理，完整實驗平均耗時2～2.5 h。本文的圖像處理算法平均耗時2～3 min，約占實驗總時長的2%，可滿足dd PCR檢測的時間要求。

5.5 算法性能對比

5.5.1 與商用儀器軟件對比

為了測試本算法的精確度，以20份人類gDNA的ddPCR實驗圖像為測試樣本，每張圖像中均包含約15 000～20 000個液滴，使用本文算法與商用儀器軟件Sight（思納福公司）對圖像中的液滴進行計數。對熒光液滴和全部液滴的識別效果如圖7所示，Sight使用十字與圓圈分別標記全部液滴和熒光液滴，本文算法用圖像掩膜分別標記全部液滴和熒光液滴。Sight漏計的液滴被本文算法準確識別。對于因算法誤差導致液滴被遺漏的情況，計算每個液滴的平均熒光強度，確定液滴的類別并統計兩種方法漏計的液滴數量，計算得出本文算法與商用儀器軟件的計數平均誤差率分別為0.64%，2.88%。對比實驗中的3組數據如表1所示。

圖7 本文算法與Sight的液滴識別效果Fig.7 Droplet recognition effect of proposed algorithm and Sight

表1 本文算法與Sight漏計液滴數Tab.1 Quantity of missing droplets of proposed algorithm and Sight

5.5.2 與同類算法對比

使用機器學習算法（Trainable Weka Segmentation，TWS）［20］與文獻［21］算法對上述20份實驗圖像的液滴進行分類和計數。統計算法漏計的液滴數量，其中3組數據如表2所示。計算得出這兩種算法對液滴計數的平均誤差率分別為3.04%，3.29%。

表2 文獻［21］算法與TWS漏計液滴數Tab.2 Quantity of missing droplets of reference［21］algorithm and TWS

本文算法的平均誤差率為0.64%，平均準確率為99.36%，與商用儀器軟件Sight相比，準確率提高了2.24%；與TWS算法和文獻［21］算法相比，檢測準確率平均提高了2.53%。本文提出的圖像處理方法對ddPCR液滴的檢測精度更高。

6 結 論

本文提出了一種ddPCR圖像熒光信息提取方法。該方法通過基于卷積運算的去噪法與灰度分布均衡化去除了圖像噪聲并增強了圖像對比度；基于灰度遍歷法采集圖像信息，通過微分分析法對液滴進行分類和計數，解決了由于液滴數量多、尺寸小、排列緊密、熒光強度不均勻所造成的液滴難以分類和計數的問題。在人類gDNA檢測實驗中，該方法的平均檢測準確率為99.36%，與商用儀器算法和同類算法相比平均檢測準確率分別提高了2.24%，2.53%。該方法為ddPCR檢測實驗提供了精確的檢測結果，為ddPCR圖像處理提供了一種可行的方法。