1.本发明涉及计算机视觉图像处理，染色体计数等技术领域，具体涉及一种染色体极性识别方法和系统。
背景技术：
2.染色体核型分析是发现染色体病的重要手段，染色体数目或者结构异常都能通过染色体核型分析的技术手段发现。为了便于临床医生根据染色体的形态结构给出相应诊断结果，分割出来的所有染色体会按照顺序排列，同时保证呈竖直状态且染色体短臂朝上，长臂朝下，即调整染色体极性，最终形成准确清晰的染色体核型图。
3.目前，染色体分析系统一般严重依赖人工调整染色体极性，医生需要通过点击或拖拽的方式调整，耗时耗力。目前已存在的技术方法往往会同时完成染色体类型识别和染色体极性识别。这些方法同时利用计算机图形学方法和深度卷积神经网络模型提取与染色体类型识别和极性识别相关的特征，然后利用两个简单的分类器来分别判断染色体的类别和极性。这些方法的优点在于可以同时完成染色体类型识别和极性识别，节省了运算时间和运算所需资源；而缺点在于，染色体类型识别应与染色体极性无关，即判断染色体类别不应受染色体极性的影响。此外，目前的方法都没有阐明如何得到处于竖直状态的染色体，这是后续染色体进行极性翻转的前提步骤。
4.目前人工调整染色体极性是广泛使用的方法，当前的染色体分析系统提供点击拖拽功能，方便临床医生对观察发现的染色体极性异常进行调整，即上下翻转极性异常的染色体，保证短臂朝上，长臂朝下的状态，需要耗费人力，效率较低，尤其有的染色体长短臂不易识别，给诊断造成困难。
5.基于深度学习的染色体识别方法充分利用深度学习技术强大的特征提取能力，同时也补充计算机图形学提取的特征，最后利用两个独立的分类器分别完成染色体类型识别和染色极性识别，但主要目的还是为了完成染色体类型识别，因此首先需要保证染色体类型识别与染色体极性无关，即无论染色体短臂是否朝上都需要模型准确的判断染色体类别，因此两个任务之前存在一定的冲突，所学到的特征也无法使两个任务同时达到最优效果。此外，该方法默认输入的染色体处于竖直状态，没有考虑从细胞分裂中期图像中分割出来的单条染色体并不必然处于竖直状态，是一个比较重大的技术缺失。
6.人工调整染色体极性效率较低，而基于深度学习的染色体识别方法的缺陷则是由于染色体类型识别和染色体极性识别所需的特征存在冲突而且难以阐明如何从分割出的、处于任意角度状态下的染色体得到处于竖直状态的染色体，以便染色体极性识别模型判断和调整染色体极性。本发明的目的旨在解决染色体极性识别过程中人工效率低下以及基于深度学习的染色体识别算法任务冲突和缺少染色体旋转关键步骤技术方案的问题，提出了一种染色体极性识别方法及系统。
技术实现要素：
7.基于上述目的，本发明拟通过计算染色体在水平面的投影距离和专注于染色体极性识别任务的深度卷积神经网络来解决现有技术存在的问题。本发明首先提供了一种基于深度学习的染色体极性识别方法，所述方法包括以下步骤：（1）收集数据集：收集从细胞分裂中期图像中通过染色体分割方法分割出来的单条染色体，旋转收集到的染色体呈竖直状态并提取染色体的极性特征，以短臂朝上或短臂朝下进行极性标注；（2）构建训练集和测试集：将步骤（1）获得的数据集划分为训练集和测试集，构建以分类网络结构为主干网络的染色体极性识别模型；（3）基于训练集对染色体极性识别模型进行学习训练，其中，softmax函数对输出的2维向量归一化，然后使用交叉熵损失函数度量预测结果与真实结果之间的距离，然后通过随机梯度下降算法对网络参数进行学习；（4）将测试集输入到染色体极性识别模型进行测试，输出待预测染色体的极性结果。
8.本发明需要收集单条染色体，这些染色体是从细胞分裂中期图像中通过染色体分割方法分割出来的单条染色体，其中轴仍会保持其在分裂中期图像中的角度。因此，为了使每条染色体呈竖直状态，本发明将沿逆时针旋转染色体并实时计算染色体在水平方向的投影距离，直至投影距离最短为止。
9.在一个优选的实施方案中，步骤（1）所述的旋转收集到的染色体呈竖直状态所需的角度计算为式（i）所示：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（i）其中，投影距离指两个投影点之间的水平距离，当前旋转角度下染色体区域最左端点在水平面上的投影点为，染色体区域最右端点在水平面上的投影点为 ，将所述染色体图像根据旋转角度θ进行逆时针旋转至竖直状态。
10.对经旋转后的呈竖直状态的染色体，进行极性标注，其中，若短臂朝下标注为“down”，短臂朝上标注为“up”。本发明充分利用卷积神经网络优异的特征提取能力，利用包括但不限于vggnet、resnet、densenet等分类模型提取与染色体极性相关特征，并最终进行极性二分类，得到染色体极性识别模型。
11.在一个优选的实施方案中，步骤（2）所述分类网络结构为resnet101。
12.更为优选地，以7:3划分训练集和测试集，并使用式（ii）对训练集和测试集进行图像归一化：（ii）其中i
o
代表原始图片且已规范至0.0到1.0范围，m代表训练集数据在各个通道上的均值，δ代表训练集数据在各个通道上的方差，i
n
代表图像归一化后的图像数据。
[0013] 在一个优选的实施方案中，步骤（3）所述交叉熵损失函数为式（iv）所示：
ꢀꢀꢀꢀ
（iv）其中，l(z,y)是指样本的预测类别与其真实类别的交叉熵损失，z表示最后一个全连接层输出的2维向量， y代表该样本的真实类别，为了简单起见，记“down”类为0，“up”类为1，则,n
cls
代表类别数目，i代表类别索引，z
y
代表输出向量在真实类别位置处的值，z
i
代表第i类位置处的值，e为自然底数。
[0014]
其次，本发明还提供了一种基于深度学习的染色体极性识别系统，所述系统包括以下模块：（1）数据收集模块，所述模块用于收集从细胞分裂中期图像中通过染色体分割方法分割出来的单条染色体，旋转收集到的染色体呈竖直状态并提取染色体的极性特征，以短臂朝上或短臂朝下进行极性标注；（2）训练集和测试集构建模块，所述模块用于将数据收集模块提供的数据集划分为训练集和测试集，并对数据集进行预处理和数据扩增；（3）训练模块，所述模块用于将经过训练集和测试集构建模块提供的训练集中的染色体图像输入到以分类网络结构为主干网络的染色体极性识别模型，并对染色体极性识别模型进行学习训练，其中，softmax函数对输出的2维向量归一化，然后使用交叉熵损失函数度量预测结果与真实结果之间的距离，然后通过随机梯度下降算法对网络参数进行学习；（4）测试模块：所述测试模块将测试集输入到染色体极性识别模型进行测试，输出待预测染色体的极性结果。
[0015]
在一个优选的实施方案中，数据收集模块所述的旋转收集到的染色体呈竖直状态所需的角度计算为式（i）所示：
ꢀꢀꢀ
（i）其中，投影距离是指两个投影点之间的水平距离，当前旋转角度下染色体区域最左端点在水平面上的投影点为，染色体区域最右端点在水平面上的投影点为，将所述染色体图像根据旋转角度θ进行逆时针旋转至竖直状态。
[0016]
在另一个优选的实施方案中，训练集和测试集构建模块所述分类网络结构为resnet101。
[0017]
更为优选地，以7:3划分训练集和测试集，并使用式（ii）对训练集和测试集进行图像归一化：
ꢀꢀꢀ
（ii）其中，i
o
代表原始图片且已规范至0.0到1.0范围，m代表训练集数据在各个通道上的均值，δ代表训练集数据在各个通道上的方差，i
n
代表图像归一化后的图像数据。
[0018]
在一个优选的实施方案中，训练模块所述交叉熵损失函数为式（iv）所示：
ꢀꢀꢀ
（iv）
其中，l(z,y)是指样本的预测类别与其真实类别的交叉熵损失， z表示最后一个全连接层输出的2维向量， y代表该样本的真实类别，为了简单起见，记“down”类为0，“up”类为1，则,n
cls
代表类别数目，i代表类别索引。z
y
代表输出向量在真实类别位置处的值，z
i
代表第i类位置处的值，e为自然底数。
[0019]
本发明首次提出了旋转分割后的染色体识别方法，染色体极性识别模型判断染色体极性，最终通过极性结果完成染色体极性翻转，得到竖直的、短臂朝上长臂朝下的染色体。本发明技术效果如下：1. 数据来源简单，因为本发明包含了基于投影距离的染色体旋转方法，因此采集的数据是从细胞分裂中期的图像分割得到的单条染色体，且无需人工调整其中轴线方位即可得到处于竖直状态的染色体，使得染色体分析自动化程度更高，流程更加简洁，可广泛推广应用。
[0020]
2. 模型设计基于深度学习分类算法，能够准确地判断当前染色体极性类别，并据此完成染色体极性调整，使得染色体均保持短臂朝上的状态。发明人用4490例标注的独立染色体图训练模型，其中训练集与测试集比例为7:3。统计显示，“up”类染色体图有2289例，“down”类染色体图有2201例，分为训练集3143例和测试集1347例，最终测试结果：真阳性（tp，true positive）=686；假阳性（fp ，false positive）=17；真阴性（tn，true negative）=612；假阴性（fn，false negative）=32；准确率（accuracy）=96.36%。
附图说明
[0021]
图1基于深度学习的染色体极性识别方法流程示意图；图2 分割后单条染色体示意图；图3 基于投影距离的染色体旋转方法原理示意图;图4 resnet101网络结构示意图。
具体实施方式
[0022]
下面结合具体实施例来进一步描述本发明，本发明的优点和特点将会随着描述而更为清楚。但这些实施例仅是范例性的，并不对本发明的权利要求所限定的保护范围构成任何限制。
[0023]
本发明是一种基于深度学习的染色体极性识别方法，主要目的是为了解决分割后的染色体旋转和极性识别问题。本发明基于深度学习resnet模型，发明人对收集到的分割后的染色体进行基于投影距离的染色体旋转，得到处于竖直状态的染色体。在对这些竖直状态的染色体进行标注后，构建训练集和测试集。模型训练完成后可以直接预测给出染色体的极性结果，据此完成染色体极性翻转。
[0024]
下面将结合附图和实施例对本发明作进一步的详细说明。
实施例
[0025]
本发明步骤流程如图1所示，主要包含以下四个步骤：1. 收集数据集：收集从细胞分裂中期图像中通过染色体分割方法分割出来的单条染色体，旋转收集到的染色体呈竖直状态并提取染色体的极性特征，以短臂朝上或短臂
朝下进行极性标注；1.1收集染色体采集徕卡cytovision 自动细胞遗传学平台记录的高分辨率显微镜视野下的处于有丝分裂中期细胞染色体真实图像样本共计4490例。收集分割得到的一例单条染色体图像如图2所示，a为分割前图像，b为分割后图像。分割后的染色体图像为一矩形图像，该矩形边缘紧密包围该染色体，且其中轴线所处角度与原图所处的角度一致。
[0026]
1.2旋转染色体（本操作可使用python语言开发并受到opencv开源计算机视觉库和numpy开源科学计算库支持），构建基于投影距离的染色体旋转方法，根据最小投影距离旋转染色体。
[0027]
本发明发现，从细胞分裂中期图像中分割出来的单条染色体，其中轴线仍会保持其在分裂中期图像中的角度，因此第一步需要找到合适的旋转角度使得旋转后的染色体处于竖直状态。然后使用大量经过人工标注的染色体极性数据训练染色体极性识别模型，该模型为图像二分类模型，仅用来判断染色体短臂朝上或朝下。最后，根据染色体极性识别模型的判断结果，决定是否需要上下翻转染色体。
[0028]
本步骤中所述旋转染色体呈竖直状态所需的角度计算为式（i）所示：（i）其中，投影距离是指两个投影点之间的水平距离，当前旋转角度下染色体区域最左端点在水平面上的投影点为，染色体区域最右端点在水平面上的投影点为，将所述染色体图像根据旋转角度θ进行逆时针旋转至竖直状态。
[0029]
本方法通过度量染色体区域最左端点到最右端点在水平上的投影距离来度量染色体是否处于竖直状态。对于任意分割得到的染色体，以图像几何中心为原点，以角度1度为步长，以角度180度为上限，逆时针旋转图像并同步计算染色体区域最左端点和最右端点在水平方向上的投影距离，当该投影距离达到最小时，记录旋转的角度，并将染色体按照该角度旋转，得到处于竖直状态的染色体。具体地，如图3所示为基于投影距离的染色体旋转方法示意图，将染色体旋转的角度记作，当前旋转角度下染色体区域最左端点记作，其在水平面上的投影为，染色体区域最右端点记作 ，其在水平面上的投影为 ，因此，投影距离是两个投影点之间的水平距离。通过最小化投影距离找到旋转角度，作为该染色体旋转至竖直状态所需角度。最后将该染色体图像根据旋转角度θ进行逆时针旋转，并在空白区域填充（255,255,255），最终得到白色背景的呈竖直状态的染色体。
[0030]
1.3标注染色体将旋转后的染色体进行数据标注，其中短臂朝下标记为“down”类型，共有2201例，短臂朝上标记为“up”类型，共有2289例。
[0031]
2．划分训练集和测试集
2.1划分数据集本发明按照比例7:3划分训练集和测试集，得到训练集3143例和测试集1347例。
[0032]
2.2进行数据预处理和数据扩增（本操作可使用使用python语言编写，使用pytorch开源机器学习库）。训练集和测试集图像均缩放至224
×
224，然后使用式（ii）对训练集和测试集进行图像归一化： （ii）其中，i
o
代表原始图片且已规范至0.0到1.0范围，m代表训练集数据在各个通道上的均值，本实施例中各通道均值分别为[0.772, 0.772, 0.772]，δ代表训练集数据在各个通道上的方差，本实施例中各通道方差分别为[0.245, 0.245, 0.245]， i
n
代表归一化后的图像数据。归一化后的数据集使得输入像素分布一致，归一化后，数据集像素值符合均值为0，标准差为1的高斯分布，模型训练更容易收敛。训练模型前，还对训练集数据进行数据扩增操作，本实施例中选择以0.5的概率对图像进行原地水平翻转，因此不会影响图像真实类别（每轮迭代时以0.5的概率原地水平翻转，不会影响最终数据量）。
[0033]
3. 将经过步骤（2）获得的训练集中的染色体图像输入到以分类网络结构为主干网络的染色体极性识别模型，并对染色体极性识别模型进行学习训练，其中，使用softmax函数对输出的2维向量归一化，然后使用交叉熵损失函数度量预测结果与真实结果之间的距离，然后通过随机梯度下降算法对网络参数进行学习3.1构建基于分类网络的染色体极性识别模型（本操作可使用python语言编写，并受pytorch开源机器学习库支持）。
[0034]
染色体极性识别模型可以采用经典的分类网络结构如vggnet、resnet和densenet等作为主要结构，本实施例以resnet101（resnet系列网络中深度为101层的神经网络）作为主干网络（本实施例中resnet101结构示意图见图4）。其中，resnet作为残差网络的重要代表用于本发明的分类任务，有效地解决了由于神经网络不断加深而带来的网络退化问题，可以训练更深的网络，增强网络的表达能力。残差网络是指网络的每层的输出h由原始输入x与卷积和非线性激活函数映射输出f(x)的和组成，即，h=f(x)+x；除resnet以外还有诸如resnext、densenet等残差网络，都可以应用于本发明技术方案。
[0035]
本实施例中，选择101层的resnet网络作为基本网络，该101层结构中包含100层卷积层和一层全连接层，并对最后一层全连接层的输出维度进行修改（原始的resnet101最后一层全连接层权重维度为2048
×
1000，为了适应本任务，将其修改为2048
×
2）。该残差网络基于bottleneck模块构造，每个bottleneck模块由三个卷积层构成，分别是两个“1
×
1”的卷积层和一个“3
×
3”的卷积层，其中两个“1
×
1”的卷积核负责减少和增加通道维度，通过bottleneck模块可以减少计算量，因而减少训练时间，降低计算和存储消耗。为了形成残差结构，与三个卷积层并行的还有一个旁路连接（shortcut），若前后两个bottleneck模块通道维度一致，则采用恒等映射将输入逐元素与三个卷积层输出相加，若前后两个bottleneck模块通道维度不一致，则额外采用一个“1
×
1”的卷积层增加输入通道维度再与三个卷积层的输出逐元素相加。本实施例主干网络采用的resnet101网络按照输出特征图尺度可以分为5个不同阶段，分别为“conv1”、“conv2”、“conv3”、“conv4”和“conv5”。“conv1”阶段由一个7
×
7的卷积核组成，移动步长为2，该阶段输出特征图，其
中c代表特征图通道维度，h和w分别代表特征图高和宽，并且h和w共同组成空间尺度。“conv2”阶段首先由一个3
×
3大小的最大池化层进行步长为2的池化操作，进一步缩小空间尺度，然后由若干个bottleneck模块组成。剩余的“conv3”、“conv4”和“conv5”同样由若干个bottleneck模块组成，其中仅每个阶段第一个bottleneck的第一个1
×
1的卷积核步长为2，以此不断缩小特征空间尺度。具体的，每个阶段bottleneck模块配置各不相同。“conv2”阶段由3个bottleneck模块组成，三个卷积层通道维度为（64,64,256）。“conv3”阶段由4个bottleneck模块组成，三个卷积层通道维度为（128,128,512）。“conv4”阶段由23个bottleneck模块组成，三个卷积层通道维度为（256，256,1024）。“conv5”阶段由3个bottleneck模块组成，三个卷积层通道维度为（512,512,2048）。“conv5”阶段后为全局平均池化层（global average pooling），该层沿着h和w维度计算均值，得到1
×1×
2048的特征向量。染色体极性识别模型的最后一层（见图4）是一个2048
×
2的全连接层（fully connected layer），该全连接层作为极性分类器，能够与主干网络形成一个端到端（end
‑
to
‑
end）的网络，可以进行联合训练，使得主干网络提取的特征更加适合染色体极性识别所需。
[0036]
3.2设置归一化参数、损失函数和超参数（超参数包括学习率和批次大小等在内的参数均为超参数），训练模型。
[0037]
3.2.12维向量归一化（本操作可使用python语言并受到pytorch开源机器学习库支持）染色体极性识别任务可以视为二分类任务（
ꢀ“
up”和“down”），因此首先使用softmax函数对输出的2维向量归一化，softmax又称归一化指数函数，在本实施例中，其能将一个含任意实数的2维向量z压缩到另一个2维向量σ(z)中，使得每一个元素的范围都在(0,1)之间，并且所有元素的和为1。softmax函数如式（iii）所示： （iii）其中， z表示最后一个全连接层输出的2维向量， y代表该样本的真实类别，为了简单起见，记“down”类为0，“up”类为1，则。 ,n
cls
代表类别数目，本实施例中为2，输出2维向量，i代表类别索引。 z
y 代表输出向量在真实类别位置处的值，z
i
代表第 i类位置处的值。对输出的2维向量归一化。2维向量是指的有两个分量的向量，是网络预测的输出，分别代表“up”类和“down”类，e是指自然底数。归一化指的是指softmax函数将输出进行缩放使得向量的两个分量的和为1，归一化的结果符合概率定义，这样两个分量就可以作为网络预测的类别的概率。
[0038]
3.2.2交叉熵损失函数度量预测结果与真实结果之间的距离（本操作可使用python语言实现并受到pytorch开源机器学习库支持）使用分类任务常用的交叉熵损失函数（cross entropy loss）来度量预测结果与真实结果之间的距离，具体地，交叉熵损失函数形式如公式（iv）所示：（iv）
其中， z表示最后一个全连接层输出的2维向量，y代表该样本的真实类别，为了简单起见，记“down”类为0，“up”类为1，则。n
cls
代表类别数目，本实施例中为2， i代表类别索引。z
y 代表输出向量在真实类别位置处的的值， z
i
代表第i类位置处的值。该函数首先对输出的2维向量使用softmax归一化，使得所有类别的输出之和等于1，随后通过交叉熵损失函数度量差距。例如，当某个样本的真实类别为“down”,网络的预测结果为 ,则差距即：交叉熵损失为: 。
[0039]
本实施例中通过交叉熵损失函数度量预测结果和真实结果之间的差距：结果显示训练集开始的损失为0.448，损失是一个动态下降过程，当训练结束时可降低至0.001。
[0040]
3.3随机梯度下降系列算法对网络参数进行学习通过随机梯度下降（stochastic gradient descent）系列算法对网络参数进行学习。网络参数是指resnet101的含参的神经元，主要是卷积层参数：包括卷积核权重参数和偏置参数，全连接层权重参数和偏置参数等，前述参数不需要具体设置，将在imagenet数据集上（krizhevsky a, sutskever i, hinton g e. imagenet classification with deep convolutional neural networks[j]. advances in neural information processing systems, 2012, 25: 1097
‑
1105.）训练得到的resnet101参数作为本实施例中的resnet101的初始化参数，然后通过基于动量的随机梯度下降算法根据损失函数相对参数的梯度更新参数。本操作可使用python语言编写并受到pytorch开源机器学习库支持。
[0041]
本实施例中选择基于动量的随机梯度下降算法 (sutskever, ilya, et al. "on the importance of initialization and momentum in deep learning." international conference on machine learning. 2013.)如（v）所示：
ꢀꢀ
（v）其中， ε表示学习率， f(θ)表示损失函数，表示损失函数相对θ
t
的梯度， v表示更新速率，t表示迭代次数，动量因子momentum(μ)设置为0.9，权重衰减因子设置为 5
×
10
‑4，学习率初始化为0.01，批次大小为32，总共训练100个epoch，其中第60个epoch和第90个epoch学习率下降至原来的1/10。通过随机梯度下降系列算法对网络参数进行学习，最终使得模型在训练集上的平均损失从0.448降低至0.001。
[0042]
以上包括图像归一化函数、resnet101这样的网络结构、softmax函数、交叉熵损失函数、基于动量的随机梯度下降算法在内的算法均可由pytorch开源机器学习库提供支持。
[0043]
4．将测试集输入到染色体极性识别模型进行测试，输出待预测染色体的极性结果。
[0044]
模型训练完毕后，将测试集经过相同的数据预处理，输入到模型中将获得每个样本的输出向量，选择值最大的索引所对应的类别作为对该染色体预测的极性结果。若预测为第0类，即“down”类，则代表染色体极性识别模型预测该染色体短臂朝下，因此通过上下翻转使该染色体短臂朝上。若预测为第1类，即“up”类，则代表染色体极性识别模型预测该染色体短臂朝上，因此无需进行上下翻转操作。
[0045]
本发明的模型设计基于深度学习分类算法，能够自动准确的识别染色体极性。发明人用4490例标注的独立染色体图训练模型，其中训练集与测试集比例为7:3。
[0046]
统计显示，“up”类染色体图有2289例，“down”类染色体图有2201例，分为训练集3143例和测试集1347例，最终测试结果真阳性（tp，true positive）=686；假阳性（fp ，false positive）=17；真阴性（tn，true negative）=612；假阴性（fn，false negative）=32；准确率（accuracy）=96.36%。