Automated segmentation and analysis of retinal microglia within ImageJ

Abstract

Microglia are immune cells of the central nervous system capable of distinct phenotypic changes and migration in response to injury. These changes most notably include the retraction of fine dendritic structures and adoption of a globular, phagocytic morphology. Due to their characteristic responses, microglia frequently act as histological indicators of injury progression. While algorithms seeking to automate microglia counts and morphological analysis are becoming increasingly popular, few exist that are adequate for use within the retina and manual analysis remains prevalent. To address this, we propose a novel segmentation routine, implemented within FIJI-ImageJ, to perform automated segmentation and cell counting of retinal microglia. We show that our routine could perform cell counts with accuracy similar to manual observers using the I307N Rho model. Tracking cell position relative to retinal vasculature, we observed population migration towards the photoreceptor layer beginning 12 h post light damage. Using feature selection with Chi² and principal component analysis, we resolved cells along a morphological gradient, demonstrating that extracted features were sufficiently descriptive to capture subtle morphological changes within cell populations in I307N Rho and Balb/c TLR2^-/- retinal degeneration models. Taken together, we introduce a novel automated routine capable of efficient image processing and segmentation. Using data retrieved following segmentation, we perform morphological analysis simultaneously on whole populations of cells, rather than individually. Our algorithm was built entirely with open-source software, for use on retinal microglia.

Keywords: Automated analysis; Cell counting; ImageJ; Kernel principal component analysis; Light damage; Retinal degeneration; Retinal microglia.

Figure 1.. Cell segmentation workflow for I307N Rho microglia.

Raw input images with an original magnification of 10x (a) were first preprocessed by normalizing stack intensities followed by rolling ball, median filter, and Gaussian smoothing routines (b). Binary candidate cell masks were produced by Otsu’s method (Otsu, 1979) and were thinned by 3D erosion. Skeletal structures (red lines) of un-eroded cells were used to repair small breaks introduced during erosion (c). Cell somas were identified by morphology using an existing algorithm (Davis et al., 2017), and their overlay with the candidate cell masks is shown in (d). Marker-controlled watershed was then used to identify distinct cells, as indicated by false color, (e) and morphological features of the cells and soma were saved. Overlay of the borders of detected cells (red) with the raw grayscale image (f) indicated high overall fidelity with occasional undetected cells and dendrites. Typical algorithm runtime was <1min. The shown example image is from the three days post LD condition.

Figure 2.. Vasculature segmentation workflow for I307N Rho mouse.

Raw 10x input images were preprocessed under rolling ball filter and grayscale attribute filtering to reduce background noise and enhance the apparent brightness of the vasculature (a). A global intensity threshold was passed as two standard deviations above the mean pixel value, effectively removing image background (b). Small breaks introduced to the continuous vasculature by this threshold were then repaired by morphological closing with an octagonal structuring element. Beginning at three days post-LD, microglia associated with the vasculature began to co-stain positive for IB4. To remove these cells, a strict circularity threshold was passed, removing particles with circularities greater than 0.1 for each slice of the image stack. While some of the vasculature was lost to this threshold, it did not affect the overall vasculature structure that was needed to characterize microglia position. The structure retained by this threshold was overlaid with the removed cells and vasculature shown in red (c). The vasculature was then skeletonized (d), and the node positions, highlighted in red, were saved to track cell position within IPython. The shown example image is from the three days post LD condition.

Figure 3.. Example vasculature mesh at 24 hours post-LD.

Produced point meshes describing the superficial (red) and deep (black) vasculature plexuses are shown for one time point (A). Microglia cell positions for this time point are also shown (blue). Cell position relative to either plexus is found by locating the nearest neighbor to the cell in either point mesh. The convex hull used to generate this point mesh is shown in (B) by its edge vectors. Axes show position within the image stack, in microns.

Figure 4.. Regression analyses for manual and automated cell counts against the mean observer.

The algorithm was run on a set of 11 images sampled at baseline and 2, 4, 12, 24, and 72 hours post-LD in the I307N Rho mouse. Relative to the mean observer, automated counts were significantly lower for one baseline image set (shown in red), thought to be caused by poor SNR. The regression was repeated with this image removed from the set. Pearson R² coefficients and regression slopes were calculated for the whole set (shown in black) and with this image removed (red) and suggested that the algorithm counts were dependent on image quality but were not significantly different from manual counts. The intercept for each regression was set at zero, and a slope of one is shown by the dashed line.

Figure 5.. Cell counts against height above the deep vasculature plexus.

Plotted data include a histogram of raw cell counts with a kernel density estimate (black). The two subpopulations produced persistent bimodality against relative height, and so were fit to the sum of two Gaussian distributions (red and blue). Time is displayed in hours post-LD, and height decreases towards the RPE at the Y-axis. Based on qualitative observations of the histological data, the subpopulation originating from the OPL at baseline underwent recruitment and migration towards the ONL beginning as early as 4 hours post-LD. Widespread activation was further apparent by 12 hours post-LD. At this time the lower subset, modelled by the blue curve, contained a majority of the total population, and its mean lay beneath the deep vasculature plexus. Activated cells were estimated as those lying outside the vasculature plexuses or within 1 standard deviation from the lower population mean beginning 12 hours post-LD. Lower population means are shown as vertical black lines, and 1 standard deviation selection thresholds as vertical blue lines. The two vasculature plexuses are shown as dashed red vertical lines, where the deep and superficial plexuses are on the left and right, respectively.

Figure 6.. kPCA projection for I307N Rho mouse.

Kernel principal component analysis produces a characteristic inverted v-curve with small local variance normal to the curve and high variance along its length, primarily captured by PC 1 within the I307N Rho training set (A). PC 1 is sufficient to distinguish early (0 – 4 hrs post-LD) and late (12 – 72 hrs post-LD) time points. Early and late time course cells are shown individually in (B) and (C) A micrograph of 7 cells reveals progression from small ameboid to large and highly ramified cells (D), demonstrating the curve’s morphological relevance. A region of high cell density is apparent near the origin. To visualize projected cell density, a contour map was generated by Gaussian kernel density estimate. Contours demonstrating the percentage population enclosed are shown in (E) and confirm the presence of a higher density of cells localized at the midpoint of the projected curve. This region is a possible indication that kPCA was unable to distinguish a substantial portion of the total cell population and prompted us to investigate limitations to discrimination performance about the curve’s midpoint. A series of Chi² tests between early and late time point cells were calculated. Between each test, a small number of cells furthest from the middle of the curve were removed. A plot of the resulting p-values against the population percentage is shown (F) and suggests that these cell populations remain significantly differentiable down to around 26.6% of the total population (p < 0.05), or a range of about −1 to 0.6 along PC 1.

Figure 7.. kPCA projection for Balb/c Tlr2^−/− mouse.

kPCA results for Balb/c TLR2^−/− cells (A) demonstrate a v-curve similar to I307N Rho cells. A density map revealing the percentage of the enclosed population is shown (B), demonstrating a density profile similar to the I307N Rho model, but with a marginally greater density cluster near the origin, indicative of higher homogeneity in the population. Because histological data for the vasculature was unavailable, it was not possible to perform an analysis of resolving power similar to that shown in Figure 6F for Balb/c TLR2^−/− cells. Cells selected at random from 7 intervals evenly spaced along the curve once again demonstrate morphological distinction in the output (C). These results act to validate the algorithm created with the I307N Rho mouse and suggest kPCA is an effective means of reducing high-dimensional data to a low-dimensional morphological spectrum.