一种微反应腔阵列式芯片的全自动快速荧光分析方法

摘要

本发明属于微流控图像处理领域，特别是一种微反应腔阵列式芯片的全自动快速荧光分析方法。该方法利用图像区域划分——图像拼接——信号采集与分析的思想来提取并分析荧光信号，拼接图片时，先将第一行第二行图片拼接起来，然后将第三行图片拼接起来，再与第一二行所拼接成的图片拼接，再进行下一行的拼接，提高了后期荧光检测准确率；该方法实现了微反应腔阵列式芯片的全自动快速荧光分析过程中的全自动化，并且快速提取微反应腔阵列式芯片上的荧光信息，避免杂散光信号的干扰，不仅提高了检测结果的精确度和检测效率，而且为整个系统的快速化、高效化、集成化发展提供了思路。

说明书

技术领域

本发明属于微流控图像处理领域，特别是一种微反应腔阵列式芯片的全自动快速荧光分析方法。

背景技术

荧光信号的检测和微流体芯片的结合已经成为诸多跨学科领域交叉解决问题的典型事件，将生物信息转化成可分析的荧光信号已广泛应用在多个领域，以数字PCR为例：数字PCR技术的概念从提出到现在，在短短不到二十年的时间里有了飞速的发展。整体来说，数字PCR技术包括前期扩增和后期检测两大部分。我们观察到，数字PCR技术的研究并没有全方面的展开，多数研究集中在芯片本身的优化和改善，但对于整个数字PCR系统中的后期荧光检测问题却几乎无人问津。

文献[1](Heyries K.A.et al.,Megapixel digital PCR,Nature Methods,2011,8(8):649:651)曾提到过对热循环反应结束后荧光信号的采集与分析的说明。该文献对dPCR芯片进行区域划分，通过手动标记来进行区域识别，然后以0.5pixels/μm的精度对所划分区域进行逐行扫描，并统计荧光强度，但是每一块区域的边缘都忽略了。虽然得到了最终结果，但是手动标记增加了操作的复杂性，不利于系统的集成化发展，而且逐行扫描降低了检测效率。

申请号为201710789184.8的专利提出了一种基于dPCR的新型全自动荧光信号采集分析方法，该专利对dPCR芯片进行分区域采集图片，自动识别拼接并统计荧光强度的方法，实现了全自动化，但是采集图片过程中每一步做到理论位置与实际位置重合，耗时长；自动识别拼接时采用先行后列的方式，将行与行之间的图片单独分隔开；在提取荧光强度时首先根据设计尺寸建立了一个同样的矩阵，其次根据一一对应关系提取照片上相应位置的数据。

对于类似于数字PCR这样的反应来说，是一个终端检测的过程，因此后期荧光检测的速度与准确性会决定整个系统的精度，甚至整个反应过程的成败，由于现有技术中针对微反应腔阵列式芯片的荧光信号采集与分析算法极少，而且共聚焦显微镜的物镜范围无法覆盖整个芯片，所以这就限制了整个系统的高效化和集成化的发展，因此在微反应腔阵列式芯片上对荧光信号的采集与分析做到快速高效化是一个亟待解决的问题。

发明内容

本发明针对现有技术研究中所存在的缺陷，提出一种微反应腔阵列式芯片的全自动快速荧光分析方法，利用图像区域划分——图像拼接——信号采集与分析的思想来提取并分析荧光信号，实现全自动化，并且快速提取微反应腔阵列式芯片上的荧光信息，避免杂散光信号的干扰，不仅提高了检测结果的精确度和检测效率，而且为整个系统的快速化、高效化、集成化发展提供了思路。

本发明实现上述目的的技术方案是，一种微反应腔阵列式芯片的全自动快速荧光分析方法，步骤如下：

步骤1：微反应腔阵列式芯片的制作。微反应腔阵列式芯片为我们利用微纳米加工技术制作出的带有特殊标记的微孔型芯片，参照图2，其整体尺寸为N*N㎜，区域A为微反应腔的部分，边缘尺寸为M*M㎜，区域A分为六块，区域A分为T块，分别用A1、A2、……、AT表示，T块Ai结构大小(大小为m*n)完全相同，并且相邻两块之间的间隔为1～22μm之间，在微反应腔阵列式芯片边缘和微反应腔A边缘构成的环形区域内有六个特殊的十字标记，参照图2中P1，P2，P3，P4，P5，P6，规定微反应腔阵列式芯片上与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中的微反应腔阵列式芯片填充荧光溶液，并拍摄得到各式各样(倍率大小、清晰度、十字标记完整度不完全相同)荧光模式下的阳性图片和阴性图片，其中，阳性图片是指含有十字标记的图片，阴性图片是指不含有十字标记的图片。通过设定所有阳性图片的ROIs(Region Of Interest的缩写，指含有十字标记的矩形区域)，创建十字标记的识别数据库。

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

子步骤1：参照图2，以P1、P2的中心点连线的方向作为水平方向，以P1、P5的中心点连线的方向作为竖直方向；将图像信息转换成矩阵，根据数字影像记录仪(例如CCD照相机)下一张图片的像素大小e*f与相应的微反应腔阵列式芯片上微反应腔的尺寸d，找到像素与单位微米之间的换算关系；利用数字影像记录仪(例如CCD照相机)下一张图片的像素大小e*f、微反应腔阵列式芯片的尺寸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：初始化当前位置为初始位置，初始位置坐标由当前位置坐标表示，从初始位置开始遍历当前微反应腔阵列式芯片的所有区域，遵循先水平方向后竖直方向的原则逐行扫描，扫描过程中，根据初始位置坐标及可编程微动控制平台的步进量确定第j步移动后的理论位置坐标(Xj0,Yj0)，同时可编程微动控制平台根据上述计算的步进量stepSize完成第j步移动，移动后的位置坐标记为(Xj0’,Yj0’)；

子步骤4：扫描全部完成，得到图片O11、O12、……、OsNr(2)sNr(1)，进行步骤4。

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

子步骤1：取第t(t＝1,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)张图片全部拼接完成；

子步骤2：记Ot1为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，

子步骤3：取第t(t＝3,……,sNr(2))行的前两张图片Ot1、Oti(i＝2,3,……，sNr(1))，按照子步骤1所述内容分别将每一行图片拼接起来，得到该行图片，分别记为O31，O41，……，Ot1，并跳转子步骤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：取c1，c2纵坐标的最小值作为首行，取c5，c6纵坐标的最小值作为末行，取c1，c3纵坐标的最小值作为首列，取c4，c6纵坐标的最小值作为最后一列，裁剪得到新的图片G，并将图片G缩放得到新的图片记为OO，使OO的尺寸为(N-2a)*(N-2a)。

步骤9：在图片OO上，根据设计尺寸，以微反应腔阵列式芯片上十字标记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的最大值Fmax，计算fi＝Fi/Fmax，则fi即为微反应腔阵列式芯片上第i个微反应腔的相对荧光强度，统计相同fi的个数并记为Numi,从而得到实验结果，实验结果由荧光强度fi及其个数Numi表示。

本发明的有益效果是：(1)在前期微反应腔阵列式芯片的设计上，将微反应腔部分划分为T块，中间有微小距离的间隔，这样在极微小的反应腔的情况下，可以用块状分析方法对每一块的荧光强度统计；(2)在采集图片过程中，微控制平台按顺序移动并拍照，可在十分钟内采集完所有图片，即在十分钟内获得微反应腔阵列式芯片上所呈现的结果；(3)拼接图片时，先将第一行第二行图片拼接起来，然后将第三行图片拼接起来，再与第一二行所拼接成的图片拼接，再进行下一行的拼接，提高了后期荧光检测准确率；(4)在将微反应腔部分缩放至设计尺寸时，采取第1,3,4,6个十字标记的中心坐标来作为新图片行列的区间，使之更精确。该方法在十分钟内就可以得到结果，大大提高了结果获取的快速性与准确性，并且更加快捷方便，为快速数字微流体诊断系统奠定了基础。

附图说明

图1为本发明的流程图。

图2为微反应腔阵列式芯片的结构图。

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

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

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

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

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

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

图9为微反应腔为20μm的微反应腔阵列式芯片的结构图。

图10(a)为实验中全部63张图片的拼接结果；(b)为(a)中黑框部分的放大图。

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

图12为实验中旋转后的拼接图。

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

图14为实验检测结果。

具体实施方式

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

一种微反应腔阵列式芯片的全自动荧光分析方法的步骤依次为：

步骤1：微反应腔阵列式芯片的制作。微反应腔阵列式芯片为我们利用微纳米加工技术制作出的带有特殊标记的微孔型芯片，参照图10，其整体尺寸为9.1*9.1㎜，区域A为微反应腔的部分，边缘尺寸为8380*8380μm，区域A分为六块，分别用A1、A2、A3、A4、A5、A6表示，六块Ai结构大小(大小为m*n)完全相同，并且相邻两块之间的间隔为5～22μm之间，在微反应腔阵列式芯片边缘和微反应腔A边缘构成的环形区域内有六个特殊的十字标记，参照图2中P1，P2，P3，P4，P5，P6，规定微反应腔阵列式芯片上与P1最接近的顶点为坐标原点，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，每一个微反应腔的尺寸为20μm，相邻两个微反应腔之间的距离(EF)为24μm，相邻两行中相邻三个微反应腔的圆心(A、B、C三点)构成等边三角形，以等六块的形式分布于微反应腔阵列式芯片上，每一块中每一行有174个微反应腔，有67*2行，则该芯片上共有139896个微反应腔。

步骤2：创建十字标记识别数据库。对步骤1中的微反应腔阵列式芯片填充荧光溶液，并拍摄得到82张各式各样(倍率大小、清晰度、十字标记完整度不完全相同)荧光模式下的阳性图片和阴性图片，其中，阳性图片是指含有十字标记的图片，有32张，阴性图片是指不含有十字标记的图片，有50张。通过设定所有阳性图片的ROIs(Region Of Interest的缩写，指含有十字标记的矩形区域)，创建十字标记的识别数据库。

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

子步骤1：参照图2，以P1、P2的中心点连线的方向作为水平方向，以P1、P5的中心点连线的方向作为竖直方向；将图像信息转换成矩阵，根据CCD照相机下一张图片的像素大小1938*1460与相应的微反应腔阵列式芯片微反应腔的尺寸50，找到像素与单位微米之间的换算关系；利用CCD照相机下一张图片的像素大小、微反应腔阵列式芯片的尺寸以及像素和微米之间的换算关系计算得到划分采集区域的总数目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：初始化当前位置为初始位置，初始位置坐标由当前位置坐标表示，从初始位置开始遍历当前微反应腔阵列式芯片的所有区域，遵循先水平方向后竖直方向的原则逐行扫描，扫描过程中，根据初始位置坐标及可编程微动控制平台的步进量确定第j步移动后的理论位置坐标(Xj0,Yj0)，同时可编程微动控制平台根据上述计算的步进量stepSize完成第j步移动，移动后的位置坐标记为(Xj0’,Yj0’)；

子步骤4：扫描全部完成，得到图片O11、O12、……、O97，进行步骤4。

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

子步骤1：取第t(t＝1,2)行的前两张图片Ot1、Oti(i＝2,3,……，9)，设置具有10个像素点宽度的比较块，在图片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)列，拼接成一张新的图片并删除原来的图片Ot1，将新的图片作为第t行第一张图片，记为Ot1，再将Ot1与下一张图片Oti按照上述步骤进行拼接，依次类推，直至该行所有7张图片全部拼接完成；

子步骤2：记Ot1为Ot，取O1、Oj(j＝2，3，……，t)两张图片，设置具有10个像素点宽度的比较块，在图片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)列，拼接成一张新的图片并删除原来的图片O1，将新的图片作为当前第一张图片，重新记为O1，

子步骤3：取第t(t＝3,……,sNr(2))行的前两张图片Ot1、Oti(i＝2,3,……，sNr(1))，按照子步骤1所述内容分别将每一行图片拼接起来，得到该行图片，分别记为O31，O41，……，Ot1，并跳转子步骤2进行列的拼接，直至所有9张图片全部拼接完成，得到最新的完整图片O1，参照图10((a)为拼接后的整张图片，包含63张图片，(b)为(a)中黑框部分放大后的情况)，其大小记为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中的十字标记(参照图11)，并利用canny算子寻找十字标记的边缘，通过霍夫转换检测并画出相关直线，从而得到Mi中十字标记的四条边缘直线，记为li1，li2，li3，li4；最后根据解析几何的算法计算出四条直线的四个交叉点，记为vi1，vi2，vi3，vi4，从而得到六个标记的中心坐标，记为c1(Xc1，Yc1)，c2(Xc2，Yc2)，c3(Xc3，Yc3)，c4(Xc4，Yc4)，c5(Xc5，Yc5)，c6(Xc6，Yc6)。

步骤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按顺时针(或逆时针)方向旋转角度α，得到新的图片(参照图12)，重新记为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：取c1，c2纵坐标的最小值作为首行，取c5，c6纵坐标的最小值作为末行，取c1，c3纵坐标的最小值作为首列，取c4，c6纵坐标的最小值作为最后一列，裁剪得到新的图片G，并将图片G缩放得到新的图片记为OO，使OO的尺寸为(N-2a)*(N-2a)。

步骤9：在图片OO上，根据设计尺寸，以微反应腔阵列式芯片上十字标记P1的中心坐标作为基点坐标，记为(0,0)，第i(i＝1,2，……，n-1，n，n为微反应腔的总数目)个微反应腔的圆心为中心坐标，作内直径为d1，外直径为d2的圆(d1<d<d2)，记为圆Ci(参照图13)，将圆Ci内(不包括内外边缘以及两边边缘所构成的区域)所有点的值提取出来相加，并记为Fi，寻找Fi的最大值Fmax，计算fi＝Fi/Fmax，则fi即为微反应腔阵列式芯片上第i个微反应腔的相对荧光强度，统计相同fi的个数并记为Numi,从而得到实验结果，实验结果由荧光强度fi及其个数Numi表示，参照图14。