Spatiotemporal dissection of the cell cycle with single-cell proteogenomics

a, The specificity of several dozen antibodies (45 in total with 10 examples presented here; Supplementary Table 3) targeting proteins with CCD expression was validated by co-localization with mNG-tagged protein. Scale bars, 10 μm. b, The specificity of several dozen antibodies (79 in total with 10 examples presented here; Supplementary Table 3) targeting proteins with non-CCD variable expression was validated by co-localization with mNG-tagged protein. Scale bars, 10 μm. c, The specificity of 29 antibodies (three presented here; Figs. 1c, 4b, Supplementary Table 3) was validated with siRNA-mediated gene silencing, which resulted in significantly lower staining intensity. Scale bars, 10 μm. For box plots: centre line, median; box, Q1 and Q3; whiskers, 1.5× IQR below Q1 and above Q3; points, outliers. TTC21B, P = 6.3 × 10⁻³¹; DUSP19, P = 0.0025; NET1, P = 9.7 × 10⁻⁸ by binomial one-sided tests. TTC21B, n = 620, 358; DUSP19, n = 125, 147; NET1, n = 162,1 49 independent cells for control and siRNA, respectively. d, The specificity of 321 antibodies was validated by independent antibody staining, as exemplified by SCIN. Scale bars, 10 μm.

a, The specificity of the anti-ANLN antibody was shown by co-localization with mNG-tagged ANLN. Scale bars, 10 μm. b, Time-lapse microscopy for mNG-tagged ANLN protein demonstrates an increase in nuclear expression over the course of the cell cycle with a peak in late G2. During mitosis and cytokinesis, ANLN localizes to the cell membrane and is involved in the formation of the cleavage furrow (t = 17–18 h and t = 35 h). Scale bar, 10 μm. c, The specificity of the anti-AURKB antibody was validated by co-localization with mNG-tagged AURKB. Scale bars, 10 μm. d, Time-lapse microscopy for mNG-tagged AURKB protein demonstrates an increase in protein expression over the course of the cell cycle with a peak in late G2. During mitosis, AURKB localizes to the mitotic chromosomes at hours 13 and 29. Scale bars, 10 μm. e, The specificity of the anti-GATAD1 antibody was validated by co-localization with mNG-tagged GATAD1. Scale bars, 10 μm. f, Time-lapse microscopy for mNG-tagged GATAD1 protein demonstrates a decrease in nuclear abundance throughout the cell cycle (mitosis at 35 h). Scale bars, 10 μm. g, Time-lapse microscopy for mNG-tagged KIFC1 protein demonstrates a CCD translocation from the nucleoplasm to the kinetochores in late G2 into mitosis at hours 14 and 30 and in early G1. Scale bars, 10 μm. h, Time-lapse microscopy for mNG-tagged NET1 protein demonstrates a CCD increase in protein abundance and translocation from the nucleoplasm to the nucleolus (mitosis at 35 h). Scale bars, 10 μm.

Gating strategy used to sort U2OS FUCCI cells into three clusters representing three distinct segments of the cell cycle. This strategy distinguished the cells expressing the red tag (CDT1, G1) from the cells expressing the green tag (GMNN, S and G2). Cells that expressed both markers simultaneously, marking the S-transition, were clustered in a third group. The RNA velocities calculated for these cells (Fig. 2c, Extended Data Fig. 7a, b) validate this sorting strategy, as a cell cycle reflecting interphase progression is clearly observed.

a, Scatterplot showing the three clusters generated by k-means clustering based on kurtosis and skewness as features for CCD proteins. b, Scatterplot showing the three clusters generated by k-means clustering based on kurtosis and skewness as features for non-CCD proteins. c, Violin plots and histograms showing the population distributions of the normalized mean intensity of each cell per protein for three selected CCD proteins (GATA6, CCNB1, and DEF6). d, Bimodal protein distributions were evaluated for cell cycle dependence separately in both low- and high-expressing cells if the two populations were determined to be distinct. This determination was performed using a Kruskal–Wallis test, adjusted for multiple testing, and if they had greater than a twofold difference in expression between them. e, GATA6 expression over the cell cycle. While GATA6 produces a bimodal population intensity distribution, it represents a single population of cells and was evaluated for cell cycle dependence as such, rather than as two populations of high- and low-expressing cells. f, SLC25A42 expression over the cell cycle is exhibited as two distinct populations of high- and low-expressing cells (left). The low-expressing cells (centre) have CCD expression, forming a second harmonic over the cell cycle, while the high-expressing cells (right) do not display correlation to the cell cycle. g, HPSE expression over the cell cycle has two distinct populations of high- and low-expressing cells that both display non-CCD expression, which may point to a temporal component related to a different cellular process.

a, Images from a time-lapse microscopy series of FUCCI cells show that the length of the cell cycle phases marked by the different colours of the FUCCI tags are 10.8 h for G1, 2.6 h for S-transition, and 11.9 h for G2. U2OS FUCCI cells allow monitoring of the cell cycle by expressing two fluorescently tagged cell cycle markers: CDT1 during G1 phase (red), GMNN during S and G2 phases (green), and both during the S-transition (yellow). Scale bars, 10 μm. b, The polar coordinate model transfers the FUCCI marker information into a linear model of pseudotime. The tagged CDT1 (red) and GMNN (green) log-intensities are displayed across pseudotime, which they are used to calculate. The solid line is the moving average across 100 cells (264,863 total), and light shading fills between the 25th and 75th percentiles of intensity. c, Cell division time can also be predicted from scRNA-seq data, and so we performed a comparison between the pseudotime calculated using FUCCI marker intensities in this work (x-axis, colour mapping) and the pseudotime inferred from scRNA-seq markers using the tool Peco⁵⁹ (y-axis). The most notable feature of this comparison is the more continuous nature of the pseudotime calculated using FUCCI marker intensities, indicating that the experimental FUCCI marker intensities provide a better basis than inferred pseudotime for the integration of single-cell proteomic and transcriptomic data. d, CNVs identified in single-cell transcriptomic data correlate to the cell cycle, thus validating the FUCCI model. e, BrdU incorporation (blue) into FUCCI cells. Incorporation of the synthetic nucleotide BrdU occurs during DNA replication in S-phase of proliferating cells and remains through the rest of the cell cycle. Here, G1 cells with red nuclei show no incorporation of BrdU, whereas most of the S and G2 cells with green nuclei show incorporation of BrdU (bottom left and centre; G1, n = 45; S-trans, n = 32; S and G2, n = 73 independent cells). For box plots: centre line, median; box, Q1 and Q3; whiskers, 1.5× IQR below Q1 and above Q3; points, outliers. Based on this, we can estimate the proportion of cells in each phase (bottom right). These results show the validity of the FUCCI model. Scale bar, 10 μm.

a, Box plot showing gene expression levels for all proteins that exhibit cell-to-cell heterogeneity and all proteins in the HPA Cell Atlas (mapped proteome). There is no significant difference (P < 0.01) in gene expression levels between the antibodies included in this study and all antibodies used in the HPA Cell Atlas (P = 0.03, two-sample Student’s t-test; variation, n = 1,607; mapped proteome, n = 9,806 independent proteins). b, Bar plot showing staining intensity levels for antibodies included in this study and all proteins in the HPA Cell Atlas (mapped proteome). The average immunofluorescence signal intensity measured for proteins showing cell-to-cell variation and for all proteins in the HPA Cell Atlas shows no significant difference for low signal antibodies. Low signal antibodies are not enriched among the cell-to-cell variability dataset (strong, P = 1.12 × 10⁻⁵; moderate, P = 0.0004; weak, P = 0.99 by binomial one-sided tests; asterisk denotes significance). c, Since the HPA antibodies are purified on affinity columns coupled with their corresponding antigen, antibody concentration after purification can serve as a proxy for antibody affinity (albeit not a perfect proxy). Box plot showing antibody concentration levels for all proteins that showed cell-to-cell heterogeneity and all proteins in the HPA Cell Atlas (mapped proteome). There is no significant difference between the antibodies included in this study and all antibodies used in the HPA Cell Atlas, hence we can conclude that the cell-to-cell heterogeneity is probably not due to differences in antibody affinity. The average antibody concentration for the antibodies published on the HPA Cell Atlas is 0.1710 mg ml⁻¹, and the average concentration for the antibodies used in this study is 0.1712 mg ml⁻¹ (P = 0.1084, two-sample Student’s t-test; variation, n = 1,415; mapped proteome, n = 14,942 independent proteins). d, Box plot showing the variance for all proteins that showed cell-to-cell heterogeneity and microtubules. Proteins showing heterogeneity show significantly higher variance than the variance of the microtubules measured from each well (P = 1.6 × 10⁻²⁹², one-sided Kruskal–Wallis test; n = 1,180 independent IF stains). e, Box plot showing the Gini index for all proteins that showed cell-to-cell heterogeneity and for microtubules. Proteins displaying heterogeneity show significantly higher Gini indexes than microtubules (P < 5 × 10⁻³²⁴, one-sided Kruskal–Wallis test; n = 1,180 independent IF stains). f, Gene ontology enrichment analysis for CCD proteins shows significantly enriched terms for the biological processes domain. Each circle represents a GO term, and line width corresponds to the number of genes that overlap between the two connected gene sets. Similar terms are grouped and labelled. g, GO enrichment analysis for non-CCD enzymes shows enrichment for basic metabolic functions, while CCD enzymes are enriched for cell cycle functions. h, Subcellular localizations of CCD and non-CCD proteins. Asterisk denotes enrichment relative to the HPA-mapped proteome. For CCD proteins: cleavage furrow, P = 0.0003; cytokinetic bridge, P = 4.3 × 10⁻⁹³; kinetochores, P = 3.5 × 10⁻⁵; midbody ring, P = 2.3 × 10⁻¹⁵; midbody, P = 9 × 10⁻³⁶; mitotic chromosome, P = 9.8 × 10⁻⁶⁶; mitotic spindle, P = 6.3 × 10⁻⁴¹; nucleoli, P = 0.005. For non-CCD proteins: intermediate filaments, P = 1.9 × 10⁻⁸; mitochondria, P = 1.9 × 10⁻¹¹; nuclear bodies, P = 3.01 × 10⁻⁵; nucleoli, P = 8.8 × 10⁻⁷. i, Protein–protein interaction (PPI) network of CCD proteins and CCD transcripts using the STRING⁷⁰ database. Proteins with known associations to the cell cycle (by GO term, teactome pathway, or cyclebase phenotype) are represented as green squares (left). Transcriptionally regulated CCD proteins are tightly clustered in the centre of the network, with an extended network of novel CCD genes and particularly ones that are not transcriptionally regulated (right). For box plots: centre line, median; box, Q1 and Q3; whiskers, 1.5× IQR below Q1 and above Q3; points, outliers.

a, Dimensionality reduction of scRNA-seq data by UMAP with RNA velocities overlaid in arrows and coloured according to clustering. Each point represents a cell, and cluster 5 (brown) appears to be cells that are sequestered from cycling after G1. b, Testing cutoffs of between 1% and 10% additional percentage variance explained by the cell cycle in scRNA-seq data over random shows that 8% is an appropriate cutoff. At cutoffs 1–3%, the cycle is completely removed from non-CCD proteins, but it appears that there are enough false-positive calls to impact the structure of the UMAP for CCD-only genes. At 4–8%, there is a small eddy in the UMAP for non-CCD genes but no cohesive cycle, and at 9–10%, there is a distinct cycle encompassing all cells. We chose to use a conservative final cutoff of 8% for calling both CCD transcripts and proteins to err on the side of having fewer false positives while balancing precision and recall for these analyses. c, Confirmation that the UMAP structure for scRNA-seq data is stable across various parameter selections. d, BIRC5 and UBE2C provide anecdotal examples of cycling genes (orange) with differential cycling between one or more transcript isoforms (green). Line, moving mean; darker shade, 25th to 75th percentile range; lighter shade, 10th to 90th percentile range; points, individual cell data.

a, Percentage variance explained for RNA is shown against the Gini index for the gene (blue, CCD; red, non-CCD). b, Randomization analysis of the protein IF data (left) and scRNA-seq data (right) for each gene was used to determine whether a protein or RNA was CCD (blue) or non-CCD (red). The significance scores, adjusted for multiple correction, on the vertical axis show that nearly all proteins and RNAs are significantly different from random, and so requiring 8% additional percentage variance explained by the cell cycle over random was the predominant cutoff. c, Examples of NFAT5 protein IF data (blue points and trace) and randomizations of cell order in pseudotime (red points and trace). These examples provided (from left to right) produced the minimum, first quartile, median, third quartile, and maximum percentage variance explained by the random fluctuations. d, LASSO analysis of marker genes was overly conservative compared to the pseudotime analysis in this work. A higher false-negative rate for calling CCD genes (top) and proteins (bottom) leaves a cyclic pattern in the UMAP dimensionality reduction expected of CCD genes and proteins (left) in the non-CCD ones (right). e, Principal component analysis showed no discernible batch effect between the three plates with 384 cells each, and instead the cell cycle phases roughly assigned by FACS sorting provide clear separation in the first two components (PCs). Results are shown before and after filtering out the non-cycling cells. f, Comparing the individual batches to the combined data for RNA-seq again confirms that no batch effects were present in the RNA-seq data. Each plot contains relative RNA expression (0 to 1, y-axes) versus cell division time (0 to 25.3 h, x-axes). Line, moving mean; darker shade, 25th to 75th percentile range; lighter shade, 10th to 90th percentile range; points, individual cell data.

Bottom, examples of temporal expression profiles for single-cell protein (blue) and RNA expression (orange). Line, moving mean; darker shade, 25th to 75th percentile range; lighter shade, 10th to 90th percentile range; points, individual cell data. The compartment for which the protein abundance was measured is denoted above: nuc, nucleus; cyto, cytoplasm; cell, entire cell. The proteins of interest are shown in green, and microtubules are shown in red. Similar temporal profiles between proteins indicate compartmentalization and association in time. For the G1 group, well-known CCD proteins such as ORC6, which is required for cell entry into S phase⁸⁷, and MCM10, which is required for DNA replication⁸⁸, were found to have similar patterns to the novel CCD protein ZNF32, whose overexpression has been associated with a shorter survival time in lung adenocarcinoma cells⁸⁹. The group that has peak expression at the end of G1 includes known proteins such as CCNE1, along with the novel CCD protein DUSP19 (Fig. 3), a phosphatase whose depletion results in increased mitotic defects⁹⁰. The G2 group includes known proteins, such as PC (Fig. 3), CCNB1, AURKB (Extended Data Fig. 2c, d) and BUB1B. Novel CCD proteins in this group include the phosphatase DUSP18 (Fig. 3); the transcription factor NFAT5 (Fig. 1F); a retrograde intraflagellar transport (IFT) protein TTC21B (Fig. 3), which has been implicated in several ciliopathies and may be implicated in additional roles beyond the primary cilium⁹¹; the oestrogen-sulfating enzyme PAPSS1 (Fig. 3), for which overexpression has been reported to affect proliferation⁹²; the methyltransferase N6AMT1 (Fig. 3); the uncharacterized protein PHLDB1 (Fig. 3); the enzyme DPH2; and transcription factor FLI1 (Fig. 3). Scale bar, 10 μm.

a, Hierarchical clustering of bulk transcript expression (log-transformed TPM) in various normal and cancer tissue types for known CCD proteins. The expression levels of the proliferation markers MCM6, CDK1, PCNA, MCM2 and KI67 are highlighted at the top as a general measure of the proliferative activity of the tissues. Four clusters are identified: c1 and c2 contain mostly normal tissues with midrange expression of the proliferation markers; c3 contains tissues with high expression of the proliferation markers, including tumours; and c4 contains normal tissues with low proliferative activity. b, Box plots displaying the average transcript level for known CCD proteins in the four different clusters from a. There is a significant difference in gene expression levels between the different clusters (c4–c1, P = 0.04; c4–c3, P = 0.006; c3–c2, P < 2 × 10⁻¹⁶; c3–c1, P = 3.5 × 10⁻⁸; c4–c2, P = 4 × 10⁻¹⁴; c2–c1, P = 6.9 × 10⁻¹⁴) by one-sided Kruskal–Wallis test. n = 238 independent tissues for each cluster. c, Hierarchical clustering of bulk transcript expression (log-transformed TPM) in various normal and cancer tissue types for novel CCD proteins. The expression levels of the proliferation markers MCM6, CDK1, PCNA, MCM2 and KI67 are highlighted at the top as a general measure of the proliferative activity of the tissues. Four clusters are identified: c1 contains normal tissues with midrange expression of the proliferation markers and tissues with high expression of the proliferation markers, including tumours; c2 contains bone marrow and bone marrow cancer tissues; c3 contains cerebral tissues with testis; and c4 contains normal tissues with low proliferative activity. d, Box plots displaying the average transcript level for novel CCD proteins in the four different clusters from c. There is a significant difference between gene expression in the different clusters, except between clusters c2 and c4 (c4–c1, P = 4.3 × 10⁻¹⁴; c4–c2, P = 0.6; c4–c3, P < 2 × 10⁻¹⁶; c1–c2, P = 2.9 × 10⁻¹²; c1–c3, P = 0.001; c3–c2, P = 3.7 × 10⁻¹⁴) by one-sided Kruskal–Wallis test. n = 301 independent tissues for each cluster. For box plots: centre line, median; box, Q1 and Q3; whiskers, 1.5× IQR below Q1 and above Q3; points, outliers.