A schematic overview of all steps from sample preparation to phenotyping is presented in Fig. 1, where the different steps can be seen grouped according to their stage. First, the steps related to the preparation of the samples are mentioned, followed by the processing and analysis of the images. Afterwards the image clustering to finally obtain the phenotypic information is presented.
Sample preparation
For each of the ten genotypes, four cuboids from individual beet roots were embedded in paraffin. At least three sections per beet root were used for staining and one representative section was chosen for image acquisition and further processing as described in material and methods.
Preprocessing of images and definition of the region of interest (ROI) by identification of non-tissue containing areas and objects touching the edges of the image
As a preprocessing step, every original image (Fig. 2A) was converted to grayscale. Furthermore, interluding, noise reduction, and opening were carried out as recommended and described in [13]. However, at the filtering and opening steps, the parameters of each function were slightly adapted to our images, using a value of sigma = 2 for the standard deviation of the Gaussian filter and by selecting the brush shape = disc, with size 18 and 21 respectively for the first and second opening cycle. Also, automatic thresholding (Otsu method, [14]) was applied. At this point, the images were converted to black and white (Fig. 2B).
Elimination of irregular edge of the peridermal area
Due to the rounded shape of the periderm, each image shows areas of the image that could not contain any biological target tissue. For the correct functioning of the workflow, the periderm should be oriented towards the left edge of the image. Figure 2A* shows an example of an irregular edge. This area is prone to hold air bubbles and other materials that are potential artifacts. Therefore, it is necessary to mask and remove this area before starting the process inside the target tissue material.
For identification and removal of this area, only the first 5% from the left image border towards the right are considered to keep the workflow efficient. A high opening (erosion followed by dilatation; brush = 50) and hole filling identified the objects of interest. Holes are defined as areas of dark pixels surrounded by lighter pixels [15]. This area is removed from the processed black and white image.
Removal of cells touching the edges of the image
At the edge of the images, incomplete cells can be seen because parts of the cells are lying outside of the photographic frame. To assess cell sizes, it is necessary to consider only complete cells; therefore, partial cells touching the edge must be eliminated.
Objects whose pixel reach the edge of the image were identified using a process proposed by the Bioconductor project [16]. The process created a mask of the objects belonging to the border and another mask with the objects not belonging to the boundary. The number and morphological features of the identified edge-touching objects were documented. The information of the left edge was merged with the incomplete edge cells to perform the edge feature calculation (Fig. 2C).
The region of interest (ROI) is the area after removing non-tissue-containing sites and edge-touching objects from the original image. Finally, the ROI was reapplied to the grayscale image to continue processing (Fig. 2D).
Closing damaged cells and final individual features extraction
Due to the delicate nature of cells, cells might be damaged during the preparation process of the tissue sections, and their surface becomes connected. The watershed algorithm allows cell walls to be rejoined. Therefore, it is necessary to generate a matrix with each pixel's distance inside the cells from the nearest cell wall pixel before determining the watershed. In our matrix, the cell interior is considered the foreground and the cell wall the background. The distance was computed with the distance map function (Distmap) in the Euclidean way.
The watershed algorithm segments the ROIs, and each object's morphological features were calculated. The area corresponding to the cell membranes and the site of the cell interior was determined in this step. In the tissue density image (Fig. 2D), high values (black) correspond to cell membranes, and low values (white) to the background. The amount of cell membrane per image was determined. With the function "ComputeFeatures" of the EBImage package, each segmented object's morphological information or features in the images were extracted into a table. This table consists of data (columns) for each detected object (lines) in the image. The following information about the shape was recorded: area, perimeter, mean radius, the standard deviation of the mean radius, max radius, and min radius. The information about the object image moments contained the center of mass x, mass y, elliptical fit major axis, elliptical eccentricity, and object angle.
The diameter of each cell was calculated as 2*mean radius; the conversion from pixel to micrometer was calculated as diameter multiplied by the factor 1.015228, representing the pixel size for the image on the 10x microscope used.
Identification of number and position of cambium rings
The cell diameter is shown on a histogram as a distance from the left edge (Fig. 3). The beetroot tissue is composed of cells of different sizes. The maximal cell size increased from the outer part of the beet on the left towards the inner part on the right. The distribution of cell sizes along the x-axis represents the alternating occurrence of small-celled areas where the cambium rings are located and big-celled regions storing parenchyma. Using part of the code and algorithm suggested by [17] using the S3 Infrastructure for Regular and Irregular Time Series from ZOO package [18], the position of the peaks (storing parenchyma) and valleys (cambium rings) in the histogram were determined. Peaks were detected after the application of smoothing to find the local maximum (using as window width value of w = 150 and as a span argument for the loss function of span = 0.2). Valleys were detected using the same method after inversion of the y-axis. In addition to the position, the number of cambium rings was also computed by numbering them from left (outermost ring) to the right starting with 1.
Detected valleys and peaks are listed and numbered from the left edge inwards. Thus, it is possible to classify the cells between which the number of valleys or peaks is found. This information was added to the table obtained by ComputedFeature. This additional information allows the differentiation of cells into smaller groups during clustering. For example, in the case of the area closest to the epidermis, where small cells are typically present, the presence of a cell of a higher proportion than these would not be identified as being different if its size is considered with the rest of the cells in the sample.
Distinction between real cells, intercellular spaces, and artifacts
After segmentation, several objects did not necessarily represent cells. This was due to biological reasons like voids occurring between adjacent cells (intercellular spaces) as well as artifacts generated due to the nature of the image itself or the filtering process by over-segmentation. For the identification of as many of these non-real objects or spaces as possible and their subsequent elimination from the analysis, different strategies were used, as explained below.
First, all objects of small size were dropped for the differentiation of cells from non-cellular objects. To accomplish this, the threshold valleys function of the benmack/threshold package was used. First, a frequency histogram of the size of the objects was determined based on the log10 (diameter) of the features (Fig. 4). If the frequency deviates from a normal distribution by showing two peaks, the benmack/threshold package determines the intermediate point between the two peaks. This point is the valley of the curve and the newly determined cut-off point to drop the lower outlier that escapes the normal distribution. Objects with a log10 (diameter) below the threshold value were eliminated. If no threshold point was found, a value of 0.5 was used, standing for a diameter of 3 µm. If valleys were located above 0.8, 0.8 was considered the new value threshold to avoid dropping viable cells.
The remaining objects can either be cells or intercellular spaces. Careful microscopic examination of the tissue revealed that the intercellular spaces mainly occur between large cells. They might be bigger than small cells in other tissue areas, and they can have quite different shapes. However, the intercellular spaces of beetroots are always smaller than the surrounding cells (Fig. 5).
For each object a subsequent analysis of the neighboring objects was performed in a two-step-process to identify the intercellular spaces. First, outlier detection analysis was conducted based only on the diameter parameter of neighboring cells. Therefore, based on the x and y coordinates of the objects, the p nearest neighbors were identified by the k-dimensional tree (nn2 function of the RANN package). The maximum number of nearest neighbors to compute was set as default to k = 6. Based on the distance of the neighboring objects to the object of interest, quantiles (Q) and Interquartile Range (IQR) were calculated. Objects with a gap above Q2 + 1.5 * IQR were not considered as direct neighbors of the object of interest. Quantiles and IQRs were calculated for log10 of the diameter of the object of interest and its immediate neighbors. Objects with a log10 (diameter) below the Q1–1.5 * IQR were considered as outlier and probable intercellular space. They were removed from the dataset, while all remaining objects were initially defined as cells.
The second, most resource-demanding strategy, is the clustering of the cells based not only on the diameter, but on all the morphological features obtained with “ComputedFeature” (except for data concerning the cell´s position; x and y, and the information about the ring number each cell is assigned to). Based on the possibility that there may remain 3 possible types of objects, respectively large cells, small cells and intercellular spaces, a clustering with 3 clusters was performed. All input data was centered and scaled with the R base function “scale”, followed by a k-means clustering from the Stats package was used with a random set of 25 (nstart), maximum iteration of 1000 (iter.max), and 3 different clusters (centers). Since the cluster value assignment (1, 2 or 3) in this function is random, a reclassification is needed to achieve that different images always receive the same order of cluster values. It was determined that the smallest cells should continuously be assigned a value of 1, and the largest cells should always be assigned a value of 3. After obtaining the clusters, the cells were repositioned in the matrix, and the neighboring cells were analyzed again, using the Kd-tree, with the difference that, in this case, the neighborhood of the different types of clusters was analyzed. Considering that when the cells surrounding the cell under analysis are large, and this analyzed cell is small, the "cell" is considered as intercellular space. The number of cells that can surround an intercellular space is variable, and many possibilities and combinations are found. Therefore, two different criteria were used for the identification of intercellular spaces: (1) Small objects, assigned to cluster 1 surrounded by big cells assigned to cluster 3 (Fig. 6A). (2) Small objects, assigned to cluster 1 surrounded by a majority of cells assigned as cluster 3 but some cells as cluster 2, for instance, three as cluster 3 and one as cluster 2 as displayed in Fig. 6B. The decision if the central cell is an intercellular space requires the definition of "majority". Therefore, the following parameters are required: Number of surrounding cells, sum (SUM) and mean (MEAN) of the assigned cluster values of the surrounding cells; and the possible maximum sum (MAX) that could be obtained if all surrounding cells were assigned to cluster 3 (number of surrounding cells * 3). If SUM is greater or equal than MAX minus MEAN, the central cell can be considered an intercellular space and be deleted from the dataset. The entire table is continuously analyzed until no outlier is found.
Final clustering and tissue discrimination
After removal of the identified intercellular spaces, in theoretical perspective, two groups of cells remain in the dataset, the small and the large cells, representing the two main tissue types of beet roots, vascular tissue and storage parenchyma. Therefore, a new k-means clustering with the remaining data was performed, with a random set of 25 and maximum iteration of 1000, but this time with two clusters. As during 3-cluster clustering, the cluster number assignment is random, and a reclassification based on size is required, always assigning 1 to the smallest and 2 to the largest cells. The two classes of cells identified by digital image processing are marked in distinct colors in Fig. 7. The two clusters discriminate the two biological tissue types, as cluster one (green) mainly holds vascular cells, while cluster two (blue) mainly holds big cells of the storage parenchyma. To facilitate understanding and differentiate the cluster names with the 3-cluster clustering, small cells reclassified with cluster 1 will be assigned the letters VT for vascular tissue and cluster 2 will be assigned the letters SP for storage parenchyma.
Gathering and storing data
After completion of the clustering process, a table with the generated data was created; additionally, the name of the image file and extra information, such as the amount of the cell walls, were recorded in the new table. The amount of cell walls in the images (variable name: cellwallcount) was determined as the sum of all black pixels within the ROI of the black and white image, thus representing all pixels within the ROI excluding the cytoplasm and intercellular spaces. Afterward, the whole process was restarted with the next image in the list of images to analyze. The tables generated for each image are compiled into a final master table that contains all the accumulated information of all the images contained in the image list (Table 1).
Table 1
Variables contained in the final master table of the digital image analysis. The final master table comprises all relevant data generated at different points in the process. Morphometric variables were measured by the ComputeFeature function. These variables were used to calculate further variables during preprocessing, processing and clustering. For each variable, the name and a description were included.
From computeFeature
|
From preprocessing, processing, and clustering
|
variable name
|
description
|
variable name
|
description
|
s.area
|
area size
|
meandiameterµm
|
diameter average calculated in µm
|
s.perimeter
|
perimeter
|
diameter_log
|
log10 of diameter
|
s.radius.mean
|
mean radius
|
Q1
|
first quantile
|
s.radius.sd
|
standard deviation of the mean radius
|
Q3
|
third quantile
|
s.radius.max
|
max radius
|
IQR
|
interquartile range
|
s.radius.min
|
min radius
|
nmbpeak
|
number of peaks
|
m.cx
|
center of mass x
|
nmdeep
|
number of valleys
|
m.cy
|
center of mass y
|
Cluster_drei
|
classification of three clusters
|
m.majoraxis
|
elliptical fit major axis
|
Cluster_zw
|
classification of two clusters
|
m.eccentricity
|
elliptical eccentricity
|
imagetotalpix
|
total image pixels
|
m.theta
|
object angle
|
innerpixarea
|
total pixel ROI
|
order
|
original order of cell sequence
|
cellwallcount
|
total pixel of cell wall
|
|
|
intercellcount
|
total pixel of cell spaces
|
|
|
dataname
|
image name
|
Validation
The correlation between manual and automated counting of sugar beet root cells is illustrated in Fig. 8. In the set of reference images (one image for each of the ten genotypes), 15100 cells were identified in total by manual identification, while 16790 were identified by automated image processing. The coefficient of determination was 0.98.
The applicability of the final clustering with two clusters was also assessed by expert interpretation of the groups regarding the known histology of the sample.
Phenotyping of different Beta genotypes
The first 6 mm from the periderm of ten Beta genotypes, nine sugar beet and one fodder beet, have been phenotyped based on the final data table created within the workflow. An average of 4081 (genotype 10) to 7259 cells (genotype 7) were identified from the four samples per genotype. The maximal cell diameter in fodder beet was 129 ± 14 µm, and for sugar beet, it ranged from 72 ± 19 µm in genotype 9 to 112 ± 23 µm in genotype 4. The mean cell diameter was highest in fodder beet (21.8 µm, genotype 10). In sugar beet, it ranged from 16.6 µm (genotype 7) to 19.5 µm (genotype 3). The number of cambium rings in the 6 mm root tissue ranged from 4.5 (genotype 1) to 5.75 (genotypes 2, 3, and 5). Statistically, significant differences occurred for the number of cells and the mean cell diameter, but not for the number of cambium rings (Table 2).
Table 2
Cell characteristics of the distal 6 mm of root tissue of ten Beta genotypes; 1–9: sugar beet, 10: fodder beet. Automatic determination from microscopic images of cell tissue with four samples per genotype. Mean and SD of the four images per genotype are indicated. MSD: minimum significant difference calculated by HSD test, alpha = 0.05. Different letters indicate significant differences between genotypes.
genotypes
|
number of cells
|
max cell diameter [µm]
|
mean cell diameter [µm]
|
number of cambium rings
|
amount of cell wall material
[% of ROI pixels]
|
1
|
6046 ± 747ab
|
82 ± 6ab
|
18.3 ± 1.1ab
|
4.5 ± 0.5a
|
0.35 ± 0.0076ab
|
2
|
6036 ± 1326ab
|
87 ± 13ab
|
18.3 ± 2ab
|
5.75 ± 1.2a
|
0.35 ± 0.0140ab
|
3
|
5319 ± 513ab
|
90 ± 10bc
|
19.5 ± 1.11ab
|
5.75 ± 0.5a
|
0.35 ± 0.0097ab
|
4
|
5174 ± 243ab
|
112 ± 23ab
|
19.2 ± 0.62ab
|
5.5 ± 0.58a
|
0.35 ± 0.0066ab
|
5
|
6288 ± 986ab
|
90 ± 11ab
|
17.7 ± 1.67ab
|
5.75 ± 0.96a
|
0.36 ± 0.0169a
|
6
|
6949 ± 1777ab
|
81 ± 19ab
|
17.0 ± 2.38b
|
5.25 ± 1.2a
|
0.37 ± 0.0162a
|
7
|
7259 ± 433 a
|
76 ± 14c
|
16.6 ± 0.68b
|
5.25 ± 1.2a
|
0.36 ± 0.0037a
|
8
|
5844 ± 1405ab
|
88 ± 23ab
|
18.6 ± 2.28ab
|
6.0 ± 1.6a
|
0.35 ± 0.023ab
|
9
|
7052 ± 1649 a
|
72 ± 19c
|
16.8 ± 2.44b
|
6.25 ± 0.9a
|
0.37 ± 0.0213a
|
10
|
4081 ± 998 b
|
129 ± 14a
|
21.8 ± 3.13a
|
5.5 ± 1a
|
0.32 ± 0.0201b
|
MSD
|
2713
|
31.8
|
4.6
|
2.5
|
0.037
|
The mean cell diameter was significantly negatively correlated to the number of cells per genotype (R2 = 0.98, p ≤ 0.001) and to the amount of cell wall material identified in the microscopic images (R2 = 0.96, p ≤ 0.001; data not shown).
The distribution of cell sizes in the 6 mm distal root tissue per genotype is shown as Kernel density estimate in Fig. 9. All genotypes showed a peak at cells with a diameter of approximately 12 µm. The comparison between genotypes showed that this cell diameter had the highest abundance in genotypes 6, 7, and 9, and the lowest in genotypes 3, 4, and 10. At larger cell diameters, this relation changed, and genotypes 3, 4, and 10 had the highest abundance of cells with a diameter above 30 µm, and genotypes 6, 7, and 9 had the lowest.
For each genotype, most cells were assigned to cluster VT. Only a small variation between genotypes was found in the ratio of cluster VT to cluster SP, 73% (genotype 9) to 79% (genotypes 4 and 10) of the cells assigned to cluster VT (Fig. 10A). Due to the smaller cell sizes, however, cluster VT only took up 28% (genotype 10) to 37% (genotype 7) of the image area (Fig. 10B). The mean cell size was calculated separately for each cluster and genotype. Genotypic differences in cell size were more pronounced in cluster SP than in cluster VT (Fig. 10C). The ratio of cluster VT to cluster SP was relatively stable across genotypes (number of cells, 10 A, and tissue area, 10 B). Genotypic differences in cell diameter were more pronounced in the storage parenchyma (cluster SP; 10 C).
For the similarity analysis of the genotypes in a dendrogram, the values of each image belonging to the final master table (Table 1), with exception of the centers of mass, Q1, Q3, IQR, file name and the original order of cell appearance, which were not considered for analysis. For the other variables, the mean, standard deviation, and maximum and minimum value were calculated. The distance-based dendrogram revealed two branches at a distance of 7.98, where genotype 10 (fodder beet), differs from the other nine sugar beet genotypes. For sugar beet, two large groups could be distinguished in which genotypes 9, 6, and 7 belong to one group, and the remaining genotypes to the other group (Fig. 11).