article: Protein subcellular localization has been systematically characterized in budding yeast using fluorescently tagged proteins. Based on the fluorescence microscopy images, subcellular localization of many proteins can be classified automatically using supervised machine learning approaches that have been trained to recognize predefined image classes based on statistical features. Here, we present an unsupervised analysis of protein expression patterns in a set of high-resolution, high-throughput microscope images. Our analysis is based on 7 biologically interpretable features which are evaluated on automatically identified cells, and whose cell-stage dependency is captured by a continuous model for cell growth. We show that it is possible to identify most previously identified localization patterns in a cluster analysis based on these features and that similarities between the inferred expression patterns contain more information about protein function than can be explained by a previous manual categorization of subcellular localization. Furthermore, the inferred cell-stage associated to each fluorescence measurement allows us to visualize large groups of proteins entering the bud at specific stages of bud growth. These correspond to proteins localized to organelles, revealing that the organelles must be entering the bud in a stereotypical order. We also identify and organize a smaller group of proteins that show subtle differences in the way they move around the bud during growth. Our results suggest that biologically interpretable features based on explicit models of cell morphology will yield unprecedented power for pattern discovery in high-resolution, high-throughput microscopy images.
  High-content screening of fluorescently tagged proteins has been widely applied to systematically characterize subcellular localizations of proteins in a variety of settings [1]. Because they employ automated liquid handling and high-throughput microscopy, these experiments result in large numbers of digital images. Previous work has demonstrated that automated image analysis approaches based on machine-learning can classify these images into groups with shared subcellular localization patterns [2]. These approaches are typically ‘supervised’ in that they rely on predefined sets of example ‘training’ images for each pattern of localization to learn specific discriminative information that defines each class [3]. In contrast, unsupervised methods offer a more exploratory approach to high-throughput data analysis in which it is not necessary to predefine patterns of interest, and therefore can discover new patterns. This also enables the analysis of patterns that are very rarely observed, which typically are hard to capture in supervised analysis as a suitable training set for classification is difficult to construct [1]. Unsupervised analysis also has the advantage that it is unbiased by prior ‘expert’ knowledge, such as the arbitrary discretization of protein expression patterns into easily recognizable classes. For these reasons, unsupervised cluster analysis has become a vital tool of computational biology through its application to genome-wide mRNA expression measurements [4]–[7], and protein-protein interaction data [8]. It has also been applied in automated microscopy image analysis [9]–[13] where it has been shown to provide complementary capabilities to supervised approaches. Here we apply unsupervised analysis to a set of high-resolution images of 4004 yeast strains, where each strain contains a different fluorescently tagged protein [14]. Because localization classes are not defined in advance, one difficulty is to identify a set of image features that reliably distinguish classes of protein expression [10]. Further, in order to allow identified statistical patterns to be directly related to our understanding of cell biology, we sought to define a small set of simple biologically interpretable measurements. This is in contrast to many automated image analysis approaches that use a large number of image features, which are typically used for object recognition in photographs [15], [16]. Although these features can be used to build powerful classifiers, the nature of the discriminative information does not need to be intelligible to allow class label recovery [3]. Recent work has demonstrated the power of incorporating cell-cycle stage into proteomics analysis (e. g. , [12], [13], [17], [18]). Several studies have identified proteins whose abundance and localization change over the cell-cycle in mammalian cells. Furthermore, unsupervised analysis has been applied to identify novel, unexpected patterns. In general, these approaches have been applied to time lapse movies of mammalian cells, although it is also possible to acquire dynamic data from still images of mammalian cells [18]. One advantage of budding yeast as a model organism is that it shows stereotypical cell-cycle dependent morphological changes, which can be used to infer cell-stage based on cell morphology in still images of asynchronous cells. Previous work has demonstrated the feasibility of uncovering and analyzing yeast morphology using automated image analysis methods [19], [20]. Although the identification of cell boundaries in images has been shown to be unnecessary for subcellular localization classification [2], [21], [22], in order to extract dynamic protein expression profiles based on changes in cell morphology, in this work we sought to accurately identify individual cells. Here, we use an explicit model of yeast cell shape in order to (1) rapidly identify cells in high-resolution images, even when they occur in clumps, (2) obtain a probabilistic confidence measure for the identified cells and (3) define biologically interpretable measurements that describe protein expression in each cell over space and time. We show that many previously defined subcellular localization patterns can be recognized in an unsupervised hierarchical cluster analysis. We find that protein complexes and small functional protein classes, which are not typically associated with their own subcellular localizations, cluster together in this analysis. Based on these observations, we show that the resolution of the hierarchical clustering is significantly higher than previous manual subcellular location assignments to discrete classes [14]. Further, we gain global insight into the cell stage dependence of protein localization; for example, we find a large cluster of nuclear proteins that seem to appear in the bud at a clearly defined time, which we believe corresponds to the inclusion of the nucleus in the daughter cell. Finally, we identify groups of proteins that show complex, dynamic patterns of localization that can not easily be predefined or described using simple localization classes; for example, many of the subunits of the exocyst complex are seen to localize to the bud periphery while the bud is small, but then move to the bud neck as the bud grows.
  Starting with a collection of 4004 strains where each protein has been systematically tagged with green fluorescent protein (GFP) [14], a red-fluorescent protein (RFP) which appears everywhere in the cell was introduced into each strain using SGA [23]. These strains were then imaged in quadruplicate at high resolution to generate two-channel fluorescent images (see Methods). The RFP was introduced to facilitate automated analysis, as it provides both a signal for cell segmentation, as well as an internal control for methodological variation in fluorescence measurements. Encouraged by the consistency and interpretability of our measurements relative to previous knowledge about yeast subcellular localization, we performed global unsupervised analysis of our time profiles of interpretable features for the mother and bud cells (Figure 6). To identify groups of proteins with similar patterns, we use agglomerative hierarchical clustering based on a maximum likelihood criterion [37] (see ‘Maximum likelihood agglomerative hierarchical clustering’ in Methods) because it does require the size (or number) of clusters to be specified, and we expect hierarchical relationships between functional classes, and a wide range in the number of proteins in each class. The hierarchical clustering results may be browsed online using the Java Treeview [38] applet at http: //www. moseslab. csb. utoronto. ca/louis-f/unsupervised/.
  Previous studies have demonstrated the feasibility of uncovering cell stage from images of unsynchronized cell populations, either from time lapse movies [17] or from still images [18]. We apply this approach to high-throughput still images of budding yeast. To do so, we devised a segmentation method to identify and separate the bud and mother cells, and uncover the cell stage based on measurement of the bud size. Our method depends critically on our estimates of bud size, and we show that the automatically estimated sizes were comparable to those obtained from manually identified cells. Several parts of the analysis may be improved. For example, since the bud-site selection is predetermined by the position of the preceding daughter cell [50], it could be used to help determine the correct mother-bud assignments. Similarly, a better model for the relation between daughter cell size and the cell cycle could be used to infer a more accurate estimate of cell stage. We presented a cell identification pipeline that includes a confidence measure which summarizes the probability that an object identified in our images is actually a correctly identified cell. To do so, we characterized the deviation of real cells from an elliptical model using several quality measures whose distribution for real cells we inferred from ellipses that had been manually fit to cells by eye. Our confidence measure allows us to distinguish correctly identified cells from artifacts and misidentified objects, without specifying what the nature of artifacts might be (Suppl. Figure S2). We believe that this type of approach for measuring the confidence of automatically identified objects in image analysis will be generally useful, because artifacts tend to vary between microscope, experiments and computational methods, whereas cell shapes are expected to be much more consistent. In addition, this confidence measure is explicitly defined as a posterior probability of an identified object to be a properly identified cell. This allows us to weight probabilistically data points according to the posterior probability. For classes of cells where our model does not fit as well, such as very early non-ellipsoidal buds, we expect to downweight all the data points, but we can still include information from these data points in our analysis. This is in contrast to the situation where we used a hard threshold to exclude artifacts. In that case, certain classes of cells are preferentially excluded (Suppl Figure 1B), and the statistical significance of downstream analyses is reduced (Suppl. Tables S2, S3, S4, S5). Typically, spatial patterns of protein expression are described by assigning labels [14] or functional annotations [51]. Such discrete classes are not sufficient to fully describe a protein' s expression if it is present in quantitatively different localizations or abundances at different cell stages, or if a protein is simultaneously present in several locations with quantitatively different fractions [52]; because our approach assigns a quantitative expression profile to each protein, we can characterize protein expression at a finer scale than the resolution currently achieved by discrete classes. Approximating protein expression patterns as discrete classes has also led to challenges for computational analysis. For example, in previous work based on discrete classes [16], [21] many proteins are often filtered out because they have either been annotated as ‘ambiguous’ or are reported to be located in several localization classes. Because we treat expression patterns quantitatively, our analysis identifies clusters of proteins that are significantly enriched in ‘ambiguous’ proteins and proteins that were manually annotated [14] as localized in multiple compartments. Furthermore, our analysis identified and organized a group of proteins that show complex patterns relating to the growth of the bud, that were not consistently annotated previously using discrete categorizations (Suppl. Table S6). To our knowledge no previous genome-scale analysis of still microscopy images has identified groups of proteins with subcellular localization patterns that change as a function of cell-stage, such as the MCM and exocyst complex subunits discussed above, although recent work on smaller collections of time-lapse images has demonstrated that functionally related proteins can be identified in unsupervised analysis of dynamic protein expression profiles [13]. One limitation of cluster analysis is that the members of each cluster identified are not always consistent between different parameter settings, or different clustering methods. Indeed, the remarkably specific groupings corresponding to specific regulatory mechanisms (such as the clustering of all 4 of MCM complex subunits and of all 3 DNA replication factor A complex subunits) were not always observed when we varied the distance metric or clustering method used (Suppl. Table S4). Despite these limitations, our analyses consistently identified clusters that were enriched in functional groups of proteins (Ribosome, Proteasome, DNA-damage pathway, exocyst complex, etc. ; see Suppl. Table S4, S5) that are not usually associated with their own subcellular compartments. Because we used hierarchical clustering of interpretable features, we could see that these functional groups of proteins showed patterns of localization similar to those localized in the same compartment, but in each case showed subtle differences in pattern that allowed them to be distinguished. These results suggest that high-resolution images could be used directly for functional discovery as has been reported for mammalian cells [12]. This work demonstrates that accurately identifying large numbers of cells for each protein allows quantitative characterization of spatial and temporal characteristics of protein expression patterns and permits direct interpretation of image-based measurements without requiring human inspection of large numbers of images to train classifiers. Our analysis gives new insight into the relationship between protein function and protein expression patterns inferred from high resolution microscope images.
  Using yeast synthetic genetic array technology [23], a new GFP collection was generated from the existing collection [14]. In this new collection, a highly expressed RFP (a tdTomato [53] fluorescent protein from the constitutive RPL39 promoter), integrated at the HO locus, was introduced into the GFP collection to mark the cell in order to facilitate automated image analysis. Micrographs were acquired using a confocal microscope (Opera, PerkinElmer). Eight micrographs were imaged (at 1331×1017,12 bit resolution) from each strain, 4 in the red channel and 4 in the green channel, yielding a dataset of 44 Gb of image data. It was noted that the background noise had a mean and variance that was not uniform across the image. Therefore, we defined a background image that was subtracted from each image. This background image was obtained by averaging all the images. The background image intensity accounts at most for a third of the RFP signal expected in mother cells, except for several defective CCD pixels which systematically report the same value. For each image, we modeled the background and foreground (cell) RFP intensity levels with Normal distributions. In order to account for punctuate noise, we used a Pseudo-2D hidden Markov model (P2DHMM) [54] to model the dependence of neighboring pixels. In order to recover the maximum likelihood parameters for the Normal distributions and state transition probabilities efficiently, we performed expectation maximization (EM) on both the image under the assumption that image rows are independent, and on the same image where columns are now assumed to be independent. Finally, we infer the probability for each pixel to belong to the foreground, as the average of the two probabilities that were calculated when we assumed rows and columns were independent. Given an image for which we know the probability of each pixel to be from the background, we want to define a map of geometric distance to background for each foreground pixel. We estimate this quantity using an iterative motion on the image grid (which includes diagonals and knight moves), where transitions from a point deterministically select the neighbor through which the shortest path to background is expected. We then compute the expected path length under the assumption that pixels reached along paths have background/foreground state transitions described by a HMM with the parameters inferred from the segmentation. The transition probabilities for diagonal and knight moves are obtained by exponentiation of the transition matrix by the distance between the two points. Since it is enforced that transitions are only allowed from point of higher expected distance to lower ones, distances can be computed directly by dynamic programming, in linear time of the number of pixels in the image. The Edge Distance map () has several uses in our pipeline: to generate the clump contours, as a quality measure for identified objects and to evaluate the distance of a protein to the periphery. We used robust regression for matching ellipsoidal shapes to the contour of the segmented area. An ellipse is characterized to be the set of points for which the algebraic error [55] is zero: (5) where is the coordinate of the ellipse centre and r is an additional parameter, proportional to the radius of a circle for a fixed matrix A. The matrix A may make the set of points with zero algebraic error correspond to a hyperbole or a line, and a superfluous scale parameter is observed in this parameterization. We therefore constrain the form the of matrix A: (6) where is the angle corresponding to the orientation of the major axis, and is a parameter which determines the eccentricity of the ellipse. This choice of this matrix to ensure that for any value for the set of 5 parameters (in equation 5) generates an ellipse with minor to major axis length ratio larger than, as they are both determined by the eigen values of the matrix A and then scale with the parameter ‘r’ [55]. Contour pixels are first identified by finding foreground pixels which are pixels away from some background pixel (using the Edge Distance Map described above). Initial guesses for ellipses are generated by first fitting a circle to 3 randomly sampled contour points (that circle is unique). Initial guesses are rejected if the circle does not fit within the rectangle clamping the contour points, or if the center is a background pixel. The initial guess ellipse will be set to match width (diameter) and center of an accepted circle. A small eccentricity corresponding to and a random angle (drawn uniformly from 0 to) is used to define its remaining parameters. If the set of contour pixels matches a single ellipse, we could directly update the ellipse coordinates by minimizing the sum of the algebraic error of all contour pixels. However, if the set of contour pixels is best explained by several ellipses, the sum of algebraic errors is likely to have local minima that are not close to any of the true ellipse parameters. Therefore, we use robust regression [56] and minimize the objective function: (7) where C is the collection of coordinates of contour pixels and is the expected error, which is chosen to be 5, matching the thickness of the contour. This effectively weights down the importance of contour points with large deviations to the current ellipse, so that the many local minima can correspond to actual ellipses. Upon convergence, we discard ellipses that are not bounded by the clamping rectangle, or that have a background pixel at the center. Since a large number of local minima are expected, we generate about 10 fold more sets of ellipse parameters than the number of expected ellipses (based on number of contour pixels) and select the ellipse with the best fit. Once we have identified the best ellipse, we remove all contour pixels that have an error smaller than, and find the next ellipse using the remaining contour pixels using the same procedure. Since some missed lone pixels may remain, we reject the ellipse and remove the corresponding pixels if the ellipse width is less than 3 pixels or if the number of removed contour pixels accounts for less than 10% of the amount expected from the ellipse parameters and known contour width. This process is repeated iteratively, until no more contour pixels can be removed. The running time of the segmentation is linear in the number of pixels in images, and the running time of cell-finding is linear with the number of randomly sampled circles for the initialization of geometric ellipse fit. On a single 2. 83 GHz Intel core, 98 seconds were required to analyze a single 1331×1017 image, which on average contained 82 cells and 31 artifacts. We want to precisely recover the cell shape, as we will be considering the size of the bud object as a cell-stage indicator, and the position of the bud neck as a point of interest for uncovering cell-stage dependent changes of protein localization. Because cells are not exactly ellipsoidal in our images, we sought to capture foreground pixels which partition the cell clumps into non-overlapping cell areas (which we refer to as ‘shapes’). In our images, cells are often separated by dim pixels, so we force boundaries to match these dim areas. We first use the watershed [57] transform to identify regions of the clump that are entirely contained within single cells. For each pixel which brighter than any of its neighboring pixels, we find the set of pixels (catchment basin) which can be reached by a path of monotonically decreasing intensities [57]. Secondly, we assign each basin to a cell, by finding the ellipse closest to each maximum intensity pixel. The proximity of a point to an ellipse is evaluated using the algebraic error (Equation 5). This procedure ensures that if two neighboring basins are assigned to different cells, we are guaranteed that the boundary pixels are all dimmer than the nearby inner pixels found inside one of the two basins. Each such basin is then assigned to the closest ellipse, such that the union of these regions forms the ‘shape’. The resulting shapes may be highly non-ellipsoidal (Figure 1iii); for example, if a cell has not be properly fitted by an ellipse, a portion of its area may appended to the area of a neighbor cell instead. In addition to the mean RFP intensity in the object, we define three shape measurements based on geometrical properties of ellipses and circles. First, we compute the best fit of an ellipse to an arbitrary shape ‘’ by evaluating the following 6 statistics on the coordinates of pixels in the shape (eq. 8). (8) where is the coordinate for a pixel from the shape ‘’, and is the number of pixels in the shape (cell size). A function defined on which takes the value ‘D’ within the area of an arbitrary ellipse has 6 degrees of freedom as well: (9) We can derive that there is a closed form for the parameters of the above function that makes the corresponding statistics defined on a continuous space match the statistics from the pixel coordinate of any shape. For instance, the center of the fitted ellipse will correspond to the center of mass of the provided shape: (10) The major and minor axes length (‘’ and ‘’) are the square root of the two solutions to a quadratic equation: (11) Finally, the recovered density ‘’ is the ratio of number of pixels to fitted ellipse area. Since the coordinates are drawn from a bitmap, we observe that the measured densities typically are bounded above by 1, except for the smallest objects. Any shape whose density is above or equal to 1 is assigned to the artifact class, otherwise we use the following first quality measure: (12) The second quality measure is based on the relationship between the perimeter and the area of an ellipse. We compute the perimeter of the shape by counting the number of pixels that have 3 or more background pixels among their 8 neighboring pixels. The theoretical relationship between the perimeter length of an ellipse and its parameters has no simple form, but may be approximated using the Ramanujan first approximation [58]: (13) The log ratio for the number of contour pixel to Ramanujan ellipse perimeter length approximation is our second confidence measure. A third quality measure captures the deviation of the shape to a circle, by reporting the log coefficient of variation of the sum of the distance to the ellipse center and the distance to the edge for each pixel in the area (eq 14). In a theoretical circle, there should not be any variance, since the two quantities are to sum up to be exactly the radius of the circle. (14) The last quality measure is the mean RFP intensity. We model each of the quality measures using a Normal distribution. We observe that the quality measure spread displays a non-trivial dependency on cell size. For this reason, we define 7 Normal distributions, for each of the 4 quality measures, which correspond to the distribution of quality measure for 7 bins of cell sizes. The quality measure vector is then modeled by the linear interpolation of a pair characterized random variables (eq. 15). (15) where and the are diagonal covariance matrices. We used the automatically identified shapes that were mapped to the 4305 manually identified cells in order to infer the parameters of the normal distribution at the 7 sizes (7 means and 7 standard deviations). In order to define the posterior probability of cell, it remains to characterize the uniform distribution for the non-cell objects and the mixing parameter (eq. 1). The uniform distributions were chosen to correspond to the extremum in quality measure obtained from the complete collection of identified objects that have not been labeled as artifacts. Finally, we used soft expected maximization (soft-EM) [59] on the complete collection to infer the mixing parameter, which rapidly converged to 9. 9% as all other parameters are already predefined. In order to evaluate the accuracy of our cell identification method, we first compared the automatically identified ellipses to a set of 4305 ellipses that had been drawn around cells manually. We assigned each manually identified ellipse to the automatically identified ellipse with closest center. We found that for 94. 2% of manually identified ellipses, there is an automatically identified one with center occurring within 10 pixels. In these cases, the average distance between the centers was 1. 86 pixels (). The correlation between the areas of the automatically identified and manually identified matched ellipse pairs was 0. 882. We next compared the center and area of the automatically identified ‘shapes’ to the set of manually drawn ellipses. Here, 92. 3% of the manually drawn ellipses have a corresponding recovered shape that has a center within 10 pixels of the manually drawn ellipse center. For these, the mean distance between the shape and the manually drawn ellipse center was 1. 41 pixels (). The area of the ‘shapes’ have a correlation of 0. 953 with area of the automatically identified ellipses, and 0. 928 with the area of manually drawn ellipses. Hence, identifying the ‘shapes’ (the hybrid operation of assigning the watershed regions to their closest ellipse) produces cells that are on average closer both in location (1. 41 vs. 1. 86 pixels) and size (correlation 0. 928 vs. 0. 882) to the manually drawn ellipses than the automatically identified ellipses. We note that the ‘shape’ -based analysis led to a slight reduction in the fraction of cells identified (92. 3% from 94. 2%) but this was acceptable to us in the context of the improvement in cell size estimation (0. 928 vs. 0. 882 correlation) because we use the cell size as an indicator of cell stage. In order to compare the accuracy of the simple cell-finding method described above with an established method for cell identification, we compared our results to Cell profiler [60]. For background correction, we used the polynomial fit to the ensemble of images, and subtracted the resulting amount from each image. We identified the primary objects under Otsu global threshold method, and used the ‘Shape’ method for defining boundaries between objects and to distinguish the clumped objects. We chose this method because the Cell profiler documentation suggests it as proper to recover round objects in clumps. Using the same method described above for our pipeline, we compared the cells identified by Cell profiler to the manually drawn ellipses. We found that 89. 0% of the manually drawn ellipses have a corresponding identified cell within 10 pixels of their area centre. The mean distance in the paired centers was 2. 23 () and the correlation in object sizes 0. 876. Although these statistics are slightly lower than for our simple methods, it did perform significantly faster, identifying the cells in a typical image in seconds. In addition, 139 artifacts were manually identified. We used this set to compute the false positive rate by pairing automatically identified cell areas to the manually identified cells and artifacts (Suppl. Figure S2). We also computed the false-positive rate as a function of cell probability threshold. For example, filtering all cells that have a cell probability below 0. 8 reduces the false positive rate. This is in agreement with previously reported results using post processing [61]. Since we we found that the number of cells is critical for the robustness of the time profile estimates (Suppl. Figure S3) and that small buds have systematically lower cell probability estimates (Suppl. Figure S1B), we prefered not to choose a hard threshold. Indeed, we found that using a 0. 8 cell probability threshold reduces the robustness of the time profiles (Suppl. Figure S3B). We also found that applying this threshold would discard of the small buds which were used to define the first four of our ten cell-stage time points. We characterize the protein expression phenotype within each cell object using the absolute intensity of the GFP, as well as geometrical distances between proteins to identified points of interest. In both cases, we use the RFP signal to normalize the observations made for the GFP signal. The RFP intensity was found to be dependent on the object size, so we characterized the expected RFP, , and used to normalize the GFP signal by the fold difference to the expectation of the mean RFP intensity (eq. 2). We defined using three linear function segments which fits the mean level of RFP in the 1. 4M automatically identified cells: (16) Some of the morphological distances require us to identify the coordinates of a point of interest; the cell center, protein mass-center and bud neck position are obtained by averaging the coordinates of the cell pixels, of GFP-tagged proteins and Mother-bud separation contour pixels, respectively. Assuming GFP intensities are proportional to protein amount, we derive the expected value for geometrical distances with respect to the position of a randomly selected protein. The position of cell center, protein mass center and bud neck are given by: (17) where is the sum of GFP intensities and ‘’ is the set of contour pixels which separates the bud from the mother cell. The other 2 distances have a slightly different form: first, the distance to the perimeter for any coordinate has been computed using Edge Map distance, so that: (18) Second, we derive the equation for the expected distance between proteins: (19) Once again we use the RFP marker to normalize these distances. In the case of distance between proteins, the distance is normalized by the expected distance between a protein and a RFP marker. For that case, the reported log ratio representing a morphological distance would be: (20) First, we model cell stage as a function of the bud size. Under the assumption that the bud volume increases at a constant rate, we expect that time scales linearly with. Because we have a number of identified cells and distribution of object size that varies throughout the collection of 4004 yeast strains, a common basis is required to enable comparisons between the expression of different proteins. For each strain, time series are defined as expected feature values for objects observed at 10 equidistant cell stage keypoints. We use local regression (LOESS) to infer the mean and variance at each keypoint (eq. 21), where the ‘’ is Gaussian kernel function with bandwidth parameter equal to 1700. In addition, because we have developed a probabilistic cell confidence, which assigns to each identified cell a posterior probability of being a properly identified cell, we use the cell confidence to compute a weighted average, which is the expected profile conditioned on each identified object being drawn from the cell class: (21) where is feature value that is expected at the cell stage keypoint from feature values, which are measured for the identified object. are the quality measures for each shape and are cell sizes for bud objects. Each protein profile is a vector of means and variances of observations. We use the Maximum likelihood clustering criterion [37] (eq. 22) in order to agglorameratively join pairs of protein profiles, proteins to cluster profiles, or pairs of cluster profiles: (22) where is the determinant of a covariance matrix. This criterion is the log likelihood ratio for two cluster of size to have their protein profiles modeled as two multivariate Normal distributions (with their corresponding parameters), to a single multivariate Normal model explaining both expression groups. Initial cluster profiles are build from individual protein profile, which corresponds to 12 concatenated time series of feature values. As such, initial covariance matrices are diagonal matrices whose values were estimated from the LOESS (see in eq. 21). New cluster profiles are characterized by a multivariate normal distribution whose parameters are obtained from merging two previous cluster profiles (eq. 23). (23) where ‘’ are cluster sizes and are normal distribution parameters for merged cluster profiles. The 4004 proteins were grouped based on exact correspondence of subcellular location annotation, as defined by Huh et al. [14] (Suppl. Table S1). 22 classes correspond to unique subcellular locations. We merged member profiles into a class profile () using the operation defined above (eq. 23). The Bhattacharyya metric (eq. 24) was used to compare each class profile (Suppl. Table S1) and class profiles were clustered using euclidean distance (Suppl. Figure S5A). (24) Protein subcellular location was characterized by Huh et al. [14] by assigning one or many annotations to each protein. We report about the enrichment of either separate annotations and/or exact localization set correspondence. For example, proteins that were only annotated to be nuclear obtained the label ‘pure nucleus’ and new labels, such as ‘nucleus AND cytoplasm’ were reported. The GO and PFAM annotation were obtained from Uniprot/SwissProt [51]. In the reported hierarchical clustering results (Figure 6), clusters were manually selected and the hypergeometric distribution was used to model the occurrence of annotations of proteins within them. Bonferroni correction was applied to the P-values (1990 hypotheses, accounting for 3. 3 in log scale). To assess the significance of the hierarchical clustering, we performed permuation tests. For each protein annotation, we find the cluster that yields the smallest P-value for annotation enrichment. We then assess the statistical significance of the sum of the smallest log P-values, ‘S’, by defining two background distributions for ‘S’. In the first, we preserved the structure of the tree, but chose random proteins to assign to each leaf. In the second, we preserved the structure of the tree, but randomly replaced the proteins with other proteins that had exactly the same set of annotations of subcellular localization. In other words, for this second ‘localization constrained’ permutation, we only allow proteins of identical characterization in subcellular location terms to be permuted, so that any enrichments of subcellular location (as displayed in Figure 6) will be preserved for any permutation. We found that the statistic ‘S’ was systematically higher in the 10000 permutations than for original hierarchical cluster. Therefore, we report the corresponding Z-score, but we note that background distribution for ‘S’ is not necessarily a normal distribution (Suppl. Table S3). In order to evaluate the resolution of functional enrichments in the hierarchical clusters, we computed the significance for subsets of annotations. We show in supplementary table S3 that complexes characterized by GO annotation are found significantly enriched, and that the ribosomal and proteasomal proteins, which typically show high protein abundance, have a limited contribution in the sum. In addition, we applied the statistical tests on 14 subsets of GO annotations based on the number of annotated proteins. This analysis was also performed on 5 alternative hierarchical clustering results: This allows us to evaluate the robustness of the results to a change of clustering algorithm (Maximum likelihood clustering, Euclidean metric with complete linkage, Correlation metric with complete linkage), and the usage of the cell confidence (as a weight or using 0. 8 as a filtering threshold) (Suppl. Table S2, S3, S4, S5). Hierarchical clusters are available to be browsed online at: http: //www. moseslab. csb. utoronto. ca/louis-f/unsupervised/. In addition, the source code for the cell identification and feature measurements, the set of 17 images in which 4305 ellipses corresponding to cells and 139 ellipses corresponding to artifacts were manually drawn, as well as a table of feature measurements for all 400 K mother-buds pairs are available.
summary: The location of a particular protein in the cell is one of the most important pieces of information that cell biologists use to understand its function. Fluorescent tags are a powerful way to determine the location of a protein in living cells. Nearly a decade ago, a collection of yeast strains was introduced, where in each strain a single protein was tagged with green fluorescent protein (GFP). Here, we show that by training a computer to accurately identify the buds of growing yeast cells, and then making simple fluorescence measurements in context of cell shape and cell stage, the computer could automatically discover most of the localization patterns (nucleus, cytoplasm, mitochondria, etc.) without any prior knowledge of what the patterns might be. Because we made the same, simple measurements for each yeast cell, we could compare and visualize the patterns of fluorescence for the entire collection of strains. This allowed us to identify large groups of proteins moving around the cell in a coordinated fashion, and to identify new, complex patterns that had previously been difficult to describe.
keywords: systems biology
  model organisms
  yeast and fungal models
  biology
  computational biology
  proteomics
title: Unsupervised Clustering of Subcellular Protein Expression Patterns in High-Throughput Microscopy Images Reveals Protein Complexes and Functional Relationships between Proteins