dPCR微阵列图像信息处理方法与流程
本发明属于生物芯片技术领域，具体涉及一种dpcr微阵列图像信息处理方法。
背景技术：
20世纪末，vogelstein等首次提出了“数字pcr(digitalpcr，dpcr)”的概念，该技术实现了核酸拷贝数绝对定量的突破。数字pcr的主要原理是将单个dna分子置于独立的反应室中，并对其进行pcr扩增，利用
微阵列技术已经广泛应用于生物信息的同步检测，完整的微阵列生物芯片分析过程包括样本采集、芯片制备、扫描成像、图像处理和数据分析等几个部分。图像处理是定位图像中的微阵列点(样点)，将每个样点对应的形状和强度量化，主要包含图像预处理、网格定位、样点分割以及分割效果的评价、信号提取等；其中，样点分割是图像处理的难点，样点分割的好坏直接影响最终信号提取的结果。现有技术一“imagegriddingalgorithmfordnamicroarrayanalyser(maziidahmukhtar,etal.2016,3rdinternationalconferenceonelectronicdesign(iced)：452-457.)”公开了一种dna微阵列图像自动网格化分析方法，包括网格化、分割和强度提取三个步骤，但该网格化处理只能针对正交垂直分布的芯片图像，使用otsu算法进行样点分割也需图像的亮度分布均匀。现有技术二“cdnamicroarrayimagesegmentationwithanimprovedmovingk-meansclusteringmethod(shaog,etal.2015,ieeeinternationalconferenceonsemanticcomputing.)”公开了一种基于抽象聚类的图像分割策略，首先提出了一种自动对比度增强方法来提高图像质量，其次进行最大类间方差网格划分，将斑点分割成单独的区域，然后将k均值聚类算法与移动k均值聚类方法相结合，以获得更高的分割精度，但是该策略仍需要经过网格化这一步骤，不利于对其它排布样式的芯片图像，而且，对每个样本采用移动k均值聚类算法进行分割，计算成本非常高，不利于高通量的生物芯片的使用。
技术实现要素：
本发明要解决的技术问题在于克服现有的生物芯片图像处理方法无法对非垂直正交排列的微阵列芯片图像做准确、自动的样点信息提取和分析的缺陷，从而提供一种全新的dpcr微阵列图像信息处理方法，以寻址定位方法代替传统的网格化步骤，以达到对高通量dpcr微阵列图像的准确快速分割、提取和分析的目的。
为解决上述技术问题，本发明采用的技术方案是：
本发明提供的一种dpcr微阵列图像信息处理方法，包括以下步骤：
(1)输入至少三个荧光通道的dpcr微阵列图像，包括通道1图像、通道2图像、通道3图像，对所述通道1图像、所述通道2图像、所述通道3图像中任意两个进行图像配准处理；
(2)对步骤(1)中的三个所述dpcr微阵列图像分别进行对比度增强处理；
(3)将步骤(2)中的经对比度增强处理的所述通道1图像、所述通道2图像、所述通道3图像融合为rgb图像，去除所述rgb图像的光照不均匀影响，并进行二值化修正，得到面积大小均一的连通域；
(4)提取步骤(3)中所述连通域的中心点坐标，并以所述中心点为形心选取每个样点的roi区域；
(5)提取步骤(4)中每个所述roi区域的信号值作为所述样点的信号值，输出所述dpcr微阵列图像上每个所述样点的信号结果。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(1)中，所述图像配准处理的方法为：以其中一个所述dpcr微阵列图像作为参考图像，以另外两个所述dpcr微阵列图像利用互信息进行图像配准处理。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(1)中，三个所述dpcr微阵列图像所对应的荧光染料分别选自a类荧光染料、b类荧光染料、c类荧光染料、d类荧光染料中的任意一种，且三个所述dpcr微阵列图像所对应的荧光染料类不同；
其中，所述a类荧光染料选自fam、sybrgreeni中的一种；
所述b类荧光染料选自vic、hex、joe、tamra、tet、cy3中的一种；
所述c类荧光染料选自texas-red、rox中的一种；
所述d类荧光染料选自cy5。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(1)中，以所述fam荧光染料、所述vic荧光染料、所述rox荧光染料分别对应的所述dpcr微阵列图像依次作为所述通道1图像、所述通道2图像、所述通道3图像，对所述通道2图像、所述通道3图像进行图像配准处理。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(2)中，所述对比度增强处理的方法为：对所述通道1图像、所述通道2图像、所述通道3图像分别做顶帽变换和底帽变换，将每组所述顶帽变换的值以及对应的所述底帽变换的值分别乘以系数并求二者差值，从而分别得到所述通道1图像、所述通道2图像、所述通道3图像的对比度增强结果。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(3)中，将所述通道1图像、所述通道2图像、所述通道3图像分别放入r通道、g通道、b通道以融合得到所述rgb图像，并采用基于二维伽马函数的光照不均匀图像自适应校正算法对所述rgb图像去除光照不均匀的影响。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(3)中，所述二值化修正包括对二值图像填充、开运算以及去除面积过小和/或过大的连通域。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(4)中，采用regionprops函数获得所有所述连通域的中心点坐标。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(4)中，以所述连通域的中心点坐标为圆心，选取所述样点大小的多边形区域作为所述样点的roi区域。
进一步优选地，该dpcr微阵列图像信息处理方法，步骤(1)和步骤(2)之间还包括对所述通道1图像以及经图像配准处理的所述通道2图像、所述通道3图像分别进行中值滤波的处理；
在步骤(2)和步骤(3)之间还包括对经对比度增强处理过的所述通道1图像、所述通道2图像、所述通道3图像分别进行均值滤波的处理。
本发明技术方案，具有如下优点：
1.本发明提供的dpcr微阵列图像信息处理方法，以一种寻址定位方法，代替传统的网格化步骤，能够对非垂直正交排列的微阵列芯片图像做准确、自动的样点信息提取和分析。
2.本发明提供的dpcr微阵列图像信息处理方法，改变了现有处理微阵列图像的步骤都需要对整幅图像做网格化操作的方式，自动化程度高，省时省力。
3.本发明提供的dpcr微阵列图像信息处理方法，相对于现有技术不能对六边形生物芯片做网格化处理的弊端，本发明提供的方法能够对六边形芯片做网格化处理。
4.本发明提供的dpcr微阵列图像信息处理方法，能够解决现有数据库的单张微阵列芯片图像样点数少、通量低的问题。
5.本发明提供的dpcr微阵列图像信息处理方法，改变了现有不能一次获取完整的芯片图从而导致后期处理需要进行图像拼接的难题。
6.本发明提供的dpcr微阵列图像信息处理方法，提高了图像分割的精确度，减少了对后续的生物学分析的准确性的负面影响。
7.本发明提供的dpcr微阵列图像信息处理方法，也可以应用到常规排列的生物芯片图像的处理，基于实验室的96通道、2万孔/单元的微阵列芯片，每幅图像包含超大量的生物信息数据，使用本发明提供的方法可以同时处理整幅图像的样点，实现对高通量芯片的处理。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例1提供的dpcr微阵列图像信息处理方法流程图；
图2为本发明实施例1提供的通道1图像、通道2图像、通道3图像；
图3为本发明实施例1提供的基于互信息配准方法流程图；
图4为本发明实施例1提供的通道1图像中值滤波前与中值滤波后的对照图；
图5为本发明实施例1提供的通道2图像经中值滤波处理后效果图；
图6为本发明实施例1提供的通道3图像经中值滤波处理后效果图；
图7为本发明实施例1提供的通道1图像对比度增强效果图；
图8为本发明实施例1提供的通道1图像经均值滤波处理后效果图；
图9为本发明实施例1提供的通道2图像经均值滤波处理后效果图；
图10为本发明实施例1提供的通道3图像经均值滤波处理后效果图；
图11为本发明实施例1提供的通道1图像、通道2图像、通道3图像融合为rgb图像自动配准结果图；
图12为本发明实施例1提供的基于二维伽马函数的光照不均匀图像自适应校正算法流程图；
图13为本发明实施例1提供的去除光照不均匀影响前后对比图；
图14为本发明实施例1提供的二值化图像修正前后对比图；
图15为本发明实施例1提供的中心点坐标定位结果图；
图16为本发明实施例1提供的圈取roi示意图；
图17为本发明实施例1提供的样点信息统计分类结果图；
图18为本发明实施例1提供的通道1图像样点荧光信号强度分布图；
图19为本发明实施例1提供的通道2图像样点荧光信号强度分布图。
具体实施方式
为了便于理解本发明的目的、技术方案和要点，下面将对本发明的实施方式作进一步详细描述。本发明可以多种不同的形式实施，而不应该被理解为仅限于在此阐述的实施例。相反，提供此实施例，使得本发明将是彻底的和完整的，并且将把本发明的构思充分传达给本领域技术人员，本发明将仅由权利要求来限定。
实施例1
本实施例提供一种的dpcr微阵列图像信息处理方法，如图1所示，包括以下步骤：
步骤一，输入三个荧光通道(通道1、2和3)的dpcr微阵列图像：
该三个dpcr微阵列图像所对应的荧光染料分别选自a类荧光染料、b类荧光染料、c类荧光染料、d类荧光染料中的任意一种，且三个dpcr微阵列图像所对应的荧光染料类不同：a类荧光染料选自fam、sybrgreeni中的一种，b类荧光染料选自vic、hex、joe、tamra、tet、cy3中的一种，c类荧光染料选自texas-red、rox中的一种，d类荧光染料选自cy5。
具体地，如图2所示，输入三个荧光通道的dpcr微阵列图像：以fam荧光染料、vic荧光染料、rox荧光染料分别对应的dpcr微阵列图像依次作为通道1图像(图2中(a))、通道2图像(图2中(b))、通道3图像(图2中(c))，且上述三张图像的规格大小相同。在荧光通道直接叠加时，为避免不同的通道存在错位现象，需将dpcr微阵列的相对位置设置细微差别。
步骤二，图像配准
为了使得不同条件环境下所获得的图像上特征点在同一视觉系统中更容易达到一对一映射，需对不同的图像进行校正相对位置平移、角度旋转和缩放尺度，以达到两张或多张图像在空间上对齐的目的。
基于互信息的医学图像配准方法，在配准之前设置配准参数：生长因子设为1.01；ε设为1.5×10-6；初始半径设为0.001；最大迭代次数设为300。使用matlab自带的imregister函数对步骤一中的通道2图像、图像3图像进行图像配准处理，包括4个基本模块：几何变换、插值算法、相似性测度和优化方法。具体地，如图3所示，以通道1图像作为参考图像r，以通道2图像、通道3图像作为浮动图像f，浮动图像f先依次经过几何变换的初始变换t、差值算法，得到变换后的浮动图像t(f)，然后对变换后的浮动图像t(f)和参考图像r进行优化算法处理以判断互信息是否最大，若不是，则将几何变换进行更新变换t，重复插值算法、优化算法的过程，直至互信息最大，输出最优的图像配准参数。基于互信息对通道2图像、通道3图像进行图像配准处理无需预处理、配准精度高、鲁棒性好。
步骤三，中值滤波
在图像的生成、扫描、传输等过程中，由于许多因素的影响，会产生各种类型的噪声，造成图像质量下降，因此减弱图像噪声、增强图像质量就显得很必要。dpcr微阵列图像中通常存在的噪声类型有：均匀分布噪声、高斯噪声、伽马噪声、椒盐噪声等。
具体地，分别对步骤一中的通道1图像以及步骤二中的经图像配准处理的通道2图像、通道3图像分别进行中值滤波处理。通道1图像经中值滤波处理后如图4所示，其中，图4(a)为中值滤波前，图4(b)为中值滤波后，对比图4(a)和图4(b)发现，椒盐噪声很明显地被消除；通道2图像经中值滤波后如图5所示，通道3图像经中值滤波后如图6所示，椒盐噪声已明显被消除。
步骤四，对比度增强
对比度增强是将图像中的亮度值范围拉伸或压缩成显示系统指定的亮度显示范围，从而提高图像全部或局部的对比度，增强原图中各部分的差异。
具体地，对步骤三中经过中值滤波处理的通道1图像、通道2图像、通道3图像分别做顶帽变换和底帽变换，得到通道1图像、通道2图像、通道3图像三组对应的顶帽变换的值和底帽变换的值，求每组顶帽变换的值和底帽变换的值乘以系数后的差值，从而分别得到通道1图像、通道2图像、通道3图像的对比度增强结果。通道1图像的对比度增强效果如图7所示，图7(a)为顶帽变换后，图7(b)为底帽变换后，图7(c)为差值效果，可见，通道1图像、通道2图像、通道3图像各自的对比度得到提高，原图中各部分的差异得到增强。
步骤五，均值滤波
在图像的生成、扫描、传输等过程中，由于许多因素的影响，会产生各种类型的噪声，造成图像质量下降，因此减弱图像噪声、增强图像质量就显得很必要。dpcr微阵列图像中通常存在的噪声类型有：均匀分布噪声、高斯噪声、伽马噪声、椒盐噪声等。
具体地，分别对步骤四中经对比度增强的通道1图像、通道2图像、通道3图像分别进行均值滤波处理，通道1图像经均值滤波处理后如图8所示，通道2图像经均值滤波处理后如图9所示，通道3图像经均值滤波处理后如图10所示，经均值滤波处理后，通道1图像、通道2图像、通道3图像整体更加平滑。
步骤六，图像融合
为提高图像信息的利用率、提升原始图像的空间分辨率和光谱分辨率，需对均值滤波处理的通道1图像、通道2图像、通道3图像融合为一张新的图像。
具体地，将通道1图像、通道2图像、通道3图像分别放入r通道、g通道、b通道经过粗配准、精配后以融合得到rgb图像。其中，粗配准采用默认参数：生长因子设为1.05；ε设为1.5×10-6；初始半径设为0.0063；最大迭代次数设为100；精配准是对上述初始参数进行了调整：生长因子设为1.01；ε设为1.5×10-6；初始半径设为0.001；最大迭代次数设为300。通道1图像、通道2图像、通道3图像自动配准的过程如图11所示，图11(a)为配准前状态，图11(b)为粗配准后状态，图11(c)为精配准后状态。经粗配准和精配准后使得通道1图像、通道2图像、通道3图像中对应于空间同一位置的点一一对应起来，从而达到信息融合的目的。
步骤七，去除光照不均匀影响
基于二维伽马函数的光照不均匀图像自适应校正算法对步骤六得到的rgb图像进行去除光照不均匀的影响。
如图12所示，具体流程如下：
1.输入rgb图像坐标f(x,y)，由rgb空间转到hsv空间；
2.保持色调h、饱和度s不变，对亮度v进行处理；
2.1对亮度v做多尺度高斯函数卷积，设置3个不同尺度，每个权重也不相同：
将f(x,y)代入公式(1)：
得到亮度的高斯函数值，将g(x,y)代入公式(2)：
得到亮度的光照分量值，其中，ω是权重，n＝3，λ是归一化常数，c是尺度因子；
2.2根据公式gamma校正处理
将i(x,y)代入公式(3)：
γ＝(1/2)(m-i(x,y))/m(3)
γ代入公式(4)：
其中，m是i(x,y)的均值(代码：m＝mean(i(：))；)，o(x,y)是新的亮度v’；
3.用新的亮度v’和原来的h、s合成新的图像，再转到rgb空间，得到最终的去除光照不均影响的结果f’(x,y)，其结果如图13所示，(a)为rgb原图，(b)为otsu阈值分割结果图，(c)为光照分量图，(d)为去除光照不均影响后的分割结果图，通过(a)、(b)、(c)、(d)，尤其(b)和(d)各选取两处对比，rgb图像经过处理后，光照不均匀问题得到了消除，图像质量得到提高。
步骤八，二值化修正
将步骤七得到的图像进行二值化修正处理，得到二值图，对二值图填充，开运算，去除面积过小过大的连通域，对于二值图上过大和过小的光斑，我们认为是噪声和其它非样点的信号，留下面积范围较为均匀的连通域作为样点的荧光信号光斑。将处理前的图像(图14(a)所示)与处理后的图像(图14(b)所示)进行对比，经二值化修正处理后，样点的荧光信号光斑分布更加均匀。
步骤九，提取中心点坐标
采用regionprops函数获得步骤八中面积大小较为均匀的连通域的中心点坐标，如图15所示，并以中心点为形心选取每个样点大小的多边形区域为roi区域，如方形、六边形、八边形等等。
具体地，以中心点为形心选取每个样点大小的六边形区域作为roi区域，如图16所示。
步骤十，提取样点的荧光信号值
提取样点的荧光信号值，得到dpcr微阵列图像所有样点的荧光信号值，由于通道1、通道2、通道3三个荧光通道的信息是相互独立的，因此在rgb图像上的任意某个区域都可以选取三个荧光通道各自的信息，本实施例中选取通道1和通道2荧光通道的信息，也即rgb图像的r通道和g通道的信息。
对于dpcr微阵列上通道1和通道2荧光通道的信息进行阴阳性分类，采用k-均值的聚类算法，通道1、通道2阴阳性两两组合，结果如图17所示的散点图，各类的样点数显示在图中右上角方框内，由图17看出(为显示清楚，在说明书附图中将图17左旋90°)，所有样点分为四类，左下角是通道1、通道2均为阴性的样点(即图17中“·”)，右下角是通道1为阳性、通道2为阴性的样点(即图17中“×”)，右上角是通道1、通道2均为阳性的样点(即图17中“+”)，左上角是通道2为阳性、通道1为阴性的样点(即图17中“○”)。
其中，k均值算法流程为：
(1)选取数据空间中的k个对象作为初始中心，每个对象代表一个聚类中心；
(2)对于样本中的数据对象，根据它们与这些聚类中心的欧氏距离，按距离最近的准则将它们分到距离它们最近的聚类中心(最相似)所对应的类；
(3)更新聚类中心：将每个类别中所有对象所对应的均值作为该类别的聚类中心，计算目标函数的值；
(4)判断聚类中心和目标函数的值是否发生改变，若不变，则输出结果，若改变，则返回步骤(2)。
本实施例中，将通道1、通道2分别二分类，再两两组合。
拷贝数(m)：
泊松分布的均值和方差：
μ＝σ2＝λ(8)
e＝p(0)＝e-λ(9)
λ＝-ln(e)(10)
m＝nλ＝-nln(e)(12)
其中，n为分区数，m为样本数，λ为每个分区的平均目标数量，c为样本浓度，vd为分区体积，p(k)为一个分区恰好包含k个样本的概率，e＝p(0)为空分区域概率，
m＝-nln(p(0))/14.55(13)
通道1、通道2的拷贝数显示在x轴下方，分别如图18和图19所示。
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。
- 还没有人留言评论。精彩留言会获得点赞！