HeLaM cells (Tiwari et al., 1987) were cultured at 37°C under 5% CO₂, in Dulbecco’s Modified Eagle’s Medium (DMEM) without Arginine, Glutamine, Lysine or Sodium Pyruvate (Gibco, #A14431-01), supplemented with 10% (vol/vol) dialyzed foetal calf serum (PAA, #A11-107), 1 mM Sodium Pyruvate (Sigma, #58636), 1 x GlutaMAX (Gibco, #35050-061) and, either; 42 mg/L ¹³C₆,¹⁵N₄-L-Arginine HCl (Silantes, #201604302) together with 73 mg/L ¹³C₆,¹⁵N₂-L-Lysine HCl (Silantes, #211604302), or 42 mg/L Arginine HCl and 73 mg/L Lysine HCl with standard isotopic constituents (Sigma, #A6969 and #L8662). Cells were passaged >7 times in these media before experiments began.

SILAC labelled HeLa cells were grown in 15 cm cell culture dishes containing 25 mL of medium, to 50–80% confluency. Recombinant Epidermal Growth Factor (Sigma, #E9644) was reconstituted in PBS (Gibco, #14190–094) and added to dishes at a final concentration of 20 ng/mL. Dishes were returned to 37°C for 20 min.

HeLa cells were grown in 35 mm cell culture dishes with embedded glass coverslips (#P35G-1.5-14-C; MatTek, Ashland, US), to approximately 30–50% confluency. The growth medium was exchanged for DMEM supplemented with 5 mM HEPES (Gibco, #15630–056), and 25 nM LysoTracker Red DND-99 (ThermoFisher Scientific, #L-7528, from 1 mM stock in DMSO), and cells were incubated at 37°C for approximately 20 min. Cells were then chilled on ice for several minutes. The growth medium was replaced with pre-chilled DMEM/HEPES supplemented with fluorescently labelled EGF (biotinylated, complexed to Alexa Fluor 488 Streptavidin; ThermoFisher Scientific, #E13345) at 2 μg/ml, as recommended by the manufacturer. Cells were incubated on ice for 30 min to pre-bind the EGF to plasma membrane EGF receptors. The medium was then replaced with cold DMEM/HEPES. One batch of cells (control) was fixed immediately, with 3% formaldehyde in PBS. The other batch was incubated at 37°C, for 30 min, to allow EGF uptake. Cells were fixed with formaldehyde, as above.

Microscopy was performed at the Imaging Facility of Max Planck Institute of Biochemistry, Martinsried, using a ZEISS (Jena, Germany) LSM780 confocal laser scanning microscope equipped with a ZEISS Plan-APO 63x/NA1.46 oil immersion objective. ZEISS Zen software was used to acquire images, and Adobe Photoshop CS6 was used for cropping and global brightness/contrast adjustments.

Protein concentrations were determined using a micro-BCA assay (Thermo, #23225). Protein concentration standards of 0, 25, 125, 250, 500, 750, 1000, and 1500 µg/mL were prepared using a stock of 2000 µg/mL BSA (Thermo, #23209) and SDS buffer. Reagent A (containing Bicinchoninic Acid solution) was mixed 50:1 with reagent B (containing Copper(II)Sulphate), and 200 µL were plated into wells of a 96-well plate. 5 µL of standard solutions and samples of unknown concentration were plated in triplicate. Plates were incubated for 35 min at 37°C before measuring the fluorescence following excitation at 595 nM using a plate reader (Infinite M200, Tecan). Sample concentrations were calculated by comparison to a standard curve.

30 µg of each SILAC light-labelled fraction were independently mixed with an equal amount of SILAC heavy reference organellar fraction. Protein was precipitated by addition of 5 volumes of ice-cold acetone, before incubation at −20°C for 30 min and subsequent centrifugation at 4°C for 5 min at 10,000 g (Centrifuge 5415R, Eppendorf). In case of nuclear, organellar, and cytosol fractions, 60 µg of SILAC heavy labelled sample only were subjected to the same precipitation regime. All subsequent steps were performed at room temperature. Supernatants were removed and pellets allowed to air-dry for 5 min. Pellets were re-suspended in digestion buffer (50 mM Tris pH 8.1, 8 M Urea, 1 mM DTT) and incubated for 15 min. Cysteines were alkylated by addition of 5 mM Iodoacetamide with incubation for 20 min. Proteins were enzymatically digested by addition of 1 µg LysC per 50 µg of protein, with incubation for 3 hr. Digests were then diluted four-fold with 50 mM Tris pH 8.1 before addition of 1 µg Trypsin per 50 µg of protein and overnight incubation.

Peptides were fractionated as previously described (Kulak et al., 2014). Briefly, SDB-RPS (Sigma, #66886-U) stage tips were activated with 100% Acetonitrile, followed by an aqueous solvent containing 30% (v/v) Methanol and 1% (v/v) TFA. Peptide mixtures were acidified with 1% (v/v) TFA, 15 µg (for single shot) or 25 µg (for fractionation) was loaded onto activated stage-tips, and washed with 0.1% (v/v) TFA. Peptides were eluted with 60 µL buffer X (80% Acetonitrile, 5% Ammonium Hydroxide) for single shot analyses (the ‘fast maps’ protocol). For fractionation (the ‘deep maps’ protocol), peptides were eluted using 20 µL SDB-RPSx1 (100 mM Ammonium formate, 40% (v/v) Acetonitrile, 0.5% (v/v) Formic acid), then 20 µL SDB-RPSx2 (150 mM Ammonium formate, 60% (v/v) Acetonitrile 0.5% (v/v) Formic acid), then 30 µL buffer X. Tryptic peptides were dried almost to completion in a centrifugal vacuum concentrator (Concentrator 5301, Eppendorf), and volumes were adjusted to 10 µL with buffer A* (0.1% (v/v) TFA, 2% (v/v) Acetonitrile).

4µL of peptides were loaded on a 50-cm column with 75-µm inner diameter, packed in-house with 1.8-µm C18 particles (Dr Maisch GmbH, Germany). Reverse phase chromatography was performed using the Thermo EASY-nLC 1000 with a binary buffer system consisting of 0.1% formic acid (buffer A) and 80% acetonitrile in 0.1% formic acid (buffer B). The peptides were separated by a linear gradient of buffer B from 2 to 30% in 130 min followed by washout (ramping to 95% B in 5 min, constant at 95% B for 5 min, ramping down to 2% B in 5 min, constant at 2% B for 5 min). The flow rate was 250 nl/min. The column was operated at a constant temperature of 55°C. The LC was coupled to a Q Exactive HF Hybrid Quadrupole-Orbitrap mass spectrometer (Scheltema et al., 2014) (Thermo Fisher Scientific, Germany) via the nanoelectrospray source (Thermo Fisher Scientifc, Germany). MS data were acquired using a data-dependent top-15 method, dynamically choosing the most abundant not-yet-sequenced precursor ions from the survey scans (maximum injection time 25 ms, 300–1650 Th). Sequencing was performed via higher energy collisional dissociation fragmentation with a target value of 1e5 ions determined with predictive automatic gain control. Isolation of precursors was performed with a window of 1.4 Th. Survey scans were acquired at a resolution of 120,000. Resolution for HCD spectra was set to 15,000 with maximum ion injection time of 55 ms. The 'underfill ratio,' specifying the minimum percentage of the target ion value likely to be reached at the maximum fill time, was defined as 20%. Repeat sequencing of peptides was kept to a minimum by dynamic exclusion of the sequenced peptides for 30 s.

In this study, two independent sets of maps were generated:

1. Six replicate maps from untreated HeLa cells, to establish the method and to generate a database of protein subcellular localisation. These maps are referred to as the ‘Static Maps’ here.

2. Three untreated (control) and three EGF-treated (+EGF) maps, as part of a comparative experiment to demonstrate the applicability of the method to detect organellar translocation events. These six maps will be referred to as the ‘Dynamic Maps’ here.

Common procedures are listed without specific reference to either method. In some cases, processing of the Dynamic Maps required extra steps; these are described at the end of each section.

Mass spectrometry raw files were processed in MaxQuant (Cox and Mann, 2008) (v1.5.3.29), using the human SwissProt canonical and isoform protein database, retrieved from UniProt (www.uniprot.org). Raw files were organized into four parameter groups (PGs): double labelled SILAC data (PG1; all organellar subfractions); label-free (LFQ) data (PGs 2, 3, and 4; corresponding to single-labelled nuclear, organellar and cytosolic fractions, kept in separate PGs). Multiplicity was set to 2 for PG1, with Lys8 and Arg10 selected as heavy labels; Re-quantify was enabled for PG1; minimum number of quantification events was set to 1. Multiplicity was set to 1 for PGs2-4, with Lys8 and Arg10 selected as labels. LFQ was enabled for PGs 2–4; minimum number of LFQ peptides was set to 1; matching between runs (within PGs only) was activated, with a match time window size of 0.7 min, and an alignment time window of 20 min. Default parameters were used for all other settings.

The primary output from MaxQuant is the ‘protein groups’ file, which was used as the basis for all subsequent analyses. Identifications were filtered by removing matches to the reverse database, matches only identified by site, and common contaminants. The remaining identifications were split into three groups, for individual processing: Group 1, the SILAC ratio data (ie the subfractions obtained by differential centrifugation); Group 2, the intensity data for the nuclear, organellar and cytosolic fractions; and Group 3, the LFQ data for nuclear, organellar and cytosolic fractions.

Group 1 data were used to derive organellar maps. The raw data for one map are five SILAC ratios (called 3 K, 5.4 K, 12.2 K, 24 K, and 78.4 K, referring to the centrifugation speeds used to derive the corresponding fractions). For each map, we applied a stringency filter to ensure high accuracy of quantification, which increases with the number of peptide quantification events. First, we only considered proteins with complete profiles (ie a set of five SILAC ratios); any proteins with missing values were rejected. Second, we only retained proteins with three or more quantifications in each subfraction. In addition, proteins with only 2 quantification events in one or more subfraction were retained, if their ratio variability for ratios obtained with 2 quantification events was below 30%.

Group 2 data were used to derive protein abundance in the organellar, nuclear and cytosolic (‘ONC’) fractions, as well as overall protein copy numbers and concentrations. The raw data for each ONC triplet were three intensity values. As a stringency filter, we included only those proteins which were identified with at least two MS/MS events (ie without matching between runs) in at least one of the subfractions. All fractions were normalized to the same total intensity. Fractions were then weighted according to protein yields. Since the ONC split was highly reproducible, the same set of weighting factors was used for all maps generated as part of one experiment (0.20 for organellar, 0.32 for nuclear, 0.48 for cytosolic fractions - see Figure 1—figure supplement 1). For each ONC triplet, weighted intensities were summed for each protein, to derive whole-cell relative intensities. These were used to normalize the individual fraction intensities, to obtain for each protein a distribution profile over organellar, nuclear and cytosolic fractions (with the three values from a triplet adding up to one). Furthermore, whole-cell relative intensities were later used to estimate protein copy numbers and concentrations (see below, proteomic ruler). Please note that we chose to perform this analysis on intensity data, not LFQ data, as the three subcellular fractions ONC have very different composition. The LFQ algorithm works optimally on fractions with similar composition, so intensity data are more appropriate here.

Group 3 data were used to investigate changes in organellar, nuclear and cytosolic fractions between experiments. Only proteins with at least one MS/MS event in at least one fraction were kept. Please note that in principle, the Group 2 data could have been used here; however, LFQ offers slightly higher accuracy, includes a normalization step, and is conveniently generated by MaxQuant. As the Group 3 data are solely used for comparing equivalent fractions (ie nuclear with nuclear, cytosolic with cytosolic, organellar with organellar, see below, ‘Global protein distribution changes’), LFQ is the preferred choice here.

Please note that Group 2 and 3 data were derived from (single-label) SILAC heavy fractions. This is a consequence of the chosen labelling strategy (Figure 1A); it would be equally possible to invert the labelling for the entire workflow (ie SILAC light reference, organellar, nuclear and cytosolic fractions; SILAC heavy organellar subfractions), without affecting the experiment or interpretation of the data.

We have implemented and tested two formats for making organellar maps, which differ in the required MS measuring time and achieved depth of coverage.

The Static Maps shown in Figure 2, Figure 2—figure supplements 1 and 2 were prepared using the ‘deep’ protocol; the Dynamic Maps shown in Figure 4 and Figure 4—figure supplement 2 were prepared using the ‘fast’ protocol.

The quality of machine-learning based classification depends strongly on the markers available to train the algorithm. Generation of a reliable organellar marker set is not trivial, as the subcellular localization of many proteins has not been unambiguously determined. Protein localization is often cell type specific, or dependent on the physiological context, further complicating the selection. Moreover, marker sets need to be extensive, to allow accurate cluster delineation. To our knowledge, there is currently no published widely accepted canonical organellar marker set. Hence, we manually curated such a set, optimized for the HeLa cell line used in this study (Supplementary file 1; 1076 proteins). Our selection criteria for a suitable marker were as follows. A) The protein must be expressed in HeLa cells and be identifiable by proteomics. Hence, all markers were picked from the list of proteins identified in the course of this study. B) The protein must have a well-documented and, as far as possible, unambiguous/predominantly unimodal steady state localization. Annotation of subcellular localization was gleaned from the various annotations in UniProt (www.uniprot.org), as well as from the protein’s primary literature. C) The proteins chosen to represent a particular organelle must, more or less, cluster with other markers of this organelle. The function of the markers is to provide a reductionist model of the structure underlying the data; the most important property of a good marker set is hence to provide clusters that are stable over all future replicate maps. To achieve this, we selected suitable markers from several pilot PCA maps (which were not included in this study). We subsequently found their clustering confirmed in the maps presented here. D) Since classification by support vector machines is a ‘boundary’ method, marker proteins near the edge of a cluster are particularly important for cluster discrimination. Hence, once a core set of markers for organelles was established, we specifically augmented the set by identifying suitable candidates (based on criteria a-c) in the cluster periphery. In many cases, we even chose markers that visibly ‘interdigitate’ neighbouring clusters in PCA maps. The aim was to obtain the best boundaries, ie the ones that afford the best predictive value for further proteins, and to minimise ‘unchartered’ space in between clusters. The marker set we provide here should serve as an excellent basis for the generation of maps in other cell types, although some fine-tuning may be required.

The choice of classes was influenced by multiple factors. First, well-established organelles (Plasma membrane, ER, Mitochondria, Golgi, Lysosomes, Endosomes, and Peroxisomes) were included. In addition, the maps themselves suggested deeper resolution in places; thus, known markers of the ERGIC formed a discrete cluster, which also contained many markers of the cis-Golgi. As cis-Golgi is formed from coalescing ERGIC, this is not surprising. Interestingly, the mid- and trans-Golgi markers formed a separate cluster, and hence we decided to define a combined ERGIC/cisGolgi compartment, separate from ER and Golgi. We also observed a small but compact cluster of proteins known to associate with high-curvature ER (REEPs, Atlastins, and Reticulons). This cluster’s map position is distinct from the much larger main ER cluster. Its vicinity to the ERGIC cluster may suggest that this cluster corresponds mostly to tubular ER structures; but since high-curvature ER markers are also found at the edges of sheet ER, we decided to call this cluster ‘ER, high curvature’.

There are several types of endosomal compartments, ranging from ‘early’, ‘late’, to multivesicular bodies (MVBs), and endolysosome (the hybrid compartment produced from fusion of MVBs and lysosomes). As lysosomes are then re-formed from endolysosomes, they share many components. Any definition of endosomal compartments through marker proteins is therefore somewhat arbitrary. In our maps, the endosomal and lysosomal clusters are adjacent. We observe a relatively clear division of established pure lysosomal markers (such as lysosomal hydrolases) from known endosomal markers (eg Retromer, Transferrin receptor). Within the endosomal cluster, typical late endosomal markers (ESCRT-0, vacuolar ATPase subunits) are more concentrated near the lysosome-proximal side. This difference in distribution between late and early endosomes is suggestive, but not to the point where it warrants separate clusters for SVM analysis.

Finally, we observed a clear division of membrane-bound organelles from large protein complexes (such as ribosomes, proteasomes, and many others, Figure 2—figure supplement 1). Many of these complexes are very well resolved. However, since the primary aim of this method is to classify organellar proteins, we considered all large protein complexes as one large cluster. The contents of the large protein complex cluster are further discussed below.

Each organellar map consist of a matrix of protein IDs associated with five SILAC ratio columns (Group1 data, see above) indicating for each protein the relative enrichment and depletion over the five subcellular fractions. For compact visualization of the data, principal component analysis was performed to achieve dimensionality reduction, using SIMCA 14 software (Umetrics/MKS - a free trial version of this software is available from their website: www.umetrics.com). Ratios were first log transformed, then scaled to unit variance. For each map, the first three principle components were calculated. We found that the best visual separation of organellar clusters was achieved by scatter plotting the scores of PC1 vs PC3.

In addition to the six individual maps (with 5 SILAC ratios each) analysed in Figure 2—figure supplement 2, a combined set with proteins common to all six maps (and thus 30 SILAC ratios) was used for generating Figure 2, as a conceptual ‘consensus’ of the six individual maps. All maps can be interactively explored online (www.MapOfTheCell.org). Similarly, the Dynamic maps in Figure 4b show maps combining the common sets from the three control maps (left), and the three EGF treatment maps (right), and are thus based on 15 SILAC ratios each. Individual control and treatment maps (5 SILAC ratios each) are shown in in Figure 4—figure supplement 2.

For the rigorous assignment of proteins to organellar clusters, a support vector machine (SVM) -based machine learning approach was used. Briefly, SVMs are a boundary separation method, which allows non-linear separation of clusters in multidimensional space. First, optimal boundaries between organellar clusters are determined using the marker proteins, with cross-validation to prevent over-fitting. Non-marker proteins falling within the boundaries of a particular cluster are then assigned to that organelle. Distance to the boundary reflects classification confidence (reviewed in Varmuza and Filzmoser [2009]).

Classification was performed on the fractionation profiles (SILAC data, Group 1, see above). Data were prepared as tab delimited text files, containing protein gene names and UniProt identifiers, as well as the five SILAC ratios constituting a single organellar map. Organellar marker proteins were annotated with their corresponding compartment (listed in Supplementary file 1). The analysis was performed in Perseus software (version 1.5.2.11; www.perseus-framework.org). The marker compartments were imported as categorical annotations. Data were first log-transformed. The support vector machine (SVM) algorithm was then trained on the marker proteins, to determine optimal classification parameters. A radial basis function (RBF) was used as the kernel. Parameters Sigma (curvature) and C (misclassification penalty) were optimized by systematic scanning, to achieve the minimal overall classification error. 16-fold cross-validation (CV) was used to prevent over-fitting. Classification was then performed with optimal settings on both marker and non-marker proteins. This time, classification performance for the marker proteins was gauged by full (leave-one-out) cross validation, as a reflection of prediction error for each compartment. For the six replicate Static Maps, the results of the full CV are shown in Supplementary file 2. The average classification performance was 94.7% correctly assigned markers. Optimal values for Sigma ranged from 0.22 to 0.48, and from 8 to 35 for C.

For each protein, the analysis provides a prediction score for each compartment. Positive scores indicate class membership (ie the protein falls within the boundaries of the compartment model derived from the corresponding marker proteins); several class memberships can be predicted for one protein. The class with the highest score is the top prediction (even if the top score is negative, ie the protein is outside the boundaries of the top prediction; in this case, the assignment is a nearest cluster ‘best guess’).

To combine the results of several SVM analyses into one joint prediction, we explored two approaches. First, we analysed the set of proteins common to all six static maps, using all (6 x 5) = 30 SILAC ratios to build an SVM model. The prediction performance on the marker set was better than that of each single map (97.5% correct); however, the depth of this ‘common’ map was substantially lower (3766 proteins vs 4500 in each single map). We therefore employed an alternative approach to combine maps. We trimmed the results from each individual SVM analysis to include only positive scores. Scores were then summed for each compartment and protein across the six different analyses. The top summed score represents a protein’s top cumulative compartment prediction. The advantages of this approach are two-fold: firstly, the full information from all replicate maps is used, not just a reduced common set; and secondly, it focusses on clear positive classification decisions. Thus, if two clusters are well separated in one map, but poorly separated in another, only the information from the better resolved dataset will be propagated. In addition, the procedure takes into account the consistency and number of times with which a protein is classified (eg a protein predicted with a medium score in six replicate maps achieves a higher cumulative score than a protein with only a single strong classification). The prediction performance of this approach is very good (96.5% correctly predicted), clearly exceeding the performance of each of the six individual maps. Furthermore, while it has similar performance to the map combining the proteins common to all six maps, it increased the number of covered proteins by 40% (from 3766 to 5265). Therefore, the results of this ‘additive’ analysis were used for construction of the HeLa spatial proteome database (Supplementary file 4).

As the numerical values of SVM scores are not immediately interpretable, we scaled the scores to percentiles. We ranked the cumulative SVM scores from all 1038 correctly predicted organellar marker proteins. Each score was then converted into a percent-rank. The SVM scores of non-marker proteins were converted into percentile scores, using this scale. For example, a protein with percentile prediction score of 30 has a cumulative SVM score that is greater than 30% of the correctly classified marker proteins.

To further put the scores into perspective, we investigated map concordance (ie prediction agreement between replicate maps, see below) as a function of scores. When the predictions from two individual maps are compared, the likelihood of concordance increases with prediction score. We defined prediction confidence classes accordingly:

For example, for proteins classified in two maps with a score > 20, both maps will make the same predictions for 99% of them. The concordance values were taken from the pairwise comparison of predictions from the six static maps; we chose the values derived from the inter-day comparisons as the most pessimistic estimate (see below, and Figure 2—figure supplement 2C).

SVMs were also used to predict subcompartments for ER (lumen and membrane) and mitochondria (matrix, outer and inner membrane). Starting with the combined set of predictions from the six static maps, all marker proteins of ER and Mitochondria which had been correctly predicted were further annotated for subcompartment localization (provided by UniProt (www.uniprot.org) and primary literature). SVMs were then applied separately to ER and Mitochondrial markers to define boundaries for the subcompartments. After optimization of Sigma and C as above, full cross-validation was used to judge prediction accuracy. Combination of the output from the six static maps showed the following levels of prediction performance for subcompartments:

SVMs were then applied to the predicted non-marker ER and mitochondrial proteins, to provide subcompartment predictions.

Processing of the six Dynamic Maps was performed likewise, but without subcompartment annotation. In addition, control maps and +EGF maps were treated separately, and combined into two prediction outputs, ie organellar predictions before and after EGF treatment. Supplementary file 2 lists the results of the prediction accuracy with full cross-validation. Since the Dynamic Maps were prepared with the ‘fast’ mass spectrometry protocol, which requires only one third of the measuring time of the ‘deep’ protocol used for the Static Maps, their performance is slightly lower.

There is no universally suitable and accepted set of organellar marker proteins (see above, Organellar marker set generation). The true prediction accuracy of the classification achieved in this study is hence difficult to gauge. Nevertheless, prediction accuracy of the training set, estimated with leave-one-out cross validation, can be used as a proxy. Thus, for a set of 1000 marker proteins, support-vector machine (SVM) models are built and predictions made 1000 times, using a different combination of 999 marker proteins each time. This essentially simulates the prediction of proteins that are not in the marker set, and ensures that the SVM models have maximum predictive value for non-maker proteins. The marker protein prediction accuracy hence provides an estimate of the overall performance for those new predictions that fall within the boundaries defined by the markers. External validation of our predictions supports this notion: the overall concordance with the recent LOPIT study (Christoforou et al., 2016) is 91.6%; the concordance with experimentally validated mitochondrial proteins (Calvo et al., 2016) is 97%, and hence close to the prediction accuracy of our mitochondrial marker set (98.8%, Table 1). Prediction confidence varies between individual proteins (see above, Organellar maps through Support Vector Machine classification), and also varies between organelles. As a conservative estimate of global prediction accuracy achieved in this study, we consider the average accuracy per membrane bound organelle, 92.7% (Table 1).

To estimate the reproducibility of our organellar maps, we determined map concordance, defined as the percentage of identical predictions between two independent map replicates. We compiled the individual top predictions and SVM scores for the six replicate Static Maps into a table (with 12 classified compartments). We then performed 15 pairwise comparisons (all combinations of the six maps), including 12 inter-day, and 3 intra-day comparisons (maps 1 vs 2, 3 vs 4, and 5 vs 6). Concordance was calculated as follows. First, the set of proteins present in both maps was determined. Next, the proportion of proteins with identical top compartment predictions was calculated, corresponding to the map concordance. The average concordance across all inter- and intra-day comparisons was then determined. As concordance is a function of prediction confidence, we next introduced SVM scores as a further criterion. Scores for individual maps were converted to approximate single map percentile scores, by dividing the summed percentile scores (see above) by 6. Average map concordances were then calculated after filtering for predictions present in both maps with a given minimum score.

Furthermore, we compared the combined output from the six Static Maps prepared with the ‘deep’ mass spectrometry protocol (60 hr) with the combined output from the three controls of the Dynamics maps, which had been prepared with the ‘fast’ protocol (20 hr). Map concordance between these two sets was calculated as a function of prediction scores. Concordance was very similar to that observed for intra-experimental comparisons of single static maps (Figure 3), at 98% with a score cut-off at 7, demonstrating a high level of agreement between maps generated with the ‘deep’ and ‘fast’ protocols.

The intensity data for nuclear, organellar and cytosolic fractions (Group 2 data, see above) were used to evaluate a protein’s global cellular distribution. It is very important to separate this analysis from the organellar assignment achieved via the SILAC based maps. Many proteins have a cytosolic pool as well as a (peripheral) membrane associated pool. If the cytosolic fraction was added to the organellar maps, it would show these proteins more or less near the cluster of purely cytosolic proteins, but not reveal the organellar association of their membrane associated pool. Hence, we decided on a two-tiered classification, which first quantifies a protein’s global distribution, and then reveals its membrane-bound compartment. The same logic applies for the numerous proteins that shuttle between nucleus and cytosol. Our intensity analysis identifies which proteins are predominantly nuclear or predominantly cytosolic, and which proteins are likely to shuttle between the two (pronounced nuclear and cytosolic pools). Again, inclusion of these data with the organellar maps would not allow the identification of nuclear/cytosolic shuttling proteins.

The intensity analysis shows for each protein a distribution across nuclear, cytosolic and organellar fractions, normalized to 1; profiles from different proteins are thus directly comparable. We broadly classified proteins based on the relative contributions of each pool, using useful arbitrary cut-offs:

This classification scheme is also incorporated into the HeLa spatial proteome database (Supplementary file 4).

Absolute protein copy numbers were calculated using the ‘Proteomic Ruler’ approach (Wisniewski et al., 2014). Briefly, this method quantifies the mass spectrometric signal of histones, which is proportional to cellular DNA content. Knowledge of the cell’s chromosomal ploidy in turn allows the conversion of the histone signal into histone protein mass per cell. The relative intensities of all proteins can then be converted into copy numbers per cell, as well as cellular concentrations. The proteomic ruler approach is remarkably accurate; comparison to ‘gold standard’ quantification with isotope labelled peptides (PrEST approach [Zeiler et al., 2012]) shows similar precision, and typical absolute deviations of less than two-fold (Wisniewski et al., 2014).

The proteomic ruler has been implemented as a convenient plug-in for the Perseus platform (www.perseus-framework.org). To derive protein copy numbers, whole-cell protein intensities (ie pooled and weighted Group 2 data, see above) were imported into Perseus 1.5.1.16. For the Static Maps, we analysed six biological replicates. For the Dynamic Maps, three replicates for control and three for EGF treatment were analysed (separately from the Static Maps). Signals were normalized to average molecular mass. Ploidy was set to three (a reasonable estimate for HeLa cells –the exact ploidy is not known and may vary). Normalization was the same across all columns. A total cellular protein content of 200 g/l was assumed (default). Medians and standard deviations were calculated from repeat experiments to provide copy number, concentrations, and error estimates.

To assess the quality of our analysis, we first compared the prediction consistency from the individual Static Maps (six biological replicates). The average Pearson correlation of copy numbers between replicates was 0.98 (calculated from 15 pairwise comparisons, log data), demonstrating very high precision (ie prediction reproducibility). To gauge prediction accuracy, we compared median copy numbers predicted from our maps with absolute copy numbers of 27 proteins derived in an independent study using the PrEST approach (Zeiler et al., 2012). Pearson correlation of the two datasets (in log space) was 0.96, and thus remarkably high, showing a near-linear relationship of our data with the PrEST data. The median ratio of copy numbers (static maps/PrEST) was 2.3, suggesting a small systematic shift in our predictions relative to the PrEST data. The two studies used different strains of HeLa cells (which may differ in cell size), and this may account for the shift. In any case, the proposed overall accuracy of the proteomic ruler is also in that range (Wisniewski et al., 2014). In sum, the copy number data obtained here show excellent internal consistency, and are in good agreement with independently obtained data.

The proteomic ruler approach was also used to estimate copy number changes in the Dynamic EGF treatment experiments. Copy numbers were log-transformed, and medians were calculated from the three control and the three EGF treatment maps. In addition, for each protein the copy number difference between control and treatment maps was analysed with a paired t-test (two-tailed). A p-value <0.01 was considered as significant (no correction for multiple hypothesis testing was applied, to maintain statistical power; in this context, false positives are not too detrimental, as the actual movement of proteins between compartments is assessed in separate statistical tests (with full FDR control). For proteins with no significant change, the overall median from all six maps (three controls and three treatment maps) was used as the basis for further calculations. For proteins with significant changes, the copy numbers for control and treatment were kept separate for subsequent calculations.

Copy numbers were then multiplied with the proportions of nuclear, organellar and cytosolic fractions (see above, intensity analysis). Average differences between corresponding fractions before and after treatment were calculated. For proteins with insignificant overall copy number changes, all compartment-specific changes balance to a net overall change of zero.

Organellar integrity is an important concern for cell lysis conditions. We assessed preparation-induced leakage of mitochondria, lysosomes, and ER, by analysing the cytosolic pool of known lumenal marker proteins of these organelles (see Supplementary file 1) for the six Static Maps (Figure 1—figure supplement 1). The underlying assumption is that these proteins should have almost no cytosolic pool if organelles stayed intact during cell lysis; conversely, any cytosolic appearance of markers should indicate the magnitude of leakage. In all three cases, we observed a large proportion of proteins with no cytosolic pool (<1%); a distribution of proteins with small cytosolic pools (1–5% for ER and mitochondria, and 1–15% for ER); and a small number of proteins with large cytosolic pools (>20%). We interpret these three groups as follows. 1) Some lumenal proteins are attached to the membrane of the compartment, and will not leak if the compartment is damaged. 2) Genuine protein leakage will vary from protein to protein, possibly with size, solubility, binding partners, etc. 3) Some proteins have splice variants with cytosolic localization, giving rise to apparent large cytosolic pools even in the absence of leakage (and this is in fact documented for many of the proteins we observe in group 3). Only the second group of proteins is hence likely to reflect genuine leakage. We therefore calculated the averages from proteins with cytosolic pools between 1–20% only, indicating leakage of 3.9% for mitochondria, 2.3% for lysosomes, and 8.3% for ER. These values suggest a very high level of organellar integrity achieved by our lysis procedure. The slightly higher value for ER is also expected, as the reticular morphology and interconnectedness with the nuclear membrane makes it impossible to extract this organelle intact.

The HeLa spatial proteome database (Supplementary file 4) is an interactive Excel-file that allows users to enter query proteins and retrieve predictions for subcellular localization, relative protein abundance, absolute copy number and concentration, as well as the distribution over organellar, cytosolic and nuclear fractions. In addition, a ‘neighbourhood analysis’ is performed, revealing which proteins share the query protein’s fractionation behaviour. Using the database is quick and easy; you only need to enter the gene name or UniProt ID of the protein(s) of interest. The database automatically retrieves all associated information. Multiple simultaneous queries can also be submitted as a list of gene names/UniProt IDs (‘batch submission). The database includes a quick start user’s guide, as well as guidelines for interpreting the output.

The database construction is straightforward. Mostly, it retrieves information from a central data table, prepared and compiled as described in the various sections of the Methods. An exception is the neighbourhood analysis, which is calculated within the database. Briefly, the squared Euclidian distance to the query protein’s profile is calculated for every protein in the database. Proteins are then sorted by ascending distance to the query. Those with low distances have fractionation profiles similar to the query. In many cases, this identifies members of stable protein complexes, but may also reveal other important/interesting connections. Useful distance cut-offs to guide the interpretation of profile similarity were derived by analysing the distances of between members of known stable complexes, as described (Borner et al., 2014).

The Dynamic EGF spatial proteome database (Supplementary file 6) is similar to Supplementary file 4, but includes information on both control and EGF treatment maps. As an additional feature, it highlights differences in protein compartment prediction, distribution, as well as copy number, before and after EGF treatment.

Calculating the mass of each subcellular compartment required the ‘Intensity data analysis of global protein distribution’ (Group 2 data), copy number estimations, and the organellar classifications. Initially the median copy numbers were multiplied with the MW (kDA) to get the mass of each protein in the cell. These masses were then multiplied by the average proportion of the intensity from the three fractions; Organellar, Nuclear and Cytosolic. This gives a mass for each protein in each of these three compartments. However, the contamination of the nuclear fraction with non-nuclear organellar material would lead to incorrect assignment of mass to the nucleus. We therefore used the 6 classifications described in ‘Intensity data analysis of global protein distribution’ to redistribute this mis-assigned nuclear mass as follows:

For example, for a protein that is mostly organellar, its calculated cytosolic pool will be as determined; its organellar pool will be the sum of its nuclear and organellar pools.

Several assumptions were made in order to define these re-distributions.

The first assumption is that everything that is found in the cytosol genuinely has a cytosolic pool. Although there will be some protein amount coming from the nucleus and leakage from organelles, it is not possible to determine which nuclear proteins genuinely have a cytosolic pool and which are contaminants due to leakage, therefore the mass of the cytosol is calculated as everything detected in the cytosolic fraction.

The second assumption is that nuclear proteins found in the organellar fraction are contaminants from the fractionation procedure. Not all nuclear material will pellet at 1 K, especially broken nuclei. Therefore the mass of the organellar fraction of proteins defined as nuclear is given to the nucleus. In addition, Nuclear/Cytosolic proteins have their organellar pool given to the nucleus for the same reason.

Conversely, organellar proteins found in the nuclear fraction are also deemed contaminants, this is particularly important since many organelles are not efficiently released upon cell lysis, due to attachment to the cytoskeleton, cells remaining partially intact, incomplete stripping of membranes that are physically attached to the nucleus, such as the ER. Therefore any protein that is classified as organellar is assigned both the organellar and nuclear pool of that protein. Likewise, for organellar/cytosolic proteins.

Proteins that are broadly distributed are assumed to be genuinely broadly distributed and are assigned to the nucleus and organellar fraction in the proportions determined in the intensity analysis.

Proteins classed as cytosolic have their organellar or nuclear mass assigned to the nucleus or the organelle fractions.

The mass of the organellar fraction can then be divided to the respective organelles using the organellar classification determined by SVM analysis of the organellar maps. From this classification, large protein complex proteins were assigned to cytosol, because these mostly pellet due to their mass, and not necessarily because of membrane association, as shown by their enrichment in the 80K fraction, where all membranous organelles are de-enriched.

Any protein assigned to Actin-binding protein cluster in the organellar maps, regardless of its split among organellar, nucleus or cytosol, is considered as part of the cytoskeleton. The cytosolic pool of actin binding proteins may truly reflect a soluble pool, hence the mass of Organellar and Nuclear pools constitute the actin cytoskeleton and the Cytosolic pool represents the soluble actin cytoskeleton components.

Translocation events are detected at two levels: 1) changes in organellar association; 2) changes in global profile distribution (organellar, nuclear, cytosolic). Each analysis is also a two-step procedure: First, a statistical test is used to identify reproducible translocations. Next, the nature of the translocation is evaluated in the context of all available spatial information (eg comparing organellar predictions, PCA maps, cytosolic vs membrane distributions, copy number changes, etc.). The rationale for decoupling detection and interpretation of movement is explained at the end of this section.

We have developed a statistical test for the detection of protein movement between organellar maps. The test combines two simple metrics, protein translocation magnitude (M) and translocation reproducibility (R), which can be combined into an easily interpreted ‘MR Plot’. The test is optimally performed on three biological repeats of control and treatment each, as demonstrated in this study. It is still useful for analyses with two repeats, but less stringent. If a comparative experiment is performed only once (eg as a pilot experiment), only the magnitude of translocation metric can be calculated. This still indicates all relevant changes, but with a much higher background (ie false positive rate). The test can conceptually be divided into several data transformation steps, and the subsequent calculation of the two metrics. The following is a step-by-step guide, assuming a simple comparative experiment has been performed in triplicate (three treated and three control maps).

Once significant outliers have been detected in the MR plots, the underlying translocation events may be evaluated through PCA based maps, SVM organellar predictions, or nearest neighbour analysis. We recommend a combination of all three approaches, as the information gleaned is complementary.

PCA maps show the movement graphically, and are a good starting point for orientation; they reveal the approximate compartment before and after the translocation, and also indicate if the translocation event is likely to be only partial. Maps also show the relative position of a protein within a cluster (central vs near the edge; in the latter case, it may be of note which other cluster is close). As an example from this study, the EGFR is very central to the plasma membrane cluster in untreated cells, but is on the edge of the endosomal cluster after EGF treatment, in close proximity to the lysosomal cluster, hinting perhaps at a late endosomal compartment, or beginning lysosomal localization. (Please note that the endolysosomal compartment may be seen as a continuous distribution of compartments, with early endosomes and lysosomes marking the extreme end points; the distinction between late endosomes, multivesicular bodies, and fused endolysosomes is somewhat fluid anyway). It is also advisable to check the movements on all available replicate maps (eg Figure 4—figure supplement 2).

Next, the organellar association as predicted by SVMs before and after treatment can be queried. To illustrate this, we have prepared an interactive database for the dynamic maps (control vs EGF treatment, Supplementary file 6). Proteins may undergo complete organellar translocations, as in the case of EGFR receptor; however, it is also possible that only a small proportion of the protein moves, and this may not be sufficient to change the overall SVM predictions. In other words, the MR plot is potentially more sensitive in detecting movement than a ‘before-after’ comparison of SVM predictions. Nevertheless, in cases of complete translocations, the SVM reveal the old and new compartments.

A further potentially revealing aspect is a protein’s ‘neighbourhood’, ie proteins with similar fractionation profiles. In case of genuine translocations, the neighbourhood should change almost completely. The database shows the identity of the old and new neighbours, before and after EGF treatment, and this can be informative. For example, GRB2 as the query protein reveals SHC1 as a close neighbour after EGF treatment (but not before), consistent with their known complex formation triggered by EGFR activation.

A central purpose of the organellar maps is to detect induced changes in protein subcellular localization. The global distribution profiles greatly add to this, as many translocations do not include two organellar compartments, but changes in organellar/cytosolic balance, or nuclear/cytosolic distribution. To determine any such global changes upon EGF treatment, we used the intensity data processed with the label free quantification (LFQ) algorithm (Group 3 data, see above). Data analysis was performed in Perseus (1.5.2.11; www.perseus-framework.org), separately for each compartment (organellar, nuclear cytosolic fractions). LFQ intensities were grouped (triplicate control vs triplicate EGF treatment experiments), log transformed, and filtered for a minimum of three valid values in at least one group. Missing values were then imputed from a normal distribution (downshift =2, width = 0.3). A (two-tailed) Student’s t-test was then performed on the difference between control and EGF treatment for each protein. To identify stringent cut-offs for significant changes, we required a minimum 2-fold change in abundance, where the protein must constitute at least 10% of the total pool, either before or after EGF stimulation, in the compartment where it is shown to be changing (the total pool distribution was obtained from the intensity (Group 2) data, see above). An FDR controlled p-value cut-off was determined by treating the six ‘Static’ control maps as a mock comparative experiment, where three maps were assigned as mock-treated and three as control. The LFQ-analysis was then performed as for the control and EGF treatment maps. For each compartment, a p-value cut-off was determined such that no false positives would be detected in the mock experiment, but changes could still be detected in the EGF experiment (corresponding to an FDR of 0). This was possible for cytosol and organellar fractions; in the case of the nuclear fraction, with optimal cut-offs, 2 false positives are expected among the 13 detected significant changes in Figure 4—figure supplement 2 (FDR=15%).

The LFQ-based statistical analysis of abundance changes in the nuclear, cytosolic and organellar fractions is very sensitive. However, it only reveals a fold change in fraction association, which then needs to be placed in the context of the protein’s global distribution. (For example: The LFQ analysis shows that the cytosolic pool of GRB2 is decreased 2 fold upon EGF treatment (Figure 4—figure supplement 2); however, only the global profile reveals that this corresponds to a drop in the cytosolic pool from 90% to 40%, Figure 4—figure supplement 3). For the EGF experiment, all profile changes can be viewed in the interactive Dynamic database (Supplementary file 6). In most cases, the interpretation of changes in global profiles is straightforward (for example, a cytosolic depletion with concomitant increase in the organellar fraction, as for GRB2; or a nuclear exit with concomitant cytosolic increase, as for NR4A3). The exception is an intra-organellar translocation event involving organelles with very different levels of enrichment in the organellar fraction (see above, “Interpretation of global intensity distribution profiles and their changes”). A good example is the EGFR. Endosomal proteins have a much higher organellar enrichment in the organellar fraction than plasma membrane proteins; a translocation from the plasma membrane to endosomes is thus accompanied by a highly characteristic shift from the nuclear to the organellar fraction (Figure 4—figure supplement 3). This case illustrates the importance of combining all levels of available data to interpret the observed changes (see also next paragraph).

Combining the results from the organellar and global translocation detection may reveal cognate changes that strengthen the individual observations, and allow further interpretation. For example, the EGF-induced translocation of GRB2 to the endo-lysosomal system is accompanied by strong recruitment to the organellar fraction, and depletion from the cytosol; the number of GRB2 copies recruited to the organellar fraction is similar to the overall number of EGFR receptor proteins in the cell. In combination these observations are consistent with recruitment of GRB2 to EGFR, with approximately 1:1 stoichiometry, and movement to endosomes/lysosomes. In this way, a systematic, integrated analysis of changes can be performed (Figure 5).

The two-step procedure described above has several advantages over direct detection of translocation events, for example by comparing SVM organellar predictions before and after treatment. Firstly, our outlier test is much more sensitive and accurate, as it works on a lower data level, the quantitative mass spectrometry data (rather than on organellar predictions based on these data). SVM predictions, although very accurate and highly concordant between maps, necessarily introduce more noise into the data. This will lead to more false positives, especially for proteins with poor predictions in both sets, and proteins whose prediction may sway between two compartments. In contrast, with the outlier test presented here, it was possible to identify translocation events with an experimentally determined FDR of 0 for the EGF experiment. Furthermore, detectable changes in organellar predictions are restricted to complete translocation events (ie a protein fully translocating from one compartment to another). The outlier test, however, will also identify partial translocation events, ie small but consistent movements on the map, which do not lead to changes in the predicted compartment. Similarly, the SVM-based analysis is limited to proteins which have clear organellar predictions before and/or after treatment; proteins that are not assigned to an organelle cannot be identified as changing. In contrast, the outlier test will identify these changes, even if they cannot be fully interpreted afterwards.

The HeLa spatial proteome database (Supplementary file 4 interactive; Supplementary file 1 compact) and the interactive webpage (www.MapOfTheCell.org) provide the classification analysis in an easily accessible format. For the best interpretation of the predictions, we recommend that users consider all types of information available for a protein.

1. The highest-level prediction is the organellar classification. The prediction score provides guidance as to the confidence of the prediction; especially low scoring predictions warrant further investigation. Low scores can occur for at least three reasons: First, the protein may not be very abundant, resulting in fewer profiles, or lower quantification accuracy. Second, the protein may not have a clear single steady state organellar association, but multiple equally dominant associations, resulting in a mixed profile that does not closely match any compartment. Third, a protein may not at all be associated with any of the compartments mapped here.

2. The PCA-based maps (as shown on the webpage) provide complementary spatial information, and show where within a cluster the protein is located. For low-scoring proteins near the edge of a cluster, or proteins with an unexpected localization prediction, it may be useful to consider what other clusters are close by. In case the protein was mis-classified, the true association is most likely one of the neighbouring clusters. Similarly, some compartments are closely related, or may even be considered as ‘snapshots’ along a continuum (eg endosomes, late endosomes, and lysosomes; ER, ER-Golgi-intermediate compartment, and Golgi). Some proteins may be present in two ‘adjacent’ compartments, assume an intermediate position, and therefore achieve only low prediction scores for either one. Again, the PCA maps may reveal if this is likely to be the case. As we provide six replicate maps, we recommend comparing a protein’s position across all maps, to check for variability.

3. Furthermore, the neighbourhood analysis may provide very useful additional information. First, it shows which proteins (if any) have so closely related profiles that they may be forming a complex with the query; in this case, the organellar predictions for these proteins may also be of interest. Second, it provides an overview of the local map environment of the protein; if many of these proteins are established markers of the query protein’s predicted class, it further supports its classification.

4. Finally, we advise users to consider the information on abundance and global protein distribution. Especially very abundant cytosolic proteins may be associated with various intracellular functions and locations (specific as well as non-specific), which may lead to broad (and difficult to interpret) profiles. Large nuclear and cytosolic pools suggest that the organellar fraction of these proteins is only minor, and the main function lies elsewhere in the cell.

The organellar fraction, which is the basis for compartment assignment, contains numerous abundant non-membranous protein particles, such as ribosomes, proteasomes, the TRiC chaperone, the tripeptidyl peptidase (TPP2), and many more. These are sufficiently large and dense to pellet even at the relatively low speeds used to pellet membrane-bound organelles. This group is necessarily very heterogeneous, as the profiling behaviour of complexes is diverse, depending on size, density, association with other structures, etc. Nevertheless, as our PCA maps reveal, the group forms a cluster that is clearly distinguishable from membrane-bound compartments; we have named this the Large Protein Complex (LPC) cluster. Consistent with the notion that the LPCs contain no membranous compartments, there are almost no transmembrane domain proteins in this cluster (Figure 2—figure supplement 1). Within the cluster, individual protein complexes are clearly resolved in many cases (with some examples shown in Figure 2—figure supplement 1). For three reasons, we deliberately did not include any of these complexes in the SVM-based compartment predictions. First, many complexes have only few members, and provide insufficient cluster size for SVM classification. Second, the composition of many complexes is already established; any ‘new’ assignments are thus likely to have a very high FDR. Third, our HeLa spatial proteome database (Supplementary file 4) includes a ‘neighbourhood analysis’, which reveals which proteins are fractionating most closely with a query; this very convenient tool can (and has been) used to identify stable protein complexes, as well as their approximate stoichiometry (Borner et al., 2014). So although the LPC cluster is resolved into smaller clusters, we de-coupled this analysis from the SVM based compartment prediction.

Many of the complexes in the LPC cluster can transiently interact with membranes, including for example vesicle coat proteins such as clathrin, COPI, COPII. Organellar maps through fractionation profiling reveal only the majority steady state distribution of a protein; hence, if a complex is predominantly cytosolic, and only a relatively small proportion membrane-bound, the complex will be in the LPC class. Ribosomes are another case in point. Based on quantitative electron microscopy (Attardi et al., 1969), only 15–20% of ribosomes are attached to the ER membranes in HeLa cells, the remainder are ‘free’. Consistent with these data, our analysis shows ribosomes in the LPC cluster, and not in the ER cluster. Furthermore, proteins may be affiliated with more than one complex, and show an ‘in between’ distribution that is not directly interpretable (for example SEC13, which has a known function in the nuclear pore complex, and in the COPII coat). We also noticed the presence of many nuclear proteins in the LPC cluster (eg histones). These proteins tend to have large nuclear pools in the membrane/nucleus/cytosol split, indicating that they are indeed predominantly nuclear in HeLa cells. Some of these proteins may well have a genuine small cytosolic pool; in addition, their presence may be due to a limited amount of nuclear breakage. We also noticed that many LPC proteins are RNA binding proteins, and it is very likely that they pellet as part of large assemblies on RNA. Finally, there are some proteins in the LPC class that are not known to be clearly associated with a particular complex, and the reasons for their presence in the LPC are currently unclear; a potential explanation is that some (abundant) soluble proteins may non-specifically stick to membranes or other complexes, giving rise to ‘rogue’ profiles. However, non-specific aggregation should not lead to consistent predictions of protein ‘neighbourhoods’, and we suggest that proteins with close neighbours are more likely to be engaged in specific complexes. The bottom line is that any protein in the LPC fraction must be associated with a particle/structure large and dense enough to pellet at ~80,000 g (the maximum speed of the fractionation protocol).

• Georg HH Borner

• Daniel N Itzhak

• Daniel N Itzhak
• Stefka Tyanova
• Jürgen Cox
• Georg HH Borner