Multiplexed protein profiling reveals spatial subcellular signaling networks.

Protein-protein interaction networks are altered in multi-gene dysregulations in many disorders. Image-based protein multiplexing sheds light on signaling pathways to dissect cell-to-cell heterogeneity, previously masked by the bulk assays. Herein, we present a rapid multiplexed immunofluorescence (RapMIF) method measuring up to 25-plex spatial protein maps from cultures and tissues at subcellular resolution, providing combinatorial 272 pairwise and 1,360 tri-protein signaling states across 33 multiplexed pixel-level clusters. The RapMIF pipeline automated staining, bleaching, and imaging of the biospecimens in a single platform. RapMIF showed that WNT/β-catenin signaling upregulated upon the inhibition of the AKT/mTOR pathway. Subcellular protein images demonstrated translocation patterns, spatial receptor discontinuity, and subcellular signaling clusters in single cells. Signaling networks exhibited spatial redistribution of signaling proteins in drug-responsive cultures. Machine learning analysis predicted the phosphorylated β-catenin expression from interconnected signaling protein images. RapMIF is an ideal signaling discovery approach for precision therapy design.

Introduction

Non-small-cell lung cancer (NSCLC) comprises 80%–85% of lung cancer, the second leading cause of death after cardiovascular disease (Stewart, 2014). The survival rates remain unimproved due to cancer relapse despite the increasing number of signaling treatments of NSCLCs (Hommura et al., 2002). Metastatic NSCLCs exhibit intrinsic resistance to chemotherapy through compensatory pathways (Stewart, 2010). Deciphering the pathogenesis of signaling pathways contributing to drug resistance can help design efficient drug combinations for enhancing targeted therapies. Protein-protein interaction (PPI) networks have shown the intricate role of multiple signaling pathways in NSCLCs (Bu et al., 2020). However, these measurements have been limited to bulk-level molecular assays. Besides, affinity purification followed by mass spectrometry has been used to study the multi-protein complexes in a high-throughput manner, but visual potential interaction maps were lacking (De Las Rivas and Fontanillo, 2012; Westermarck et al., 2013). There is a crucial need for an image-based signaling assay to study subcellular PPIs and distributions in NSCLCs at the single-cell level.

WNT-1 is overexpressed in NSCLC to overcome the cancer-associated defects in apoptosis through WNT pathways (Chen et al., 2001; He et al., 2004). One of the WNT pathways, canonical or β-catenin-dependent WNT signaling, functions in mediating cell adhesion and cell growth. In the absence of WNT ligands, a β-catenin destruction complex is formed consisting of axis inhibition protein (AXIN), adenomatous polyposis coli (APC), and glycogen synthase kinase 3β (GSK-3β). This destruction complex can lead to the degradation of β-catenin by phosphorylating it and blocking the downstream signaling transduction. However, when the WNT ligand is present, it binds to members of the Frizzled (FZD) receptors, forming a complex between WNT, FZD, lipoprotein receptor-related protein (LRP), disheveled (DVL), and AXIN, leaving β-catenin unphosphorylated. The unphosphorylated β-catenin can migrate to the nucleus and activates the transcription of target genes, including cyclin D1, cyclin E, c-MYC, etc (Stewart, 2014). DKK1 and DKK2 as negative regulators of the WNT/β-catenin pathway play complex roles in tumor development. The mechanism of how DKK1 affects the β-catenin-independent WNT pathway in NSCLC remains unclear. It has been found that decreased DKK1 can inhibit cell migration and invasion, confirming its role in promoting NSCLC cell EMT via the Wnt/β-catenin pathway. The stabilization of β-catenin upon WNT activation can promote the epithelial-mesenchymal transition (EMT) process (Zhang et al., 2017). DKK2 as another antagonist of the WNT pathway by binding to LRP5/6 could promote tumor progression by suppressing cytotoxic immune cell activation in APC-mutated cells in colorectal cancer and NSCLC (Shen et al., 2019; Xiao et al., 2018). However, the location of DKK2 is still unclear. RNF43 also functions as a negative regulator of the WNT pathway by ubiquitinating and degrading WNT receptors (Hao et al., 2016); however, it may also promote EMT by inhibiting the function of E-cadherin in NSCLC (Zhang et al., 2019). EMMPRIN has been discovered as a tumor promoter in stimulating cell migration and invasion of tumor cells by increasing the expression of MMPs, and EMMPRIN may upregulate the WNT/β-catenin pathway by destabilizing E-cadherin/β-catenin complex to increase free-β-catenin in the cytosol. However, the role of EMMPRIN in cell-cell contacts through the AKT/mTOR pathway remains elusive (Sidhu et al., 2010). In the WNT/β-catenin pathway, the phosphorylation of GSK3β can lead to the phosphorylation of β-catenin and subsequent degradation. In the AKT/mTOR pathway, GSK3β behaves differently. The phosphorylated AKT can inhibit GSK3β kinase activity by phosphorylating it at Ser9 (Duda et al., 2020; McCubrey et al., 2014a; Zheng et al., 2020). Even though some studies showed no direct relationship between the two pathways, the simultaneous inhibition of GSK3beta activity induced by the two pathways may indicate the role of both pathways in promoting β-catenin accumulation.

Another commonly activated pathway in NSCLC is the PI3K/AKT/mTOR pathway. In 10%–20% NSCLC, the epidermal growth factor receptor (EGFR) mutations can lead to constituent activation of the downstream AKT/mTOR cascade, resulting in promoting oncogenesis through enhancing cell proliferation, growth, and survival. Besides, this pathway mediates the resistance to EGFR tyrosine kinase inhibitors (TKI) (Cheng et al., 2014). There are two main dominant downstream pathways of EGFR, including AKT/mTOR controlling apoptosis and the MEK/ERK pathway controlling cell proliferation (Siegelin and Borczuk, 2014). The EGFR-mutated cell lines have shown resistance to chemotherapeutic agents; however, the response of EGFR-mutated cells to EGFR TKIs became sensitive by inhibiting the AKT-mediated pathway, indicating the role of the AKT pathway in regulating the drug resistance of EGFR-mutated cell cultures (Raz et al., 2006). WNT/β-catenin and PI3K/AKT/mTOR pathways share some common effectors, such as GSK3β (Duda et al., 2020). The function of both pathways is widely described; however, the interaction between two signaling pathways is still unclear (Prossomariti et al., 2020). The downstream signaling proteins of both pathways play a heterogeneous role, further complicating signaling cascades. Thus, unraveling the complexity of the PPIs or crosstalk between pathways at the subcellular level is key to identifying the signaling networks that regulate cellular function.

Osimertinib (AZD9291), a third generation of irreversible TKIs, has become an approved drug for both first-line and second-line treatment. Osimertinib has shown more efficacy than the previous generation EGFR-TKIs in treating EGFR mutant NSCLC. Patients receiving osimertinib exhibited significantly longer progression-free survival, resulting in a trend of longer overall survival and better safety profiles than patients treated with either gefitinib or erlotinib (Ramalingam et al., 2019). However, a subset of patients receiving osimertinib has developed relapse and drug resistance. MEK/ERK signaling pathway as the downstream of EGFR has been studied to delay or prevent the acquired resistance to osimertinib by using western blot, flow cytometry, cell viability, and apoptosis assays (Gu et al., 2020a). Other inhibitors, such as activated Cdc42-associated kinase 1 (ACK1), have complemented therapies to potentially overcome the acquired resistance (Gu et al., 2020b). However, the bulk level protein analysis fails to demonstrate the cell-to-cell heterogeneity in responding to osimertinib or acquiring the resistance. The role of signaling crosstalk in cell proliferation and survival under osimertinib treatment is unknown, mainly because of the crosstalk between the WNT/β-catenin and AKT/mTOR pathway as a key downstream of EGFR. Therefore, multiplexing proteins in signaling pathways are needed to evaluate the spatial dynamics between the two pathways in response to osimertinib.

Highly multiplexed protein profiling methods are emerging for deciphering the complexity of the tumor microenvironment at the single-cell level. Such technologies include multiplexed ion beam imaging (MIBI) and imaging mass cytometry (IMC) utilizing a time-of-flight (TOF) instrument to detect multiple markers simultaneously by measuring metal and rare earth isotopes tagged to antibodies. However, these tissue-level imaging techniques are limited to lower spatial resolution and require specialized equipment with a low scanning speed (Baharlou et al., 2019; Giesen et al., 2014). Similarly, CO-Detection by indEXing (CODEX) is another multiplexed imaging technique that employs cyclic addition and removal of fluorescently labeled DNA-barcoded antibodies to enable deep protein phenotyping using 50 or more markers in tissues (Schürch et al., 2020). However, CODEX exhibits high complexity in conjugating DNA to the antibodies of interest, making it challenging to develop custom protein panels for signaling targets. A fast and straightforward multiplexed protein assay is lacking for dissecting the signaling transduction and crosstalk at the single-cell level.

To provide a solution to this important need, we demonstrate a rapid multiplexed immunofluorescence (RapMIF) technique that can quickly screen a large number of protein targets in single cancer cells using an automated stainer, yielding a subcellular protein network of signaling pathways. The RapMIF method uses an iterative labeling process to expand the molecular targeting palette of single-cell imaging approaches, akin to sequential sample processing manually performed in cyclic immunofluorescence (CycIF) (Lin et al., 2015, 2016, 2018). The RapMIF method overcomes the limitation of using 4–6 channels to study the crosstalk between fluorophores by conventional immunofluorescence (IF) imaging. In each cycle, the fluorescence images of the proteins of interest are repeatedly collected from the same sample. The processed RapMIF data visualized and quantified up to 25-plex signaling markers involved in the WNT/β-catenin and AKT/mTOR pathways, together with additional epigenetic markers to assess histone modifications (Barlési et al., 2007; Leng et al., 2020; Zhu et al., 2018). The RapMIF visualized the physical, regulatory, and functional protein associations of WNT/β-catenin and AKT/mTOR pathways and provided the signaling protein’s subcellular location and expression level in NSCLC cultures and lung adenocarcinoma tissue samples. RapMIF experiments of cultures treated with EGFR TKIs demonstrated the signaling responses of two compensatory pathways. Spatially resolved signaling networks have yielded alterations of AKT and WNT pathways in lung cancer cells when treated with osimertinib. RapMIF then was used to profile formalin-fixed, paraffin-embedded (FFPE) samples to study the signaling proteins’ subcellular localization, translocation, and phosphorylation states in a large tissue microarray. Signaling pathways coordinate many protein functions that involve phosphorylation or activation of proteins, cell-cycle regulation, transcription factor promotion or suppression, downstream transductive activation, and translocation of effectors. Therefore, multiplexed protein maps are imperative to identify the distinct cell states, cell-to-cell heterogeneity, and the complexity of signaling dynamics. Image-based protein methods are ideally suited to extract the subcellular protein distributions and translocation dynamics of the signaling networks.

Results

Rapid multiplexed immunofluorescence for subcellular spatial protein analysis

To resolve the subcellular spatial protein networks in the WNT/β-catenin and AKT/mTOR pathways, the presented RapMIF method was performed on A549 cells, EGFR-mutant PC9 cells, and lung adenocarcinoma tissues. RapMIF technology has been built upon the validated protocols of the CycIF method that typically utilizes robust fluorophore inactivation and repeated labeling of dyes to overcome the limitations of the microscopes constrained with 4–6 color channels in conventional IF (Lin et al., 2016). The RapMIF approach incorporates an autostainer setup and pre-conjugated antibodies to expedite the sample processing for highly multiplexed imaging of cells or tissues located on coverslips (Figure 1A). Specifically, the RapMIF reduces the sample staining time from 12–16 h to 1 h. The RapMIF achieved imaging of up to 25 markers and revealed the signaling proteins as part of the WNT/β-catenin and AKT/mTOR pathways using quantitative expression profiles and spatial distributions of target biomarkers (Figure 1B). The spatial signaling networks were altered upon osimertinib treatment in PC9 cells (Figure 1C). The RapMIF approach on tissue samples yielded signaling alterations in complex tumor microenvironments (Figure 1D).

Figure 1

Rapid multiplexed immunofluorescence (RapMIF) for subcellular spatial protein analysis

(A) A schematic setup of a Rapid, automated, and multiplexed imaging system was presented. The samples were stained with multiple cycles of IF using an autostainer. Created with BioRender.com.

(B) The left panel indicates the activated Wnt and AKT signaling pathways; the markers highlighted in red were included in our panel. The right diagram demonstrates the predicted PPIs among these proteins. Created with BioRender.com.

(C) Predicted expression of signaling markers in PC9 cells was provided from the multiplexed protein data. Created with BioRender.com.

(D) Schematic of the tissue samples stained with pan-cytokeratin used for tumor and stroma analysis of spatial signaling networks. Created with BioRender.com.

(E) An illustration of the data processing and analysis pipeline was shown. Multiplexed images provided the nuclear and cytosolic masks in A549 cells, correlation plot, protein prediction, the p-EGFR spatial intensity plot, and the pixel clustering maps of an individual PC9 cell. Created with BioRender.com.

For the multiplexing experiment, RapMIF utilized a panel of 11 cycles, resulting in a total of up to 25 markers, including 5 structural, 17 signaling, and 3 epigenetic markers. For robust segmentation, the panel included Concanavalin A (Endoplasmic reticulum), Phalloidin (Actin), and wheat germ agglutinin/WGA (Golgi and plasma membrane) as cell painting markers (Bray et al., 2016) and β-actin (Actin filaments) and β-tubulin (Microtubules) to accurately delineate nucleus and cell body using an automated algorithm. To determine the best staining conditions, we performed titration of each antibody in the multiplexed panel. IF staining evaluated antibody staining using two different dilution rates (Figure S1). We compared the signal-to-noise ratios (SNRs) of these IF images, yielding the condition with higher SNR used for RapMIF experiments (Table S1). The staining dilutions were consistent with the manufacturers’ suggested dilution range. To reduce the time of sample incubation in indirect IF, most antibodies were conjugated with one of the 488, 555, and 647 Alexa dyes using Lightning-Link Rapid Alexa Fluor Antibody Labeling Kit (Novus/Abcam). These markers were evaluated in four conditions based on conjugated antibodies, unconjugated antibodies, overnight incubation at 4°C, and 1-h incubation at room temperature (RT) (Figure S2). The comparisons of antibody staining experiments at overnight and RT temperature incubation yielded similar signaling patterns, denoted in the cross-sections of plots across pixel histograms (Figure S3). For each round of staining, the sample was incubated with antibodies conjugated to distinct dyes for 1 h at RT after blocking with cell staining media (CSM). Following imaging of labeled cells, the dyes were bleached in the NaOH-H2O2 mixture under white light for 1 h (Figures S4). Although the base-hydrogen peroxide mixture bleaches Alexa 488, 555, and 647 at different rates, a common incubation for 30–45 min within a fluorophore inactivation buffer was sufficient (Lin et al., 2018). Post-bleaching images were acquired after the inactivation of the dyes before the next round of staining (Figure S5). To confirm the bleaching efficacy, quantitative analysis was performed on post-bleach images of the A549 cells (Figures S6 and S7). Negligible cell loss was noticed in the first and second cycles in the sample, corresponding to the 2%–5% and 0%–2% of cells per cycle, subsequently (Lin et al., 2018). Following fluorophore inactivation, the coverslip was re-blocked with CSM to make the cells ready for the next round of staining. The details of customized parameter settings for fluid transfer and antibody incubation are included in Table S2. The 11-cycle RapMIF experiment yielded up to 25 protein-type images on the stained slides at the single-cell level. The multiplexed single-cell measurements revealed the subcellular distribution of signaling protein networks. Compared with CycIF, RapMIF has reduced the cycle time from 12 h to 5 h (Table S3). This prevents cell degradation, a common issue for the multiplexed experiments performed over long-time periods. To maximize the efficiency of epitope detection in tissue, overnight incubation for staining was performed, consistent with the settings in t-CycIF. Current experiments rely on the commercially available tissue on 1-mm thick glass slides; however, tissue can be sectioned on the thin coverslips to be compatible with the acrylic holder used for RapMIF automation (Figure S4D).

To reduce the need for conjugated antibodies for each marker, we combined both indirect labeling and direct labeling in our panel. To minimize the effect of unspecific binding of secondary antibodies to conjugated primary antibodies, the indirect labeling of p-EGFR was designed in the first cycle. Indirect labeling of β-tubulin and β-actin allowed cell segmentation, and pan-cytokeratin distinguished the cancerous tissue regions in the last cycle. The images were taken at 40X with a step size of 0.4 μm to reveal the slowly varying signaling molecules across the z axis of a cell (Figure S8). The best focus images were generated from +/− 2 z-stacks, and they were used for post-processing of the single cells (Figures S9A and S9B). To avoid the effect of the dye inactivation solution on Phalloidin-labeled structures, p-EGFR targeting antibodies were directly incubated between the first and second cycles for the cell and tissue experiments without an intermediate bleaching step.

To study the spatial distributions of proteins, we performed a single-pixel analysis to demonstrate the nuclear and cytosolic expression levels and translocations of signaling proteins in WNT/β-catenin and AKT/mTOR pathways. The measurement of the expression level of each marker in two regions, cytosol and nucleus, was quantified separately (Figure 1E). Pixel-level colocalization of signaling proteins was then performed by unsupervised clustering of the protein images in subcellular regions, revealing the spatially distinct signaling neighborhood regions (Figure 1E). Machine learning (e.g., Random forest) analysis of subcellular maps and expression profiles of the signaling targets yielded the phosphorylation activity of four markers, p-AKT, p-mTOR, p-EGFR, and p-β-catenin, demonstrating the interdependency of spatial signaling networks that are altered in wild-type and drug-perturbed cells (Figure 1E). Together, the RapMIF approach resolved the cell-to-cell heterogeneity in expression profiles and spatial signaling maps and deciphered the phosphorylation activities of the kinases and their direct role in regulating signaling pathways.

Highly multiplex spatial decoding of WNT/β-catenin and AKT/mTOR in A549 cells

To optimize the staining protocols and data analysis pipeline, RapMIF was performed to quantify the expression and spatial distribution of signaling markers involved in WNT/β-catenin and AKT/mTOR pathways in the A549 NSCLC cell line (Figure 2A). The raw IF images (Figure 2B) were processed for cell segmentation and quantification of the protein expression levels (Figures 2C, 2D, S9C and S10). To consider the cell morphology difference, the mean intensity was used for correlation analysis (Figures 2E and 2F), and total intensity correlation also showed the location-dependent intensity comparison (Figures S11 and S12). The analysis of WNT/β-catenin and AKT/mTOR pathways in stabilizing β-catenin demonstrated that non-p-β-catenin was positively associated with the expression of cyclin E and cyclin D1 (R = 0.6–0.7). This finding indicated non-p-β-catenin’s role in the progression of the cell cycle (Figures 2E, S11 and S12). The co-expression of cytosolic protein distributions was detected in non-p-β-catenin, AKT, and p-AKT signaling markers (Figures 2F, S11B and S12B). We speculate that the p-AKT can promote the activity of β-catenin by inhibiting GSK3β, one member of the β-catenin degradation complex (McCubrey et al., 2014b). The stabilization of β-catenin by p-AKT could then facilitate the activity of cyclin E and cyclin D1 in the nucleus. The pairwise analysis between two proteins determined their correlation in expression levels, yielding 272 pairs in cytosol and nucleus of individual cells. We also showed the possibility of using 3-way intensity correlations to illustrate the 3-protein neighborhood in single cells (Figure S13). The expression level of cytosolic p-mTOR proteins yielded a higher correlation with cytosolic p-β-catenin than with cytosolic non-p-β-catenin in the A549 cell line (Figure S12B). This pattern may be due to the role of mTORC1 in suppressing the WNT/β-catenin pathway by inhibiting LRP6 phosphorylation and FZD activation (Zeng et al., 2018). To confirm the removal of fluorescence signals, we demonstrated the marker correlations of single cells from the post-bleach images (Figure S14).

Figure 2

Expression profiles and correlation analysis of multiplexed signaling proteins in A549 cells

(A) The correlation analysis pipeline of protein intensity at the single-cell level. Created with BioRender.com.

(B) Raw staining images in A549 cells from 23 protein markers and DAPI from 11 cycles were shown.

(C) A workflow of data analysis was presented. The raw images were segmented and masked using WGA, concanavalin A, WNT1, APC, and p-EGFR. The signaling intensity was quantified based on the cytoplasmic and nuclear masks.

(D) Normalized mean intensity per cell of 25 markers for one ROI. Sixty-three regions of interest and a total of 3,457 numbers of A549 cells were analyzed. The normalized total intensity was shown in Figure S10.

(E) Pairwise marker correlation based on single-cell mean intensity in A549 cells was evaluated using Pearson Correlation and average linkage method on Euclidean distance. The correlation of total intensity per marker per cell was shown in Figure S11, and the correlation of total intensity was shown in Figure S12. The pairwise Pearson correlations were corrected by the Holm-Šídák method, yielding significance with a p value of p<=0.0001 (∗∗∗∗).

(F) Pairwise marker correlation based on single-cell mean intensity in the cytoplasm and nucleus-dependent manner was evaluated using Pearson Correlation using average linkage method on Euclidean distance from the same cell as (D). The pairwise Pearson correlations were corrected by the Holm-Šídák method, yielding significance with a p value of p<=0.0001 (∗∗∗∗).

To capture the spatial information among the two pathways, signaling protein maps of the multiplexed panel were converted into 33-pixel clusters and uniquely colored back onto the original images of the single cancer cells (Figures 3A–3D). Each pixel cluster corresponded to a set of signaling markers co-expressing in a pixel distribution based on similarities of multiplexed intensity profiles (Gut et al., 2018). Background pixels were filtered out based on semi-automatically gated parameters (Figure S15). The UMAP visualizes pixel-wise clusters in A549 cells (Figure 3B). We observed that cluster 13 has relatively higher expression levels of WNT1 and p-AKT, indicating the potential neighborhood between the two pathways (Figure 3D). The pixel clustering map demonstrated that non-p-β-catenin was distributed in the nucleus in cluster 2, whereas another portion of the signal was detected in the cytosol and along the cell membrane in cluster 16 (Figure 3D). Such subcellular signaling patterns may indicate non-p-β-catenin functions as a component of cell-cell adhesion structures by interacting with E-cadherin and actin (Kase et al., 2000). Besides, p-β-catenin (Ser45) was highly expressed in the cytosol in cluster 21 of A549 NSCLC cells (Figure 3B), supporting the degradation of the phosphorylated β-catenin in the cytosol, but not in the nucleus. Cytosolic p-β-catenin degradation is dependent on the phosphorylation at Ser33/Ser37/Thr41 by GSK3β (Li et al., 2012; Parker et al., 2020). The p-β-catenin signaling has previously been found almost exclusively expressed in the nucleus in melanoma and colorectal cancer cells, whereas p-β-catenin has been detected in both cytoplasm and nucleus in invasive breast carcinomas (Chung et al., 2001; Kielhorn et al., 2003; Sinnberg et al., 2011). The presence of p-β-catenin in the nucleus may be due to the overexpression of p-β-catenin. The overexpressed p-β-catenin weighs much larger than the proteins undergoing degradation, thereby leading to the translocation to the nucleus (Nakopoulou et al., 2006). The p-EGFR (Tyr1086) signaling proteins were mainly distributed in cytosol shown in clusters 14 and 32, whereas some activity was detected in the nucleus denoted in cluster 2 (Figure 3D). This translocation pattern is possibly due to the activation of EGFR at different phosphorylation sites. This can facilitate the translocation of EGFR to the nucleus, whereas nuclear EGFR functions differently from the EGFR on the cell membrane, for instance, to promote the transcription of cyclin D1 (Brand et al., 2011, 2013).

Figure 3

Spatial signaling network analysis in single pixels and prediction of phosphoproteins in A549 cells

(A) Pixel clustering analysis pipeline was demonstrated. Created with BioRender.com.

(B) Spatial signaling maps of 33 clusters at the single-pixel level were provided. The UMAP representation of the distribution of these clusters in A549 cells was shown. The legend showed the cluster name in a frequency-descending order for each cell.

(C) A heatmap of 33 clusters on a Z score scale was shown. Each cluster represented one distinct expression profile of 21 protein markers in A549 cells. Data are represented as mean expression per cluster.

(D) Pixel-level clustering of signaling proteins in one ROI in A549 cells was demonstrated. The highly expressed signaling proteins of different clusters were indicated in each cancer cell.

(E) The importance of 20 protein markers in A549 cells was evaluated based on the random forest algorithm model mean decrease in impurity (MDI). The MDI was defined as the total decrease in node impurity (weighted by the probability of reaching that node averaged over all trees of the ensemble in the random forest algorithm.)

(F) The expression level of p-β-catenin in A549 cells was predicted based on 20 features (R = 0.75, p value <0.001).

(G) Comparison between the intensity level of the raw p- β-catenin and the predicted intensity in one ROI was obtained based on random forest analysis in A549 cells.

(H) The quantification of cyclin D1 expression across DAPI expression level per cell in A549 cells. Most of the cells were harvested in the G1 phase. The right images show the comparisons of the cyclin D1 expression level in three cells in three different phases. Cell number 25 has a relatively higher cyclin D1 intensity than cell 23 and cell 28.

We also observed the colocalization of WNT/β-catenin and p-AKT pathway components in the spatial signaling network maps, demonstrating the physical neighborhood between two pathways in the same pixel area (Figure 3D), providing the common downstream effectors. Both pathways can cause the stabilization and accumulation of β-catenin and increase the activity of cyclin D1 and cyclin E to facilitate the progression of the cell cycle at the single-pixel level. In addition to the indirect stabilization of β-catenin by p-AKT via GSK3β, p-AKT and non-p-β-catenin may also have potential direct neighborhood due to the co-expression and colocalization of the two proteins (Figures 2E and 3D). AKT can phosphorylate β-catenin at Ser552 directly, causing β-catenin to dissociate from cell contacts, and promote its translocation from cytosol to nucleus to increase transcriptional activity and tumor cell growth (Fang et al., 2007).

Next, we designed a machine learning model to predict the β-catenin phosphorylation activity based on the 20 biomarkers using a random forest (RF) algorithm (Figures 3E–3G). The activity of β-catenin is highly correlated with cell growth, proliferation, and migration. There is no reliable biomarker to predict the activity of β-catenin. In our prediction model, the subset of signaling proteins including WNT1, RNF43, and AKT exhibited relatively high importance in predicting the phosphorylation of β-catenin, yielding an r = 0.76 fit (p < 0.001) in prediction accuracy. The same RF model provided predictions of other phosphorylated markers, including p-AKT, p-mTOR, and p-EGFR (Figure S16). The multiplexed dataset revealed potential biomarkers as features to predict the phosphorylation activity of β-catenin.

We also conducted cell-cycle analysis based on the quantification of nuclear cyclin D1 expression level, classifying the cells into G1, S, and G2 phases (Figure 3H). Most of the cells were harvested in the G1 phase. This indicated the regulatory role of cyclin D1 in hyper-phosphorylating the retinoblastoma tumor suppressor protein (Rb1) by forming cyclin D1/CDK4 complex, necessary for the cell to enter the S phase in the cell cycle (Gautschi et al., 2007). The number of cells expressing cyclin D1 was low during the S and G2 phases due to the exportation of cyclin D1 from the nucleus to the cytoplasm for degradation by the proteasome (Guo et al., 2005). Similarly, cell-cycle analysis was assessed based on cyclin E expression levels. Cyclin E demonstrated similar results that most cells expressing cyclin E were harvested in G1, and some in S phases (Figure S17), indicating its role in promoting G1 and S phases for cell proliferation and growth by binding to CDK2 (Hwang and Clurman, 2005). In our study, cyclin E accumulated around the nuclear membrane and was mainly distributed in the cytosol (Figure 3D). The accumulation of cyclin E in the cytoplasm may be the result of the formation of low-molecular-weight cyclin E (LMW-E) with the loss of the canonical NH2-terminal nuclear localization sequence. The overexpression of LMW-E is associated with breast cancer recurrence and poor overall survival; it has also been observed in NSCLC (Hunt et al., 2017; Koutsami et al., 2006). Besides, we also looked into the role of DKK1 in NSCLC. DKK1 plays multiple roles in different diseases; it functions as an oncogene in hepatocellular carcinoma, liver cancer, prostate cancer, and myeloma (Kagey and He, 2017; Zhang et al., 2017). However, it can also promote epithelial-mesenchymal transition (EMT) by stabilizing β-catenin and inducing its nuclear localization (Zhang et al., 2017). In the RapMIF data, the nuclear DKK1 proteins exhibited a higher correlation to non-p-β-catenin than p-β-catenin in both cytosol and nucleus (Figure S12B). This reveals that the nuclear DKK1 may act as a promoter in the WNT/β-catenin pathway. However, cytosolic DKK1 yielded a similar correlation to cytosolic p-β-catenin and non-p-β-catenin (Figure S12B). The modulation of DKK1 in cancer cell development may not be limited to WNT signaling pathways. Cytoskeleton-associated protein 4 (CKAP4) has been found as a key receptor of DKK1. The binding of DKK1 to CKAP4 could activate PI3K/AKT pathways, increasing cell proliferation (Kagey and He, 2017). Therefore, the role of DKK1 in the WNT/β-catenin pathway may be location dependent, and further studies are needed to decipher the function of DKK1 in NSCLC.

Spatial networks in EGFR mutant cells for analysis of EGFR, WNT, and AKT-mTOR crosstalks

To reconstruct the signaling networks in EGFR mutant cells under drug perturbations, the RapMIF was used to profile multiple signaling molecules in EGFR mutant NSCLC-sensitive PC9 cell lines. Osimertinib has shown a prominent effect in inhibiting PC9 cells, and the cell survival rate has decreased, as the concentration has increased from 0 to 10 nM under a 3-day treatment (Gu et al., 2020a, 2020b). PC9 cells responding to osimertinib treatment have been shown to develop acquired resistance. However, in previous work, the western blot experiments of signaling proteins in PC9 cells treated with the combination of osimertinib and either ACK1 inhibitor or MEK inhibitor have demonstrated effectiveness to overcome the acquired resistance to osimertinib (Gu et al., 2020b). Besides, it has been shown that osimertinib can effectively inhibit the p-AKT (S473) level in PC9 cells but not AKT itself. Similar results have been obtained from the signaling studies in H1975 cells, a drug-resistant cell line (Gu et al., 2020a). Osimertinib has inhibited the phosphorylation of EGFR in a phosphorylation-site-dependent manner and decreased the p-EGFR (Y1173) level but failed to alter the p-EGFR (Y1068) in PC9 cells (Gu et al., 2020a, 2020b). However, osimertinib could affect the downstream effector of EGFR, Grb2, by suppressing its binding to p-EGFR (Tyr1086) in H1975 cells (Lee et al., 2018).

In our study, we first asked whether the EGFR pathway would be inhibited by osimertinib in PC9 cells, eventually affecting the proliferation and growth of the PC9 cells. A 48-h osimertinib treatment reduced the number of adherent PC9 cells compared with untreated ones. Higher drug concentrations further decreased the cell viability compared with lower drug treatments (Figures S18-S20). Following drug treatments, the RapMIF was used to generate the multiplexed imaging maps of up to 25-plex signaling proteins in PC9 with and without the treatment of 40nM osimertinib (Figures 4A, 4B, and S21). The drug treatment disrupted spatial continuity in the p-EGFR distributions on the cell membrane (Figures 4A, 4B, S22 and S23). This finding may be due to the effect of osimertinib on inhibiting the dimerization and phosphorylation of EGFR. Furthermore, the p-EGFR signaling proteins were highly expressed in both cell membrane and nucleus in the absence of osimertinib. However, p-EGFR was mainly localized in the cell membrane, and the nuclear p-EGFR expression was relatively lower under the drug treatment (Figures 4A and 4B). Similar trends were observed in p-EGFR validations across four different drug concentrations (Figures S22 and S23). As a positive control, β-tubulin exhibited uniform distributions across the cytoplasm of the cell with and without drug treatments. In addition to staining in the cytosol, the β-tubulin expression level was enriched in the nucleolus (Figure S23) due to the role of β-tubulin in the mitotic spindle in helping nucleolus division at mitosis or nonspecific staining in the nucleolus (Jouhilahti et al., 2008; Walsh, 2012). The nuclear p-EGFR proteins are typically shed from the cell membrane in the sample without osimertinib. The translocation of p-EGFR occurs via a nuclear importing system, and the nuclear p-EGFR can promote cell proliferation in an EGF-dependent manner to control the cyclin D1 gene expression (Lin et al., 2001; Nishimura et al., 2008). Also, previous literature suggested that the p-EGFR yields a high correlation to the cells with high proliferation in native tissues (Hoshino et al., 2007; Lin et al., 2001). Thus, the osimertinib treatment not only inhibited the phosphorylation of p-EGFR at the cell membrane but also reduced the translocation of p-EGFR to the nucleus, impeding cell proliferation.

Figure 4

The effect of osimertinib on signaling expression profiles and correlation analysis in PC9 cells

(A) Raw staining images in PC9 cells from 23 protein markers and DAPI from 12 cycles were presented. The cells were incubated in cell media for 48 h and seeded on coverslips for multiplexing IF. n = 531 number of cells were utilized.

(B) Raw staining images in PC9 cells from 23 protein markers and DAPI from 12 cycles were shown. The cells were treated with 40nM osimertinib for 48 h and seeded on coverslips for multiplexing IF. n = 451 number of cells were used.

(C) The correlation of the single-cell mean intensity of 23 markers for control PC9 was evaluated using Pearson Correlation using the average linkage method based on Euclidean distance. The pairwise Pearson correlations were corrected by the Holm-Šídák method, yielding significance with a p value of p<=0.0001 (∗∗∗∗).

(D) The correlation of the single-cell mean intensity of 23 markers for PC9 cells treated with 40nM osimertinib for 48 h was evaluated using Pearson Correlation based on the average linkage method implemented on Euclidean distance. The pairwise Pearson correlations were corrected by the Holm-Šídák method, yielding significance with a p value of p<=0.0001 (∗∗∗∗).

(E) The p-EGFR (Tyr1086) intensity in PC9 cells across four-drug concentrations, 0, 20, 40, and 60 nM was analyzed and normalized to β-tubulin. The cells were treated with drug or cell media for 48 h and seeded on coverslips. The cells were stained with p-EGFR, and β-tubulin overnight at 4°C, followed by secondary antibody staining at RT 1 h, 15 min DAPI staining. Asterisk indicates the statistical significance for pairwise comparison; p-value calculated using Wilcoxon rank-sum test with Bonferroni correction (ns: 0.05 < p, ∗: 0.01 < p ≤ 0.05, ∗∗: 0.001 < p ≤ 0.01 ∗∗∗: 0.0001 < p ≤ 0.001, ∗∗∗∗: p<=0.0001). Boxplot showing the distribution of the data with minimum, first quartile (Q1), median, third quartile (Q3), and maximum.

(F) The comparison of expression profiles of 24 markers between control and 40nM osimertinib samples was demonstrated. Boxplot showing the distribution of the data with minimum, first quartile (Q1), median, third quartile (Q3), and maximum.

(G) The p-AKT(S473) intensity comparison between PC9 control and PC9 treated with 40-nM osimertinib multiplexing samples was provided. Asterisk indicates the statistical significance for pairwise comparison; p-value calculated using Wilcoxon rank-sum test with Bonferroni correction (ns: 0.05 < p, ∗: 0.01 < p ≤ 0.05, ∗∗: 0.001 < p ≤ 0.01 ∗∗∗: 0.0001 < p ≤ 0.001, ∗∗∗∗: p<=0.0001).

To evaluate the effect of osimertinib on EGFR’s downstream proteins, we investigated whether p-EGFR would correlate with p-AKT expression level. The p-EGFR signaling proteins demonstrated a high correlation to p-AKT after drug treatment based on the mean-intensity correlation of single cells (Figures 4C and 4D), suggesting the role of osimertinib in inhibiting both EGFR and the downstream targets of AKT/mTOR pathways. Furthermore, the p-EGFR expression decreased as the drug concentration increased from 0 nM to 20–40nM (Figure 4E), whereas p-EGFR exhibited relatively higher expression under 60nM than 40nM. This may be due to concentration-dependent drug response, causing the discontinuity in receptor distributions and downregulation of the p-EGFR expression. Interestingly, the p-EGFR (Tyr1086) expression level remained minimally altered under the osimertinib treatment (Figure 4F). This may be due to the difference in phosphorylation sites, as a different site of p-EGFR (Y1173) was downregulated upon osimertinib treatment (Gu et al., 2020b). Therefore, multiple phosphorylation sites of p-EGFR can contribute to the effect of osimertinib in inhibiting p-EGFR. More importantly, p-AKT was consistently downregulated (Figure 4G), demonstrating the function of osimertinib in decreasing the downstream pathway of p-EGFR (Gu et al., 2020b).

Next, we asked whether osimertinib would alter the WNT pathway in PC9 cells. Prior studies on breast tumors and HC11 mammary epithelial cells yielded constituent activation of WNT1 to transactivate the EGFR pathway by increasing the ligand availability (Civenni et al., 2003; Faivre and Lange, 2007; Schlange et al., 2007). However, in the RapMIF results, WNT1 showed a low correlation with p-EGFR in both conditions, whereas a high correlation with p-AKT in the presence of the drug was obtained (Figures 4C and 4D). This observation indicated that the activation of WNT1 would not act on EGFR directly but could transactivate EGFR downstream targets. Also, the drug-treated cells exhibited a higher positive correlation between p-EGFR and non-p-β-catenin compared with the control (Figures 4C, 4D, and S24), and the osimertinib treatment increased the WNT1 signaling (Figure 4F). Thus, the inhibited EGFR pathway could promote the activity of WNT signaling and β-catenin under the osimertinib treatment.

We then analyzed the effect of osimertinib in spatial distributions and networks in PC9 cells. The spatial signaling network maps of the PC9 cells after drug treatment exhibited more dispersed distribution in the cytosol, noted as clusters 0, 1, and 2, compared with the control sample (Figures 5A–5C). Subsets of spatially resolved clusters exhibited high expression levels of WNT1, indicating that inhibition of the AKT pathway could increase the activity of WNT1 as a compensatory pathway. Subcellular spatial signaling networks were represented as a physical neighborhood map for the assessment of protein-protein proximity at the pixel level (Figure 5D). Each node corresponds to one of the most expressed signaling markers in that cluster of pixels (Figures 5E and S25). The color of the node matches the cluster color assigned to the pixel’s identity. The length of the connected line represents the pixel-to-pixel closeness between two clusters. The color of the line demonstrates the probability of neighboring with another cluster. The size of the node shows the percentage of the cluster out of all clusters. For instance, cluster 2 (WNT1) and cluster 6 (non-p-β-catenin) occupied a large proportion (frequency) after drug treatment, indicating the osimertinib inhibited the AKT pathway while upregulating the WNT pathway. Also, cluster 6 (non-p-β-catenin) was mainly distributed in the nucleus, indicating the location of non-p-β-catenin. Cyclin D1 was highly expressed in cluster 14, indicating the activation of cyclin D1/CDK4 complex in the cytoplasm to promote the G1 phase and cell proliferation (Gladden and Diehl, 2005). Another interesting finding was that the degree of the node of WNT1 in cluster 2 reduced from 10 to 4 under the drug treatment. This phenomenon may be due to the effect of osimertinib in reducing the neighborhoods and clustering between the WNT and AKT pathways. The p-β-catenin in cluster 12 and cluster 19 exhibited a higher degree of a node and more neighborhoods with p-EGFR clusters under the drug treatment. These observations indicated that the production and degradation of most p-β-catenin resulted from the EGFR pathway. The p-AKT displayed consistent downregulation in the replicates of multiplexed signaling experiments (Figures S26-S30). Therefore, RapMIF data supported the role of the WNT1 pathway as a compensatory pathway, contributing to drug resistance.

Figure 5

The effect of osimertinib in spatial subcellular distribution and network in PC9 cells

(A) Pseudo-colored cells by pixel clustering of 20 signaling protein markers in PC9 cells without osimertinib treatment. The legend indicates the cluster name in a frequency-descending order for each cell.

(B) Colored cells by pixel clustering of 20 signaling protein markers in PC9 cells with 48-h 40-nM osimertinib treatment. The legend provides the cluster name in a frequency-descending order for each cell.

(C) Heatmap of 20 clusters on a Z score scale was shown. Each cluster represented one distinct expression profile of 20 protein markers in PCI cells. Data are represented as mean expression per cluster.

(D) Analysis pipeline for spatial networks was outlined. Created with BioRender.com.

(E) Spatial signaling networks of 19 clusters in control and 40nM osimertinib-treated samples were provided. Each node was labeled with one of the most expressed signaling markers (excluding the epigenetic markers) in that cluster of pixels shown in Figure S25D. The color of the node matched with the cluster color shown in (C), the color of the line represented the probability of neighboring with another cluster, and the size of the node provided the pixel distribution across all clusters.

Subpopulation analysis on NSCLC cells

To visualize the heterogeneity of signaling expression levels, we performed subpopulation analysis on both EGFR-wild type (A549) and EGFR-mutant (PC9) NSCLC cell lines. We firstly clustered A549 into 12 subpopulations based on their single-cell mean intensity levels (Figure S31). Each subpopulation represented a distinct signaling profile across 17 signaling markers. A549 cells in subpopulation 3 displayed relatively lower expression levels of all signaling proteins. In contrast, cells in subpopulation 0 and subpopulation 6 showed higher signaling expression levels and highly expressed WNT1, p-EGFR, cyclin E, mTOR, p-AKT, and non-p-β-catenin, indicating these cells may have higher activation of WNT/β-catenin and AKT/mTOR signaling. To decipher the cell heterogeneity, we performed a similar analysis on PC9 cells in the presence and absence of drug treatment (Figure S32). The single cells were organized into 12 subpopulations based on their single-cell mean intensity. PC9 40-nM osimertinib-treated cells were mainly distributed in subpopulation 3 and subpopulation 0. However, the two populations displayed distinct signaling profiles. Subpopulation 3 demonstrated a signaling profile with high expression levels of 17 signaling markers, whereas the cells in subpopulation 0 have relatively low signaling expression levels. This population structure suggested that cells can exhibit various signaling profiles even after drug treatment. Therefore, dissecting the heterogeneity in the signaling distributions of cells in response to drug perturbations is necessary for combinational therapy design.

The signaling network concept and subpopulation analysis can be extended to other cell types. Thus, RapMIF examined a 7-plex protein panel in UC-mesenchymal stem cells (UC-MSC) (Figures S33 and S34). The more stretched cells provided better visualization of the signaling networks in the cytoplasm. Because WNT and AKT pathways also play an important role in cell differentiation and development, it would be worthwhile to study the signaling networks in MSCs for stem-cell-based therapeutic designs (He et al., 2021; Takam Kamga et al., 2021).

Cross-scale analysis of spatial signaling networks in lung tissues

To further verify the feasibility of signaling measurements in situ, we investigated RapMIF on FFPE tissues. The RapMIF on cell culture resolved the subcellular distribution among target proteins. However, pixel-level signaling analysis in tissues was challenging due to the densely packed cells with a limited spatial barcoding space in the cytosol of each cell. Thus, single-cell signaling maps were generated in tissues to decipher the architecture of normal and diseased cellular distributions. Hematoxylin and eosin (H&E) staining and immunohistochemistry (IHC) have been widely used to diagnose disease (Lin et al., 2018). However, the information from H&E is limited, and IHC is a single-channel imaging method. Tissue imaging typically requires intricate spatial analysis for feature extraction and cell classification algorithms due to the native tissue architecture. The RapMIF made it possible to study potential signaling imaging biomarkers to predict drug response using the baseline tissues obtained from the patients before targeted therapies.

RapMIF screened up to 25-plex signaling proteins in a tissue lung adenocarcinoma microarray (BS04081a, Biomax) that contained 63 cores, 21 patients’ cases, varying from normal, stage I, II, and III tumors. The images were first processed to correct the line artifacts caused by the uneven illumination from the flat field and dark field using the BASIC function in ImageJ (Figure S35) (Peng et al., 2017). Background corrections improved the downstream image analysis. Single cells were then segmented in the tissue cores (Figure S36) and classified into pan-cytokeratin positive and negative regions (Figure 6A and S37). Patients with higher tumor stages demonstrated upregulated WNT and AKT pathways with larger variation in the pan-cytokeratin-positive region (Figures 6B and S38). The constituent activation of the WNT and AKT pathway promotes the stabilization of non-p-β-catenin and the degradation of p-β-catenin (Jin et al., 2017). However, we observed upregulated p-β-catenin in tumor regions (Figures 6B and S39). The crosstalk mechanisms between the WNT pathway and other pathways, such as JAK-STAT and NF-κB signaling pathways, may contribute to the upregulation of p-β-catenin (Bai et al., 2017; Ma and Hottiger, 2016). Also, the activity of β-catenin is phosphorylation site dependent. Of note, β-catenin is phosphorylated at Thr41, Ser37, and Ser33 residues by GSK3 and then at Ser45 by CKI to undergo ubiquitination and degradation by the proteasome (Li et al., 2012).

Figure 6

Expression and spatial signaling analysis on tissue microarray from diverse lung cancers

(A) Multiplexed protein images of lung tissue samples at normal, malignant stages IB, IIA, and IIIA were presented. The tissue was masked based on pan-cytokeratin-positive staining.

(B) The normalized total intensity of each marker per cell was provided in pan-cytokeratin-negative (false) and -positive (true) regions. Asterisk indicates the statistical significance for pairwise comparison; p-value calculated using Wilcoxon rank-sum test (ns: 0.05 < p, ∗: 0.01 < p ≤ 0.05, ∗∗: 0.001 < p ≤ 0.01 ∗∗∗: 0.0001 < p ≤ 0.001, ∗∗∗∗: p<=0.0001). Boxplot showing the distribution of the data with minimum, first quartile (Q1), median, third quartile (Q3), and maximum.

(C) The comparison of signaling expression profiles of 24 markers was provided. The tissue cores were classified by stages and pan-cytokeratin expression. Classification of signaling maps was demonstrated for 21 patients and 55 tissue cores using the average linkage method based on Euclidean distance to cluster cores and markers along the x and y axis. Data are represented as mean expression per core.

(D) Heatmap of 19 clusters on a Z score scale was shown. Each cluster represented one distinct expression profile of 17 protein markers in lung microarray, excluding segmentation and epigenetic markers. Data are represented as mean expression per cluster.

(E) Multiplexed signaling protein images from tissue cores were analyzed by single-cell level clustering of 17 protein markers in four stages of tumor. The expression profile was clustered in pan-cytokeratin-positive regions. Each stage contained three images from the same patient.

Next, we sought to assess the signaling proteins across various cancer stages. Multiplexed signaling maps across 55 cores were classified by cancer stages and pan-cytokeratin expression at the single-cell level. The signaling markers in pan-cytokeratin-positive regions exhibited higher expression levels than normal tissue regions (Figures 6C and S40A). The majority of signaling markers, RNF43, EGFR, H3K27me3, WNT1, p-EGFR, AKT, mTOR, and cyclin D1, demonstrated upregulation in the tumor regions (Figure S40B). The variances of signaling expression levels were larger in pan-cytokeratin-positive regions, especially in stage IB tissues (Figure S40C). However, there was no obvious trend in the signaling levels across different stages of the tumor (Figure S40D). For patients 8 and 18, the expression levels of EMMPRIN, non-p-β-catenin, p-EGFR, APC, and cyclin E were higher in both pan-cytokeratin-positive and -negative regions.

To investigate the spatial distribution of signaling proteins, we performed a single-cell clustering analysis in the pan-cytokeratin-positive areas. The malignant tissue exhibited more variance in signaling expression from the same patient (Figures 6D, 6E, and S41). Five distinct clusters were highly distributed in the malignant IIIA tumors, whereas only one cluster was mainly expressed in the normal tissues (Figure 6E). In malignant stage IIIA, the sample was dominated by clusters 15 and 16, corresponding to the high expression levels of non-p-β-catenin, cyclin E, p-AKT, and p-mTOR. This finding indicated the upregulation of both WNT and Akt pathways. The signaling maps were similar in the pan-cytokeratin positive regions of normal and malignant stage IIA tissue (cluster 12 and 14, respectively); however, the tissues presented higher expression of p-EGFR and non-p-β-catenin in malignant stage IIA, revealing the upregulation of AKT pathways (Figures 6D and 6E).

The next step of RapMIF in tissues would be to study how AKT/mTOR and WNT/β-catenin pathways coordinate in wild-type, EGFR mutant FFPE tissues and fresh frozen tissues. Fresh frozen tissues may provide varying staining results, dependent on the conditions of fixation and antigen retrieval (Shi et al., 2008). Signaling networks in tissues of EGFR mutant patients in response to EGFR TKI treatments would help predict the drug response and isolate signaling mechanisms responsible for drug resistance. Considering the complexity of tumor samples, combinations of additional signaling targets would enhance the prediction accuracy of patient classification using deep learning methods.

Discussion

Signaling pathways mediate cell communication, proliferation, differentiation, and migration. The pathways are usually regulated through protein complexes, controlled via PPIs. To avoid the need of modifying samples genetically, which is not feasible with patients’ samples, proximity ligation assay has shown the capability of detecting the PPI in situ (distance <40nm) at endogenous protein levels (Alam, 2018). Förster resonance energy transfer (FRET) has also been used for observing dynamics and reversible PPIs in vivo by transferring the fluorophores from a donor to an acceptor in its closeness, proximity. Time-resolved FRET advances the current FRET methods with a faster, more sensitive, and less complex detection of PPIs (Maurel et al., 2008; Rajapakse et al., 2010). Due to spectral bleed and limited energy transfer efficiency, it is difficult to implement FRET with super-resolution imaging methods. Biomolecular fluorescence complementation (BiFC) has shown the feasibility of detecting PPIs based on fluorescent complementation reporters. Implementing BiFC on a super-resolution imaging approach, photoactivated localization microscopy (PALM), achieves imaging of the subcellular dynamics of PPIs at high spatial-temporal resolution (Liu et al., 2014).

Phosphorylation can be used to modulate PPIs, thereby regulating signaling transduction through either modifying proteins post-translationally or producing secondary messengers. Protein kinases play an important role in activating cellular processes, and the target proteins can be phosphorylated at different sites resulting in heterogeneous functions (Xue et al., 2012). Conventional methods for identifying kinase substrates suffer from several limitations. Kinase assay is one of the methods to identify kinase-substrate pairs by using analog-sensitive alleles. The entire in vitro substrate can be assayed for a particular kinase. Various synthesized analogs provided a panel for substrate identification and protein kinase function (Elphick et al., 2007). However, this method is not applicable when the introduced analog does not affect the kinase’s function (Mok et al., 2009). Protein microarray analysis can also reveal the potential phosphorylation reactions as well as important functions. To identify potential substrates, the protein kinases are added to arrays of immobilized proteins and generate signals upon detection (Meng et al., 2008). However, the relationships are identified in vitro and may not represent the relationships in vivo (Mok et al., 2009). Kinase substrate tracking and elucidation (KESTREL) has also been used to identify protein kinase substrates. The high-background phosphorylation in the cell extracts can obscure the detection of putative phosphorylation; however, KESTREL achieves the identification of physiologically relevant substrates of proteins in crude cell extracts with an improved signal-to-noise ratio by limiting the phosphorylation reaction time and increasing concentrations of protein kinase of interest. Performing chromatography steps before the analysis could also greatly improve the SNR ratio. Also, to establish high physiological relevance, two or three closely related protein kinases with similar binding specificities were used for one substrate detection. However, the temporal and spatial information is largely ignored in the cell lysate (Knebel et al., 2001). However, the assays have experienced challenges with the background in cell extracts and medium-to-low specific activity to kinases (Knebel et al., 2001; Troiani et al., 2005).

High-throughput techniques such as mass cytometry have been widely used to identify the phosphorylation sites of the target protein. Although the bulk level correlation between the protein expression levels of markers was previously mapped out, the spatial correlations and the precise association between the kinases and their direct substrates are still unclear (Xue et al., 2012). Also, the cell-to-cell heterogeneity in the signaling transduction resulting from factors crucial for cancer development could not be systematically profiled due to technical limitations (Lun and Bodenmiller, 2020). Our study could advance the current single-cell mass cytometry (CyTOF) in detecting kinases involved in signaling pathways using additional antibodies in the multiplexed imaging data (Han et al., 2015; Leelatian et al., 2015; Mingueneau et al., 2014). RapMIF provides spatial information on the protein distributions compared with traditional protein analysis assays. Because kinase assay requires the synthesis and introduction of ATP-analogs and is only applicable when the introduced analogs do not interfere with kinase’ function (Mok et al., 2009), RapMIF circumvents the complex process of synthesizing ATP analogs by directly detecting protein activity from multiple cycles of IF. Also, unlike KESTREL, RapMIF conserves the subcellular structure by avoiding lysate in the cells.

The multiplexed single-pixel analyses have quantified the subcellular protein correlations, indicating potential direct physical neighborhood among signaling protein subsets. Distinct multiplexed intensity profiles of signaling markers revealed the crosstalk between heterogeneous signaling pathways. The presented RapMIF method simplified multiplexed protein analysis by pre-conjugating the antibodies with commercially available dyes. Also, we observed that the bleaching solution interfered with Phalloidin signals in A549 cells. Therefore, in the A549 data analysis, a combination of p-EGFR, concanavalin A, B-actin, WGA, APC, and WNT1 markers was used to outline the cell boundary for segmentation. Also, to avoid the bleaching effect, cycle 2 (Phalloidin and WGA) staining was done directly after cycle 1 (p-EGFR) without an intermediate bleaching step in PC9 multiplexing samples. A549 cells showed a larger cytosolic area than that of PC9 cells, providing a more detailed spatial network model for studying the protein distribution in the cytosol. The spatial maps demonstrated the translocation patterns of signaling markers corresponding to their functions in both cell types. The combination of ERK or MEK inhibition with osimertinib could delay or overcome the acquired osimertinib resistance (Gu et al., 2020a). Thus, the spatial signaling patterns in the MEK/ERK pathway, another downstream pathway of EGFR, regulating cell proliferation and survival could further be determined using additional antibody targets.

Even in the same population of cells, the cells may display heterogeneous signaling profiles due to the difference in signaling-acting pathways. Autocrine signaling involves the production and secretion of an extracellular mediator, acting on the same cell, whereas paracrine signaling occurs between different types of cells. Both signaling pathways mediate tumor growth, invasion, and metastasis. During cancer development, the growth factors are usually produced in an autocrine manner, allowing the cell to simulate itself in positive feedback (Caicedo, 2013; Ungefroren, 2021). Current methods have been developed to decipher the cell-cell communication in single cells; however, they failed to demonstrate the heterogeneity in autocrine and paracrine signaling interactions. CellChat as a predictive model for cell-cell communication has been achieved to study the autocrine-acting versus paracrine-acting pathways in the same cell type or same subpopulation of cells. However, this approach highly relies on the accuracy of the ligand-receptor databases, and it produces fewer interactions than other methods when detecting communications from distant cells (Jin et al., 2021). Therefore, combining predictive models and proteomics could optimize cell-cell communication methods. Also, the complex crosstalk among signaling pathways could contribute to the heterogeneity in signaling profiles. Fluorescent reporters can be introduced and utilized to track protein activity in the subpopulation of cells following drug treatment (Kurppa et al., 2020). However, the current study focuses on the design of signaling pathway analysis using molecular snapshots, instead of the subpopulation validation in real time.

RapMIF multiplexed 17 signaling markers and converted their distributions into 33-pixel spatial clusters. Physical proximity analysis between two distinct proteins yielded a total of 272 (17 choose 2 combinations) possible pairwise neighborhoods, yielding their correlations in cytosolic and nuclear expression levels at the single-pixel level. Pairwise analysis of protein neighborhoods could reveal potential critical signaling proteins or networks (Ramana, 2019). Besides, coordination of the multi-protein neighborhood can be extracted from the multiplexed signaling data. For instance, the tri-protein interactions in cytosol and nucleus are assessed with 1,360 (17 choose 3 combinations) possible three-way neighborhoods. Combinatorial signaling analysis would predict the dynamics of signaling networks when multiple proteins physically interact, a task that would be challenging to achieve using bulk signaling assays.

Limitations of study

Our system design achieved a fully automated staining system including blocking, staining, imaging, and bleaching. The fluidic automation takes the antibody diluted in CSM media in the multi-well in a 2.5-fold dilution. Therefore, to save the antibody volume, the antibody loading step by step could be done by manually pipetting the antibodies to the sample. If needed, the automated system could be further optimized to reduce the dilution rate of the antibody during the staining step. Aspiration settings and washing speed could be adjusted to minimize cell loss across cycles. The p-EGFR signaling receptors in PC9 cells were dependent on the sample quality due to the step of drying the coverslip before mounting it to the holder. The freshness of the sample may affect the morphology of p-EGFR. Thus, the multiplexing experiments in PC9 cells were initiated on the same day of cell fixation to make sure the quality of the samples. Thus, the RapMIF could provide a precision oncology framework for defining the signaling dynamics in both intensities and spatial distribution of human cancers.

For multiplexing experiments, two factors are concerned to ensure cell quality, the efficiency of the bleaching approach, and the total spent time for one set experiment. RapMIF utilizes fluorophore-conjugated antibodies for direct labeling. Fluorophore-conjugated antibodies generally produce relatively weaker signals than indirect labeling. However, indirect labeling raised some problems in the context of multiplexed experiments. It is constrained by the availability of the reactivity of secondary antibodies to limited species (e.g., mouse, human, donkey). Also, multiplexing indirect IF requires antibody stripping rather than a regular bleaching approach. The harsh process may cause cell loss and tissue detachment. Compared with the antibody stripping approach removing the whole secondary antibody, the RapMIF bleaching approach is effective in removing dyes while preserving the cells. One concern of the RapMIF is whether the direct IF staining at RT for 1 h for each cycle can label all the antigens in the cell. This trade-off can be addressed by improving the antibody concentration to generate similar signals in 1-h RT staining compared with overnight incubation. Our approach aims to achieve rapid profiling of signaling pathways while maximizing the efficiency of protein detection. The staining settings were kept consistent across different samples, making it reliable for studying drug efficiency. Also, staining the sample for 1 h instead of overnight reduces the 12-cycle multiplexing experiments from 24 days to 6 days and prevents cell degradation. Another drawback of the RapMIF on FFPE tissues is that FFPE tissue generates a high autofluorescence background and is affected by shading or vignetting due to imaging artifacts. Signaling proteins are expected to express nearly everywhere, making it harder to distinguish between background and real signals. Therefore, we performed background correction steps using the BASIC algorithm to reduce the contribution of background on real signals (Peng et al., 2017). Also, post-bleach imaging data serve as a control of native autofluorescence signals and confirm the background corrected images using the stained raw tissue image that is subtracted from the post-bleach image.

In terms of antibody specificity validation, a positive control (e.g., a cell line expressing the target) and negative control (e.g., a cell line not expressing the target) can be utilized to confirm the selective binding of the antibody. Sometimes, it is difficult to find cell lines that do not express the target. Therefore, knockout (KO) models such as CRISPR–Cas9 and siRNA can serve as negative controls and be applied to various approaches, including mass spectrometry, western blot, IHC, and IF (Kurppa et al., 2020; Uhlen et al., 2016). In previous studies, independent antibody strategies are also suggested to confirm antibody specificity. When two antibodies bind to different epitopes of the same target, they can yield correlated signals, suggesting both antibodies recognize the same target (Uhlen et al., 2016). However, this is beyond the scope of current work, and independent validations of antibody specificity are lacking in the current study. In our study, the antibody staining has been evaluated using conventional IF protocol before performing the multiplexing to improve the staining quality of the antibodies.

Another important consideration is the physical dimensions of signaling protein visualizations. The diffraction limit of light microscopy is around 200 nm, whereas the antibody size is around 10 nm (Erickson, 2009; Tan et al., 2008; Xing et al., 2016). Due to the low spatial resolution, standard light microscopy fails to resolve the protein complex interactions based on the fluorescent signals. To determine the colocalization of two proteins, their spatial distributions are needed. When two proteins are colocalized with one another, the overlay of the two spatial distributions should show a correlation of intensity at the pixel level. In our case, the pixel size is around 188 nm using an oil immersion 40X objective lens with a 1.4 numerical aperture and the protein molecule is around 2–9 nm (Erickson, 2009). For a single cell, the number of protein molecules is 2 × 109, considering the number of ribosomes, tRNA, synthesis rate, and translation rate, among others (Princiotta et al., 2003). Each pixel would contain around up to 15,700 proteins. Thus, due to the limitation on the resolution of the microscope, it is unable to resolve the single protein interaction. Additional color deconvolution methods might be needed to reach single protein detection. The pixel-level correlations can distinguish the random color overlap due to protein compartmentalization from protein colocalization and indicate that true colocalization only occupied 3% of total colocalization (Costes et al., 2004). To remove nonspecific interactions in data analysis, a probabilistic protein interaction network has previously been developed to identify protein interactions based on the normalized spectral abundance factors using Bayes' approach. The probabilistic model predicts each pair of proteins by providing quantitative information on the preference of each interaction. It takes advantage of quantitative proteomics by measuring spectral counts to generate a probabilistic network of protein interactions, providing an alternative framework to quantify probabilistic pixel-level protein interactions in the RapMIF data (Sardiu et al., 2008). Although it is not implemented, it can provide an alternative data analysis approach to quantify probabilistic pixel-level protein distributions in the same pixel area in the RapMIF data. Current studies focus on the development of the RapMIF pipelines for studying the crosstalk of multiple signaling pathways and the protein neighborhood likelihood using spatial networks by studying the neighboring protein distribution in the same pixel area. However, super-resolution can be implemented in this technology to further resolve protein complexes. Protein interaction detection assays, such as proximity ligation assays, can be multiplexed to achieve multiparameter detection of PPIs.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ahmet F. Coskun (ahmet.coskun@bme.gatech.edu).

Materials availability

This study did not generate new unique reagents.

Data and code availability

Supplemental videos data have been deposited at Zenodo and are publicly available as of the date of publication. DOIs are listed in the Key resources table. The imaging data reported in this study cannot be deposited in a public repository because the data size is larger than 10GB that is beyond the limit of Mendeley Data uploads. All data reported in this paper will be shared by the lead contact upon request.

All original code has been deposited at Github and is publicly available as of the date of publication. DOIs are listed in the Key resources table.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Experimental model and subject details

Our study does not use experimental models (e.g., organisms) typical in the life sciences.

Method details

Cells

A549 and PC9 cell lines were seeded on coverslips treated with 0.01% poly-L-lysine in a 6-well plate overnight, and cultured for 48 h respectively at a 37°C incubator. A549 was purchased from ATCC. The NSCLC-sensitive cell line, PC9, was provided by Dr. Sun Shi-yong (Emory University).

Drugs

The Osimertinib was provided by Dr. Sun Shi-yong (Emory University) at 10 mM concentration. The cells were treated with 0, 20, 40, and 60nM osimertinib for 48 h in a 37°C incubator. Following that, the cells were then fixed and permeabilized with 1.6% formaldehyde in 1X PBS for 10 min at RT and 0.5% Triton X-100 for 10 min at RT, respectively (Figure S4A). The cells were then blocked using a cell staining medium (CSM) containing 5days. The effect of osimertinib on p-EGFR (Tyr1086) and cell viability was tested. PC9 cells were stained with p-EGFR and β-tubulin overnight followed by 1-h incubation of secondary antibodies at RT, and 15 min DAPI staining. The cells were then imaged at 40X magnification, and the cell number was analyzed using ImageJ. For the multiplexing experiment, PC9 cells were treated with 40nM osimertinib for 48 h in a 37°C incubator, and then fixed, permeabilized, and blocked. The coverslip was mounted to an acrylic holder for multiple cycles of staining.

Tissues

FFPE lung cancer microarray tissue blocks were purchased from Biomax with 61 cores and 21 cases (BS04081a). The slide was baked in a 60°C oven for 2 h followed by deparaffinization and heat-induced epitope retrieval (HIER). Sectioned slices were ready for multiple rounds of staining, imaging, and bleaching.

Antibodies

In our study, a total of 26 antibodies were used. Cell segmentation included Phalloidin (A34055, Thermofisher Scientific), wheat germ agglutinin (WGA, W32466, Thermofisher Scientific), Concanavalin A (C11252, Thermofisher Scientific), β-tubulin (sc-5274, Santa Cruz Biotechnology), and β-actin (4970BF, Cell Signaling Technology). For tumor tissue region identification: pan-cytokeratin (914,204, BioLegend) was used. 11 Markers targeting the WNT/β-catenin pathway included non-β-catenin (Ser45) (70034S, Cell Signaling Technology), APC (ab239828, Abcam), DKK1 (sc-374574, Santa Cruz Biotechnology), Cyclin E (sc-247, Santa Cruz Biotechnology), EMMPRIN (sc-21746, Santa Cruz Biotechnology), WNT1 (sc-514531, Santa Cruz Biotechnology), p-β-catenin (Ser45) (44–208G, Thermofisher Scientific), RNF43 (ab84125, Abcam), DKK2 (PA5-11616, Thermofisher Scientific), AXIN1 (2087BF, Cell Signaling Technology), Cyclin D1 (ab156448, Abcam). 6 Markers involved AKT/mTOR pathway included p-AKT (Ser473) (4060BF, Cell Signaling Technology), AKT (4691BF, Cell Signaling Technology), p-mTOR (Ser2448) (5536BF, Cell Signaling Technology), mTOR (2983BF, Cell Signaling Technology), p-EGFR (Tyr1086) (36–9700, Thermofisher Scientific), EGFR (sc-373746, Santa Cruz Biotechnology). 3 epigenetic markers included H3K27me3 (9733BF, Cell Signaling Technology), H3K9ac (ab203951, Abcam), H3k4me3 (ab227060, Abcam). For rapid staining of the sample, the BSA-free primary antibodies were pre-conjugated with Alexa Fluor 488, 555, or 647 using the Lightning-Link Rapid Conjugation kit (ab236553, ab269820, ab269823). Before conjugation, the concentrations of the antibodies were measured using Nanodrop to ensure the required antibody amount is in the range. The overall time for the conjugation process is around 20 min. 15 antibodies involved in signaling pathways were first tested in IF of A549 cells. For conjugated antibodies, IF was tested under four conditions, indirect/direct labeling, overnight 4°C/RT 1 h, respectively. For antibodies conjugated to Alexa Fluor dyes, the direct labeling was tested at overnight 4°C/RT for 1 h (Figure S2). To improve the staining quality, titrations were performed on antibodies, and IF was evaluated at two different dilution rates (Figure S1). The improved staining condition of each antibody was determined based on the SNR enhancement. The SNR of each condition was measured based on calculations of the max and min pixel values of five random cell regions in ImageJ (Table S1). For IF staining on cell culture, the antibodies and DAPI (62,248, Thermofisher Scientific) were diluted in CSM containing 0.5% BSA, and 0.02% sodium azide in PBS. For IF staining on tissue, the antibodies and Hoechst 33,342 (H3570, Thermofisher Scientific) were diluted in a protein block buffer (DAKO, X0909).

RapMIF

The RapMIF setup consisted of an autosampler that transferred antibodies and performed washes of the sample automatically and controlled by a custom-written python code (Figure S4). The buffers used during the procedure are 1X PBS (D8537, Sigma), CSM, and Fluorophore bleaching solution. PBS 1X was the wash buffer for the stabilization of the cells and the tubes. CSM was the imaging buffer that helped in the annealing of antibodies with the cells. Fluorophore bleaching solution was the strip buffer used to remove staining after each cycle. The walls on the 96-well plate for the antibodies were pre-washed with 1x PBS. The well plate was covered with foil stickers to ensure a dark condition. The microscope stage setup included tubes for transporting buffers and an acrylic holder for mounting cells on coverslips. The stage was designed to be compatible with the Keyence BZ-X810 microscope. The flow of the buffers was controlled by a valve. The stage included three 3D-printed ports, connecting fluidic tubes to needle tips. One port is used for delivering the fluids to the solution, while the other two function as aspirating fluids and preventing overflow using a vacuum (Figure S4B). To mount the cells, the coverslip containing cells was placed on an adhesive spacer. It was attached to the acrylic holder using instant adhesive such that the cells were on the side facing the well. The acrylic holder was then secured onto the microscope stage. Once the cells were mounted and the antibodies and buffers were ready for transfer, automatic staining was done by running Python code. The code used the following as inputs: 1) position of the first well, 2) the number of imaging cycles, 3) time to wait after staining, 4) whether to bleach before staining, and 5) time to wait after bleaching. The default time for staining was 60 min. The default time for the bleaching step was 60 min. In between subsequent cycles, the dyes attached to the antibodies were bleached out using the fluorophore bleaching solution, followed by the re-labeling of another set of three antibody-dye conjugates (Figure S4C). The parameter settings including fluid transfer and antibody incubation for the RapMIF setup are shown in Table S2. The total time taken for one staining cycle is approximately 5 h including imaging and bleaching steps. RapMIF advances the current CycIF approach by saving time spent on each cycle. RapMIF takes around 5 h per cycle, while CycIF takes more than 12 h per cycle (Table S3.). Although another 1 h for each cycle was spent for cell reblocking to avoid the nonspecific binding, the overall time of RapMIF is still shorter.

Cell multiplexing

Each coverslip prepared with cells was stained with 25 markers using a total of 11–12 cycles (Figure S4C). The staining process followed the steps below: 1) block the sample with CSM for 1 h at RT; 2) stain the sample with a set of conjugated-antibodies and DAPI diluted in CSM for 1 h at RT followed by 3 times 1X PBS wash; 3) image the sample using Keyence fluorescence microscopy; 4) deactivate the fluorophores using 3% H2O2 and 20mM NaOH made up in PBS for 1 h at RT in the presence of white light, followed by 3 times 1X PBS wash; 5) reimage the same region of interest (ROI) of the sample to check the residuals after bleaching, following the step 1–5 for next round of staining. To quantify the frequencies belonging to the signal that disappears after bleaching, a Fourier analysis of before and after bleaching was performed to indicate the difference in frequencies belonging to the signal (Figure S6).

To examine the efficiency of bleaching, the signaling levels of p-AKT-488, Cyclin E−555, and EGFR-647 before and after bleaching were quantified (Figure S7). The A549 cells were stained with conjugated antibodies for 1 h at room temperature. The sample was imaged in ten different regions. The signaling levels at three conditions are visualized in the histogram. For all three markers, the signaling levels of them returned back to the background level after 1-h incubation of the bleaching solution. Bleaching for 60 min should be sufficient to reduce the intensity to 102 to 103-fold (Lin et al., 2018). The signals before staining, after staining, and after bleaching were quantified after background subtraction and applying a cell mask (Figure S7B). Also, quantification of single-cell marker intensity with comparison to after bleach intensity for PC9 cells was performed to examine the efficiency of the bleaching step (Figure S7C).

Tissue multiplexing

RapMIF is compatible with tissues mounted on coverslips and the tissue holder would be processed on a cooled stage holder (Figure S4D). For antibody incubation on tissue, 4°C can be performed without removing the sample from the microscope to ensure the same position of the sample for imaging. The temperature is controlled using a cooling and heating stage, a temperature controller, and a BTC-W heat exchange unit. The heat exchange unit uses liquid water to bring down the temperature of the stage during cooling. A tubing system is used for a closed-loop of water circulation. However, the commercial microarrays are available in the tissue slide format, necessitating a modified pathological sample preparation pipeline to make microarrays on coverslips. While microarray on coverslip is feasible, here we worked with available tissue microarrays. Thus, for tissue specimens mounted on slides, to maximize the efficiency of detecting epitopes in tissue, instead of incubating at room temperature for 1-h, overnight incubation was followed similar to the t-CycIF protocols (Lin et al., 2018). More specifically, the iterative process was modified in blocking and bleaching steps: after antigen retrieval, the sample follows the steps below: 1) block the sample with protein block buffer for 30 min at RT; 2) stain the sample with a set of conjugated-antibodies and Hoechst diluted in protein block buffer overnight at 4°C, followed by 3 times 1X PBS wash; 3) mount the sample in 10% glycerol made in 1X PBS, and cover the slide using 24 × 50mm No. 1 coverslip (3322, Thermo Scientific) to prevent evaporation during imaging; 4) image the sample using Keyence fluorescence microscopy; 5) de-coverslip the sample by placing the slide in a vertical jar containing 1X PBS for around 10 min, and the coverslip was released due to gravity. 6) deactivate the fluorophores using 4.5% H2O2 and 24 mM NaOH made up in PBS for 1 h at RT in the presence of white light, followed by 3 times 1X PBS wash; 7) followed by mounting, reimage the same ROI of the sample to check the residuals after bleaching, following the step 2–7 for next round of staining. After IF staining, we performed H&E staining at the last cycle.

Imaging

A wide-field microscope, Keyence BZ-X700, was used for fluorescence and brightfield imaging. For each cycle of imaging, four channels were captured with 40% excitation light: Channel 1 detects DAPI/Hoechst at 360 nm, Channel 2–4 detects fluorophores at Alexa Fluor 488 nm (FITC), 555 nm (TRITC), and 647 nm (Cy5), respectively. The exposure time was varied, but the exposure time for each marker across control and drug-treated PC9 cells was consistent. The exposure time for Channel 1, DAPI/Hoechst was around 1/400s, and Channel 2–4 varied from 1s to ½0s. The sample for cells grown on the coverslip was imaged at a 40X oil lens for 48–62 region of interest (ROI), and each ROI was imaged across 25–30 z-stacks with 0.4 μm/stack. To determine the optimal step size for imaging, A549 cells stained with WNT1 marker were imaged at four different step sizes, 0.1-μm, 0.4-μm, 0.8-μm, and 1.2-μm (Figure S8). The step size of 0.4-μm was chosen for the imaging system to reveal more continuity of the dim signaling molecules. Even though the acquisition time was around 1 h for each cycle, and slower compared to images at a larger step size, we could better resolve the signaling molecules. The 40X oil lens achieved a high-resolution acquisition at 0.18872 μm/pixel. The tissue microarray was imaged at a 20X dry lens, each ROI was imaged using autofocus. The 20X dry lens provided a resolution of 0.37742 μm/pixel. The best focus z level was determined by calculating the blur score of each DAPI channel image. The blur score of each image was calculated by taking the fast Fourier transform (FFT) transform of the image and filtering out low frequencies, the resulting inverse Fourier transform image’s mean magnitude was calculated. All z-level images of one location were then ranked based on the mean magnitude as a blurriness score. The best focus images around +/− 2 z-stacks were used for signaling analysis. 3D subcellular spatial signaling levels were also assigned into 26-pixel clusters in 3D across 21 z-stacks, demonstrating the 3D distributions of signaling proteins as a proof of concept.

Background correction

To correct the uneven illumination caused by the flat and dark field, a BASIC algorithm in the ImageJ plugin was used to remove the line artifacts in each tissue core (Peng et al., 2017). Since the patterns of the uneven illumination across different channels were different, and the shape and size of tissue cores varied, the correction was applied per tissue core for each channel separately across 12 IF cycles. The software estimated the flat and dark fields and subtracted the raw images from the two fields to generate corrected images. Then, images were post-processed for quantification of signaling networks and intensity levels.

Segmentation

For cell culture, single-cell nuclei regions were segmented using Cellpose (Stringer et al., 2021) by cell painting markers (Bray et al., 2016). The best focus images around +/− 2 z-stacks were used for segmentation of the cells. DAPI or Hoechst marker was used to segment cell nuclei boundaries. In the A549 experiment, the Phalloidin channel showed strong background noise due to the effect of the bleaching solution, therefore we combined p-EGFR, Concanavalin A, B-actin, WGA, APC, and WNT1 to outline the cell boundary for segmentation. For PC9 data, Phalloidin marker was used for cell boundary segmentation. For FFPE lung cancer microarray tissue, DAPI was used to segment the nuclei region. The cytosol region was calculated by expanding the nuclei segmented region by 10 pixels. The 10 pixels were determined empirically by looking at nuclei segmentation and cytosol segmentation markers (Figure S36).

Image intensity

From nuclei and cytosol cell masks, the mean intensity of markers in each cell area was calculated. Each cell nuclei and cytosol region were separated and mean intensity for each marker was calculated and cell masks with a colormap showing the intensity were generated to visualize the variation across signaling markers within and outside the nuclei area.

Correlation analysis

Correlation between pairs of cells was performed based on the mean intensity level for a pair of markers across all the cells (Figure 2A). We investigate the correlation between pairwise markers at the cell level to identify co-expression patterns. For the correlation analysis, the Pearson correlation of pairwise markers at the single-cell level was quantified by using the mean intensity distribution for pairs of markers across the cells (Figures 2E, 2F, 4C, and 4D). We first extracted the per cell mean intensity for all markers to obtain a vector of mean intensities per cell. We then separated the nuclei and cytosol region for each cell and extracted one vector of mean intensities in the cytosol and another vector of mean intensities in the nuclei per cell. To confirm the bleaching efficiency, we compared the mean intensity correlations of the before and after bleach images (Figure S14). The 3-way correlation was calculated with the adjusted-Rˆ2 value of the linear regression model of one marker predicted by the two others. Specifically, we predicted the expression of the maker on the z axis by the two markers on the x and y axis using a 3D scatterplot (Figure S13).

Marker intensity prediction

Per cell intensity of individual markers was predicted from other available markers using a Random Forest regressor with a 5-fold cross-validation training. The training set was composed of 2885 single-cell data points with a testing set of size 572. The regression was evaluated by the r-score between true intensity per cell compared to predicted intensity per cell. The R value found was 0.76 with a p value<0.01. Feature importance was plotted using the impurity index.

Pixel clustering

Pixel phenotypes were clustered in a two-step clustering pipeline (Figure 3A). From each pixel location within the cell segmented region, the intensity value of each marker expression was extracted. The resulting feature matrix consisted of n rows of a total number of pixels and p columns equal to the number of markers. Each column of the feature matrix was min-max normalized. Then to filter out the background, each pixel location (row in the matrix) with all markers intensity lower than 0.3 was considered as background and dropped from the feature matrix. The cut-off of 0.3 was determined empirically so noise was filtered out of the cell segmented region (Figure S15). The feature matrix was first clustered using the K-Means algorithm with a high number of clusters (k = 60). Then, the mean expression per cluster was calculated and the cluster was merged based on hierarchical clustering on the cosine similarity of the mean expression using the average linkage method with a tolerance of 0.2 of the maximum distance between clusters. The latent feature space was obtained using Uniform manifold approximation and projection (UMAP) and cluster labels were visualized by assigning each cluster color and reprojected on the UMAP embedding. A Heatmap of the mean expression level of marker per cluster was generated to show the variation of cluster phenotype.

Spatial signaling network

The network representation of spatial signaling was represented using the spatial information of pixel clustering. Each cluster was represented by a node with node size in the network proportional to the number of pixels in the corresponding cluster. For each cluster, the distribution of the neighbor cluster was calculated by looking at the 2-hop neighbor’s cluster of all pixels in the cluster. The edge color represented low (blue) to high (red) neighboring probability between two clusters.

Subpopulation clustering

A549 and PC9 dataset subpopulations were defined based on signaling phenotype using unsupervised clustering with the Leiden algorithm from multiplex protein mean expression level per cell (Traag et al., 2019). Both A549 and PC9 datasets were clustered into 12 clusters.

Receptor discontinuity

The centroid of nuclei of each cell was calculated from each cell mask. Then each cell was divided into regions defined by the absolute distance from the centroid. In each region, the mean intensity of the marker was quantified as a function of distance from the nuclei centroid. Distances from nuclei centroid in each cell were shown as a radius map and overlaid with marker intensity with a custom colormap. The colormap for marker intensity was transformed from matplotlib package ‘greens’ with alpha set using the NumPy function linspace going from 0 to 1 with N equal to the number of RGB quantization levels in the colormap.

Cell survival

The effect of osimertinib on PC9 cell proliferation and survival was measured using ImageJ. The cell survival was measured based on the cell number on the coverslip. The cells were seeded on coverslips at the same starting concentration, 0.833 × 105, followed by 48-h treatment of osimertinib at 0, 20, 40, and 60 nM. The cells were stained with p-EGFR, and β-tubulin and three ROIs were imaged for each concentration at 40X (Figures S19 and S20). The cell number on each ROI was counted in ImageJ following the steps: import the Channel image, make binary, convert it to masks, analyze particles, and then choose “masks” under the “show” option.

Quantification and statistical analysis

Pearson correlations with Holm-Šídák correction method were used to correct pairwise marker correlation at the single-cell levels. Wilcoxon rank sum tests with Bonferroni correction method were used for comparison of intensity distribution between groups. The level of significance was set at: not significant (ns): 0.05 < p, ∗: 0.01 < p ≤ 0.05, ∗∗: 0.001 < p ≤ 0.01 ∗∗∗: 0.0001 < p ≤ 0.001, ∗∗∗∗: p<=0.0001. All data visualization and statistical analyses were carried out using in-house Python scripts.

Additional resources

Our study does not have additional resources.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
p-EGFR (Tyr1086)	Thermofisher Scientific	Cat# 36–9700; RRID: AB_2533287
Phalloidin	Thermofisher Scientific	Cat# A34055
WGA	Thermofisher Scientific	Cat# W32466
non-p-β-catenin (Ser45)	Cell Signaling Technology	Cat# 70034S; RRID:AB_2799770
APC	Abcam	Cat# ab239828
DKK1	Santa Cruz Biotechnology	Cat# sc-374574; RRID:AB_10989416
Cyclin E	Santa Cruz Biotechnology	Cat# sc-247; RRID:AB_627357
EMMPRIN	Santa Cruz Biotechnology	Cat# sc-21746; RRID:AB_626911
Wnt1	Santa Cruz Biotechnology	Cat# sc-514531
p-AKT(Ser473)	Cell Signaling Technology	Cat# 4060BF; RRID:AB_2315049
p-β-catenin (Ser45)	Thermofisher Scientific	Cat# 44–208G; RRID:AB_2533604
RNF 43	Abcam	Cat# ab84125; RRID:AB_2181252
AKT	Cell Signaling Technology	Cat# 4691BF; RRID:AB_915783
p-mTOR (Ser2448)	Cell Signaling Technology	Cat# 5536BF; RRID:AB_10691552
mTOR	Cell Signaling Technology	Cat# 2983BF; RRID:AB_2105622
DKK2	Thermofisher Scientific	Cat# PA5-11616; RRID: AB_2091482
AXIN1	Cell Signaling Technology	Cat# 2087BF: RRID:AB_2274550
EGFR	Santa Cruz Biotechnology	Cat# sc-373746; RRID:AB_10920395
Cyclin D1	Abcam	Cat# ab156448
H3k27me3	Cell Signaling Technology	Cat# 9733BF; RRID:AB_2616029
H3k9Ac	Abcam	Cat# ab203951
H3k4me3	Abcam	Cat# ab227060
Concanavalin A	Thermofisher Scientific	Cat# C11252
β-actin	Cell Signaling Technology	Cat# 4970BF; RRID:AB_2223172
β-tubulin	Santa Cruz Biotechnology	Cat# sc-5274; RRID:AB_2288090
Pan-cytokeratin	BioLegend	Cat# 914204; RRID:AB_2616960
Biological samples
FFPE tissue	Biomax	BS04081a
Experimental models: Cell lines
A549	ATCC	https://www.atcc.org/products/ccl-185
PC9 cell lines	Dr. Sun Shi-yong (Emory University)	N/A
Software and algorithms
ImageJ	(Schindelin et al., 2012)	https://imagej.nih.gov/ij/
Cellpose	(Stringer et al., 2021)	https://github.com/MouseLand/cellpose
Dask	N/A	https://ml.dask.org/clustering.html
Statsannotations	N/A	https://github.com/trevismd/statannotations
Pandas	N/A	https://github.com/pandas-dev/pandas
Seaborn	N/A	https://github.com/mwaskom/seaborn
Networkx	N/A	https://github.com/networkx/networkx
Code related to this study	This study	https://github.com/coskunlab/Spatial-signaling
Other
Supplemental videos	This study	https://zenodo.org/record/5636798#.Yo2F5FTMJhE