Phospho-seq: integrated, multi-modal profiling of intracellular protein dynamics in single cells

Abstract

Cell signaling plays a critical role in neurodevelopment, regulating cellular behavior and fate. While multimodal single-cell sequencing technologies are rapidly advancing, scalable and flexible profiling of cell signaling states alongside other molecular modalities remains challenging. Here we present Phospho-seq, an integrated approach that aims to quantify cytoplasmic and nuclear proteins, including those with post-translational modifications, and to connect their activity with cis-regulatory elements and transcriptional targets. We utilize a simplified benchtop antibody conjugation method to create large custom neuro-focused antibody panels for simultaneous protein and scATAC-seq profiling on whole cells, alongside both experimental and computational strategies to incorporate transcriptomic measurements. We apply our workflow to cell lines, induced pluripotent stem cells, and months-old retinal and brain organoids to demonstrate its broad applicability. We show that Phospho-seq can provide insights into cellular states and trajectories, shed light on gene regulatory relationships, and help explore the causes and effects of diverse cell signaling in neurodevelopment.

Fig. 1. Phospho-seq experimental workflow and pilot experiment.

A Schematic of Phospho-seq workflow. B Schematic of antibody conjugation procedure. C Protein gel results of two antibody (Ab) purification methods using the Abcam BSA Removal Kit (left panel) and Promega Magne Protein G beads (right panel). This experiment was performed once. D Protein gel of mTz-PEG4-NHS labeled antibodies incubated with different quantities of TCO-PEG4-NHS labeled ssDNA tags. This experiment was performed once. E Flow cytometry plots of K562 and iPS cells stained with unconjugated and conjugated SOX2 antibodies (left panel) and cells stained with unconjugated and conjugated + single-stranded DNA binding protein SOX2 antibodies (right panel) with unstained controls. F UMAP representation from scATAC-Seq of K562 and iPS cells colored by demultiplexed HTOs assigned to each cell. G UMAP representation of K562 and iPS cells colored by normalized ADT values for OCT4. H UMAP representation of K562 and iPS cells colored by normalized ADT values for GATA1. I Coverage plot of chromatin accessibility of K562 and iPS cells at the POU5F1 (OCT4) genomic locus. J Violin plot of chromVAR scores for the OCT4 TF binding motif (MA1115.1) in K562 and iPS cells. K Coverage plot of chromatin accessibility of K562 and iPS cells at the GATA1 genomic locus. L Violin plot of chromVAR scores for the GATA1 TF binding motif (MA0140.2) in K562 and iPS cells. M Schematic of PI3K/AKT/mTORC1 pathway activation and repression paradigm used in this experiment. N Scatter plot of pseudobulked chromatin accessibility data in 5 kb windows across the genome comparing inhibited K562 cells with stimulated K562 cells. Red line indicates perfect correlation between the two conditions. O Violin plot of normalized pRPS6 values in stimulated (Stim) and inhibited (Inhib) K562 cells. P Violin plot of normalized RPS6 values in stimulated (Stim) and inhibited (Inhib) K562 cells.

Fig. 2. Phospho-seq on retinal organoids.

A Schematic of retinal organoid differentiation. B UMAP representation of cells and cell type assignments based on the ATAC-seq modality in Phospho-seq. RPC retinal precursor cell, MG müller glia, PC precursor cell, PR-PC photoreceptor precursor, RPE retinal pigment epithelium. C Coverage and violin plots of the gene promoter and protein level respectively of nuclear protein SOX9. Color and order are the same as in (B). D Coverage and violin plots of the gene promoter and protein level respectively of cytoplasmic protein RCVRN. Color and order are the same as in (B). E Violin Plot of pRPS6 ADT quantification in Rods and Cones. F Confocal images of 33-week-old retinal organoids harvested and immunostained with antibodies against the cone marker ARR3 (green), Rod marker RHO (blue), and pRPS6 (red) with examples of high-pRPS6 cones indicated by arrowheads. Scale bars are 50 μm. This experiment was performed once. G Quantification of ARR3, RHO, and pRPS6 in 240 Rods and 64 Cones from (F).

Fig. 3. Extending Phospho-seq with transcriptomic measurements.

A Schematic of Phospho-seq using the 10x multiome kit. B Weight-nearest neighbors UMAP of Phospho-seq-Multi experiment labeled by cell type. AC/HC amacrine cells/horizontal cells. C Heatmap comparison of 11 ADTs between Phospho-seq experiments using either the Multiome kit or the ATAC-seq kit. D Schematic of bridge integration procedure. E UMAP representation of retinal organoid cells based on ATAC-seq modality with cell type assignments from bridge integration. BC bipolar cells. IN intermediate. F Alluvial plot demonstrating cell label transfer when using bridge integration.

Fig. 4. Phospho-seq on brain organoids.

A Schematic of antibody panel used in brain organoid Phospho-seq experiment and the cellular compartment of the target protein. B Schematic of brain organoid differentiation. C UMAP representation of cells and cell type assignments based on the ATAC-seq modality in brain organoid Phospho-seq. GE ganglionic eminence. D Coverage and violin plots of the gene promoters and protein levels respectively of nuclear protein OTX2 and cytoplasmic protein S100B, color and order are same as in (C). E Diffusion MAP plot of diencephalic and telencephalic differentiation trajectories. Cell-type colors correspond to legend in (C). F Heatmap of ADT expression of proteins and phospho-proteins associated with MAPK/ERK and mTOR signaling. G Schematic of MAPK/ERK signaling. H Coverage plot of RPS6KA1 split into quantiles of pRPS6 levels across telencephalic and diencephalic differentiation trajectories with OTX2 binding motifs indicated. I Coverage plot of RPS6KB1 split into quantiles of pRPS6 levels across telencephalic and diencephalic differentiation trajectories. J Rank-correlation plot showing correlation between pRPS6 levels and motif accessibility across the whole Phospho-seq dataset. The top hits are indicated in red. K Rank-correlation plot showing correlation between pMAPK levels and motif accessibility across the whole Phospho-seq dataset.

Fig. 5. Identification of GLI3 transcriptional co-factors.

A GLI3 protein levels correlate with the accessibility of HD/2 motifs (blue), but not the canonical GLI3 motif (red). Shown is a rank-correlation plot of all ADT/motif correlations. B GLI3 protein levels correlate with GLI3 RNA levels. Rank-correlation plot shows GLI3 ADT correlation with all genes. C Upon GLI3 perturbation, differentially accessible regions (DAR) are enriched for HD/2 motifs. Shown is the Rank-Fold Enrichment plot of motifs enriched in differentially accessible GLI3-bound regions in GLI3 knockout vs WT (left) and motifs enriched in non-differentially accessible regions (right). D Upon GLI3 perturbation, only genes linked to DAR are enriched for expression changes, which are predominantly negative as expected. Shown is a density plot of the fold change of genes linked to DAR and non-DAR. ** = p = 3.89 × 10⁻²³ using a two-sided t-test. Phospho-seq data in WT cells predicts which regions will respond to functional perturbation. E, F show two exemplary loci, both bound by GLI3 in CUT&Tag data. The region in (E) contains HD/2 and GLI3 motifs. In the GLI3 KO data, the region is differentially accessible (middle), and in the Phospho-seq data, the accessibility is correlated with GLI3 protein expression. The region in (F) is also bound by GLI3 but contains no HD/2 motifs, is not associated with GLI3 ADT levels in Phospho-seq, and is not DAR upon GLI3 perturbation. Density plot in (G) shows that this pattern holds genome-wide: the correlation between peak accessibility and GLI3 ADT levels in Phospho-seq data discriminates responsive (DAR) and non-responsive (nonDAR) peaks in GLI3 KO cells ** = p = 2.2 × 10⁻¹⁶ using a two-sided t-test.

Fig. 6. Gene regulation inference with Phospho-seq.

A Schematic approach to use metacelling to discover cis-regulatory elements associated with individual proteins. B Tn5 cut-site footprinting between cells with high OTX2 expression and low OTX2 expression in peaks that are highly correlated with OTX2 (left panel) and uncorrelated with OTX2 (right panel). C Example of an inferred transcription activating peak associated with indicated proximal gene expression for OTX2. Coverage plots and violin plots are ordered by quantile ADT expression. Red lines indicated the location of an OTX2 binding motif. D Example of an inferred transcription repressing peak associated with indicated proximal gene expression for OTX2. Coverage plots and violin plots are ordered by quantile ADT expression. Red lines indicated the location of an OTX2 binding motif. E Coverage and violin plots of GAS7 split by quantiles of pRPS6 levels across the whole dataset with TEAD1 motifs indicated by red lines. F Density plot showing ChIP-Seq depth at all peaks with a TEAD1 motif compared to pRPS6:TEAD1 peaks identified in this study. The significance of p = 1.9 × 10⁻¹⁵ is determined by a Wilcoxon rank-sum test (two-sided). G Violin plot of normalized pSTAT3 ADT levels in a subset of progenitor and glial cells. H Violin Plot of STAT3 motif chromVAR scores in a subset of progenitor and glial cells. I Violin plot of STAT3 activated gene module scores in a subset of progenitor and glial cells.