基于dPCR的新型全自动荧光信号采集分析方法与流程

本发明属于图像处理领域，特别是一种基于dpcr(digitalpolymerasechainreaction数字聚合酶链式反应)的新型全自动荧光信号采集分析方法。

背景技术：

dpcr技术的概念从提出到现在，在短短不到二十年的时间里有了飞速的发展。这项技术在生物体体外对基因进行扩增，是用于放大扩增特定的dna片段的分子生物学技术。目前，dpcr技术作为靶基因的绝对定量化检测、转基因检测、基因缺陷筛查、癌症预防、基因突变检测等应用的有力工具，在生物医学以及精准医疗方面有着极其广泛的应用。但是现在大多数研究集中在针对dpcr芯片本身进行优化与改善，而对于整个dpcr系统中后期荧光信号检测的问题却成为这一重大研究的“灰区”。

经过对现有文献资料的检索发现，文献[1](qiangyuanzhu,linqinandetal,digitalpcronanintegratedself-primingcompartmentalizationchip,labchip,labchip,2014,1176:1185)曾提到浙江大学的qiangyuanzhu等人设计了一种用pdms制作，真空封装的无阀，无泵的微芯片式dpcr，并对其进行实验验证，但是文章在阳性反应计数，样品的均匀性以及重复性上并未做进一步讨论。

文献[2](heyriesk.a.etal.,megapixeldigitalpcr,naturemethods,2011,8(8):649:651)曾提到过对热循环反应结束后荧光信号的采集与分析的说明。该文献对dpcr芯片进行区域划分，通过手动标记来进行区域识别，然后以0.5pixels/μm的精度对所划分区域进行逐行扫描，并统计荧光强度。虽然得到了最终结果，但是手动标记增加了操作的复杂性，不利于系统的集成化发展，而且逐行扫描降低了检测效率。

后期荧光信号采集与分析的过程决定着系统的检测精度，甚至整个反应过程的成败，由于现有技术中针对dpcr系统的荧光信号采集与分析算法极少，而且共聚焦显微镜的物镜范围无法覆盖整个芯片，所以这就限制了整个dpcr系统的高效化和集成化的发展，因此在pcr反应结束后对荧光信号的采集与分析做到自动化与快速高效化是一个亟待解决的问题。

技术实现要素：

本发明针对现有dpcr技术的研究中所存在的缺陷，提出一种基于dpcr的新型全自动荧光信号采集分析方法，利用图像区域划分——图像拼接——信号采集与分析的思想来提取并分析荧光信号，实现全自动化，避免杂散光信号的干扰，不仅提高了检测结果的精确度和检测效率，而且为dpcr技术的高效化和集成化发展提供了思路。

本发明实现上述目的的技术方案是，一种基于dpcr的新型全自动荧光信号采集分析方法，其创新点在于：该方法的步骤如下：

步骤1：dpcr芯片的制作。dpcr芯片为我们利用微纳米加工技术制作出的带有特殊标记的高密度微孔型dpcr芯片，参照图2，其整体尺寸为n*n㎜，区域a为微反应腔的部分，边缘尺寸为m*m㎜，在dpcr芯片边缘和微反应腔a边缘构成的环形区域内有六个特殊的十字标记，参照图2中p1，p2，p3，p4，p5，p6，规定dpcr芯片上与p1最接近的顶点为坐标原点，p1、p2中心点连线并指向p2的方向为x正方向，p1、p5中心点连线并指向p5的方向为y正方向，若设p1的中心坐标为p1(a,b)(0<a,b<(n-m)/2)，则其他五点的中心坐标分别为p2(n-a,b)，p3(a，n/2)，p4(n-a，n/2)，p5(a，n-b)，p6(n-a，n-b)；该六个十字标记的中心所构成的矩形区域的大小为(n-2a)*(n-2a)。

步骤2：创建十字标记识别数据库。对步骤1中的dpcr芯片拍摄，得到各式各样(倍率大小、清晰度、十字标记完整度不完全相同)的阳性图片和阴性图片，其中，阳性图片是指含有十字标记的图片，阴性图片是指不含有十字标记的图片。通过设定所有阳性图片的rois(regionofinterest的缩写，指含有十字标记的矩形区域)，创建十字标记的识别数据库。

步骤3：初始化，设定各方向和各顶点，采集图像。该步骤由以下子步骤完成。

子步骤1：参照图2，以p1、p2的中心点连线的方向作为水平方向，以p1、p5的中心点连线的方向作为竖直方向；将图像信息转换成矩阵，根据ccd照相机下一张图片的像素大小e*f与相应的dpcr芯片微反应腔的尺寸d，找到像素与单位微米之间的换算关系；利用ccd照相机下一张图片的像素大小e*f、dpcr芯片的尺寸n*n以及像素和微米之间的换算关系计算得到划分采集区域的总数目snr(1)*snr(2)(分别代表在水平方向和竖直方向所需采集照片的数目)和可编程微动控制平台的步进量stepsize(1)和stepsize(2)(分别代表在水平方向和竖直方向上的步进量)；由e-stepsize(1)和f-stepsize(2)得到相邻两张图片之间的重叠距离overlap(1)和overlap(2)(分别代表在水平方向和竖直方向上相邻两张图片的重叠距离)，判断overlap(1)>e/4～e/3并且overlap(2)>f/4～f/3是否成立，如果成立，进行子步骤2，如果不成立，增大snr并按照子步骤1所述方法重新计算stepsize和overlap，直至满足条件为止，进行子步骤2；

子步骤2：从芯片的四个顶点中任选一个作为初始顶点，与该初始顶点最相近的一个十字标记记为pz(z的取值为1,2,5,6)，将含有初始顶点、pz和部分微反应腔等三部分的图片作为第一张拍摄图片，该图片的位置即当前位置，当前位置坐标由初始顶点坐标来表示；

子步骤3：初始化当前位置为初始位置，初始位置坐标由当前位置坐标表示，从初始位置开始遍历当前dpcr芯片的所有区域，遵循先水平方向后竖直方向的原则逐行扫描，扫描过程中，根据初始位置坐标及可编程微动控制平台的步进量确定第j步移动后的理论位置坐标(xj0,yj0)，同时可编程微动控制平台根据上述计算的步进量stepsize完成第j步移动，移动后的位置坐标记为(xj0’,yj0’)，判断(xj0,yj0)与(xj0’,yj0’)是否相同，若相同，进行拍照等候下一步移动；若不同，可编程微动控制平台连续五次尝试移动到(xj0,yj0)位置，在五次移动过程中，每次移动后的位置坐标记为(xji,yji)(i＝1,2,3,4,5),第i次移动完成后，判断(xji,yji)与(xj0,yj0)是否相同，并得到实际位置(xji,yji)与理论位置(xj0,yj0)的偏差er，若(xji,yji)与(xj0,yj0)相同，则进行拍照并进行第j+1步移动，若五次移动完成后(xji,yji)与(xj0,yj0)仍不同，则根据偏差er自动调整至理论位置(xj0,yj0)，再进行拍照并进行第j+1步移动；

子步骤4：扫描全部完成，得到图片o11、o12、……、osnr(2)snr(1)，进行步骤4。

步骤4：在水平方向和竖直方向拼接图片，每个方向上两张待拼接图片的位置分别有两种情况，参照图3(a)(b)和图4(a)(b)，由以下子步骤来完成：

子步骤1：取第t(t＝1,2,……,snr(2))行的前两张图片ot1、oti(i＝2,3,……，snr(1))，设置具有1/30～1/50*overlap个像素点宽度的比较块，在图片ot1、oti之间进行相关性比较，将k行q列比较区域内的相关性记为r1、r2、……、rn，找到最大的相关性rmax＝max{r1、r2、……、rn}和重叠位置第k’行第q’列(参照图3a1，ot1的大小为e’*f’)，则取ot1的(k’+1：f’)行，(1，q’)列，取oti的(1：f-k’)行，(1，e)列，拼接成一张新的图片(参照图3a2)并删除原来的图片ot1，将新的图片作为第t行第一张图片，记为ot1，再将ot1与下一张图片oti按照上述步骤进行拼接，依次类推，直至该行所有snr(1)张图片全部拼接完成；得到每行拼接完成的最新图片o11，o21，……，ot1，然后转子步骤2。

子步骤2：记o11，o21，……，ot1分别为o1，o2，……，ot，取o1、oj(j＝2，3，……，t)两张图片，设置具有1/30～1/50*overlap个像素点宽度的比较块，在图片o1、oj之间进行相关性比较，将kk行qq列比较区域内的相关性记为r11、r22、……、rnn，找到最大的相关性rmmax＝max{r11、r22、……、rnn}和重叠位置第k’’行第q’’列(参照图4a1)，则取o1的(1：k’’)行，(1，q’’)列，取oj的(1：f2)行，(q’’+1，e2)列，拼接成一张新的图片(参照图4a2)并删除原来的图片o1，将新的图片作为当前第一张图片，重新记为o1，再将o1与下一张图片oj按照子步骤2进行拼接，直至所有snr(2)张图片全部拼接完成，得到最新的完整图片o1，参照图5，其大小记为n’*m’。

步骤5：参照图5，将图片o1划分为六个区域，分别记为m1，m2，m3，m4，m5，m6，每个区域的大小均为n’/2*m’/3，该六个区域内各含有一个十字标记；取mi(i＝1，2,3,4,5,6)，根据步骤2创建的十字标记识别数据库识别mi中的十字标记，并利用canny算子寻找十字标记的边缘，通过霍夫转换检测并画出相关直线，从而得到mi中十字标记的四条边缘直线，记为li1，li2，li3，li4(参照图6)；最后根据解析几何的算法计算出四条直线的四个交叉点，记为vi1，vi2，vi3，vi4(参照图6)，从而得到六个标记的中心坐标，记为c1(xc1，yc1)，c2(xc2，yc2)，c3(xc3，yc3)，c4(xc4，yc4)，c5(xc5，yc5)，c6(xc6，yc6)(参照图7)。

步骤6：由tan^-1((ycj-yci)/(xcj-xci))(j＝2,4,6；i＝1,3,5)计算水平方向上c1和c2、c3和c4、c5和c6的相对偏转角度并记为α1、α2、α3，由tan^-1((xcj-xci)/(ycj-yci))(j＝3,4,5，6；i＝1,2，3，4)计算竖直方向上c1和c3、c1和c5、c3和c5、c2和c4、c2和c6、c4和c6的相对偏转角度并记为α4、α5、α6、α7、α8、α9，分别取其平均值，即αx＝(α1+α2+α3)/3，αy＝(α4+α5+α6+α7+α8+α9)/6，αx和αy即为图片o1在水平方向和竖直方向上的偏转角度，取α＝(αx+αy)/2，则α即为o1的偏转角度，根据yc1-yc2与零的关系确定o1的偏转方向，从而将o1按顺时针(或逆时针)方向旋转角度α，得到新的图片(参照图8)，重新记为o1，其大小记为n”*m”。

步骤7：重复步骤5，得到新图片o1上六个标记的中心坐标，重新记为c1(xc1，yc1)，c2(xc2，yc2)，c3(xc3，yc3)，c4(xc4，yc4)，c5(xc5，yc5)，c6(xc6，yc6)，并转至步骤8；

步骤8：参照图8，裁剪掉o1上由ci(i＝1,2，……，6)所构成的矩形区域g以外的部分，并将g以内的部分放大(或缩小)，得到新的图片记为oo，使oo的尺寸与步骤1中所述dpcr芯片同一部分的尺寸相同，为(n-2a)*(n-2a)。

步骤9：参照图9，以dpcr芯片上十字标记p1的中心坐标作为基点坐标，记为(0,0)，第i(i＝1,2，……，n-1，n，n为微反应腔的总数目)个微反应腔的圆心为中心坐标，作直径为d的圆，记为圆cei，圆cei的位置即为图片oo上第i个微反应腔的位置，将圆cei内(包括边缘)所有点的值设置为i，其余值均设置为0，建立(n-2a)*(n-2a)的矩阵b，设矩阵b中数据i的位置为b1行b2列，则提取图片oo中b1行b2列的数据fy，取fy的平均值并记为fi，fi即为dpcr芯片上第i个微反应腔的荧光强度，统计相同fi的个数并记为numi,从而得到实验结果，实验结果由荧光强度fi及其个数numi表示。

本发明的有益效果是：(1)首次以工程化的方式将荧光信号采集分析方法应用于dpcr后期的结果检测，使之成为一个完整的工程检测系统，提出一种新型的全自动荧光信号采集分析方法，目前，国内外尚未见到类似于本发明的文献报道。(2)通过合理的设计使整个系统实现了全自动化，避免了人工干涉例如手动标记、逐行扫描等等所产生的测量误差，准确性好，操作控制方便，有效提高了系统的检测精度和检测效率。(3)通过上述拼接步骤可以将所采集的零散图像完美拼接成设计时的原图像大小，高度还原原图，方便后期对荧光信号的分析与检测，提高了检测的精确度。

附图说明

图1为本发明的流程图。

图2为dpcr芯片的结构图。

图3(a)、(b)为行拼接时两张图片的相对位置情况图；(a1)、(b1)为行拼接时比较块所寻找到的重叠位置示意图；(a2)、(b2)为同一行两张图片拼接起来的示意图。

图4(a)、(b)为列拼接时两张图片的相对位置情况图；(a1)、(b1)为列拼接时比较块所寻找到的重叠位置示意图；(a2)、(b2)为同一列两张图片拼接起来的示意图。

图5为全部拼接完成后的图片及六个含有十字标记区域的划分示意图。

图6为图片mi中检测到的十字标记的四条直线和四条直线相交的四个交点。

图7为所检测的六个十字标记的中心图。

图8为拼接完成的图片旋转后的图。

图9为根据dpcr芯片设计大小创建微反应腔位置矩阵的示意图。

图10为微反应腔为50μm的dpcr芯片的结构图。

图11为实验中第一行第一张图片和第二张图片的拼接结果。

图12(a)为实验中利用十字标记识别数据库识别mi中的十字标记的图片；(b)为实验中从mi中提取出来十字标记的部分。

图13为实验中检测并标记出m1中边缘直线及交叉点的十字标记图。

图14为实验中检测出的dpcr芯片上六个标记的中心坐标图。

图15为实验中所建微反应腔的位置矩阵与实际位置吻合的部分图片。

图16为实验检测结果。

具体实施方式

以下结合附图对本发明做进一步详细阐述。

一种基于dpcr的新型全自动荧光信号采集分析方法的步骤依次为：

步骤1：dpcr芯片的制作。我们利用微纳米加工技术制作出了带有特殊标记的高密度微孔型dpcr芯片，如图10：其整体尺寸为9.1*9.1㎜，中间微反应腔的部分的边缘尺寸为8380*8380μm，在dpcr芯片边缘和微反应腔a边缘构成的环形区域内有六个特殊的十字标记，记为p1，p2，p3，p4，p5，p6，规定p1、p2中心点连线并指向p2的方向为x正方向，p1、p5中心点连线并指向p5的方向为y正方向，若设p1的中心坐标为p1(200,200)(单位为μm)，则其他五点的中心坐标分别为p2(8900,200)，p3(200，4550)，p4(8900，4550)，p5(200，8900)，p6(8900，8900)(单位为μm)；参照图10，每一个微反应腔的尺寸为50μm，相邻两个微反应腔之间的距离(ef)为55μm，相邻两行中相邻三个微反应腔的圆心(a、b、c三点)构成等边三角形，以等六块的形式分布于dpcr芯片上，每一块中每一行有76个微反应腔，有29*2行，则该芯片上共有26448(76*29*2*6)个微反应腔。

步骤2：创建十字标记识别数据库。对dpcr芯片进行拍摄，得到218张各式各样(倍率大小、清晰度、十字标记完整度不完全相同)的阳性图片和阴性图片，其中，阳性图片是指含有十字标记的图片，有109张，阴性图片是指不含有十字标记的图片，有109张。通过设定109张阳性图片的rois，创建十字标记的识别数据库。

步骤3：初始化，设定各方向和各顶点，采集图像。该步骤由以下子步骤完成。

子步骤1：参照图2，以p1、p2的中心点连线的方向作为水平方向，以p1、p5的中心点连线的方向作为竖直方向；将图像信息转换成矩阵，根据ccd照相机下一张图片的像素大小1938*1460与相应的dpcr芯片微反应腔的尺寸50，找到像素与单位微米之间的换算关系；利用ccd照相机下一张图片的像素大小、dpcr芯片的尺寸以及像素和微米之间的换算关系计算得到划分采集区域的总数目snr为5*7(分别为水平方向和竖直方向所需采集照片的数目)和可编程微动控制平台的步进量stepsize分别为1820和1300(分别代表在水平方向和竖直方向上的步进量)；由e-stepsize(1)和f-stepsize(2)得到相邻两张图片之间的重叠距离overlap为118和160(分别代表在水平方向和竖直方向上相邻两张图片的重叠距离)，判断overlap(1)>300并且overlap(2)>500是否成立，如果成立，进行子步骤2，如果不成立，增大snr并按照子步骤1所述方法重新计算stepsize、和overlap，直至满足条件为止，得到snr分别为7和9，进行子步骤2；

子步骤2：从芯片的四个顶点中选取p1作为初始顶点，与该初始顶点最相近的一个十字标记记为p1，将含有初始顶点、p1和部分微反应腔等三部分的图片作为第一张拍摄图片，该图片的位置即当前位置，当前位置坐标由初始顶点坐标来表示；

子步骤3：初始化当前位置为初始位置，初始位置坐标由当前位置坐标表示，从初始位置开始遍历当前dpcr芯片的所有区域，遵循先水平方向后竖直方向的原则逐行扫描，扫描过程中，根据初始位置坐标及可编程微动控制平台的步进量确定第j步移动后的理论位置坐标(xj0,yj0)，同时可编程微动控制平台根据上述计算的步进量stepsize完成第j步移动，移动后的位置坐标记为(xj0’,yj0’)，判断(xj0,yj0)与(xj0’,yj0’)是否相同，若相同，进行拍照等候下一步移动；若不同，可编程微动控制平台连续五次尝试移动到(xj0,yj0)位置，在五次移动过程中，每次移动后的位置坐标记为(xji,yji)(i＝1,2,3,4,5),第i次移动完成后，判断(xji,yji)与((xj0,yj0)是否相同，并得到实际位置与理论位置的偏差er，若(xji,yji)与(xj0,yj0)相同，则进行拍照并进行第j+1步移动，若五次移动完成后(xji,yji)与(xj0,yj0)仍不同，则根据偏差er自动调整至理论位置(xj0,yj0)，再进行拍照并进行第j+1步移动；

子步骤4：扫描全部完成，得到图片o11、o12、……、o97，进行步骤4。

步骤4：在水平方向和竖直方向拼接图片，每个方向上两张待拼接图片的位置分别有两种情况，参照图3和图4，由以下子步骤来完成：

子步骤1：取第t(t＝1,2,……,9)行的前两张图片ot1、oti(i＝2,3,……，7)，设置具有10个像素点宽度的比较块在图片ot1、oi之间进行相关性比较，将k行q列比较区域内的相关性记为r1、r2、……、rn，找到最大的相关性rmax＝max{r1、r2、……、rn}和重叠位置第k’行第q’列(参照图3a1，ot1的大小为e’*f’)，则取ot1的(k’+1：f’)行，(1，q’)列，取oti的(1：f-k’)行，(1，e)列，拼接成一张新的图片并删除原来的图片ot1，将新的图片作为当前第一张图片，记为ot1(参照图11，第一行前两张图片的拼接结果)，再将ot1与下一张图片oti按照上述步骤进行拼接，依次类推，直至第t行的所有7张图片全部拼接完成；得到当前最新图片o11，o21，……，ot1，然后转子步骤2。

子步骤2：记o11，o21，……，ot1分别为o1，o2，……，ot，取o1、oj(j＝2，3，……，t)两张图片，设置具有10个像素点宽度的比较块在图片o1、oj之间进行相关性比较，将kk行qq列比较区域内的相关性记为r11、r22、……、rnn，找到最大的相关性rmmax＝max{r11、r22、……、rnn}和重叠位置第k”行第q”列，则取o1的(1：k”)行，(1，q”)列，取oj的(1：f2)行，(q”+1，e2)列，拼接成一张新的图片并删除原来的图片o1，将新的图片作为当前第一张图片，重新记为o1，再将o1与下一张图片oj按照子步骤2进行拼接，直至9张图片全部拼接完成，得到最新的完整图片o1，其大小记为n’*m’。

步骤5：参照图5，将图片o1划分为六个区域，分别记为m1，m2，m3，m4，m5，m6，每个区域的大小均为n’/2*m’/3，该六个区域内各含有一个十字标记；取mi(i＝1，2,3,4,5,6)，根据步骤2创建的十字标记识别数据库识别mi中的十字标记(参照图12)，并利用canny算子寻找十字标记的边缘，通过霍夫转换检测并画出相关直线，从而得到mi中十字标记的四条边缘直线，记为li1，li2，li3，li4(参照图13,为第一个十字标记m1中检测到的直线及其交点)；最后根据解析几何的算法计算出四条直线的四个交叉点，记为vi1，vi2，vi3，vi4(参照图13)，从而得到六个标记的中心坐标，记为c1(xc1，yc1)，c2(xc2，yc2)，c3(xc3，yc3)，c4(xc4，yc4)，c5(xc5，yc5)，c6(xc6，yc6)(参照图14)。

步骤6：由tan^-1((ycj-yci)/(xcj-xci))(j＝2,4,6,i＝1,3,5)计算水平方向上c1和c2、c3和c4、c5和c6的相对偏转角度并记为α1、α2、α3，由tan^-1((xcj-xci)/(ycj-yci))(j＝3,4,5，6,i＝1,2，3，4)计算垂直方向上c1和c3、c1和c5、c3和c5、c2和c4、c2和c6、c4和c6的相对偏转角度并记为α4、α5、α6、α7、α8、α9，分别取其平均值，即αx＝(α1+α2+α3)/3，αy＝(α4+α5+α6+α7+α8+α9)/6，αx和αy即为图片o1在水平方向和竖直方向上的偏转角度，取α＝(αx+αy)/2，则α即为o1的偏转角度，根据yc1-yc2与零的关系确定o1的偏转方向，从而将o1按逆时针方向旋转角度α，得到新的图片，重新记为o1，其大小记为n”*m”。

步骤7：重复步骤5，得到新图片o1上六个标记的中心坐标，重新记为c1(xc1，yc1)，c2(xc2，yc2)，c3(xc3，yc3)，c4(xc4，yc4)，c5(xc5，yc5)，c6(xc6，yc6)，并转至步骤9；

步骤8：参照图8，裁剪掉o1上ci(i＝1,2，……，6)所构成的矩形区域g以外的部分，并将g以内的部分放大(或缩小)，得到新的图片记为oo，使oo的尺寸与步骤1中所述dpcr芯片同一部分的尺寸相同，为8700*8700。

步骤9：参照图9，以dpcr芯片上十字标记p1的中心坐标作为基点坐标，记为(0,0)，第i(i＝1,2，……，26448)个微反应腔的圆心为中心坐标，作直径为d的圆，记为圆cei，圆cei的位置即为图片oo上第i个微反应腔的位置(参照图15)，将圆cei内(包括边缘)所有点的值设置为i，其余值均设置为0，建立8700*8700的矩阵b，设矩阵b中数据i的位置为b1行b2列，则提取图片oo中b1行b2列的数据fy，取fy的平均值并记为fi，fi即为dpcr芯片上第i个微反应腔的荧光强度，统计相同fi的个数并记为numi,从而得到实验结果，实验结果由荧光强度fi及其个数numi表示，参照图16。