Automating Quantitative Confocal Microscopy Analysis

Quantitative confocal microscopy is a powerful analytical tool used to visualize the associations between cellular processes and anatomical structures. In our biological experiments, we use quantitative confocal microscopy to study the association of three cellular components: binding proteins, receptors, and organelles. We propose an automated method that will (1) reduce the time consuming effort of manual background correction and (2) compute numerical coefficients to associate cellular process with structure. The project is implemented, end-to-end, in Python. Pure Python is used for managing file access, input parameters, and initial processing of the repository of 933 images. NumPy is used to apply manual background correction, to compute the automated background corrections, and to calculate the domain specific coefficients. We visualize the raw intensity values and computed coefficient values with Tufte-style panel plots created in matplotlib. A longer term goal of this work is to explore plausible extensions of our automated methods to triplelabel coefficients.


Introduction
Light microscopes capture energy emitted from fluorescently labeled-proteins within a biological sample. Fluorescent labels are bound to molecules of interest in the sample. The corresponding pixel intensity in the captured image is proportional to the amount of molecule in the sample. Multiple molecules can be labelled simultaneously by using fluorescent labels with different excitation/emission spectra. We designed and executed a biological experiment to determine the presence of a binding protein and a receptor protein at sub-cellular structures over time. The experiment was analyzed by quantitative confocal microscopy and resulted in a set of 933 RGB (red, green, and blue) images. Colocalization of binding protein, receptor, and subcellular structure is represented by RGB intensities in a pixel. The cooccurrence of signal in multiple channels signifies interesting biological phenomena. Therefore, we employed statistical methods of colocalization to quantify co-occurrence of RGB. The following sections describe our methods of quantifying the data contained in these experiments.
Conventional light microscopes produce a two-dimensional image from a three-dimensional sample by flattening its Z-axis into one visual plane [Cro05]. Thus, the notion of depth is removed by merging deep and shallow material into a single cross-section in the XY-plane. Confocal microscopes maintain Z-axis fidelity by performing repeated scans of very thin (∼5µm) XY-sections at fixed depths. A stack of confocal images represents the original three-dimensional sample. In an RGB confocal image, the brightness of a two-dimensional pixel represents the intensity of fluorescence in each of the three RGB color channels.
Background noise is the portion of the intensity signal that does not represent true biological phenomena. Confocal microscopy inherently reduces background noise from autofluorescence of cellular material, light refractive scatter, and detection artifacts [Cro05]. It is further reduced by choosing appropriate (1) microscope hardware, (2) fluorescent labels, and (3) computer software settings [Bol06], [Cro05]. Even the best confocal microscopy technique and practice produces images that contain background noise. For a detailed description of basic confocal optics and digital imaging, see [Bol06]. Pre-processing tools decrease background noise, but images often need additional manual background correction [Bol06], [Zin07], [Gou05]. Image processing filters, deconvolution, background subtraction and threshold techniques reduce background noise using different algorithms [Rob12]. Each technique has application specific advantages and weaknesses.

Biological Context and Experimental Model
We used confocal microscopy to investigate the post-endocytosis transport of two proteins in neurons. Specifically, we assessed the localization of binding proteins and their receptors to sub-cellular structures. Post-endocytosis transport of proteins is a highly regulated, complex process [Yap12]. Briefly, the intracellular transport pathway is initiated when an extracellular protein binds to its receptor on the cell membrane. Once internalized, the proteins may be localized to three sub-cellular structures: endosomes, lysosomes, and recycling vesicles. Proteins are internalized in endosomes, degraded in lysosomes, and transported back to the cell membrane in recycling vesicles. In our model, neuroblastoma cells were treated with a binding protein over different treatment times (10, 15, 30, 60, or 120 minutes). Following binding protein treatment, we stained cells for binding protein (red), receptor (green), and sub-cellular structure (blue). In different treatments, blue represents different sub-cellular structures. We performed six replicates of each condition, resulting in 6 Series for each condition. At each experimental Time, a set of 6 image stacks were captured with 5-12 optical XY-sections comprising one stack.
In these experiments, the binding protein is brain-derived neurotrophic factor (BDNF), the receptor is the truncated trkB receptor (trkB.t1), and the sub-cellular structures are endosomes, lysosomes, and recycling vesicles. For the biological importance of this system, see [Fen12]. The co-occurrence of red, green, and blue represents the presence of BDNF and trkB.t1 at one of the sub-cellular structures.

Manual Thresholding
We applied a manual thresholding procedure to reduce background noise. For each channel (R, G, and B) within an image, we (1) visually assessed the single-channel histogram and determined a threshold intensity, (2) mapped all intensity values at or below the threshold to zero, and (3) linearly scaled the remaining values to the intensity range [1, 255]. Additionally, we recorded the range, [low, high], around the manual threshold value that resulted in equivalent expert visual perception of the thresholded image. The thresholding procedure was repeated for each channel. Consequently, all intensity values for red, green, and blue below their respective thresholds are attributed to background noise and discarded. The major drawback to manual thresholding is the large time involvement of an imaging expert. Within-and betweenexperimenter reliability, differences in color output between visual displays, and access to expensive software packages are additional drawbacks to manual thresholding.

Automating the Thresholding Procedure
Initially, we manually determined threshold values for one randomly selected stack per experimental condition called our training set. Later, we manually thresholded the entire image set. Using the training set, we developed a linear regression model of the manual thresholds. Applying this linear model, we predicted thresholds for the full image set.
To generate automated background thresholds, we first extracted the deciles of the intensity histograms after removing non-responder pixels (see Visualization of Colocalization). Then, we considered linear regression models from (1) the intensity deciles and the channel to (2) the midpoint of the expert threshold range. For model development, we used only the training set of images. Our initial model included all deciles and the channel. Only the 8th and 9th deciles (80-and 90-percentiles) and the channel had statistically significant coefficients. We retained only these features in our model with a resulting R 2 of 0.6907 with p < 2.2e − 16. We evaluated the predictive ability of the model on the full dataset. The mean absolute error against the midpoint was 6.1313; the mean distance from the [low, high] threshold range was 2.2662. While these metrics are encouraging, we are more interested in the overall effect of automated thresholding on the computed colocalization coefficients, discussed below.
Finally, we compared the images generated by applying manual and automated thresholds. Both methods produced visually similar images (Figure 1). In both cases, the greatest amount of background correction occurred in the green channel. This is expected due to natural autofluorescence of cellular material in the green channel. However, the green channel also demonstrated the greatest difference between methods: the automated method under-corrected.

Visualization of Colocalization
In total, the images contain approximately 1 billion pixels. Only a small percent of the pixels represent protein, receptor, or subcellular structure. Therefore, the majority of the image pixels have zero intensity in all channels. These pixels are non-responders and are removed from further analysis. Channels values of 255 are considered to be over-saturated and are removed because they likely represent experimental or imaging artifacts. We computed the bivariate probability distributions of intensity values for each pair of channels across Time and Organelle. Due to the very large probability mass for low intensity values, we graphed the logprobabilities to visualize the behavior of the distribution tails. We generated a Tufte-style [Tuf01] panel plot of the bivariate histograms for all conditions. The panel plot for Time=10, Or-ganelle=Endosome is shown in Figure 2.
From the panel plot, we see that the bivariate distributions under manual and automated thresholding are qualitatively similar. For example, the RG histograms show low green intensities distributed over a wide range of red, with green showing a skew towards higher red intensities. The RB histograms show more even distributions over both channels. The GB histograms show lower green intensities over a wider range of blue. The patterns are the same for both thresholding methods. Next, we discuss quantitative assessments of colocalization.

Quantification of Colocalization
In dual-and triple-label confocal microscopy, several measures of association are used to quantify the degree of colocalization among labeled molecules [Bol06], [Zin07]. The two most commonly used measures are Pearson and Manders coefficients [Man92], [Man93], [Com06], [Zin07]. Other measures of colocalization are described below. We call all of these measures the colocalization coefficients.
Here, we consider the two-dimensional grid of RGB pixels as three one-dimensional vectors of intensity values for each color channel. In analogy with the moments of a random variable (as opposed to sample statistics), we define the colocalization coefficients for vectors x and y of the same length n.
Let mean(x) = sum(x)/n, dot(x, y) = ∑ i x i y i , cov(x, y) = dot(x − mean(x), y − mean(y))/n, and var(x) = cov(x, x): Pearson(x, y) = cov(x, y)/ var(x)var(y) The split k-overlap coefficients are: Let θ xy be the angle between x and y and recall dot(x, x) is the length of x: Manders(x, y) = cos(θ xy ) = dot(x, y)/ dot(x, x)dot(y, y) When a threshold is applied manually, the background noise is minimal (E-H). Automated thresholding methods reduce background noise to similar levels compared to manual thresholding (I-L). The green channel has more background noise after automated thresholding (K), compared to manual (G). Panels A, E, and I are RGB; Panels B, F, and J are the red channel; Panels C, G, and K are the green channel; Panels D, H, and L are the blue channel. The black and white panels are detailed views of the outlined squares in the left-most column.
Generally, the colocalization coefficients have the following interpretations when applied to vectors. Pearson is the degree of linear relationship between the two vectors. Pearson 2 is the fraction of the variance in y explained by the linear relationship with x.
Manders, more broadly known as the cosine similarity, is the cosine of the angle between the two intensity vectors. m 1 is the proportion of x, summed when y is above threshold, to the sum total of all x values; m 2 is likewise for y. k 1 (equivalent to cos(θ xy )length(x)/length(y)) is the ratio of the length of x and y times the cosine similarity between them.
In colocalization analysis, the colocalization coefficients have the following semantics. Pearson describes the linear relationship between two channels. Manders describes the directional similarity between the two channels. Thus, Manders is not sensitive to variation in total intensity, which may happen with different fluorophores. m 1 describes the amount of channel one intensity when channel two is on to the total amount of channel one intensity. k 1 is similar to Manders, but weights the degree of directional similarity by the ratio of the lengths of x and y.  [Bro00], [Phe01], [Val05], [Li04], [Rei12].
We computed the set of all colocalization coefficients efficiently by noting the common mathematical components of the coefficients and computing the common values only once. In the m-coefficients, the threshold T x is taken to be zero, since the coefficients are computed after manual or automated thresholding. 1 import math 2 import numpy as np 3 from numpy.core.umath_tests import inner1d 4 # inner1d computes inner product on last dimension 5 # and broadcasts the rest

Colocalization Coefficient Results
We computed the colocalization coefficients, for the manual and automated threshold images, over each time point for the Endosome organelle after grouping image stacks (Figure 3). The coefficients were used to compare the effects of manual versus automated thresholding on the scientific interpretation of the confocal images. For this analysis, correlation coefficients were calculated for each channel pair (Table 1). In the RG channel pair, there is a similar pattern seen between automated and manually thresholded images, for all correlation coefficient calculated (Figure 3). For instance, Pearson at Endosomes, 10, Manual is 0.32±0.02 (mean ± standard error over Series) while for Endosome, 10, Automated is 0.35±0.01. The Pearson coefficient for Endosomes, 30, Manual is 0.55±0.03 and Endosomes, 30, Automated is 0.55±0.03. By Endosomes, 60, the Pearson's coefficient for Manual is 0.35±0.04 and Automated is 0.39±0.03. The scientific interpretation of the coefficient data, regardless of Manual versus Automated, suggests that binding protein (red) and receptor (green) are associated with each other at all times, but that their greatest association occurs 30 minutes post-treatment time. The same conclusions are obtained from interpreting Manders (Table  1). We can use the combined data from all channel pairs to develop a model of intracellular localization of binding protein and receptor.

Applications
The automated background correction method we used can be applied to images generated from any type of microscopy studies including wide-field, live-cell, and electron microscopy. A second biological application for background correction is microarray analysis. Microarrays are tools used to study experimental differences in DNA, protein, or RNA, which often produce very large datasets [Hell02]. Multi-channel microarray experiments have similar background noise challenges as confocal microscopy. Most microarray experimental data is captured in the form of twocolor channel images with background noise generated from nonspecific label binding or processing artifacts. A third biological application for our automated thresholding method is magnetic resonance imaging (MRI) [Bal10]. In MRI images, background correction is often needed for phase distortion and general background noise. While other methods need to be applied to correct for phase distortion, our methods could be applied to reduce general background noise. Other biological applications include 2-D protein gel electrophoresis, protein dot blots, and western blot analysis [Dow03], [Gas09]. For any of these techniques, the background noise in the resulting images must be corrected prior to quantification of biological phenomena. Non-biological applications for our background correction method include, but are not limited to, photo restoration and enhancement [Dep02]. The correlation coefficient processing can be applied in many of these applications or any generic RGB image workflow.

Conclusions
Confocal microscopy is a powerful tool to investigate physiological processes in morphological context. Quantitative analysis of confocal images is possible using optimized image capture settings, background correction, and colocalization statistics. We used confocal microscopy to quantify the intracellular colocalization of a binding protein and a receptor to a specific organelle, over time. There were two major hurdles: (1) the time and consistency required for manually thresholding a large number of images and (2) batch processing of large image sets for statistical analysis. In 2005, Goucher et al. developed an open source image analysis program, in Perl, to batch process colocalization for RGB images using an ad hoc association metric [Gou05]. The purpose of our methods was to further this type of automated process to combine automated thresholding with batch processing of colocalization coefficients using Python. The benefits of our model are: (1) reducing the time consuming effort of manual background correction and (2) batch processing of multiple correlation measures for multi-color images. While our experiments focus on applying automated quantification methods to better understand intracellular protein transport, our computational methods can be used to study a wide range of biological and non-biological phenomena. A longer term goal of this work is to explore plausible extensions of our automated methods to triple-label coefficients.
Source code, under a BSD license, for computing colocalization coefficients, panel plots, and various other utilities is available at https://github.com/mfenner1/py_coloc_utils .  Pearson, Manders, m-, and k-overlap coefficients were calculated for manual and automated threshold images. The coefficients were calculated for each channel pair. Similar patterns for correlations coefficients are seen between manual and automated threshold images. The data in this figure was taken from the experimental condition Endosomes (i.e., B represents endosome) over all Times and Series . Values in one vertical line, a strip, come from the six repeated Series in that condition. Left to right, triples of strips are from increasing Time.