一种数字PCR微阵列图像分析方法

摘要

本发明提供一种数字PCR微阵列图像分析方法，包括：输入第一通道荧光图像，并对其进行预处理；对预处理后的荧光图像进行图像分割；识别分割图像的有效单元，并提有效单元的质心坐标；根据坐标定位并提取第二与第三通道图像中有效单元，并根据其有效单元统计第二与第三通道图像中有效单元的荧光信号强度；根据其信号强度，绘制第二与第三通道图像中反应单元荧光强度的散点图与柱状图，并统计出阳性单元比例并结合泊松分布，计算出原始反应液中目标核酸分子浓度。本发明以第一通道图像建立掩膜从而间接分析第二通道与第三通道图像，可有效分析非均匀光照下的多通道数字PCR微阵列图像，提高了图像分析的准确度，实现对目标核酸分子绝对定量分析。

说明书

技术领域

本发明涉及图像处理领域，具体而言，涉及一种数字PCR微阵 列图像分析方法。

背景技术

数字PCR是一种具有高灵敏度、高特异性以及绝对定量的核酸 检测技术，利用“分而治之”的思想，将PCR反应体系细分成几万个 微反应体系，使所有反应体系之间独立的进行PCR循环扩增反应， 待反应结束后，统计出阳性单元比例，并利用泊松分布便可计算出原始反应液中目标核酸分子的浓度。近年来，随着芯片制造技术的不断 发展，通过微加工技术制造出含有数以万计的微反应单元的芯片。 2013年由Life Techmologies公司推出的QuantastudioTM3D芯片式数 字PCR系统中的芯片已包含两万个微单元结构。相比于qPCR可以 不依赖标准曲线作为参照实现绝对定量分析，已经在分子诊断、病毒 检测以及食品安全等领域得到广泛应用。然而，对PCR反应之后阳 性单元的统计准确率并不高，目前主要利用图像处理技术对数字PCR 微阵列图像中阳性单元进行识别，随着图像分析技术的日益更新，对 于荧光图像的分析主要含有：图像预处理、图像分割、图像识别以及 后期数据分析等，其中起到关键作用的是图像分割算法。当前主要利 用网格化分析法、自动小室定位法以及引入深度学习技术进行分割， 这些方法缺点是只能针对正交垂直分布的微阵列图像进行分割，不适 用非垂直正交的图像。而通过Otsu法计算分割阈值或者利用Otsu法 计算双阈值进行分割均需在图像光照均匀的情况下进行。而对于光照 不均匀以及图像模糊的图像分割准确率低。

因此，针对目前数字PCR图像分析技术存在分割不准确、识别 率低以及鲁棒性弱等问题，开发高准确度数字PCR微阵列图像分析 方法具有重要意义。

发明内容

本发明的目的在于提供一种数字PCR微阵列图像分析方法，以 解决上述存在的问题。

为实现上述目的，本发明实施例提供一种数字PCR微阵列图像分 析方法，包括

输入第一通道荧光图像，并对其进行预处理；

对预处理后的第一通道荧光图像进行图像分割，以获得分割后图 像；

识别所述分割图像的有效单元，并提取所述有效单元的质心位置 坐标；

根据所述有效单元质心位置坐标定位并提取第二通道与第三通 道图像中有效单元，并根据所述第二通道与第三通道图像中有效单元 统计所述第二通道与第三通道图像中有效单元的荧光信号强度；

根据所述第二通道与第三通道图像中有效单元荧光信号强度，绘 制第二通道与第三通道图像中反应单元荧光强度分布的散点图与柱 状图，并统计出阳性单元比例；

基于所述阳性单元比例并结合泊松分布，计算出原始反应液中目 标核酸分子浓度。

进一步的，输入第一通道荧光图像，并对其进行预处理具体为：

将所述第一通道荧光图像进行灰度处理；

将灰度处理后的图像进行双线性插值处理，以获得双线性插值结 果；

基于所述双线性插值结果，将所述所述第一通道荧光图像进行图 像增强处理；

将图像增强后的图像依次进行中值滤波与高斯滤波，以完成预处 理。

更进一步的，所述基于所述双线性插值结果，将所述第一通道荧 光图像进行图像增强具体为：

将所述第一通道荧光图像进行自适应直方图均衡化，以获得均衡 化后的图像；

将所述均衡化后的图像进行顶帽变换，以获得顶帽变换结果；

将所述均衡化后的图像进行底帽变换，以获得底帽变换结果；

将顶帽变换结果与双线性插值结果求和后再与底帽变换结果做 差，完成图像增强。

进一步的，所述对预处理后的第一通道荧光图像进行图像分割， 以获得分割后图像具体为：

对预处理后的第一通道荧光图像使用Otsu法进行分割，以获得第 一分割图像；

对所述第一分割图像进行膨胀，以获得最大连通区域图像；

对最大连通区域图像和第一分割图像进行逻辑与运算，以获得感 兴趣区域分割图；

获取所述感兴趣区域分割图中欠分割区域的灰度图，并对其进行 预处理；

根据图像信息熵改进Otsu法，对预处理后的欠分割区域进行再次 分割，以获得第二分割图像；

将所述第一分割图像与第二分割图像进行合并，完成图像分割， 其中所述合并利用逻辑运算完成。

更进一步的，所述获取所述感兴趣区域分割图中欠分割区域的灰 度图，并对其进行预处理具体为：

设置面积阈值上限，以从欠分割区域中筛选出大面积区域；

对所述大面积区域进行膨胀，以获取膨胀后连通区域灰度图；

对所述膨胀后连通区域灰度图行预处理。

进一步的，述根据图像信息熵改进Otsu法，对预处理后的欠分割 区域进行再次分割，以获得第二分割图像具体为：

设预处理后的欠分割区域中最大灰度值和最小灰度值分别为 Zmax、Zmin，计算二者均值作为初始分割阈值T0：

计算全局图像信息熵：

其中，E表示全局图像信息熵，h(i)表示灰度值为i的像素个数， N表示整幅图像像素点总数；

以步长1的速度减少图像灰度级两端灰度值及其像素，设减少之 后的最小和最大灰度值分别为α和β；

计算(α，T0-1)和(T0，β)区间的信息熵，分别为:

计算(α，T0-1)和(T0，β)区间像素数占(α，β)区间像 素总数的比例，分别为:

计算(α，β)区间平均信息熵:

E

设存在任意大于0的数，当最小灰度值为α时，最大灰度值为β 时，满足迭代条件：

ΔE＝|E-E

根据感兴趣灰度区间(α，β)，计算感兴趣灰度区间类间方差 最大时的阈值；

根据所述阈值对所述欠分割区域进行分割。

进一步的，所述识别所述分割图像的有效单元，并提取所述有效 单元的质心位置坐标具体为：

根据面积阈值上下限分别为t

设置圆形度阈值，将区域A2中不满足圆形度阈值的区域合并到 区域A1中；

获取区域A3周围的连通区域，将其合并到区域A1中，形成新的 区域A12；

将区域A3周围的连通区域从区域A2中排除，得到区域A22；

对预处理后的第一通道荧光图像使用canny算子进行边缘检测， 并将区域A12与边缘检测结果合并，得到区域A13；

对区域A12进行膨胀，与所述区域A13进行逻辑与运算形成区域 A14；

对区域A14进行填充以及开运算，得到区域A15；

设置离心率阈值，去除区域A15中不满足离心率阈值范围的连通 区域，得到区域A16；

将区域A16与区域A22合并得到图像A4；

计算区域A22中所有区域的平均灰度值；根据所述平均灰度值设 置灰度阈值上下限分别为k

计算图像A41中所有连通域的孔洞大小；设置孔洞阈值，去除图 像A41中不满足孔洞阈值的连通域，得到图像A42；

计算图像A42中所有连通区域的圆形度大小，按照圆形度阈值去 除不在阈值范围内的连通域，得到图像A43；

去除图像A43中面积过大和过小的连通域，完成有效单元的识 别；

将最终识别的连通域的质心标记在第一通道荧光图像的灰度图 上。

进一步的，所述根据所述有效单元质心位置坐标定位并提取第二 通道与第三通道图像中有效单元，并根据所述第二通道与第三通道图 像中有效单元统计所述第二通道与第三通道图像中有效单元的荧光 信号强度具体为：

统计第一通道图像中有效单元的数量与质心位置，读取第二通道 和第三通道图像中有效单元的荧光信号强度。

进一步的，所述根据所述第二通道与第三通道图像中有效单元荧 光信号强度，绘制第二通道与第三通道图像中有效单元荧光强度分布 的散点图与柱状图，并统计出阳性单元比例：

根据第二通道和第三通道图像所有有效单元的荧光信号强度绘 制直方图；

寻找所述直方图中波峰跨度最大的两个极值点坐标，取其均值作 为阳性判定阈值，若荧光强度大于阈值，则判定为阳性单元，反之为 阴性单元；

根据所述判定阈值以及有效单元荧光信号强度，统计出阳性单元 数量以及有效单元总数，并绘制出第二通道和第三通道图像荧光强度 的散点图、直方图以及双参数散点图。

进一步的，所述基于所述阳性单元比例并结合泊松分布，计算出 原始反应液中目标核酸分子浓度，计算公式为：

其中p表示阳性单元数与有效单元总数的比例，N

本发明提供一种多通道数字PCR微阵列图像分析方法，输入第一 通道荧光图像，并对其进行预处理；对预处理后的荧光图像进行图像 分割，以获得分割图像；识别所述分割图像的有效单元，并提取所述 有效单元的质心位置坐标；根据所述有效单元质心位置坐标定位并提 取第二通道与第三通道图像中有效单元，并根据所述第二通道与第三 通道图像中有效单元统计所述第二通道与第三通道图像中有效单元 的荧光信号强度；根据所述第二通道与第三通道图像中有效单元的荧 光信号强度，绘制第二通道与第三通道图像中反应单元荧光强度分布 的散点图与柱状图，并统计出阳性单元比例；基于所述阳性单元比例 并结合泊松分布，计算出原始反应液中目标核酸分子浓度。本发明以 第一通道图像建立掩膜从而间接分析第二通道与第三通道图像，可有 效分析非均匀光照下的多通道数字PCR微阵列图像，提高了图像分析 的准确度，实现对目标核酸分子绝对定量分析。

附图说明

为了更清楚地说明本发明实施方式的技术方案，下面将对实施方 式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了 本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域 普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些 附图获得其他相关的附图。

图1为本发明实施例提供的一种数字PCR微阵列图像分析方法 方法流程示意图。

图2为本发明的第一通道荧光图像；

图3为本发明的第一通道荧光图像自适应直方图均衡化后的效 果图；

图4为本发明的第一通道荧光图像图像增强后的效果图；

图5为本发明的第一通道荧光图像中值滤波后的效果图；

图6为本发明的第一通道荧光图像高斯滤波后的效果图；

图7为本发明的Otsu法分割图像；

图8为本发明的去除周边杂质后的图；

图9为本发明的Otsu法分割欠分割区域灰度图；

图10为本发明的欠分割区域分割图像；

图11为本发明的第一通道图像最终分割图；

图12为本发明的Otsu法分割后的大面积区域图像；

图13为本发明的大面积区域及其周围连通域图像；

图14为本发明的正常面积区域中排除大面积周围连通域图像；

图15为本发明的小面积区域；

图16为本发明的大面积周围区域通过边缘检测图像；

图17为本发明的离心率筛选的对比图；

图18为本发明的经过孔洞、圆形度、面积筛选的前的图像；

图19为本发明的经过孔洞、圆形度、面积筛选的后对比图；

图20为本发明的最终识别的有效单元标记图；

图21为本发明的第二通道的荧光图像；

图22为本发明的和第三通道荧光图像；

图23为本发明的根据第一通道图像提取出的第二通道图像的有 效单元灰度图；

图24为本发明的根据第一通道图像提取出的第三通道图像的有 效单元灰度图；

图25为本发明的数字PCR图像分析APP界面。

具体实施方式

为使本发明实施方式的目的、技术方案和优点更加清楚，下面将 结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行 清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方 式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通 技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式， 都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实 施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅 表示本发明的选定实施方式。

参考图1，本发明实施例提供了一种数字PCR微阵列图像分析方 法，包括：

S11:输入第一通道荧光图像，并对其进行预处理。

在本实施例中，首先将所述第一通道荧光图像进行灰度处理；然 后将灰度处理后的图像进行双线性插值处理，以获得双线性插值结 果。

基于所述双线性插值结果，将所述所述第一通道荧光图像进行图 像增强处理。图像增强处理过程具体为：

将所述第一通道荧光图像进行自适应直方图均衡化，以获得均衡 化后的图像，参考图3；

将所述均衡化后的图像进行顶帽变换，以获得顶帽变换结果；

将所述均衡化后的图像进行底帽变换，以获得底帽变换结果；

将顶帽变换结果与双线性插值结果求和后再与底帽变换结果做 差，参考图4，完成图像增强。

将图像增强后的图像依次进行中值滤波，与高斯滤波，参考图 5、图6，以完成预处理。

S12:对预处理后的第一通道荧光图像进行图像分割，以获得分 割后图像。

在本实施例中，首先对预处理后的第一通道荧光图像使用Otsu法 进行分割，以获得第一分割图像。具体地，利用Otsu法计算出使图像 类间方差最大的阈值Thre1，利用阈值Thre1对全局图像进行分割，公 式如下；

其中，f

在本实施例中，然后对所述第一分割图像进行膨胀，以获得最大 连通区域图像。具体过程是：

求所有连通域的面积，将面积按递增排序，取中间80％面积的平 均值作为正常单元的平均面积。

求正常单元的面积均值对应的等径R:

公式中S是面积。

以等径R作为膨胀结构元大小，对分割图像进行膨胀。

统计膨胀后所有连通域的面积大小，保留面积最大的连通域，去 除其他连通域。

在本实施例中，对上述最大连通区域图像和第一分割图像进行逻 辑与运算，以获得感兴趣区域分割图，参考图8。

在本实施例中，获取所述感兴趣区域分割图中欠分割区域的灰度 图，参考图9，并对其进行预处理；预处理具体过程如下：

设置面积阈值上限，以筛选出大面积区域；例如可以设置上述正 常面积的1.5倍为面积阈值上限t1、0.5倍为面积阈值下限t2；再根据面 积阈值上限筛选出欠分割的区域。当然需要说明的是面积阈值上下限 也可以是其他参数，这些方案均在本发明的保护范围。

在本实施例中，对所述大面积区域进行膨胀，获取膨胀后连通区 域灰度图；例如膨胀结构体大小可为正常单元等径的0.8倍。

在本实施例中，对所述膨胀后连通区域灰度图行预处理。所述预 处理具体包括：

对所述膨胀后连通区域灰度图进行双线性插值处理，以获得双线 性插值结果。

对膨胀后连通区域灰度图像进行自适应直方图均衡化。

对自适应直方图均衡化后的图像进行顶帽变换，以获得顶帽变换 结果。

对自适应均衡化后的图像进行底帽变换，以获得底帽变换结果。

将顶帽变换结果与双线性插值结果求和，再与底帽变换结果做 差。

对求和做差后的图像依次进行中值滤波与高斯滤波，以完成所述 膨胀后连通区域灰度图行预处理。

在本实施例中，根据图像信息熵改进Otsu法，对预处理后的欠分 割区域进行再次分割，以获得第二分割图像。具体分割过程是：

设预处理后的欠分割区域中最大灰度值和最小灰度值分别为 Zmax、Zmin，计算二者均值作为初始分割阈值T0

计算全局图像信息熵：

其中，E表示全局图像信息熵，h(i)表示灰度值为i的像素个数， N表示整幅图像像素点总数；

以步长1的速度减少图像灰度级两端灰度值及其像素，设减少之 后的最小和最大灰度值分别为α和β。

计算(α，T0-1)和(T0，β)区间的信息熵，分别为:

计算(α，T0-1)和(T0，β)区间像素数占(α，β)区间像 素总数的比例，分别为:

计算(α，β)区间平均信息熵:

E

设存在任意大于0的数，当最小灰度值为α时，最大灰度值为β 时，满足迭代条件：

ΔE＝|E-E

其中ε可取0.3E。

根据感兴趣灰度区间(α，β)，计算感兴趣灰度区间类间方差 最大时的阈值Thte

公式中f

将所述第一分割图像与第二分割图像进行合并，去除因膨胀导致 反应单元之间的不连接的现象，完成图像分割，参考图11，其中所述 合并利用逻辑运算完成。

S13:识别所述分割图像的有效单元，并提取所述有效单元的质 心位置坐标。

在本实施例中，首先根据面积阈值上限t1和面积阈值下限t2，将 感兴趣区域分成小面积区域、正常面积区域以及大面积区域三类；其 中将面积小于t

设置圆形度阈值，将区域A2中不满足圆形度阈值的合并到区域 A1中；具体地，先求出区域A2中所有单元的圆形度，圆形度的计算 公式为：

其中S和C分别是每个单元的面积和周长。当F为1时，区域为圆 形，当F大于1时，区域为其他形状。通过设置圆形度阈值为0.86，将 不在圆形度范围内的区域合并到区域A1中。

获取区域A3周围的连通区域，如图13，所示将其合并到区域A1 中，形成新的区域A12。

将区域A3周围的连通区域从区域A2中排除，得要区域A22，区 域A22，如图14所示。

对预处理后的第一通道荧光图像使用canny算子进行边缘检测， 以获得区域A13。

将区域A12(参考图15)与边缘检测结果合并，得到区域A13。

对区域A12进行膨胀，膨胀的结构体大小可为0.2倍的正常点等级 大小，膨胀后图像与得到区域A13的图像进行逻辑与运算形成区域 A14，结果参考图16。

对区域A14进行填充以及开运算，得到区域A15。

设置离心率阈值为0.6，筛除不满足离心率阈值的连通区域，去除区 域A15中不满足离心率阈值范围的连通区域，得到区域A16，如图17。

将区域A16与区域A22合并得到初步筛选后的区域A4，参考图 18。

计算区域A22中所有单元的平均灰度值，记为g；根据所述平均 灰度值设置灰度阈值上下限分别为k

计算图像A41中所有连通区域的孔洞大小，孔洞的大小是指连通 域中灰度值为0的像素个数，设置孔洞阈值可为12，去除连通区域中 孔洞大小大于阈值的连通区域，得到图像A42。

计算图像A42中所有连通区域的圆形度大小，按照圆形度阈值去 除不满足圆形度阈值的连通域，得到图像A43。

去除图像A43中面积过大和过小的连通域，如图19所示。

将最终识别的连通域的质心标记在第一通道的灰度图上，如图20 所示。

S14:根据所述有效单元质心位置坐标定位并提取第二通道与第 三通道图像中有效单元，并根据所述第二通道与第三通道图像中有效 单元统计所述第二通道与第三通道图像中有效单元的荧光信号强度。

在本实施例中，以第一通道图像作为参考图像，第二通道图像， 如图21，和第三通道图像，如图22，作为浮动图像，通过互信息作为 度量进行图像配准。

根据识别出的第一通道图像有效单元的质心位置和数量，定位第 二通道和第三通道图像中有效单元位置，提取得第二通道和第三通道 图像中有效单元，如图23和图24示；

根据第二通道和第三通道图像中的有效单元，统计出所有有效单 元的荧光信号强度。

S15:根据所述第二通道与第三通道图像中有效单元荧光信号强 度，绘制第二通道与第三通道图像中反应单元荧光强度分布的散点图 与柱状图，并统计出阳性单元比例。

在本实施例中，根据第二通道和第三通道图像所有有效样点的荧 光信号强度绘制直方图，寻找直方图中波峰跨度最大的两个极大值 点，取极大值点均值作为阳性判定阈值，若荧光强度大于阈值时，判 定为阳性单元，反之为阴性单元。

根据判定阈值以及反应单元荧光值，统计出阳性单元数量以及有 效单元总数。绘制出第二通道和第三通道图像荧光值的散点图、柱状 图、双参数散点图以及芯片效果图，通过制作APP界面将上述图像展 示，如图25中绘图区所示.

S16：基于所述阳性单元比例并结合泊松分布，计算出原始反应 液中目标核酸分子浓度。

在本实施例中，原始反应液中目标核酸分子浓度，计算公式为：

其中p表示阳性单元数与有效单元总数的比例，N

根据统计学原理计算出所求浓度的置信区间：

其中λ是泊松分布的均值，

本实施例提供的一种多通道数字PCR微阵列图像分析方法，包 括：输入第一通道荧光图像，并对其进行预处理；对预处理后的荧光 图像进行图像分割，以获得分割图像；识别所述分割图像的有效单元， 并提取所述有效单元的质心位置坐标；根据所述有效单元质心位置坐 标定位并提取第二通道与第三通道图像中有效单元，并根据所述第二 通道与第三通道图像中有效单元统计所述第二通道与第三通道图像 中有效单元的荧光信号强度；根据所述第二通道与第三通道图像中有 效单元的荧光信号强度，绘制第二通道与第三通道图像中反应单元荧 光强度分布的散点图与柱状图，并统计出阳性单元比例；基于所述阳性单元比例并结合泊松分布，计算出原始反应液中目标核酸分子浓 度。第一通道图像建立掩膜从而间接分析第二通道与第三通道图像， 可有效分析非均匀光照下的多通道数字PCR微阵列图像，提高了图像 分析的准确度，实现对目标核酸分子绝对定量分析。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的 技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神 和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发 明的保护范围。