CN102279146A - 基于激光鞘流技术的血液细胞五分类方法 - Google Patents

Abstract

本发明公开了一种基于激光鞘流技术的血液细胞五分类方法，包括以下步骤：A、收集每个细胞或粒子的至少4路光信号，并表征成与各细胞或粒子的各路光信号强度相关的四维初始向量；B、以加权欧式距离的K均值聚类法进行四类细胞集合的初次聚类；C、对各聚类排序，得到W₁、W₂、W₃和W₄四个聚类，删除W₁聚类；D、将W₂和W₄两个聚类的质心经过空间坐标变换得到用于产生分类计算使用的特征向量的旋转矩阵和旋转向量；E、用旋转向量或旋转矩阵将各细胞的初始向量转换成特征向量，计算并提取相应空间坐标内细胞特征向量的统计特征，通过决策树对细胞进行五分类。本方法无需使用特殊试剂，而且计算效率高、分类结果更为准确。

Description

基于激光鞘流技术的血液细胞五分类方法

技术领域

本发明涉及血液细胞的分类方法，具体涉及基于激光鞘流技术的血液细胞五分类方法。

背景技术

基于流式细胞技术的血液分析仪都是通过收集和分析细胞的二维或二维以上的多维特征来识别包含在被鞘液包裹的悬液中的细胞以将它们划分成不同的类别。在基于流式细胞技术的血液分析仪中，鞘流包裹着细胞悬液匀速通过一个光照区域，在光照区域中细胞逐个的受到激光的照射以及经偏振光产生通路处理的而产生不同的光信号，如前向透射光，小角度散射光，大角度散射光和大角度消偏振振光等。各路聚光透镜将这些光信号分别汇聚，检测器，通常是PMT(光电倍增管)分别收集这些来自同一个细胞的特征信号。分析系统将由这些检测器收集来的信号生成二维、三维或高维的散点图，可以在散点图上划分多个区域，根据表征细胞的多维特征向量落在散点图上的分布趋势和统计特性，再结合散点图的物理、生物医学上的实际解释含义，可以确定得到血液细胞五分类的结果。传统的做法基本上有两种，一种是在细胞信号特征形成的散点图上使用固定边界进行分类，另外一种是在经过简单投影变换以后在投影图上查找浮动阈值，然后采用浮动阈值作为边界进行分类。第一种方式对于大部分带有正常样本特征的样本能够有得到正确的分类结果，优点是实现简单、计算速度快且对计算平台的硬件要求较低，其缺点是不能针对具有特异性的样本进行边界调整，当某些样本的细胞特征显著不同于固定边界所表达的特征时会出现较大的误差，且不能真实的从宏观上反映出样本细胞特征的分布趋势及特征。例如专利US.Pat.No.4987086所披露的一种是用前向散射光和侧向散射光所形成的散点图上通过划分出固定边界的方式从全血细胞样本中区分出淋巴细胞、中性粒细胞核单核细胞的方式，该方式将特征值落在由固定边界划分的区域中的细胞识别为同一类细胞。第二种方式是对第一种方式的一种改进，较之第一种方式考虑到了不同样本细胞特征分布的差异性和趋势特征，具有较好的分类正确性，但是由于投影变换过程中受到各类细胞差异性不显著的特征值的相互干扰，不容易找到最佳的分类边界阈值，所以正确性及用于临床检验的实用性有待进一步的提高。专利US.Pat.No.4599307、4704891、4727020、4987086、6014904都披露了一些使用该方式对血液细胞样本进行识别、分类和计数的方法。

公开号为CN101226190的发明专利公开了一种流式细胞术的自动分类方法及装置，包括以下步骤：A1.根据收集地每个细胞或粒子逐个通过光照区域时产生的至少两路光信号，将每个细胞或粒子表征为一个与其各路光信号强度有关的、至少二维的向量；B1.计算所有有效细胞或粒子两两之间的距离，且距离越近，两个细胞或粒子之间的相似程度越高；C1.将相似程度高的细胞或粒子聚集成同一个类；D1.重复步骤C1，至少将所有有效细胞或粒子聚集成样本根据测量原理所应有的类别数L。该专利使用传统的均值聚类算法获得细胞或粒子的聚类数目，对于差异性较小的类别，不能准确的划分出来，导致分类不够准确，不适合应用在血液细胞五分类之中。

一些其它的血液细胞五分类技术使用附加的装置或加入特殊的试剂以获得更多的细胞特异性。如COULTER公司专利的VCS三维分析技术在流式细胞技术的基础之上使用了三个独立的能量源，以测定细胞的体积、对电磁波的传导性及光散射，通过该方式，可以获得较好的细胞分类效果。但该方法固定的对8192个细胞进行检测统计，适用于对细胞分类并不适用于进行分类的同时进行计数，且由于附加了特殊的处理部件，造成了制造和维护成本的上升。以Sysmex公司的SF-3000为代表的产品，除使用了流式细胞技术外还使用了组化染色提取、加强不同类型白细胞之间的差异，例如在嗜酸性粒细胞的检测中加入特殊的染液以改变嗜酸性粒细胞的散射光特性，同中性粒细胞区分开；在嗜碱性粒细胞的检测中加入特殊的溶血素以溶解其他的白细胞形成裸核，只保持嗜碱性粒细胞的完整性，从而使得嗜碱性粒细胞在体积大小上与其他白细胞表现出明显差异，再根据体积大小进行区分。该方法存在检测过程中需要使用特殊的试剂，用户使用成本较高的问题。

发明内容

本发明要解决的技术问题是提供一种无需使用特殊试剂，而且计算效率高、分类结果更为准确的基于激光鞘流技术的血液细胞五分类方法。

本发明所述的基于激光鞘流技术的血液细胞五分类方法，包括以下步骤：

A、收集每个细胞或粒子通过鞘流器的激光照射区域时产生的至少4路光信号，将每个细胞或粒子的表征为一个与其各路光信号强度相关的四维的初始向量，分别为前向透射光强度P0、3°至10°小角度散射光强度P10、90°大角度散射光强度P90和大角度偏正光强度P90p四个属性；

B、以加权欧式距离的K均值聚类法进行四类细胞集合的初次聚类，其中距离的权值采用变异系数赋权法得到；

C、分别提取出上述各聚类质心初始向量中的P0、P10和P90三个属性，将每个聚类质心按照P0的属性值由小到大进行排序，得到排序后的四个聚类分别编号为W₁、W₂、W₃和W₄，其中W₁是由破损细胞碎片引起的干扰信息，该信息通过第一门限判别单元G1之后被删除；

D、将W₂和W₄两个聚类的聚类质心经过空间坐标变换得到用于产生分类计算使用的特征向量的旋转矩阵和旋转向量；

E、用旋转向量或旋转矩阵将各细胞的初始向量转换成特征向量，计算并提取相应空间坐标内细胞特征向量的统计特征，通过决策树对细胞进行五分类。

上述方法中：

所述步骤D包括以下步骤：

D1、提取W₂和W₄聚类的质心Q₂(P0₂，P10₂，P90₂)与Q₄(P0₄，P10₄，P90₄)，拟合一条在由P0、P10和P90组成的三维欧式空间中通过点P₂和P₄的直线L1；

D2、将直线L1投影到由P0和P10组成的平面直角坐标系中，得到投影线L2，计算投影线L2与坐标轴P0的夹角θ；

D3、通过夹角θ计算得到P0和P10坐标系中的特征向量变换所需要的旋转矩阵 $\mathrm{Mr}1=\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right];$

D4、将直线L1投影到由P10与P90组成的平面直角坐标系中，得到投影线L3，计算投影线L3与坐标轴P90的夹角β；

D5、通过夹角β计算得到P10与P90坐标系中的特征向量变换所需要的旋转向量Vr＝(sinβ·cosβ，cosβ·cosβ)^T。

所述步骤E包括以下步骤：

E1、在P10、P90组成的坐标平面上将各细胞的初始向量用旋转向量转换成特征向量，计算过程为：设第i个细胞的初始向量中P10、P90属性特征分别为(P10_i、P90_i)，转换得到的特征向量P90_i‘＝(P10_i·sinβ·cosβ+P90_i·cosβ·cosβ)；

E2、将转换后的特征向量P90₀’，P90₁’，…，P90_i’，…，投影到坐标轴P90上，计算其投影的统计直方图H1；

E3、在统计直方图H1上查找波峰、波谷，得到波峰PK1、PK2和波谷VL1；

E4、波峰PK1、PK2和波谷VL1送第二门限判别单元G2进行评价判别，若该组峰谷值合理，则转入步骤E5，否则转入步骤E3重新识别；

E5、以波谷值VL1作为阈值，设置第一分类计算单元CL1，P90属性的特征向量值大于等于该阈值的细胞，划分为多核类细胞；小于该阈值的细胞，划分为单核类细胞；

E6、在P10、P0平面上将被划分为单核类的细胞的初始向量用旋转矩阵转换成特征向量，计算过程为：设第i个单核类细胞的初始向量的P10、P0属性特征分别为(P10_i、P0_i)，转换得到的特征向量 $({{P10}_i}^{\′},{{P0}_i}^{\′})=\left[\begin{array}{cc}{P10}_i& {P0}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right];$

E7、转化后的特征向量(P10’₁、P0’₁)，(P10’₂、P0’₂)，…，(P10’_i、P0’_i)，…，投影到坐标轴P0上，并计算其投影的统计直方图H2；

E8、在统计直方图H1上查找波峰、波谷，得到波峰PK3、PK4、PK5和波谷VL2、VL3；

E9、波峰PK3、PK4、PK5和波谷VL2、VL3送第三门限判别单元G3进行评价判别，若该组峰谷值合理，则转入步骤E10，否则转入步骤E8重新识别；

E10、以波谷值VL2、VL3作为阈值，设置第二分类计算单元CL2，将P0属性的特征向量值大于等于VL2且小于VL3的细胞，划分为淋巴细胞；而P0属性的特征向量值大于等于VL3的细胞，划分为单核细胞；

E11、提取步骤E5得到的多核类细胞中各细胞初始向量中的P90属性和P90p属性，形成分别以P90和P90p为纵、横坐标轴的散点图；

E12、对生成的散点图做变换、统计计算，得到散点数据的线性划分线L90p，用线性划分线L90p的参数设置第三分类计算单元CL3；

E13、将步骤E5得到的多核类细胞中各细胞初始向量中的各属性(即P、P10、P90和P90p属性)送入第三分类计算单元CL3进行分类计算，位于线性划分线L90p之上的细胞划分为嗜酸细胞；其余多核细胞划分为中性粒细胞；

E14、提取已分出的单核细胞聚类与中性粒细胞聚类中细胞初始向量中的P0、P10和P90属性，分别重新计算该两类细胞集合的质心Q_D(P0_D，P10_D，P90_D)和Q_Z(P0_Z，P10_Z，P90_Z)，以及两类细胞集合的标准差SD1和SD2；

E15、将Q_D、Q_Z、SD1、SD2作为参数设置第四分类计算单元CL4，并将单核细胞集合与中性粒细胞集合送入第四分类计算单元CL4进行分类计算，识别出集合中的嗜碱细胞。

在步骤E12中包括以下步骤：

F1、计算多核细胞聚类中各细胞初始向量中P90p和P90属性的均值，得到

其中m为多核细胞的数目；令SD90’＝0；

F2、在由P90和P90p属性为坐标轴构成的平面坐标系中计算得到拟合线L90a的参数

拟合线L90a即为：P90p＝A+B*P90；

F3、计算L90a与坐标轴P90的夹角γ，得到P90-P90p坐标系中的旋转矩阵 $\mathrm{Mr}2=\left[\begin{array}{cc}\cos \mathrm{\γ}& \sin \mathrm{\γ}\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right];$

F4、在P90、P90p平面上将被划分为多核类的细胞的初始向量用旋转矩阵Mr2进行旋转变换，计算过程为：设第i个多核类细胞初始向量的P90、P90p属性特征分别为(P90_i、P90p_i)，变换后得到的特征向量为 $({{P90}_i}^{\′},{{P90p}_i}^{\′})=\left[\begin{array}{cc}{P90}_i& {P90p}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\γ}& \sin \mathrm{\γ}\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right];$

F5、将变换后的各细胞的特征向量(P90’₁、P90p’₁)，(P90’₂、P90p’₂)，…，(P90’_i、P90p’_i)，…，投影到坐标轴P90p上，并计算其投影的统计直方图H3；

F6、在统计直方图H3上查找波峰、波谷，得到波峰PK6和波谷VL4；

F7、波峰PK6波谷VL4送第四门限判别单元G4进行评价判别，若该组峰谷值合理，则可以通过以下算式计算得到SD90p＝VL4-PK6；否则SD90p＝|SD/R90|；其中，

m为多核类细胞的数目；

F8、以A、B、SD90为参数，生成划分线L90p’：P90p＝A*P90+B+SD90，送第五分类计算单元CL5重新划分待分类的多核细胞集合；

F9、若SD90’＝SD90，则转入步骤F10；否则，令SD90’＝SD90，将经过步骤F8重新划分的多核细胞集合送入步骤F1，开始新的一轮迭代计算过程；

F10、以A、B、SD90为参数，生成划分线L90p：P90p＝A*P90+B+SD90，送入第三分类计算单元CL3。

本发明所述方法中，波峰、波谷的查找以常规方式进行。为了提高波峰、波谷查找的准确度，可在步骤E7和E8之间设置对直方图进行滤波的步骤。滤波可按常规方法进行，优选指数平滑滤波。

本发明所述述方法的步骤D包括以下步骤：

B1、以实验获得的应用经验值作为四个分类细胞集合的中心；

B2、计算样本集中每个细胞对象与这些中心细胞对象的加权欧式距离，并根据最小的加权欧式距离重新对相应的细胞对象进行划分；

B3、对于中心细胞对象有变化的分类，重新计算它们的中心细胞对象；B4、重复执行B2到B3的迭代循环，直到每个聚类的中心细胞对象稳定不再发生变化为止。

为了更进一步减少计算量，可在步骤A和B之间设置以下步骤：即通过设定阈值，将不符合条件的噪声或粒子的数据删除。

与现有技术相比，本发明的优点在于：1)本发明在空间投影分析计算时使用了初次聚类计算的结果，减少了迭代计算查找表征最佳分类趋势的计算量，同时提高了查找的准确性；2)采用变异系数赋权法对聚类计算时的距离进行加权，有效的减少了细胞碎片等其它干扰信息对聚类结果的影响；3)本发明将高维分布的特征向量投影到二维笛卡尔坐标系中进行处理并在各个坐标系中的计算结果之间以及计算结果与医学检验和临床试验的先验知识之间进行交叉验证，在减少计算量、提高计算效率的同时，保证了分类结果的准确性；4)无需使用特殊试剂，在不需要附加能量作用于检测装置的同时，可以准确、高效的对血液细胞进行五分类。

附图说明

图1为应用本发明所述方法的全自动血液细胞五分类分析系统结构示意图；

图2为本发明所述方法对血液细胞进行分类的流程示意图；

图3为图2所示方法中步骤S1处理过程的流程示意图；

图4为图2所示方法中经步骤S2处理后得到的聚类结果的示意图；

图5为图2所示方法中步骤S3处理过程的流程示意图；

图6为图2所示方法中步骤S4处理过程的流程示意图；

图7为嗜酸性细胞分类的示意图；

图8为利用传统方法进行分类的结果示意图；

图9为本发明所述方法对同一样本进行分类的结果示意图。

图中标号为：

1、液路控制模块；11、废液处理模块；12、流动检测室；13、进样控制模块；2、检验光束形成模块；21、激光发生器；22、检验光束形成与校准模块；3、信号检测与预处理模块；31、光信号检测器；32、信号预处理器；4、分类分析统计处理模块。

具体实施方式

本发明所述方法适用于基于流式细胞技术的全自动血液分析仪。

在本实施方式中，本发明所述方法建立在如图1所示的全自动血液细胞五分类分析系统上，当然也可建立在其它类似结构的系统上。

图1所示的全自动血液细胞五分类分析系统，包括：液路控制模块1、检验光束形成模块2、信号检测与预处理模块3和分类分析统计处理模块4。所述液路控制模块1进一步包括从真空管中采集血样并将其混匀的装置及相互连接的，所述进样控制模块13用于产生被鞘液包裹的含有待测血液细胞的样本液，并使之匀速通过流动检测室12，流动检测室12是一个透明的器件，中间有一个圆形的管道，待测试细胞在鞘液的包裹下逐个通过该管道，接受检验光束的照射。所述废液处理模块11用于排出检验过的样本废液并对进样液路进行清洗以消除样本间的交叉污染，避免因携带污染影响下一个样本检验准确性的情况。检验光束形成模块2包括一个氦-氖激光发生器21和用于将散射光整理成检验所需要光束形状及特性的检验光束形成与校准模块22，检验光束聚焦在流动检测室12的中心管道上，待测细胞在鞘液包裹下通过中心管道时产生各种角度的光信号，如前向透射光信号、小角度散射光信号，大角度散射光信号。信号检测与预处理模块3，用于收集血液细胞逐个通过激光照射区域时产生的多路光信号并对光信号进行域处理，所述信号检测与预处理模块3进一步包括用于收集光信号的光信号检测器31(通常是PMT)和用于滤除干扰信号的信号预处理器32。本发明所述方法中使用四路独立的光通道对通过检测区的细胞进行检测，分别是：前向透射光检测通道、小角度散射光检测通道、大角度散射光检测通道和大角度消偏振光检测通道，所得特征分别记为P0、P10、P90和P90p(它们即为各细胞初始向量中的四个属性)。信号预处理器32用于计算个通道光强度和使用“设门”的方式设置合理的阈值从而剔除掉不合理的信号值，减少下一级处理的运算量。分类统计处理模块4，用于根据信号检测与预处理模块3检测到的合理的光信号，将每个细胞表征为一个与其各路光信号强度相关的四维向量，并根据计算统计各个血液细胞的特征得到细胞的分布趋势从而将血液细胞进行五分类。

在基于流式细胞技术的血细胞分析仪中，每个细胞通过光检测区域时(之前经过常规的溶血操作)，都会在对应的光检测通道产生光信号，从而形成可用于表征该细胞特性的多维光信号，习惯上将这一过程称为触发了一个检测事件(Event)。假设分析仪器具有m维光信号检测通道，则当第i个细胞通过光检测区域触发事件E_i时，可得到一个m维向量E_i＝(x_i1，x_i2，Λ，x_im)，其中x_ij表示第j通道上的光信号强度，由信号预处理器32使用MIE散射算法计算得到。由于仪器工作过程中的本底噪声和使用环境叠加的一些高斯噪声，有可能在检测过程中事件会被误触发，因此通过“设门”的方式对由硬件实现的信号预处理器32进行设置，使符合设定阈值要求的事件才能够被触发，杜绝带有错误特征数据的无效事件参与后续运算，提高运算的准确性和速度。设一次测量过程结束时，通过了n个细胞，则认为触发了n次事件，从而得到细胞信息的原始数据集合

本发明所述方法对数据集合D_n×m进行分析处理的过程如图2所示，包括以下步骤：

步骤S1，采用基于加权欧式距离的K均值聚类法对细胞进行初次聚类，以便更好的提取细胞的分布趋势，具体包括以下S11～S16步骤，实施流程如图3所示：

步骤S11，以大量临床实验获得的应用经验值作为四类细胞集合的初始聚类中心，这样可以减少聚类计算收敛时的迭代次数，减少计算的时间；

步骤S12，根据每个聚类中所有细胞对象的均值(中心细胞对象)，计算样本集中每个细胞对象与这些中心细胞对象的加权欧式距离。距离的权值采用变异系数赋权法得到；其中，一组一维的细胞特征向量数据集合的变异系数是它的标准差除以均值的绝对值，即对数据集合中的n个数据x₁，x₂，……，x_n，记为均值为

标准差为

变异系数为

细胞的各维特征向量记为z₁，z₂，……，z₄，每个维度的特征向量都有各自的变异系数，用v_i表示z_i的变异系数，特征向量z_i相应的权重系数由该公式计算得到：

两个向量x_i，x_j之间的加权欧式距离表示为：

d(x_i，x_j)＝(w₁(x_i1-x_j1)²+w₂(x_i2-x_j2)²+K+w₄(x_i4-x_j4)²)^1/2；

步骤S13，根据具有最小的加权欧式距离的细胞相似度最高的原则，重新为相应的细胞划分类别；

步骤S14，重新计算每个有变化聚类的中心；

步骤S15，若聚类中心稳定，则执行步骤S16，否则重复执行步骤S12到S14的迭代循环；

步骤S16，输出聚类分析的结果W1，W2，W3和W4。

在步骤S2，进一步剔除细胞碎片等引起的干扰信息，得到参与分类计算的有效细胞数据集。虽然前述信号预处理器32通过“设门”的方式滤除了一些噪声带来的干扰，但是为了提高仪器检测的灵敏度，门限值一般不会设定得太高，因此，一些如破损细胞碎片也会产生有效的事件。为保证分类计算的准确性和速度，在分类计算前必须将其滤除。提取出步骤S16输出结果中每个聚类质心初始向量的属性中的P0(前向透射光强度)、P10(3°至10°小角度散射光强度)，P90(大角度散射光强度)三个属性，并按照属性P0进行由小到大的排序，得到排序后的聚类结果分别记为W₁，W₂，W₃，W₄，如图4所示。其中W₁是由破损细胞碎片引起的干扰信息，该信息通过第一门限判别单元G1之后被删除。该区域如图4中的标记为Noise Area的区域。

在步骤S3，根据步骤S2得到的聚类结果计算最佳划分投影变换的参数，其实施流程如图5所示，包括以下S31～S35步骤：

步骤S31，根据W₂和W₄的聚类质心Q₂(P0₂，P10₂，P90₂)与Q₄(P0₄，P10₄，P90₄)，拟合一条在由P0、P10和P90组成的三维欧式空间中通过点Q₂和Q₄的直线L1；

步骤S32，将直线L1投影到由P0-P10组成的平面直角坐标系之中(由于横、纵坐标的指定与后续计算无关，且三维坐标在观察视图中可以处于不断旋转之中，因此横、纵坐标随着视角的变换是不停变化的，所以可以不进行横、纵坐标的指定，如果要指定的话，随便指定一个即可，如纵坐标P10，横坐标P0，下同)，得到投影线L2，计算投影线L2与坐标轴P0的夹角θ；步骤S33，通过夹角θ计算得到P0与P10坐标系中的旋转矩阵

步骤S34，将直线L1投影到由P10与P90组成的平面直角坐标系之中，得到投影线L3，计算投影线L3与坐标轴P90的夹角β；

步骤S35，通过夹角β计算得到P10与P90坐标系中的旋转向量Vr＝(sinβ·cosβ，cosβ·cosβ)^T。

在步骤S4中，通过使用步骤S3得到的变换参数对细胞数据进行变换投影分析，将细胞进行五分类。注意到淋巴细胞体积在所有白细胞中最小，虽然也有大淋巴细胞，但其数量相当少，体积在6至10um细胞数居多数。由于淋巴细胞的单核结构，根据MIE散射理论，可通过透射光强度和小角度散射光的强度直接从白细胞中区分出来。可定性的认为淋巴细胞形成的透射光和小角度散射光强度在所有细胞中最小。单核细胞是白细胞中体积最大的细胞，其正向的透射光强信号最大，由于其单核结构，在小角度散射光信号强度较若，根据正向透射和小角度散射信号，可定性的认为单核细胞透射光最强，小角度散射光较弱。粒细胞包括中性、嗜酸、嗜碱三类细胞，由于三类细胞的结构类似，因此分类较淋巴和单核细胞困难。嗜酸性粒细胞内含有充满粗大(直径约0.5μm)、均匀的嗜酸性颗粒，这些颗粒可以将有强偏振性质的激光散射为无偏振性质的各向均匀的无偏振光，通过在大角度光消偏振光强度就将其从粒细胞中分辨出来。中性粒细胞核内的颗粒数量较多且细小，嗜碱性粒细胞核内颗粒粗大且较中性粒细胞少，通过3个光通道的信号强度综合判断区分。该步骤的实施流程图如图6所示，包含如下S41～S415步骤：

步骤S41，在P10、P90平面上将细胞的初始向量使用旋转向量Vr进行变换，计算过程如下：设第i个细胞P10、P90属性特征分别为(P10_i、P90_i)，变换之后的特征向量P90_i‘＝(P10_i·sinβ·cosβ+P90_i·cosβ·cosβ)

步骤S42，将经过步骤S41变化过后的特征向值P90₀’，P90₁’，…，P90_i’，…，投影到坐标轴P90上，计算其投影的统计直方图H1并进行指数滤波平滑；

步骤S43，在统计直方图H1上进行波峰、波谷的识别，得到波峰PK1、PK2与波谷VL1；

步骤S44，波峰PK1、PK2与波谷VL1送第二门限判别单元G2进行评价判别，若该组峰谷值合理，则转入步骤S45，否则转入S43重新识别；

步骤S45，以波谷值VL1作为阈值，设置第一分类计算单元CL1，P90属性的特征向量值大于等于该阈值的细胞，划分为多核类细胞；小于该阈值的细胞，划分为单核类细胞。

步骤S46，在P0、P10平面上将被划分为单核类的细胞的初始向量使用旋转矩阵Mr1进行旋转变换，计算过程如下：设第i个单核类细胞P10、P0属性特征分别为(P10_i、P0_i)，变换之后得到特征向量 $({{P10}_i}^{\′},{{P0}_i}^{\′})=\left[\begin{array}{cc}{P10}_i& {P0}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right].$

步骤S47，将经过步骤D16旋转变化过后的特征向量(P10’₁、P0’₁)，(P10’₂、P0’₂)，…，(P10’_i、P0’_i)，…，投影到坐标轴P0上，并计算其投影的统计直方图H2，并进行指数滤波平滑；

步骤S48，在统计直方图H1上进行波峰、波谷的识别，得到波峰PK3、PK4、PK5与波谷VL2、VL3；

步骤S49，波峰PK3、PK4、PK5与波谷VL2、VL3送第三门限判别单元G3进行评价判别，若该组峰谷值合理，则转入步骤S410，否则转入S48重新识别；

步骤S410，以波谷值VL2、VL3作为阈值，设置第二分类计算单元CL2，P0属性的特征向量值大于等于VL2且小于VL3的细胞，划分为淋巴细胞；大于等于VL3的细胞，划分为单核细胞；

步骤S411，提取步骤S45得到的多核细胞的P90属性和P90p属性(大角度偏正光强度属性)，形成分别以P90和P90p为纵、横坐标轴的散点图；

步骤S412，对步骤S411生成的散点图做变换、统计计算之后，得到散点数据的线性划分线L90p，使用L90p的参数设置第三分类计算单元CL3；

步骤S413，将多核细胞对应的特征数据送入CL3进行分类计算，位于L90p之上的细胞划分为嗜酸性粒细胞；其余多核细胞划分为中性粒细胞；

步骤S414，提取已分出的单核细胞聚类与中性粒细胞聚类中细胞初始向量中的P0、P10和P90属性，分别重新计算得到的两类细胞集合的质心Q_D(P0_D，P10_D，P90_D)和Q_Z(P0_Z，P10_Z，P90_Z)，以及两类细胞集合的标准差SD1和SD2；

步骤S415，将Q_D、Q_Z、SD1、SD2作为参数设置第四分类计算单元CL4，并将单核细胞集合与中性粒细胞集合送入第四分类计算单元CL4进行分类计算，识别出集合中的嗜碱性粒细胞。

其中，步骤S42和S47中，为了提高波峰、波谷查找的准确程度，分别对直方图进行了指数平滑滤波，所述指数平滑滤波以S42的为例进行说明，包含如下S421～S426的步骤：

步骤S421，设Pi为待滤波直方图序列的第i项值，E′_i是一次指数滤波结果序列的第i项，另E′_i＝P₀。同时，令平滑系数a＝0.1；

步骤S422，由近似的指数平滑公式E′_i＝aP_i+(1-a)E′_i-1迭代计算得到平滑后的直方图的每一项；

步骤S423，使用均方差MSE评价该平滑指数的滤波效果，算式如下：

记录下该MSE_a及其对应的a值；

步骤S424，若a＜0.9，则a＝a+0.1，重复执行S422到S424步骤的迭代计算；否则从已计算得到的MSE_0.1，MSE_0.2……MSE_0.9集合里选择其中的最小值MSE_a，对应的a值就是最佳一次平滑滤波系数；

步骤S425，使用步骤S424得到的a值和步骤S422所示的计算公式，对直方图进行一次指数平滑滤波；

步骤S426，若0.5≤a＜0.9，则需要对直方图进行二次指数平滑滤波，算式由下式给出：

迭代初值由

给出；

步骤S47所述的滤波流程与上述过程一致。

在步骤S412之中，包含以下S4121～S41210步骤：

步骤S4121，计算属于多核细胞的细胞特征P90p与P90的均值，得到

其中m为多核细胞的数目；令SD90’＝0；

步骤S4122，在由P90与P90p特征分量为坐标轴构成的平面坐标系中计算得到拟合线L90a的参数

拟合线L90a即为：P90p＝A+B*P90；

步骤S4123，计算L90a与坐标轴P90的夹角γ，得到P90-P90p坐标系中的旋转矩阵

步骤S4124，在P90、P90p平面上将被划分为多核类的细胞的初始向量使用旋转矩阵Mr2进行旋转变换，计算过程如下：设第i个多核类细胞P90、P90p属性特征分别为(P90_i、P90p_i)，变换之后得到的特征向量为 $({{P90}_i}^{\′},{{P90p}_i}^{\′})=\left[\begin{array}{cc}{P90}_i& {P90p}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\γ}& \sin \mathrm{\γ}\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right].$

步骤S4125，将经过步骤S4124旋转变化过后的特征向量(P90’₁、P90p’₁)，(P90’₂、P90p’₂)，…，(P90’_i、P90p’_i)，…，投影到坐标轴P90p上，并计算其投影的统计直方图H3；

步骤S4126，在统计直方图H3上进行波峰、波谷的识别，得到波峰PK6与波谷VL4；

步骤S4127，波峰PK6波谷VL4送第四门限判别单元G4进行评价判别，若该组峰谷值合理，则可以通过如下算式计算得到SD90p＝VL4-PK6；否则SD90p＝|SD/R90|。其中，

m为多核类细胞的数目；

步骤S4128，以A、B、SD90为参数，生成划分线L90p’：P90p＝A*P90+B+SD90，送第五判别单元CL5重新划分待分类的多核细胞集合；

步骤S4129，若SD90’等于SD90，则转入步骤S41210；否则，令SD90’＝SD90，将经过步骤S4128重新划分的多核细胞集合送入步骤S4121，开始新一轮的迭代计算过程；

步骤S41210，以A、B、SD90为参数，生成划分线L90p：P90p＝A*P90+B+SD90，送第三判别单元CL3。

其中，在步骤S4128所述的第五分类计算单元CL5包含以下功能：使用线性划分线L90p’将待处理的多核细胞集合划分为预分类嗜酸细胞集合和待分类多核细胞集合。经过步骤S412和第五分类计算单元CL5处理过后的嗜酸、中性粒细胞的分类图示如图7所示，其中位于分类线之上的是嗜酸性粒细胞。

本发明所述方法中，G2、G3和G4均以人血的极限值作为设定值，G1只是用于把属于W1类的粒子信息删除。

以下是一个用于血液细胞五分类分析仪的实现例：

在基于流式细胞技术的血液细胞五分类仪中利用前向透射光检测通道、小角度散射光检测通道、大角度散射光检测通道和大角度消偏振光检测通道将全血中的白细胞分为五个子类，分别是淋巴细胞(Ly)、单核细胞(Mo)、中性粒细胞(Neu)、嗜酸性粒细胞(Eos，简称嗜酸细胞)和嗜碱性粒细胞(Baso，简称嗜碱细胞)，图7、图8和图9是对同一样本进行分类的结果，其中图7和图9采用了本发明所述方法进行分类，图9的散点图是经过了旋转变化后的重新投影。图8是未使用本发明所述方法对同一样本的一种分类结果。从图中可以看出，经过本发明所述方法进行分类，分类界限变得更为优化，分类结果更为合理。

Claims (7)

1.基于激光鞘流技术的血液细胞五分类方法，其特征在于包括以下步骤：

A、收集每个细胞或粒子通过鞘流器的激光照射区域时产生的至少4路光信号，将每个细胞或粒子的表征为一个与其各路光信号强度相关的四维的初始向量，分别为前向透射光强度P0、3°至10°小角度散射光强度P10、90°大角度散射光强度P90和大角度消偏振光强度P90p四个属性；

B、以加权欧式距离的K均值聚类法进行四类细胞集合的初次聚类，其中距离的权值采用变异系数赋权法得到；

C、分别提取出上述各聚类质心初始向量中的P0、P10和P90三个属性，将每个聚类质心按照P0的属性值由小到大进行排序，得到排序后的四个聚类分别编号为W₁、W₂、W₃和W₄，其中W₁是由破损细胞碎片引起的干扰信息，该信息通过第一门限判别单元G1之后被删除；

D、将W₂和W₄两个聚类的聚类质心经过空间坐标变换得到用于产生分类计算使用的特征向量的旋转矩阵和旋转向量；

E、用旋转向量或旋转矩阵将各细胞的初始向量转换成特征向量，计算并提取相应空间坐标内细胞特征向量的统计特征，通过决策树对细胞进行五分类。

2.根据权利要求1所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：所述步骤D包括以下步骤：

D1、提取W₂和W₄聚类的质心Q₂(P0₂，P10₂，P90₂)与Q₄(P0₄，P10₄，P90₄)，拟合一条在由P0、P10和P90组成的三维欧式空间中通过点P₂和P₄的直线L1；

D2、将直线L1投影到由P0和P10组成的平面直角坐标系中，得到投影线L2，计算投影线L2与坐标轴P0的夹角θ；

D3、通过夹角θ计算得到P0和P10坐标系中的特征向量变换所需要的旋转矩阵

D4、将直线L1投影到由P10与P90组成的平面直角坐标系中，得到投影线L3，计算投影线L3与坐标轴P90的夹角β；

D5、通过夹角β计算得到P10与P90坐标系中的特征向量变换所需要的旋转向量Vr＝(sinβ·cosβ，cosβ·cosβ)^T。

3.根据权利要求1所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：所述步骤E包括以下步骤：

E1、在P10、P90组成的坐标平面上将各细胞的初始向量用旋转向量转换成特征向量，计算过程为：设第i个细胞的初始向量中P10、P90属性特征分别为(P10_i、P90_i)，转换得到的特征向量P90_i‘＝(P10_i·sinβ·cosβ+P90_i·cosβ·cosβ)；

E2、将转换后的特征向量P90₀’，P90₁’，…，P90_i’，…，投影到坐标轴P90上，计算其投影的统计直方图H1；

E3、在统计直方图H1上查找波峰、波谷，得到波峰PK1、PK2和波谷VL1；

E4、波峰PK1、PK2和波谷VL1送第二门限判别单元G2进行评价判别，若该组峰谷值合理，则转入步骤E5，否则转入步骤E3重新识别；

E5、以波谷值VL1作为阈值，设置第一分类计算单元CL1，P90属性的特征向量值大于等于该阈值的细胞，划分为多核类细胞；小于该阈值的细胞，划分为单核类细胞；

E6、在P10、P0平面上将被划分为单核类的细胞的初始向量用旋转矩阵转换成特征向量，计算过程为：设第i个单核类细胞的初始向量的P10、P0属性特征分别为(P10_i、P0_i)，转换得到的特征向量 $({{P10}_i}^{\′},{{P0}_i}^{\′})=\left[\begin{array}{cc}{P10}_i& {P0}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right];$

E7、转化后的特征向量(P10’₁、P0’₁)，(P10’₂、P0’₂)，…，(P10’_i、P0’_i)，…，投影到坐标轴P0上，并计算其投影的统计直方图H2；

E8、在统计直方图H1上查找波峰、波谷，得到波峰PK3、PK4、PK5和波谷VL2、VL3；

E9、波峰PK3、PK4、PK5和波谷VL2、VL3送第三门限判别单元G3进行评价判别，若该组峰谷值合理，则转入步骤E10，否则转入步骤E8重新识别；

E10、以波谷值VL2、VL3作为阈值，设置第二分类计算单元CL2，将P0属性的特征向量值大于等于VL2且小于VL3的细胞，划分为淋巴细胞；而P0属性的特征向量值大于等于VL3的细胞，划分为单核细胞；

E11、提取步骤E5得到的多核类细胞中各细胞的初始向量中的P90属性和P90p属性，形成分别以P90和P90p为纵、横坐标轴的散点图；

E12、对生成的散点图做变换、统计计算，得到散点数据的线性划分线L90p，用线性划分线L90p的参数设置第三分类计算单元CL3；

E13、将步骤E5得到的多核类细胞中各细胞初始向量中的各属性送入第三分类计算单元CL3进行分类计算，位于线性划分线L90p之上的细胞划分为嗜酸细胞；其余多核细胞划分为中性粒细胞；

E14、提取已分出的单核细胞聚类与中性粒细胞聚类中细胞初始向量中的P0、P10和P90属性，分别重新计算得到的两类细胞集合的质心Q_D(P0_D，P10_D，P90_D)和Q_Z(P0_Z，P10_Z，P90_Z)，以及两类细胞集合的标准差SD1和SD2；

E15、将Q_D、Q_Z、SD1、SD2作为参数设置第四分类计算单元CL4，并将单核细胞集合与中性粒细胞集合送入第四分类计算单元CL4进行分类计算，识别出集合中的嗜碱性粒细胞。

4.根据权利要求3所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：步骤E12中包括以下步骤：

F1、计算多核细胞聚类中各细胞初始向量中P90p和P90属性的均值，得到

其中m为多核细胞的数目；令SD90’＝0；

F2、在由P90和P90p属性为坐标轴构成的平面坐标系中计算得到拟合线L90a的参数

拟合线L90a即为：P90p＝A+B*P90；

F3、计算L90a与坐标轴P90的夹角γ，得到P90-P90p坐标系中的旋转矩阵 $\mathrm{Mr}2=\left[\begin{array}{cc}\cos \mathrm{\γ}& \sin \mathrm{\γ}\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right];$

F4、在P90、P90p平面上将被划分为多核类的细胞的初始向量用旋转矩阵Mr2进行旋转变换，计算过程为：设第i个多核类细胞初始向量的P90、P90p属性特征分别为(P90_i、P90p_i)，变换后得到的特征向量为 $({{P90}_i}^{\′},{{P90p}_i}^{\′})=\left[\begin{array}{cc}{P90}_i& {P90p}_i\end{array}\right]\·\left[\begin{array}{cc}\cos \mathrm{\γ}& \sin \mathrm{\γ}\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right];$

F5、将变换后的各细胞的特征向量(P90’₁、P90p’₁)，(P90’₂、P90p’₂)，…，(P90’_i、P90p’_i)，…，投影到坐标轴P90p上，并计算其投影的统计直方图H3；

F6、在统计直方图H3上查找波峰、波谷，得到波峰PK6和波谷VL4；

F7、波峰PK6波谷VL4送第四门限判别单元G4进行评价判别，若该组峰谷值合理，则可以通过以下算式计算得到SD90p＝VL4-PK6；否则SD90p＝|SD/R90|；其中，

m为多核类细胞的数目；

F8、以A、B、SD90为参数，生成划分线L90p’：P90p＝A*P90+B+SD90，送第五分类计算单元CL5重新划分待分类的多核细胞集合；

F9、若SD90’＝SD90，则转入步骤F10；否则，令SD90’＝SD90，将经过步骤F8重新划分的多核细胞集合送入步骤F1，开始新的一轮迭代计算过程；

F10、以A、B、SD90为参数，生成划分线L90p：P90p＝A*P90+B+SD90，送入第三分类计算单元CL3。

5.根据权利要求3所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：步骤E7和E8之间还包括对直方图进行滤波的步骤。

6.根据权利要求1所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：所述步骤D包括以下步骤：

B1、以实验获得的应用经验值作为四个分类细胞集合的中心；

B2、计算样本集中每个细胞对象与这些中心细胞对象的加权欧式距离，并根据最小的加权欧式距离重新对相应的细胞对象进行划分；

B3、对于中心细胞对象有变化的分类，重新计算它们的中心细胞对象；

B4、重复执行B2到B3的迭代循环，直到每个聚类的中心细胞对象稳定不再发生变化为止。

7.根据权利要求1～6中任何一项所述的基于激光鞘流技术的血液细胞五分类方法，其特征在于：步骤A和B之间，还包括通过设定阈值，将不符合条件的噪声或粒子的数据删除的步骤。

Images