The overall work-flow of SCoPE2 is illustrated in Fig. 1a. Single cells are isolated in individual wells, lysed, and their proteins digested to peptides. The peptides from each single cell are covalently labeled (barcoded) with isobaric tandem-mass-tags (TMT), and therefore, labeled peptides with the same sequence (and thus mass) appear as a single mass/charge cluster in the MS1 scans. The MS instrument isolates such clusters and fragments them, see Fig. 2. In addition to generating peptide fragments which facilitate peptide identification, fragmentation generates reporter ions (RI), whose abundances reflect protein abundances in the corresponding samples (single cells) (Fig. 1a). Below, we summarize the key advances of SCoPE2 over SCoPE-MS.
Automated and miniaturized cell lysis
Instead of lysing cells by focused acoustic sonication, SCoPE2 lyses cells by Minimal ProteOmic sample Preparation (mPOP) [24]. mPOP uses a freeze-heat cycle that extracts proteins efficiently in pure water, thus obviating cleanup before MS analysis. mPOP allows sample preparation in multiwell plates, which enables simultaneous processing of many samples in parallel with inexpensive PCR thermocyclers and liquid dispensers. This advance over SCoPE-MS allows SCoPE2 to decrease lysis volumes 10-fold, from 10 to 1 μl, to reduce the cost of consumables and equipment over 100-fold, and to increase the throughput of sample preparation over 100-fold by parallel processing.
Improved data integration across single cells
SCoPE2 also introduces a reference channel composed of a reference sample used in all sets. While commonly used in bulk proteomic experiments to integrate quantitative proteomic measurements across multiple multiplexed experiments, this concept was not implemented in SCoPE-MS. In SCoPE2, the reference is about 5-fold more abundant than a single-cell proteome so that the higher abundance results in improved ion-counting statistics while remaining comparable to that of single cells, and thus likely to be within the linear range of quantification. The reference channel for the experiments involving monocytes and macrophages is composed of an equal number of cell equivalents of monocyte and macrophage cell material (2.5 cell equivalents each) from cells counted by hemocytometer, mixed, and prepared in bulk.
Increased throughput and improved quantification
SCoPE2 introduces a shorter nLC gradient, which allows analyzing more cells per unit of time. Furthermore, improved peptide separation and a narrower isolation window (0.7 Th) result in improved ion isolation and thus improved quantification, see Fig. 2 and Additional file 1: Fig. S5.
Systematic parameter optimization
The sensitivity and quantification of LC-MS/MS experiments depend on numerous instrument parameters whose tuning often depends on trial-and-error approaches [25]. For example, tuned parameters may increase the number of chromatographic peaks sampled at their apexes, thus increasing the copies of peptide ions used for quantification, see Fig. 2. To make such optimizations more systematic, SCoPE2 employs DO-MS to interactively display data required for the rational optimization of instrument parameters and performance.
Enhanced peptide sequence identification
Once the MS data are acquired, SCoPE2 can use additional features of the data to enhance their interpretation. Specifically, SCoPE2 uses a Bayesian framework (DART-ID) to incorporate retention time information for increasing the confidence of assigning peptide sequences to MS spectra while rigorously estimating the false discovery rate [26].
Optimizing and evaluating MS analysis with standards modeling SCoPE2 sets
The quality of LC-MS/MS data strongly depends on numerous interdependent parameters (e.g., chromatographic packing, LC gradient steepness, and ion accumulation time), see Fig. 2. To optimize such parameters, we applied DO-MS on replicate samples, termed “master standards.” Each 1xM standard is a 1% injection of a bulk sample (100xM) composed of eight TMT-labeled (barcoded) samples: two 5000-cell carrier samples (one for each cell type) and six 100-cell samples (three for each cell type), as shown in Fig. 1b. Each 1-μl injection of this bulk sample constitutes a 1xM standard which contains peptide input equivalent to 50 cells in each carrier channel, with the remaining six channels each containing peptide input equivalent to a single cell. Thus, the 1xM standards approximate idealized SCoPE2 standards that enabled us to focus on optimizing LC-MS/MS parameters using identical samples, i.e., independently of the biological variability between single cells.
First, we optimized our analytical column configuration and LC gradient settings. Each 1xM sample was analyzed over a 60-min active gradient since our goal was to optimize the number of proteins quantified across many cells, rather than merely the number of proteins quantified per sample [6]. By varying chromatographic parameters and benchmarking their effects with DO-MS, we minimized elution peak widths. Sharper elution peaks increase the sampled copies of each peptide per unit time and reduce the probability that multiple peptides are simultaneously isolated for MS2 analysis, see Fig. 2. Concurrent with optimizing peptide elution profiles, we optimized the data-dependent acquisition MS settings, such as minimum MS2 intensity threshold, MS2 injection time, and the number of ions sent for MS2 analysis per duty cycle (i.e., TopN), to increase the probability of sampling the apex of the elution peak of each peptide. This optimization increased the number of ion copies sampled from each peptide [25].
The 1xM standards also permitted estimating the instrument measurement noise—independently of biological and sample preparation noise—in the context of diluted bulk standards modeling SCoPE2 sets. This noise estimate was motivated by our concern that factors unique to ultra-low abundance samples, such as counting noise [1, 4, 6, 29], may undermine measurement accuracy. To estimate the measurement reproducibility, we correlated the RI intensities of lysates corresponding to the same cell type: an average correlation of 0.98 suggested good reproducibility. However, this reproducibility benchmark cannot evaluate the accuracy of relative protein quantification, i.e., the ability of SCoPE2 to quantify changes of proteins across cell types [8, 30, 31]. To benchmark accuracy, we compared the fold changes of proteins between human embryonic kidney cells and monocytes (HEK-293/U-937 protein ratios) estimated from two cell lysates diluted to single-cell levels against the corresponding ratios estimated from the bulk samples used as carriers (Fig. 1b). The high concordance of these estimates (Spearman ρ = 0.84) strongly indicates that the instrument noise in quantifying proteins from standards with isobaric carriers is small, consistent with our arguments that the abundance of proteins in single mammalian cells is high enough to minimize the sampling (counting) noise [6].
To further evaluate relative quantification, beyond the results for two samples diluted to single-cell level (Fig. 1b), we consolidated the data from 73 1xM standards and computed all pairwise correlations. This 584-dimensional matrix was projected onto its first two principal components (PC). The largest PC accounts for 64% of the total variance in the data and perfectly separates all samples corresponding to HEK-293 cells or monocytes. Crucially, the cell lysates diluted to single-cell levels cluster the same way as the corresponding carrier samples, indicating that the clustering is driven by cell type-specific protein differences rather than by reproducible artifacts. Thus, the standards with isobaric carriers (modeling SCoPE2 sets) can reliably quantify protein abundances at the single-cell level.
Model system and technical benchmarks
Having demonstrated that proteins from 1xM standards can be quantified with low noise, we next applied SCoPE2 to the analysis of single cells. As a model system, we chose monocytes differentiated into macrophage-like cells in the presence of an agonist of protein kinase C, phorbol-12-myristate-13-acetate (PMA), Fig. 3a. We chose this system since it provides a clear benchmark—the ability to identify two closely related but distinct cell types. This system also presents an open research question: are macrophage-like cells produced from this differentiation as homogeneous as the monocytes from which they originate or more heterogeneous? To answer this question independently from the heterogeneity inherent to primary monocytes, we used a homogeneous monocytic human cell line, U-937 [32].
The SCoPE2 workflow (Fig. 1a) can be used with manually picked cells, FACS-sorted cells, or cells isolated by microfluidic technologies that minimize the volume of the droplets containing cells [6]. Here, we used a BD FACSAria I cell sorter to sort single cells into 384-well plates, one cell per well, see Fig. 3a and Additional file 1: Fig. S1. The single cells from two biological replicate preparations of the differentiation protocol were sorted and prepared into 179 SCoPE2 sets, as depicted in Fig. 1a and Fig. 3a. The sorting followed a randomized layout to minimize biases, and sample preparation was automated as described in the “Methods” section. Each set has an isobaric carrier corresponding to about 200 cells and a reference corresponding to 5 cells. The size of the carrier was chosen to optimize the number of quantified proteins and the number of copies sampled per protein based on guidelines from controlled experiments [29]. The carriers of some sets were sorted individually, with 200 cells per well, while the carriers and reference channels of other sets were diluted from a larger bulk sample. We analyzed two biological replicates available at massIVE [34, 35]: the first used TMT 11-plex [34], and the second used TMT pro 16-plex [35].
First, we sought to identify successfully analyzed single cells based on low variability for the relative protein quantification derived from different peptides, Fig. 3b. Specifically, different peptides originating from the same protein should provide similar estimates for the protein fold changes across single cells, as observed in SCoPE2 data (Additional file 1: Fig. S6). The variability between estimates from different peptides was quantified by the coefficient of variation (CV) for all peptides originating from the same protein, i.e., standard deviation/mean for the RI ratios. Then, the median CV (across all proteins) of a single cell provides a measure for the consistency of relative protein quantification in that cell.
To estimate the background noise, we incorporated control wells into SCoPE2 sets, Fig. 2a. These control wells did not receive cells but were treated identically to wells with single cells, i.e., received trypsin, TMT additions, and hydroxylamine. Thus, any reporter ion intensities detected in the control wells should correspond to the coisolation and background noise; this noise should be highly variable and uncorrelated for different peptides originating from the same protein because the peptides likely elute at different points in the chromatography and have different mass-to-charge ratios. This high variability is reflected in the high median CV for the control wells, higher than for single-cell wells. The distribution of CVs for all single cells and control wells indicates that 1490 single cells have lower quantification variability than the control wells, Fig. 3b. These single cells were analyzed further, see the “Evaluating single-cell quantification” section.
As explained in Fig. 2, the quantitative accuracy of LC-MS/MS measurements depends on sampling the apex of chromatographic peaks. Thus, we used DO-MS to optimize the instrument parameters so that chromatographic peaks are sampled close to their apexes [25]. As a result, most ions were sampled for MS2 analysis within 3 s of their apexes (Fig. 3c). The combination of sampling close to the apex, sharp chromatographic peaks, and narrow isolation widows (0.7 Th) minimized the simultaneous coisolation of different peptide ions for MS2 analysis. Indeed, MaxQuant estimated over 97% median purity of the ions sent for MS2 analysis, Fig. 3d. Independent estimates of coisolation by Proteome Discoverer confirm the purity of MS2 spectra, Additional file 1: Fig. S5b. The high spectral purity and the use of retention times to bolster peptide sequence identification [26] allowed assigning a peptide sequence to about 35% of the MS2 spectra from each SCoPE2 run at 1% FDR, Fig. 3e. The high spectral purity supports accurate quantification, consistent with a linear dependence between the measured signal and the input amount shown in Additional file 1: Fig. S3a.
On average, SCoPE2 quantifies over 2500 peptides corresponding to about 1000 proteins per SCoPE2 set analyzed on 1-h chromatographic gradient, Fig. 3f. While longer gradients can increase this coverage, they will reduce the number of cells analyzed per unit time. Since most single-cell analysis requires analyzing large numbers of single cells, we focused on shorter nLC gradients that increase the number of proteins quantified across many single cells rather than merely the number of proteins per cell [6]. We applied further filtering of peptides and proteins to ensure FDR below 1% within each peptide, and within each protein, not just across all peptides and data points; see the “Methods” section. Not all proteins are quantified in all single cells because some protein values are missing. The first type of missing values is due to peptides sent for MS2 analysis but resulting in too low signal-to-noise ratio (SNR) for quantification. This is a minor contribution to missingness. Indeed, only 10% of the data for peptides identified in a SCoPE2 set were missing because the RI signal in some single cells is below the SNR threshold. This missingness can be reduced further by increasing ion sampling (accumulation) times or by narrowing chromatographic peaks, see Fig. 2 [6, 29]. The second and major type of missing values is due to peptides not sent for MS2 analysis. This is because the MS instrument does not have time to send for MS2 scans all ions detected at MS1 survey scans, and thus, it isolates for MS2 analysis different peptide subsets during different SCoPE2 runs (see Fig. 1a and Fig. 2), a well-described phenomenon for data-dependent acquisition (DDA) [11, 12]. As a consequence, we obtained quantification for a diverse set of 3042 proteins spanning a wide dynamic range, Additional file 1: Fig. S2. Thus, the average number of data points (single cells) per protein is about 305. This type of missing data, unlike the first, is a product of experimental design and not necessarily a lack of sensitivity. To increase the reproducibility of protein measurement between SCoPE2 sets, a specific list of peptides can be targeted for analysis in each set [6, 17].
Cell type classification
High-throughput single-cell measurements are commonly used to identify and classify cell types. Thus, we sought to test the ability of SCoPE2 data to perform such classification. As a first approach, we performed principal component analysis (PCA), Fig. 4a. Using all 3042 quantified proteins, PCA separates the cells into two mostly discrete clusters along PC1, which accounts for 29% of the total variance of the data. Color-coding the cells by their labels indicates that the clusters correspond to the monocytes and the macrophages, Fig. 4a. These results are qualitatively recapitulated if the PCA is performed on the raw data without imputation and batch correction, Additional file 1: Fig. S7. To evaluate whether we can assign cell types based on the abundance of monocyte and macrophage associated proteins, we color-coded each cell by the median abundance of proteins identified to be differentially abundant from analyzing bulk samples of monocytes and macrophages, Fig. 4b. The resulting color-coding of the cells coincides with their cell types, supporting the utility of SCoPE2 data for classifying cell types.
The macrophage cluster appears more spread-out, which might reflect either the larger number of macrophage cells or increased cellular heterogeneity during the differentiation process. To distinguish, between these possibilities, we computed the pairwise correlations between monocytes and those between macrophages, Fig. 7a. The lower pairwise-pairwise correlations between macrophages suggest that indeed differentiation increases heterogeneity, Fig. 7a.
While the PCA clustering is consistent with the known cell types, it is inadequate to benchmark protein quantification. As a direct benchmark, we averaged in silico the single-cell data to compute monocytes/macrophages protein fold changes and compared these estimates to the corresponding fold changes measured from bulk samples, i.e., averaging across single cells by physically mixing their lysates, Fig. 4c. The correlation (ρ = 0.88) indicates that SCoPE2 can accurately measure protein fold changes in single cells.
In principle, the cell type classification in Fig. 4a, b may be driven by relatively few abundant structural proteins while less abundant regulatory proteins, such as kinases, receptors, and transcription factors, might be poorly quantified. To evaluate this possibility, we applied the PCA analysis using only proteins functioning in signaling (Fig. 4d) or only the least abundant proteins quantified by SCoPE2 (Fig. 4e). The results indicate that these protein groups also correctly classify cell types (Fig. 4d, e), albeit the fraction of variance captured by PC1 is lower for the signaling proteins (24%) shown in Fig. 4d. The relative quantification of proteins from both sets correlates positively (ρ = 0.75, ρ = 0.72) to the corresponding bulk protein ratios, Fig. 4d,e.
Macrophages exhibit a continuum of proteome states
Next, we turn to the question of whether the homogeneous monocytes differentiated to similarly homogeneous macrophages, Fig. 3a. To this end, we performed an unsupervised spectral analysis of the cell population, which allowed us to characterize cellular heterogeneity without assuming that cells fall into discrete clusters. To do so, the cells were connected into a graph based on the correlation of their proteome profiles. We then examined the eigenvector of the graph Laplacian matrix with the smallest non-trivial eigenvalue (Fiedler vector, Eq. 1 [36]), which captures the most prominent structural axis of the cell similarity graph; see the “Methods” section. We then sorted cells according to their Fiedler vector values, thus capturing the most prominent aspect of cellular heterogeneity Fig. 5a. This cell clustering is based on hundreds of proteins with differential abundance between monocytes and macrophage-like cells, shown in Fig. 5a: Proteins with higher abundance in monocytes are enriched for proliferative functions, including chromatin organization and translation [37]. Proteins with higher abundance in macrophages are enriched for cell surface signaling and cell adhesion proteins, including intermediate filament protein vimentin, Fig. 5b. These enrichment results are consistent with the functional specialization of monocytes and macrophages and further validate the ability of SCoPE2 data to recapitulate known biology.
To explore the heterogeneity within the macrophages, we applied the same spectral analysis as in Fig. 5a, but this time only to the macrophage-like cells. The distribution of the elements of the associated Fiedler eigenvector (defined by Eq. 1) suggests that the cellular heterogeneity observed in this population is better described by a continuous spectrum rather than by discrete clusters, Fig. 5c. Indeed, the heatmap of protein levels for macrophage-like cells ordered based on the Fiedler eigenvector shows that most proteins change gradually across a continuous spectrum, Fig. 5c. Analyzing the proteins from this gradient, we observed a remarkable trend: genes previously identified as differentially expressed between M1 and M2-polarized primary human macrophages [38] are also differentially expressed between single macrophage cells. For example, the cells at the left edge of Fig. 5c show high expression of genes upregulated in M2-polarized macrophages, decreasing monotonically from the left to right of Fig. 5c. Genes upregulated in M1-polarized primary human macrophages appear to be expressed in a reciprocal fashion, with lower expression at the left edge of Fig. 5c, increasing monotonically across the figure. These results from Fig. 5a, c are qualitatively recapitulated if the spectral clustering is performed on the raw data without imputation and batch correction, Additional file 1: Fig. S7.
Joint analysis of single-cell protein and RNA data
We also analyzed single cells from two biological replicates of differentiating monocytes (Fig. 3a) by scRNA-seq using the 10× Chromium platform [39] and compared the cell clustering and differential genes measured by scRNA-seq and SCoPE2. We first consider measurement noise since it may contribute to apparent differences between the protein and RNA measurements [31]. Specifically, we can expect to confidently identify differential abundance between cell types only for genes whose biological variability exceeds the measurement noise.
A major source of measurement noise for both SCoPE2 and scRNA-seq is sampling low copy number of molecules per gene. A low copy number of sampled molecules results in significant counting noise [6]. This noise arises because scRNA-seq and SCoPE2 sample only a subset of the molecules from a single cell. This sampling process contributes to a counting error: The standard deviation for sampling n copies is sqrt(n) (from the Poisson distribution), and thus the relative sampling error, estimated as standard deviation over mean, is sqrt(n)/ n = 1/sqrt(n), Fig. 6a. Thus our optimization of SCoPE2 aimed to increase ion delivery not merely the number of identified peptides [25]. Similarly, to increase the sampling of mRNAs, we sequenced the 10× library to an average depth of 100,000 reads per cell. Furthermore, we chose for subsequent analysis only representative cells with at least 10,000 UMIs per cell, Additional file 1: Fig. S4.
To estimate our sampling error, we sought to convert the RI abundances (i.e., the barcodes from which SCoPE2 estimates peptide abundances) into ion copy numbers. To do so, we extracted the signal to noise ratios (S/N) for RIs and multiplied these ratios by the number of ions that induces a unit change in S/N. Since our Q-Exactive basic Orbitrap operated at 70,000 resolving power, a S/N ratio of 1 corresponds to about 6 ions [42, 43]. Thus, for our system, a S/N ratio of 50 corresponds to 300 ions; see the “Methods” section for details. The results in Fig. 6a indicate that SCoPE2 samples 10–100-fold more copies per gene than 10× Genomics, which corresponds to smaller sampling (counting) errors. The sampling can be increased further by increasing the ion sampling times, see Fig. 2.
To investigate the similarities and differences at the RNA and protein levels, we compared the correlation vectors of the scRNA-seq and SCoPE2 data as previously described [13, 40]. Specifically, for each dataset, we computed the matrices of pairwise Pearson correlations between genes, averaging across the single cells. Then to quantify the similarity of covariation of the ith gene, we correlated the ith vectors of RNA (ri) and protein (pi) correlations. The correlations between ri and pi for all genes are bimodally distributed, with a large positive mode, Fig. 6b. Yet, they become close to zero when the order of genes is permuted, see the null distribution in Fig. 6b. We used the two modes of the distribution of correlations between ri and pi to define genes with similar protein and RNA covariation (cluster 1) and genes with opposite covariation (cluster 2). The genes not included in either cluster are not differentially abundant between monocytes and macrophages at either the protein or the RNA level. The genes within a cluster correlate to each other similarly as indicated by the high correlations between ri and pi when these vectors are comprised only by the correlations of genes from cluster 1 or from cluster 2, Fig. 6b. To further confirm the coherence of these clusters, we displayed the RNA and protein levels of their genes as heatmaps in Fig. 6c using common PCA to order the cells [44]. Ordering the single cells and the cluster 1 genes based on the first common principal component (CPC 1) of the RNA and protein datasets confirms that the majority of genes exhibit qualitatively similar RNA and protein profiles. Gene set enrichment analysis [45] of these genes identified many biological functions, including antigen presentation, cell adhesion, cell proliferation, and protein synthesis, Additional file 1: Fig. S9. This similarity of RNA and protein profiles extends to the macrophage heterogeneity that we discovered in the SCoPE2 data Fig. 5c: The gradient is also observed with the RNA data, and it correlates to polarization marker genes albeit with smaller correlations, Additional file 1: Fig. S8. The genes from cluster 2 have opposite RNA and protein profiles and are enriched for functions including Rab GTPase activity, and protein complex assembly, Additional file 1: Fig. S10. The similarity between the RNA and protein abundances of genes from cluster 1 (Fig. 6b, c) suggests that the RNA and protein data might be projected jointly. To this end, we used Conos to generate a joint graph integrating all analyzed cells and to project them on the same set of axes, Fig. 6d–f. The joint projection results in two distinct clusters corresponding to the monocytes and macrophages, Fig. 6d. The cells analyzed by SCoPE2 are at the center of the clusters and are surrounded by more diffused clusters of the cells analyzed by scRNA-seq.
Color-coding the cells by biological and technical replicate reveals no significant clusters and thus no residual artifacts and batch effects, Fig. 6e. In contrast, color-coding the cells with monocytes and macrophage markers reveals that the clusters correspond to the cell types, Fig. 6f. These results demonstrate that, at least for our model system, the low dimensional projections of single-cell protein and RNA measurements result in similar clusters albeit the clusters based on RNA measurements are more spread out, suggesting more biological and technical variability in the RNA measurements.
Transcriptional and post-transcriptional regulation
Next, we sought to explore transcriptional and post-transcriptional regulation based on joint RNA and protein analysis. Indeed, it has been suggested that the variability between single cells can be used to infer gene regulation [46, 47]. For example, transcriptional regulation should be reflected in the correlations between transcription factors (TFs) and the abundance of mRNAs whose transcription they regulate while post-transcriptional regulation can be seen in the joint distributions of RNA and protein abundances of the same gene. For such analysis, we needed to pair cells analyzed by 10× Genomics and by SCoPE2. To pair cells in similar states, we used the observations that proteins and transcripts of many genes covary together as shown in Fig. 6. Thus, we ordered the cells based on the loadings of their common first principal components and paired cells having the same rank as shown in Fig. 7a for NCL. For this gene, the protein and RNA abundances exhibit similar trends as reflected in a correlated joint distribution, Fig. 7a. This pattern is typical for abundant genes from cluster 1 as defined in Fig. 6b. Some less abundant genes from cluster 1 also exhibit similar trends though with more discrete variability at the RNA level, as exemplified with RHOC, Fig. 7b. Other genes (particularly from cluster 2 as defined in Fig. 6b) show more divergent trends in their RNA and protein profiles as exemplified by TNFSF13B (Fig. 7c) and by the tumor suppressor protein p53 (Fig. 7d). This divergence in the levels of p53 protein and its transcript levels suggests post-transcriptional regulation, consistent with previous research demonstrating that p53 levels are regulated by protein degradation [48].
Having paired RNA and protein levels, we next explored the correlations expected from the transcriptional regulatory network. Specifically, we used the manually curated database TRRUST [49] to identify target genes whose transcription is either activated or repressed by TFs. We focused on TFs and their targets that are detected both by 10× Genomics and SCoPE2 and are known to be either activated or repressed by the TFs. We correlated the abundance of these TFs (at either mRNA or protein level) to the abundance of RNAs whose transcription they regulate. For p53 (encoded by TP53), the protein abundance correlates positively to the activated targets and negatively to the repressed targets, Fig. 7e. Other trans factors (including RELA, E2F1, and SIRT1) also correlate to their target genes as expected (negatively to repressed targets and positively to activated targets), but the number of transcripts characterized in the TRRUST database [49] and quantified in our dataset is too small to establish the statistical significance of their trends. In contrast, the correlations between the mRNA coding for p53 differ from the expectation (Fig. 7c), consistent with recent observations [50].
These results suggest that the measured protein abundances may be more informative for TF activity. Indeed, p53 is known to be regulated by proteolysis, which alters its protein levels without changing RNA abundances [48]. All correlations between the TFs and their target RNAs are weak in our dataset, which likely reflects both measurement noise and unaccounted for regulatory mechanisms, such as nuclear localization, post-translational modifications of these TFs, and additional regulators of the transcription and degradation same RNAs.