Computational model for sperm cells identification and morphological classification by automatic analysis of digital micrographs

Heidy Hernandez, Wolfram Research South America

Join the Conversation #WolframTechConf

Introduction

The analysis of the sperm cells involves concentration, motility analysis and morphological classification

Nowadays, the methods used to evaluate sperm morphology are mostly manual and lack objectivity

Criticism to the computerized tools existing nowadays: high cost, difficulty for calibration, validation of results and parameters standardization

The goal of this project was to develop a computational model to characterize the morphological aspects of the sperm cells head by analyzing digital micrographs and later on the classification of those cell’s head in normal or abnormal using machine learning techniques.

Source images

First source

Repository of digital micrographs made available in the Center of Digital Spermiogramma Assisted by Internet (CEDAI): 19 images with a total of 240 sperm cells identified

fileNames1=NotebookDirectory[]~~"images/source1/"~~ToString[#]~~".bmp"&/@Range[19];
sourceImages1=Import/@fileNames1;
sourceImages1[[1]]

Second source

World Health Organization laboratory manual for the examination and processing of human semen (5th edition): 28 images with a total of 226 sperm cells identified

fileNames2=NotebookDirectory[]~~"images/source2/"~~ToString[#]~~".PNG"&/@Range[28];
sourceImages2=Import/@fileNames2;
sourceImages2[[5]]

For both sources, sperm cells were identified in each image and were previously classified as normal or abnormal for specialists in the field.

Image Binarization

The first source contained a mask associated to each image

binarizedImages1=Import/@(NotebookDirectory[]~~"images/source1/masks/mask_"~~ToString[#]~~".bmp"&/@Range[19]);
binarizedImages2=Import/@(NotebookDirectory[]~~"images/source2/masks/mask_"~~ToString[#]~~".PNG"&/@Range[28]);
{sourceImages1[[5]],Import[NotebookDirectory[]~~"images/source1/masks/mask_5.bmp"]}

For the second source, there was needed a process to obtain a binarized version of the images to highlight the sperm cells head

Adjusting threshold

The FindThreshold function finds a global threshold value that partitions the intensity values in a given image into two intervals, in this case, using Otsu's algorithm.

sourceImages2[[5]]

threshold=FindThreshold[sourceImages2[[5]]]

Later, the function Binarize creates a binary image from by replacing all values above the specified threshold with 1 and others with 0.

Binarize[sourceImages2[[5]],{0,threshold}]

Manipulate[SelectComponents[Binarize[sourceImages2[[5]],{0,threshold+a}],{"AdjacentBorderCount","Area"},#1==0&&50<#2<500&],{a,-0.25,0.25,0.05}]

Example result:

{sourceImages2[[5]],SelectComponents[Binarize[sourceImages2[[5]],{0,0.53}],{"AdjacentBorderCount","Area"},#1==0&&50<#2<500&]}

This procedure allowed to identify 221 of the 226 cells present in the images, that is the 97.8% of them.

Feature Extraction

There is no consensus in the metrics that must be used to characterize the sperm cells head

Some of them are used repeatedly in studies, those are the head’s length, width, area and perimeter

In addition, it was also taken into account the eccentricity of the ellipse that approximates the best the sperm head and the ratio between the width and the length, to represent the idea of a regular oval shape

Component Measurements

The ComponentMeasurements function detect automatically components in a binary 2D image, and compute several characteristics associated to those components.

showComponents[headMask_]:=Block[{rects,numbers},(rects=Rectangle[#[[2,1]],#[[2,2]]]&/@ComponentMeasurements[headMask,"BoundingBox"];numbers=Inset[#[[1]],#[[2]]]&/@ComponentMeasurements[headMask,"Centroid"];
Show[headMask,Graphics[{Red,Opacity[0.5],rects,Black,Opacity[1],numbers}]])]

showComponents[binarizedImages1[[8]]]

measureComponents[image_]:=Block[{components},
(
components=ComponentMeasurements[image,{"Length","Width","Area","PerimeterLength","Eccentricity"}][[All,2]];
components=Append[#,#[[2]]/#[[1]]]&/@components;
Inner[Rule,{"Length","Width","Area","Perimeter","Eccentricity","Ratio W/L"},#,Association]&/@components
)]

measureComponents[binarizedImages1[[8]]]//Dataset

Additional steps

Mapping the identified cells with the specialist classification assigning a label to each of the records

map1=Import[NotebookDirectory[]~~"map1.csv","CSV"];
map1Fmt=Inner[Rule,Range[First[#]],Rest[#][[1;;First[#]]],List]&/@map1;
map2=Import[NotebookDirectory[]~~"map2.csv","CSV"];
map2Fmt=Inner[Rule,Range[First[#]],Rest[#][[1;;First[#]]],List]&/@map2;

Convert all this values measured in pixels to micrometers, using a different conversion factor for each of the sources

factor1=0.21;factor2=0.28;

Cleaning the data by dropping components with no classification associated

Data generation:

generateData[headMask_,labels_,factor_:1]:=Block[{components,features},
(
components=ComponentMeasurements[headMask,{"Length","Width","Area","PerimeterLength","Eccentricity"}][[All,2]];
components[[All,1]]=components[[All,1]]*factor;
components[[All,2]]=components[[All,2]]*factor;
components[[All,3]]=components[[All,3]]*factor*factor;
components[[All,4]]=components[[All,4]]*factor;
components=Append[#,#[[2]]/#[[1]]]&/@components;
Append[components[[#]],#/.labels]&/@Range[Length[components]]
)]

Generating data for a single image:

generateData[binarizedImages1[[1]],map1Fmt[[1]],factor1]

Generating data for all the images:

rawData1=Flatten[generateData[binarizedImages1[[#]],map1Fmt[[#]],factor1]&/@Range[Length[binarizedImages1]],1];
rawData2=Flatten[generateData[binarizedImages2[[#]],map2Fmt[[#]],factor2]&/@Range[Length[binarizedImages2]],1];
rawData=Join[rawData1,rawData2];
rawData=Select[rawData,Last[#]≥0&];
Dataset[<|"Total"->Length[rawData],"Normal"->Count[rawData,u_/;Last[u]==0],"Abnormal"->Count[rawData,u_/;Last[u]==1]|>]

Converting data to Machine Learning training data format:

trainingData=Most[#]->Last[#]&/@rawData;
trainingData[[15;;20]]

Classification using Support Vector Machines

Statistical approach to machine learning, where each entity (sperm cell’s head) is represented by a set of features

Support Vector Machines are widely used along bioinformatics application because of the high levels of accuracy obtained, the ability to handle multidimensional data and its flexibility

Data Standardization

Using the Standardize function of Mathematica the data was modified and scaled so they have mean 0 and variance 1

This procedure is recommended to obtain better results when using SVMs: same opportunities to all the features, regardless of their numeric value and reduction of loss of information due to numeric operations

dataStandardized=Standardize[#[[1]]]->#[[2]]&/@trainingData;
dataStandardized[[1;;5]]

Select features: PCA

This method allows to extract relevant information, reduce the problem dimensionality and discover non trivial relations between the features, redundant data, etc.

The PrincipalComponents function of Mathematica transforms a data matrix into its principal components matrix by transforming the columns on the original matrix in non-related (independent) columns sorted from left to right in descending order of variances

principalComponents=dataStandardized;
principalComponents[[All,1]]=PrincipalComponents[dataStandardized[[All,1]]];

Simple SVM Accuracy Test Function

simpleSVMAccuracyTest[data_,times_:100,kernelType_:Automatic,gamma_:Automatic,margin_:Automatic]:=Block[{trainingSet,testSet,defaultSVMClassifier},
(
ParallelTable[(trainingSet=RandomSample[data,Ceiling[0.7*Length[data]]];
testSet=Complement[data,trainingSet];
defaultSVMClassifier=Classify[trainingSet,Method->{"SupportVectorMachine","KernelType"->kernelType,"SoftMarginParameter"->margin,"GammaScalingParameter"->gamma}];
ClassifierMeasurements[defaultSVMClassifier,testSet,"Accuracy"]),times]//Mean
)]

Selecting components

PCAAccuracyTest[data_,components_]:=simpleSVMAccuracyTest[#[[1,;;components]]->#[[2]]&/@data,10];
accuracyList=((#->PCAAccuracyTest[principalComponents,#])&/@Range[6])//Association

(List@@DiscretePlot[accuracyList[x],{x,1,6},PlotRange->{{0,6},{0.5,0.9}},FrameLabel->{"Principal Components","Accuracy"},PlotTheme->"Detailed"])[[1]]

For more than three components, the accuracy percentage does not get considerably better.

data=#[[1,;;3]]->#[[2]]&/@principalComponents;

Parameter calibration

Kernel type

The selection of the kernel type depends on the amount of data, the number of features, etc.

Some of the most popular kernels are Linear, Polynomial and RadialBasisFunction

The RadialBasisFunction kernel is generally a good choice:

It handles well non-linear relations between the features and the class of the entities

The linear kernels may behave similarly for some parameter values

The difficulty in the selection of the algorithm´s parameters is lower when compared with the polynomial kernel

There were three tests performed with the Linear, Polynomial and RadialBasisFunction kernels:

accuracyList=(#->simpleSVMAccuracyTest[data,10,#])&/@{"Linear","Polynomial","RadialBasisFunction"}

Grid search

The parameters for the SVM algorithm depend on the kernel type selected

For the RBF kernel there are the penalty factor C and the kernel parameter gamma

The strategy was a grid search, considering different pairs (C,gamma)

The values of the parameters were incremented exponentially, to escape local optima

Loose grid-search

accuracyValues=Table[simpleSVMAccuracyTest[data,100,"RadialBasisFunction",(2^i),(2^j)]*100,{i,-7,3,1},{j,-8,2,1}];
cTicks=Transpose[{Range@11,Range[-7,3,1]}];
gTicks=Transpose[{Range@11,Range[-8,2,1]}];
ArrayPlot[accuracyValues,PlotLegends->Automatic,FrameTicks->{cTicks,gTicks},FrameLabel->{"SoftMarginParameter","GammaScalingParameter"},ColorFunction->"Rainbow"]

Fine grid - search

accuracyValues=Table[simpleSVMAccuracyTest[data,100,"RadialBasisFunction",(2^i),(2^j)]*100,{i,-1,0,0.1},{j,-6,-4,0.2}];
cTicks=Transpose[{Range@11,Range[-1,0,0.1]}];
gTicks=Transpose[{Range@11,Range[-6,-4,0.2]}];
ArrayPlot[accuracyValues,PlotLegends->Automatic,FrameTicks->{cTicks,gTicks},FrameLabel->{"SoftMarginParameter","GammaScalingParameter"},ColorFunction->"Rainbow"]

Average Accuracy

accuracyValues=Table[simpleSVMAccuracyTest[data,1,"RadialBasisFunction",(2^(-0.5)),(2^(-5))]*100,100];

Histogram[accuracyValues,AxesLabel->{"Accuracy (%)","Tests"}]

Mean[accuracyValues]

Testing the Model

applyModel function

Options[applyModel]={"ShowComponents"->False};

applyModel[image_,classifierFunction_,actualClass_,factor_:1,OptionsPattern[]]:=Block[{components,stdFeatures,principalComponents,predictedClass,outputData,legend={-1->"Not classified",0->"Normal",1->"Abnormal"}},
(
components=ComponentMeasurements[image,{"Length","Width","Area","PerimeterLength","Eccentricity"}][[All,2]];
components[[All,1]]=components[[All,1]]*factor;
components[[All,2]]=components[[All,2]]*factor;
components[[All,3]]=components[[All,3]]*factor*factor;
components[[All,4]]=components[[All,4]]*factor;
components=Append[#,#[[2]]/#[[1]]]&/@components;
stdFeatures=Standardize/@components;
principalComponents=PrincipalComponents[stdFeatures][[All,;;3]];
predictedClass=classifierFunction/@principalComponents;
outputData=Transpose[{components,actualClass/.legend,predictedClass/.legend}];
outputData=Flatten/@outputData;
outputData=Prepend[outputData[[#]],#]&/@Range[Length[outputData]];
outputData=Inner[Rule,{"ID","Length","Width","Area","Perimeter","Eccentricity","Ratio","ActualClass","PredictedClass"},#,Association]&/@outputData;
If[OptionValue["ShowComponents"],
{showComponents[image],outputData},
outputData]
)]

Test cases

customClassifier=Classify[data,Method->{"SupportVectorMachine","KernelType"->"RadialBasisFunction","GammaScalingParameter"->(2^(-5)),"SoftMarginParameter"->(2^(-0.5))}]

i=8;
testImage=binarizedImages1[[i]]

result=applyModel[testImage,customClassifier,map1[[i,2;;map1[[i,1]]+1]],factor1,"ShowComponents"->True];
result[[2]]//Dataset

result[[1]]

Acknowledgements

Special thanks to my advisor Phd. Cesar Beltran Castañon, profesor at the Pontificia Universidad Catolica del Peru for his support and his advices through the project

References

Image sources:

Chang, V., Saavedra, J. M., Castañeda, V., Sarabia, L., Hitschfeld, N., & Hartel, S. (2014). Gold-standard and improved framework for sperm head segmentation. Computer Methods and Programs in Biomedicine, 1-13.

World Health Organization. (2010). WHO laboratory manual for the Examination and processing of human semen (5th ed.). Geneva, Switzerland: World Health Organization.

About Support Vector Machine parameters calibration:

Chih-Wei,H.,Chih-Chung,C.,& Chih-Jen,L.(2010).A Practical Guide to Support Vector Classification.Taiwan:Department of Computer Science,National Taiwan Univeristy.

The end

Join the Conversation #WolframTechConf