数字液滴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图像处理提供了一种可行的方法。