Abstract

   The spatial proteomic profiling of complex tissues is essential for
   investigating cellular function in physiological and pathological
   states. However, the imbalance among resolution, protein coverage, and
   expense precludes their systematic application to analyze whole tissue
   sections in an unbiased manner and with high resolution. Here, we
   introduce panoramic spatial enhanced resolution proteomics (PSERP), a
   method that combines tissue expansion, automated sample segmentation,
   and tryptic digestion with high-throughput proteomic profiling. The
   PSERP approach facilitates rapid quantitative profiling of proteomic
   spatial variability in whole tissue sections at sub-millimeter
   resolution. We demonstrated the utility of this method for determining
   the streamlined large-scale spatial proteomic features of gliomas.
   Specifically, we profiled spatial proteomic features for nine glioma
   samples across three different mutation types (IDH1-WT/EGFR-mutant,
   IDH1-mutant, and IDH1/EGFR-double-WT gliomas) at sub-millimeter
   resolution (corresponding to a total of 2,230 voxels). The results
   revealed over 10,000 proteins identified in a single slide, which helps
   us to portray the diverse proteins and pathways with spatial abundance
   patterns in the context of tumor heterogeneity and cellular features.
   Our spatial proteomic data revealed distinctive proteomic features of
   malignant and non-malignant tumor regions and depicted the distribution
   of proteins from tumor centers to tumor borders and non-malignant tumor
   regions. Through integrative analysis with single-cell transcriptomic
   data, we elucidated the cellular composition and cell–cell
   communications in a spatial context. Our PSERP also includes a
   spatially resolved tumor-specific peptidome identification workflow
   that not only enables us to elucidate the spatial expression patterns
   of tumor-specific peptides in glioma samples with different genomic
   types but also provides us with opportunities to select combinations of
   tumor-specific mutational peptides whose expression could cover the
   maximum tumor regions for future immune therapies. We further
   demonstrated that combining tumor-specific peptides might enhance the
   efficacy of immunotherapy in both patient-derived cell (PDC) and
   patient-derived xenograft (PDX) models. PSERP efficiently retains
   precise spatial proteomic information within the tissue context and
   provides a deeper understanding of tissue biology and pathology at the
   molecular level.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s13045-025-01710-5.

   Keywords: Spatial proteomics, Tumor heterogeneity, Tumor
   microenvironment, Glioma

Introduction

   Tissues are composed of various microscopic features, cell types, and
   phenotypically diverse subpopulations, and the spatial localization of
   cells within the tumor and their spatial neighborhood is crucial for
   elucidating their identity and functions [[54]1]. The cellular
   composition of tissue has a substantial influence on measured
   co-expression signals within the molecular profiles of bulk tissue and
   their microenvironment, contributing to the influence of cellular
   function, signaling, and different disease outcomes [[55]2]. Spatial
   molecular profiling of complex tissues is essential for elucidating
   cellular function in physiological and pathological states, providing
   critical insights into disease mechanisms.

   Spatial techniques such as spatial transcriptomics have undergone rapid
   progress, from early spatial detection methods, such as the 10
   × Genomics Visium platform, serial microtomy sequencing (tomo-seq), to
   recent high-resolution spatially resolved transcriptomic methods, such
   as Stereo-seq, SeqScope, and PIXEL-seq [[56]3–[57]5]. Scientists can
   now profile transcript expression in an unbiased and panoramic manner.
   Similar to spatial transcriptomics, emerging methods involving spatial
   genomics, metabolomics and epigenomics are also being developed. For
   instance, spatial ATAC-seq can be performed by in situ Tn5
   transposition and probe ligation via microfluidics devices, followed by
   digestion and sequencing to profile chromatin accessibility [[58]6].
   Moreover, spatially resolved metabolic techniques, such as
   matrix-assisted laser desorption/ionization imaging mass spectrometry
   (MALDI-IMS) [[59]7], hold significant potential for revealing the
   spatial distribution of metabolic species and molecular interactions.

   Despite the rapid development of spatial-resolved techniques for
   transcriptomics and metabolomics, spatial proteomics remains
   challenging due to the limited tools available. Currently, widely used
   spatial proteomic analysis methods include laser capture
   microdissection (LCM)-based methods and antibody-based spatial
   proteomics methods, such as MACSima Imaging Cyclic Staining (MICS)
   [[60]8], tissue-based cyclic immunofluorescence (t-CycIF) [[61]9], and
   codetection by IndEXing (CODEX) [[62]10]. While LCM-based spatially
   resolved proteomics is an elegant approach that enables the microscopic
   isolation of defined regions from tissue-preserving crucial spatial
   information, it is time-consuming, requires sophisticated equipment,
   and is limited by its efficiency. Consequently, its spatial coverage is
   inevitably limited to a comparatively small sample size, typically
   hundreds of micrometers. Antibody-based methods, on the other hand,
   depend on the affinity and specificity of antibodies, which compromises
   their reproducibility and multiplexing capacity (currently limited to
   approximately 50–100 proteins), thus impeding scalability. Therefore,
   there is an urgent need to develop an unbiased, high-throughput, and
   panoramic spatial proteomic approach that facilitates systematic and
   comprehensive analysis of tissue proteomics.

   The rapid development of spatial techniques has provided an opportunity
   to better decipher the complex tumor microenvironment (TME) and enhance
   immunotherapy efficacy [[63]1, [64]6]. Currently, immune checkpoint
   inhibitor (ICI) and chimeric antigen receptor T-cell (CAR-T) therapies,
   which leverage the human immune system to target and eliminate cancer
   cells, represent the most promising therapeutic strategies for treating
   malignant cancers. However, patients with an immune"desert"in their
   microenvironment benefit little from immunotherapy [[65]11].
   Neoantigens, which are derived from nonsynonymous somatic mutations and
   are specific to tumor cells, offer the potential for targeted tumor
   destruction without toxic side effects on normal tissues [[66]12,
   [67]13]. The personalized neoantigen vaccine strategy can be a major
   step in improving the poor tumor microenvironment. However, identifying
   and selecting personalized neoantigens remains a significant challenge.
   Current methods for neoantigen prediction lack expression data [[68]13,
   [69]14], which limits their therapeutic potential and complicates the
   validation process. Developing an improved method for identifying and
   quantifying neoantigens, which is able to demonstrate their
   tumor-specific expression and spatial distribution, could enhance the
   efficacy of neoantigen-based immunotherapies.

   Spatially resolved proteomics using tissue expansion is a cutting-edge
   approach, in which samples are prepared by adhering a tissue section to
   a stretchable hydrogel [[70]15, [71]16]. Then, defined-sized tissue
   fragments can be isolated from regions of interest within the tissue
   section and subsequently subjected to mass spectrometry analysis
   [[72]15]. This approach enhances spatial resolution; however, it is
   technically complex, labor-intensive, prone to potential batch effects,
   and limited by insufficient throughput due to manual processing.
   Additionally, the limited throughput of current tissue expansion-based
   spatial proteomics techniques hinders their integration with spatial
   transcriptomics and single-cell transcriptomics data. These limitations
   hinder the large-scale implementation of unbiased spatial proteomic
   analysis in complex tumor ecosystems.

   Spatial proteomics via tissue expansion significantly enhanced the
   resolution, yet since the existing procedure for tissue expansion-based
   proteomics is very sophisticated and relies on manual labor, it is
   still technically fragile and challenging to implement and has
   insufficient throughput [[73]15]. Limitations in detection throughput
   and potential batch effects due to manual labor have prevented the
   utilization of this approach in the completion of unbiased spatial
   proteomics analysis of larger tumor ecosystems. In parallel, the
   throughput of the current tissue expansion spatial proteomics technique
   is still inadequate, making profiling results difficult to utilize for
   further integrative analysis with data from spatial transcriptomics and
   single-cell transcriptomics.

   To address these challenges, we developed panoramic spatial enhanced
   resolution proteomics (PSERP), a method that combines tissue expansion,
   automated sample segmentation using a three-dimensionally (3D) printed
   cutting toolkit, and tryptic digestion with high-throughput proteomic
   profiling. PSERP enables rapid quantitative profiling of proteomic
   spatial diversity across the entire tissue section at sub-millimeter
   resolution (Methods). In this proof-of-principle study, we validated
   the PSERP method using gliomas with different genomic alterations and
   demonstrated its capability for in-depth characterization of tumor
   ecosystems, including the cellular composition, distribution, and
   cell–cell communication within the TME. PSERP includes a spatially
   resolved tumor-specific peptidome identification workflow, enabling the
   systematic analysis of tumor-specific peptides and their impact on
   immune therapies. By utilizing PSERP, we delineated the distinctive
   proteomic signatures of malignant and nonmalignant regions,
   illustrating the dynamic expression patterns of proteins from tumor
   centers to borders and nonmalignant areas. Integrating this data with
   single-cell transcriptomics, we provided a comprehensive view of the
   tumor ecosystem and cell‒cell interactions in glioma. We also
   identified shared and distinct proteomic features across glioma samples
   with or without IDH1/EGFR mutations and analyzed the spatial expression
   of tumor-specific peptides (potential neoantigens), confirming their
   impact on the tumor microenvironment in both patient-derived cell (PDC)
   and patient-derived xenograft (PDX) models. PSERP has the potential to
   provide novel insights into the molecular landscape of tumors, thereby
   paving the way for future advancements in precision medicine and
   targeted immunotherapy.

Results

Spatially resolved proteomics of glioma using a panoramic spatial enhanced
resolution proteomics (PSERP) approach

   To achieve high-resolution spatially resolved proteomics in tissue
   samples, we developed the PSERP approach, comprising the following
   systematic steps (Fig. [74]1A and S1 A, Methods). In Step 1, we
   performed tissue expansion process, enabling the physical enlargement
   of the tissue while preserving protein localization (Methods). In Step
   2, we employed the custom-designed 3D-printed cutting toolkit, allowing
   the expanded tissue to be precisely sectioned into thousands of
   “voxels” of defined dimensions (Figure S1B-D, Methods). In Step 3,
   these divided voxels were transferred into 96-well plates, with their
   spatial coordinates meticulously recorded to maintain positional
   context (Figure S1E, Methods). Step 4 involved in-gel digestion and
   peptide extraction from each voxel, preparing the samples for
   downstream analysis. Finally, in Step 5, the extracted peptides were
   subjected to high-throughput data-independent acquisition (DIA) liquid
   chromatography-tandem mass spectrometry (LC–MS/MS) for comprehensive
   spatial proteomic profiling, phosphoproteomic profiling, and neoantigen
   profiling (Fig. [75]1A, Methods). To facilitate a comprehensive
   understanding of the PSERP workflow, we have included a detailed
   timeline from tissue preparation to LC–MS/MS analysis in Fig. [76]1B.
   This approach was optimized for time efficiency, and the majority of
   the workflow can be automated on a robotic workstation (Figure S1 A,
   Methods).

Fig. 1.

   [77]Fig. 1
   [78]Open in a new tab

   Overview of the panoramic spatial enhanced resolution proteomics
   (PSERP) workflow and its application. A A schematic representation of
   the PSERP workflow, encompassing the steps involved in neoantigen
   profiling. B Timeline illustrating the duration of each step in the
   PSERP process, from sample processing to LC–MS/MS data acquisition. C
   Hematoxylin and eosin (H&E) stained images and histological annotations
   of samples obtained from glioma patients. The samples include: G1, G4,
   and G5 (IDH1-WT/EGFR-mutant); G2, G6, and G7 (IDH1-mutant); G3, G8, and
   G9 (EGFR & IDH1 wild-type). The dotted line delineates the malignant
   tumor regions

   The primary objective of this proof-of-principle study was to validate
   the capabilities of our PSERP-based spatial proteomics method to
   capture the spatial heterogeneity of tumors and their corresponding
   microenvironments. The mutual exclusivity and microenvironment
   diversities of EGFR- and IDH1-mutant gliomas are well-documented in
   previous studies [[79]17]. EGFR mutations predominantly occur in
   IDH1-wildtype (WT) tumors. Thus, we focused on IDH1-mutant and two
   types of IDH1-wildtype (WT) gliomas: IDH1-WT/EGFR-mutant (G1, G4 and
   G5), IDH1-mutant (G2, G6 and G7), and IDH1/EGFR-double-WT (G3, G8 and
   G9) gliomas.

   For this proof-of-principle study, we sampled sections from a total of
   nine samples (G1-G9), and processed them with PSERP procedures
   (Fig. [80]1C, Methods). The isolated voxels were then processed
   following the PSERP workflow and analyzed using the DIA strategy
   (Methods). A total of 12,190 proteins were identified across 2,230
   voxels from all glioma samples. The median number of proteins
   identified was 1,650 and the standard deviation was 562 across samples.
   Herein, we obtained the panoramic spatial proteomic profiles with
   comprehensive coverage of whole sections of glioma samples.

The fidelity, quantification and reproducibility of PSERP

   In this proof-of-principle study, we consistently applied a sixfold
   expansion factor across all samples, resulting in an isotropic increase
   in length, width, and height. This corresponded to a 216-fold (6 by 6
   by 6) volumetric increase in the sample. To illustrate the changes in
   tissue size before and after expansion, we have provided representative
   images of the samples in their pre-expanded and post-expanded states
   (Figure S2 A). Additionally, we applied PSERP to consecutive replicate
   sections of glioma tissue, and provided the images of the corresponding
   replicate samples after expansion to show the consistency (Figure S2B).
   The replicate samples expanded uniformly by the same expansion factor,
   indicating a high degree of consistency between replicate samples
   before and after expansion. Furthermore, to validate the consistency of
   the expanded tissues, we performed isotropic distortion analysis to
   assess the expansion distortion. Specifically, the pre- and
   post-expanded tissue images were recorded using the Biorad imaging
   system (ChemiDoc, Biorad), a B-spline-based image registration MATLAB
   package was used to register post-PSERP to pre-PSERP image by landmark
   selection [[81]15]. The root-mean-square (RMS) length measurement error
   of feature measurement was calculated to quantify the percentages of
   the distortion in the expanded tissue. In this way, we found that in
   our study, after tissue expansion over length scales up to 1,500 μm
   with the linear expansion factor of six fold, the average of feature
   RMS errors in our samples was 5% (Figure S2 C and S2D), demonstrating
   high fidelity in the preservation of tissue architecture.

   All samples in our study were processed using the same protocol,
   ensuring comparability and reliability across different samples. To
   evaluate the reproducibility, we utilized two consecutive tissue slides
   from one sample block (G1 sample) for the replicate experiment. These
   two consecutive slides were named G1-R1 and G1-R2, respectively (Figure
   S3 A and S3B). We processed these two tissue slides following the same
   protocol, including gel embedding, expanding, cutting, transferring,
   peptide extraction, and MS detection steps. During transferring
   process, we recorded the (x,y) coordinates of all voxels for both G1-R1
   and G1-R2 samples to map them back to their original positions of the
   tissue slide (Methods). Then, we specifically selected 12 voxels within
   the malignant regions and 15 voxels within the non-malignant regions of
   the G1 and G1 replicate samples, respectively. The corresponding voxels
   with the same (x,y) coordinates in G1-R1 and G1-R2 samples were
   utilized to assess the reproducibility between samples (Figure S3 A and
   S3B). We generated proteomic data from selected voxels in G1-R1 and
   G1-R2 samples for further analysis (Figure S3 C). For the G1-R1 sample,
   we identified an average of 1,440 proteins per voxel in the
   non-malignant region, and 1,722 proteins per voxel in the malignant
   region. For the G1-R2 sample, we identified an average of 1,518
   proteins per voxel in the non-malignant region, and 1,729 proteins per
   voxel in the malignant region. Additionally, 93% of the proteins
   identified in the non-malignant voxels of G1-R1 were consistently
   identified in their corresponding replicates of G1-R2, while 97% of
   proteins overlapped between malignant voxels. This substantial overlap
   supported the reproducibility of our method. We observed high
   correlation coefficients (R > 0.92) between corresponding voxels in
   G1-R1 and G1-R2 samples in both the malignant region and non-malignant
   region (Figure S3D). Clustering analysis revealed a clear distinction
   between voxels from malignant and non-malignant regions (Figure S3E).
   This finding indicated that the differential proteomic patterns were
   consistently preserved in both G1-R1 and G1-R2 samples.

   Additionally, the evaluation of protein diffusion between voxels was
   necessary for assessing the reliability of the method. In light of
   this, we assessed protein diffusion by sampling a homogeneous glioma
   tissue. After the same processing protocol mentioned above, we
   collected proteomic data from voxels containing the tissue sample and
   adjacent blank gel voxels without any tissue. Then, we compared the
   results between the tissue gel voxels and the adjacent blank gel
   voxels. The results showed that an average of 2,641 proteins were
   identified in tissue gel voxels (n = 10), while only an average of 57
   proteins were identified in blank gel voxels (n = 10) (Figure S3 F). We
   observed that marker proteins, such as GFAP, EGFR, VIM, and PARP1, were
   significantly more abundant in glioma tissue voxels, while they were
   barely detectable in the blank gel voxels (Figure S3G). Although
   several proteins were identified in both blank gel voxels and tissue
   voxels, the abundance of proteins identified in blank voxels was
   significantly lower than that in tissue voxels (lower than 1%) (Figure
   S3H). These results further suggested that protein diffusion between
   voxels was nearly negligible. Moreover, voxels with fewer than 500
   identified proteins were excluded from further analysis in our study
   (Methods).

   Collectively, by utilizing the PSERP procedure, we could obtain the
   panoramic spatial proteomic profiles covering whole pathological
   sections with deep coverage and high quantitative reproducibility,
   accurately illustrating the diverse proteomic features of malignant and
   nonmalignant regions within one tumor section.

Spatial proteomic profiling of malignant and non-malignant regions in glioma

   To assess whether the PSERP method could capture spatial heterogeneity,
   we first conducted an unsupervised Louvain clustering analysis using
   the spatial proteomic data from the IDH1-WT/EGFR-mutant G1 sample
   (Methods). The results revealed six Louvain clusters among all voxels
   (Fig. [82]2A). These six clusters could be further partitioned into two
   major groups by applying hierarchical clustering (Fig. [83]2B and C).
   The principal component analysis (PCA) revealed clear separation
   between the malignant and non-malignant clusters (Fig. [84]2D). To
   explore the relationship between the clusters and the glioma sample’s
   spatial information, we mapped the clusters to their original tissue
   positions (Fig. [85]2E). Interestingly, when combined with H&E
   staining, the clusters corresponded precisely to the malignant and
   non-malignant regions of the sample. Specifically, cluster 2 (C2) and
   cluster 3 (C3) were identified as malignant clusters, while cluster 0
   (C0), cluster 1 (C1), cluster 4 (C4), and cluster 5 (C5) were
   designated as non-malignant clusters (Fig. [86]2E). These results
   demonstrated the capability of identifying distinct spatial regions
   within a single sample section based solely on the spatial proteome.

Fig. 2.

   [87]Fig. 2
   [88]Open in a new tab

   Spatial proteomic dissection of malignant and non-malignant regions in
   IDH1-WT/EGFR-mutant gliomas using PSERP. A Uniform Manifold
   Approximation and Projection (UMAP) plots of spatial spots from the
   IDH1-WT/EGFR-mutant G1 sample, colored according to their Louvain
   clustering identities. B Phylogenetic tree illustrating proteomic
   similarities across clusters in the IDH1-WT/EGFR-mutant G1 sample. C
   UMAP plots of spots from the IDH1-WT/EGFR-mutant G1 sample, with colors
   representing malignant and non-malignant regions based on histological
   annotations. D Principal Component Analysis (PCA) plot comparing the
   proteomic profiles of malignant and non-malignant spots from the
   IDH1-WT/EGFR-mutant G1 sample. E H&E staining (left), spatial cluster
   distribution (middle), and spatial annotations of malignant versus
   non-malignant regions (right) for the IDH1-WT/EGFR-mutant G1 sample. F
   Volcano plot depicting differentially expressed proteins between
   malignant and non-malignant clusters from the IDH1-WT/EGFR-mutant G1
   sample (FDR < 0.05). G Bar plots showing the enrichment pathways of
   upregulated proteins (FDR < 0.05, Log2 FC > 1.5) in malignant (left)
   and non-malignant (right) clusters from the IDH1-WT/EGFR-mutant G1
   sample. H Heatmaps illustrating the spatial distribution of
   differentially expressed proteins between malignant and non-malignant
   regions in the IDH1-WT/EGFR-mutant G1 sample

   Further comparative analysis and enrichment analysis revealed that
   upregulated proteins in malignant clusters were enriched in cancer
   hallmark pathways, including MYC targets, MTORC1 signaling, E2 F
   targets, and mitotic spindle pathways (Fig. [89]2F-G, Table S1). In
   contrast, proteins upregulated in non-malignant clusters were enriched
   in coagulation, adipogenesis, and metabolism pathways (Fig. [90]2F-G,
   Table S1). In addition, spatial differential expression analysis
   revealed that glioma malignancy-related proteins, such as VIM, GFAP,
   PARP1 and EGFR, were significantly upregulated in the malignant
   regions. Conversely, immune-associated proteins, such as C1QB, and
   adipogenesis-associated proteins, such as APOA2, APOA4, and APOC2, were
   significantly upregulated in the non-malignant regions of the
   IDH1-WT/EGFR-mutant sample (Fig. [91]2H). Moreover, comparative and
   enrichment analyses validated the reproducibility of spatial expression
   profiles between corresponding voxels in G1-R1 and G1-R2 samples
   (Figure S3I and S3 J; Table S1). Specifically, proteins upregulated in
   the malignant region of replicates showed enrichment in the Myc target
   V1 and G2M checkpoint pathways (Figure S3 J). In contrast, proteins
   upregulated in the non-malignant region showed enrichment in
   coagulation and complement pathways (Figure S3 J). These findings
   demonstrated that the differential expression of proteins and pathways
   observed in the IDH1-WT/EGFR-mutant G1 sample was consistent and
   reproducible across replicates. Overall, the strong quantitative
   reproducibility across replicates underscored the robustness and
   reliability of our PSERP approach.

   To minimize the impact of interpatient heterogeneity on the diverse
   proteomic features of malignant and non-malignant regions identified by
   PSERP, we further analyzed the differences between malignant and
   non-malignant regions in three IDH1-WT/EGFR-mutant samples (G1, G4 and
   G5). The results of enrichment analysis and differentially expressed
   proteins were consistent across three distinct IDH1-WT/EGFR-mutant
   samples (Figure S4 A-B). We also performed paired comparisons between
   malignant and non-malignant clusters across the three
   IDH1-WT/EGFR-mutant samples (G1, G4, and G5 samples). As a result, the
   expression levels of tumor-related pathways, such as Myc targets, E2 F
   targets, and mitotic spindle, were increased in the malignant regions,
   while the expression levels of immune or metabolism-related pathways,
   such as coagulation, inflammatory response, and adipogenesis pathways,
   were increased in the non-malignant regions (Figure S4 C), indicating a
   consistent pattern of pathway alterations between malignant and
   non-malignant regions. The consistency of these findings across
   distinct IDH1-WT/EGFR-mutant samples highlighted the robustness of our
   PSERP method in capturing tumor heterogeneity, and further supporting
   the reproducibility of spatial proteomic variations in
   IDH1-WT/EGFR-mutant gliomas.

High-resolution spatial proteomic heterogeneities in different regions of
gliomas

   To evaluate whether our PSERP approach could comprehensively
   characterize the molecular diversity of distinct tumor regions at high
   resolution, we conducted a comparative analysis of the six clusters in
   the IDH1-WT/EGFR-mutant G1 sample (C0-C5) and performed pathway
   enrichment analysis on proteins that exhibited differential expression
   patterns within each cluster (Fig. [92]3A-3E and S5 A-C, Methods). The
   results revealed that immune-related pathways were significantly
   enriched in cluster C0, heme metabolism pathways were enriched in
   cluster C1, and clusters C2 and C3 exhibited similar patterns of
   pathway enrichment, primarily involving the G2M checkpoint and mRNA
   metabolic process pathways (Fig. [93]3B). Clusters C4 and C5 displayed
   enrichment in complement and coagulation cascades, complement
   activation pathways, and other related pathways (Fig. [94]3B).

Fig. 3.

   [95]Fig. 3
   [96]Open in a new tab

   Spatially resolved proteomic profiling of gliomas reveals intra-tumor
   heterogeneity. A H&E staining images and spatial proteomic cluster
   distribution in the IDH1-WT/EGFR-mutant G1 sample (left panel; scale
   bar, 1 mm). The details of each regions in the IDH1-WT/EGFR-mutant G1
   sample (right panel; scale bar, 100 μm). B Dot plots illustrating the
   enrichment of pathways associated with significantly upregulated
   proteins across six spatial clusters in the IDH1-WT/EGFR-mutant G1
   sample. C Heatmaps and boxplots illustrating significantly upregulated
   pathways in non-malignant clusters (C0, C1, C4, and C5) compared to
   malignant clusters (C2 and C3) in the IDH1-WT/EGFR-mutant G1 sample
   (adjusted p value < 0.05). D Heatmaps and boxplots showing
   significantly upregulated pathways in cluster C2 compared to cluster C3
   in the IDH1-WT/EGFR-mutant G1 sample (adjusted p value < 0.05). E
   Heatmaps and boxplots showing significantly downregulated pathways in
   cluster C2 compared to cluster C3 in the IDH1-WT/EGFR-mutant G1 sample
   (adjusted p value < 0.05). F Heatmap demonstrating differentially
   expressed proteins related to key pathways between clusters C2 and C3
   of the IDH1-WT/EGFR-mutant G1 sample (FDR < 0.05). G Heatmaps and
   boxplots illustrating significant differences in the immune score,
   immune cell types, and stromal score between clusters C3 and C2 in the
   IDH1-WT/EGFR-mutant G1 sample (FDR < 0.05)

   Importantly, by integrating pathological information from H&E images,
   we found that the spatial distributions of the six clusters exhibited
   significant diversity: clusters C2 and C3 were enriched in malignant
   regions, whereas clusters C0, C1, C4 and C5 were enriched in the
   non-malignant regions (Fig. [97]3A). Consistent with their regional
   distributions, we observed that more tumor-related pathways were
   significantly upregulated in malignant regions (clusters C2 and C3)
   than in non-malignant regions (clusters C0, C1, C4, and C5), whereas
   the tumor microenvironment and immune-related pathways were
   significantly upregulated in non-malignant regions (clusters C0, C1,
   C4, and C5) (Fig. [98]3B, C and S5 A).

   Although both clusters C2 and C3 were enriched in malignant tumor
   tissue regions, their separation into two distinct clusters suggested
   molecular heterogeneity between spatially different malignant tumor
   regions within the same sample. To further explore the heterogeneity
   between clusters C2 and C3, we compared the pathway ssGSEA scores of
   spots within these clusters to identify differentially enriched
   pathways (Figure S5B, Methods). The analysis revealed that Myc targets
   and G2M checkpoint pathways were significantly upregulated in cluster
   C2 (Fig. [99]3D), suggesting enhanced tumor cell proliferation in this
   cluster. This was further confirmed by the elevated multi-gene
   proliferation score (MGPS) in cluster C2 (Fig. [100]3D and F). In
   contrast, cluster C3 exhibited elevated expression of proteins involved
   in fatty acid metabolism, the interferon alpha response, and antigen
   processing and presentation pathways was detected (Fig. [101]3E and F).
   This result was consistent with the higher immune score and increased
   scores of immune cell types in cluster C3 (Fig. [102]3G, Methods). In
   parallel with the results of ssGSEA analysis, the upregulation of cell
   proliferation-associated molecules such as MRTOR4, CDK6, MCM5, MCM6,
   and AURKB was observed in cluster C2, while proteins related to fatty
   acid degradation (HADH, ACADVL, and ACAT2), the tricarboxylic acid
   (TCA) cycle (LDHA, PDHB, and ACO2), the interferon alpha response
   (PSMB8/9, PSME1/2, and STAT2), and antigen presentation (HLA-DQA1, DRA,
   and DRB1) were significantly upregulated in cluster C3 compared to
   those in cluster C2 (Fig. [103]3F, Methods). These results confirmed
   the spatial heterogeneity of glioma tumor regions, suggesting that
   cluster C2 might have a stronger proliferative capacity, whereas
   cluster C3 exhibited a more pronounced tumor immune response.

   In summary, our findings revealed glioma spatial heterogeneity through
   PSERP. By integrating the H&E pathology slides, we observed distinct
   histological characteristics within different clusters. The analysis of
   significantly upregulated proteins and pathway enrichment allowed us to
   identify molecular features associated with each cluster. Notably, the
   differences observed in the malignant and non-malignant regions on the
   H&E images and the pathway enrichment patterns among clusters were
   highly consistent. Furthermore, the heterogeneity between the two
   malignant tumor clusters, C2 and C3, was characterized by differences
   in proliferation- and immune-related pathways, suggesting the presence
   of distinct molecular subtypes within spatially distinct malignant
   tumor regions. These findings provided valuable insights into the
   spatial complexity and molecular characteristics of glioma, and
   contributed to a deeper understanding of this heterogeneous disease.

Characterization of Tumor Border Characteristics and Interactions with the
Tumor Microenvironment Using Spatial Proteomics

   Additionally, to illustrate the diverse proteomic features along the
   malignant-border-nonmalignant axis in the glioma sample, our
   high-resolution PSERP approach enabled a more precise exploration of
   unique protein expression patterns and microenvironmental features
   relative to the borders between malignant and non-malignant tumor
   regions. Based on the tumor borders identified by pathologists using
   the H&E-stained image of the IDH1-WT/EGFR-mutant sample, we divided all
   the voxels of the IDH1-WT/EGFR-mutant sample into three groups: tumor
   center, border, and paratumor (Fig. [104]4A, Methods). To assess the
   alignment between pathological-annotated tumor boundaries and those
   defined by PSERP data, we also applied Cottrazm (Construction of Tumor
   Transition Zone Microenvironment) [[105]18] using the spatial proteomic
   data to further delineate the malignant region, tumor border, and
   non-malignant region of the IDH1-WT/EGFR-mutant sample (Figure S5 C and
   S5D, Methods). We then compared the tumor border identified through H&E
   histological annotation with the border determined by Cottrazm. The
   results demonstrated a high degree of consistency between the tumor
   borders predicted by Cottrazm and those annotated by pathology (Figure
   S5E and S5 F). To understand the proteomic feature related to the tumor
   core spots, we calculated pathway enrichment scores based on spatial
   proteomic data of the IDH1-WT/EGFR-mutant G1 sample, and assessed their
   correlations with inferred copy number variation (CNV) scores. The
   results indicated that the hypoxia signaling pathway was the most
   significant correlated with inferred CNV scores (p-value = 1.6e-7,
   Figure S5G). Consistently, the hypoxia pathway scores increased at the
   tumor core within the Cottrazm-identified malignant regions and
   gradually decreased with the increasing distance from the tumor core
   (Figure S5H).

Fig. 4.

   [106]Fig. 4
   [107]Open in a new tab

   The spatial characteristics along the tumor border. A Strategy to
   quantify the Euclidean distance of spatial spots from the malignant
   tumor border. The red line delineates the malignant tumor border. B
   Heatmap illustrating the distance of each spot from the malignant tumor
   border. C Line graphs depicting distinct groups of proteins, clustered
   by Mfuzz, to show relative expression changes across the spatial
   proteomic landscape. D Enrichment analysis results of proteins within
   each Mfuzz cluster. E Variation in G2M checkpoint pathway scores of
   spatial spots with respect to the distance from the tumor region to the
   para-tumor region. F Variation in immune scores of spatial spots with
   respect to the distance from the tumor region to the para-tumor region.
   G Heatmap demonstrating the expression level changes of representative
   proteins significantly correlated with distance from the tumor border

   These findings prompted us to then compute the Euclidean distances from
   tumor and paratumor spots to border spots (Fig. [108]5A and B,
   Methods). We then applied the Mfuzz algorithm [[109]19] to identify
   different proteomic expression patterns associated with the distance
   from the border, and to explore potential crosstalk among the tumor,
   border and paratumor regions (Methods). This analysis revealed six
   distinct protein groups representing diverse distances from the border
   (Fig. [110]4C). Among these groups, group 2 represented proteins
   displaying specific higher expression levels at the tumor center
   region, and proteins enriched in group 4 exhibited the highest
   expression within the closest proximity to the border, implying that
   proteins in this group represent molecular features of tumor border
   regions. We subsequently performed gene set enrichment analysis (GSEA)
   on the proteins categorized within each group to characterize the
   proteomic features of the tumor center and tumor border (Fig. [111]4D).
   Specifically, proteins in group 2 were associated with cell cycle
   features, consistent with the elevated expression of G2M checkpoint
   pathway (Fig. [112]4E and G). In contrast, group 4 showed enrichment of
   proteins involved in coagulation, innate immune system and complement
   cascade pathways, suggesting a prominent role for immune-related
   processes in shaping expression dynamics at the tumor border
   (Fig. [113]4F and G). Based on these findings, we investigated immune
   cell infiltration in the IDH1-WT/EGFR-mutant sample using the xCell
   deconvolution algorithm (Methods). The results revealed a gradient of
   immune scores across the malignant and non-malignant regions of the
   tumor (Fig. [114]4F). As the distance from the border increased, the
   immune scores gradually decreased (Fig. [115]4F), suggesting the
   regulatory mechanisms in the border region contribute to an activated
   immune response. Furthermore, our analysis revealed that the expression
   levels of immune-related proteins such as C1QB, C3, IGHG1, IGLC2, F2,
   FGA, and FGB were upregulated with the increasing distance from the
   tumor region (Fig. [116]4G). These spatially dependent alterations in
   protein expression were shown on the spatial proteomic map of the
   IDH1-WT/EGFR-mutant sample (Figure S5 J).

Fig. 5.

   [117]Fig. 5
   [118]Open in a new tab

   Spatial heterogeneities across glioma tissues with distinct mutations
   (IDH1-WT/EGFR-mutant, IDH1-mutant, and IDH1/EGFR-double-WT). A H&E
   staining images, and proteomic cluster maps in different samples from
   glioma samples with distinct mutations: G1(IDH1-WT/EGFR-mutant), G2
   (IDH1-mutant), and G3 (IDH1/EGFR-double-WT) samples. B Significantly
   upregulated marker proteins in the malignant regions compared to
   non-malignant regions for each sample (FDR < 0.05). C Boxplots and
   heatmaps showing ssGSEA scores of pathways upregulated in the malignant
   regions compared to non-malignant regions of the IDH1-mutant G2 sample.
   D Boxplots and heatmaps showing ssGSEA scores of pathways upregulated
   in the malignant regions compared to non-malignant regions of the
   IDH1/EGFR-double-WT G3 sample. E Heatmaps depicting significantly
   upregulated proteins in the malignant regions c non-malignant regions
   in the IDH1-mutant G2 sample (FDR < 0.05). F Heatmaps showing
   significantly upregulated proteins in the malignant region
   comparedompared to non-malignant regions in the IDH1/EGFR-double-WT G3
   sample (FDR < 0.05). G Heatmaps and boxplots illustrating significantly
   downregulated xCell scores for immune cell types in the malignant
   regions compared to non-malignant regions of the IDH1-mutant G2 sample.
   H Heatmaps and boxplots illustrating significantly upregulated xCell
   scores for immune cell types in the malignant region compared to
   non-malignant regions of the IDH1/EGFR-double-WT G3 sample

High-resolution spatial proteomic heterogeneities revealed intra- and
inter-tumor heterogeneity in glioma

   To evaluate whether our PSERP approach could robustly reflect the
   diverse spatial proteomic heterogeneity of different samples, we
   conducted spatial proteomic profiling of glioma samples with distinct
   mutational statuses (IDH1-WT/EGFR-mutant, IDH1-mutant, and
   IDH1/EGFR-double-WT gliomas). We performed unsupervised Louvain
   clustering analysis for the IDH1-mutant G2 sample and EGFR & IDH1-WT G3
   sample utilizing the same method as we used for the EGFR-mutant G1
   sample (Methods). This analysis resulted in 7 clusters (G2 C0 to G2 C6)
   across 405 voxels of the IDH1-mutant G2 sample, and 6 clusters (G3 C0
   to G3 C5) across 630 voxels of the IDH1/EGFR-double-WT G3 sample
   (Fig. [119]5A). Using the H&E image-based histological annotations of
   the G1, G2 and G3 samples, we found that in the IDH1-WT/EGFR-mutant G1
   sample, clusters G1 C2 and G1 C3 were enriched in malignant regions. In
   the IDH1-mutant G2 sample, clusters G2 C1 and G2 C3 were enriched in
   malignant regions. In the IDH1/EGFR-double-WT G3 sample, clusters G3
   C0, G3 C1, G3 C4, and G3 C5 were enriched in the malignant regions
   (Fig. [120]5A). These results indicated that the assigned clusters of
   the G1, G2, and G3 samples reflected the spatial diversity between
   malignant and non-malignant regions of gliomas. Previous proteogenomic
   studies have demonstrated that genomic alterations can influence
   downstream protein expression. For instance, IDH1 mutations can
   upregulate proteins related to the hypoxia and glutamate metabolism
   pathways [[121]20, [122]21]. In parallel, mutations in EGFR can
   increase the expression of its cognate protein and proteins genes
   involved in the TCA cycle [[123]22]. To explore the spatial
   distributions reflecting the distinctive or common proteomic features
   of gliomas with different mutational statuses, we examined the spatial
   expression patterns of proteins known to be altered, such as GFAP,
   EGFR, PARP1, and GLUD1 [[124]23–[125]26], among IDH1-WT/EGFR-mutant,
   IDH1-mutant and IDH1/EGFR-double-WT samples. Specifically, we observed
   significant upregulation of key proteins in the malignant regions
   across distinct samples (Fig. [126]5B), which could reflect the
   underlying mutation statuses and the tumor’s biological
   characteristics. A significant upregulation of EGFR was detected in the
   malignant regions of IDH1-WT/EGFR-mutant samples (G1, G4, and G5),
   consistent with the well-established role of EGFR mutations in glioma
   pathogenesis [[127]27] (Fig. [128]5B and Figure S6 A). In IDH1-mutant
   gliomas (G2, G6, and G7), we observed a significant upregulation of
   GLUD1 in the malignant regions (Fig. [129]5B and Figure S6B). GLUD1, a
   key enzyme in the glutamate-glutamine cycle, plays a crucial role in
   maintaining cellular metabolic homeostasis. Its upregulation in
   IDH1-mutant gliomas suggested a compensatory mechanism to maintain
   metabolic flux, further supporting the metabolic reprogramming in these
   tumors [[130]28, [131]29]. PARP1, a DNA repair enzyme involved in
   repairing single-strand DNA breaks, was significantly upregulated in
   the malignant regions of both IDH1-WT/EGFR-mutant gliomas (G1, G4, and
   G5) and IDH1/EGFR-double-WT gliomas (G3, G8, and G9) (Fig. [132]5B,
   Figure S6 A and Figure S6 C). The enhanced expression of PARP1 in these
   regions reflected the increased DNA damage and the reliance of tumor
   cells on DNA repair mechanisms to sustain their uncontrolled
   proliferation and survival [[133]30]. These findings suggested a
   possible vulnerability in these gliomas, which could be targeted for
   therapeutic intervention. Additionally, GFAP, a key protein of
   astrocytic differentiation [[134]31], was upregulated in the malignant
   regions of all three glioma mutation statuses, including
   IDH1-WT/EGFR-mutant, IDH1-mutant, and IDH1/EGFR-double-WT gliomas
   (Fig. [135]5B and Figure S6 A -S6 C). The increase in GFAP expression
   across all mutation statuses indicated its importance as a general
   marker of malignancy in gliomas, regardless of their mutational
   profiles. The upregulation of key proteins, such as EGFR, GLUD1, PARP1,
   and GFAP reflected critical aspects of glioma progression and tumor
   malignancy, including cell proliferation, metabolic reprogramming, DNA
   repair, and reactive gliosis.

   In addition, we conducted pathway enrichment analysis to more
   comprehensively depict inter-tumor heterogeneity across gliomas with
   different mutations (Table S4). We compared the differential expression
   of proteins between the malignant and non-malignant regions of the
   glioma samples characterized by distinct mutational profile.
   Specifically, we found that oxidative stress-induced senescence and
   glutathione catabolic process pathways were upregulated in the
   malignant regions of the IDH1-mutant samples (G2, G6, and G7)
   (Fig. [136]5C and Figure S7 A-B). In contrast, coagulation, interferon
   α response, and interferon γ response pathways were upregulated in the
   malignant regions of the IDH1/EGFR-double-WT samples (G3, G8, and G9)
   (Fig. [137]5D and Figure S7 C-D). Consistently, the expression levels
   of proteins, such as H2BC5, H2BC1, MINK1, H2 AC4, MAPK8, H2 AX, and
   CEACAM8, were significantly higher in the malignant regions compared
   with the non-malignant regions of the IDH1-mutant samples (G2, G6, and
   G7) (Fig. [138]5E and Figure S7E). Immune-related proteins, including
   CD9, ANXA1, C3, C9, HLA-A, HLA-B and HLA-C, were upregulated in the
   malignant regions of IDH1/EGFR-double-WT samples (G3, G8, and G9)
   (Fig. [139]5F and Figure S7 F). Furthermore, the xCell deconvolution
   results (Methods) revealed that the spatial distribution of immune
   cells between the non-malignant and malignant regions was distinct in
   IDH1/EGFR-double-WT samples compared to IDH1-mutant samples.
   Specifically, immune cells, including M2 macrophages, cDC, iDC, and CD8
   + T effector memory cells, were significantly upregulated in the
   malignant regions of IDH1/EGFR-double-WT samples (G3, G8 and G9)
   (Fig. [140]5G and Figure S8 A). In contrast, the immune score and
   scores for immune cells, including B cell and plasma cells, were
   significantly downregulated in the malignant regions of the IDH1-mutant
   samples (G2, G6 and G7) (Fig. [141]5H and Figure S8B), consistent with
   the inhibitory effect of IDH1 mutations on immune infiltration
   [[142]32].

   Given that PSERP approach can capture the proteomic characteristics of
   whole tissue sections without bias, we systematically compared the
   molecular features of both malignant and non-malignant regions in
   glioma samples with different genomic mutations. These results
   demonstrated that our PSERP approach accurately reflected the
   distinctive features of both intra- and inter-tumor heterogeneity in
   glioma samples. This approach also enabled us to analyze how different
   mutations influence the tumor and its surrounding microenvironment.

Spatial phosphoproteomic profiling of gliomas revealed spatial heterogeneity
of kinase-substrate signal transduction

   Besides elucidating the spatial proteomic features, our PSERP can also
   portray the spatial distributions of phosphoproteomic features
   (Methods). Specifically, based on phosphoproteomic data from our
   previously published article [[143]17], we built a phospho-spectral
   library using the FragPipe computational platform [[144]33]. The raw
   DIA data generated in our research by PSERP were then searched against
   the phospho-spectral library. The results revealed an average of 1,177
   phosphosites in the IDH1-WT/EGFR-mutant sample, 1,447 phosphosites in
   the IDH1-mutant sample, and 1,310 phosphosites in the
   IDH1/EGFR-double-WT sample (Figure S8 C-S8H). Intriguingly, compared to
   the non-malignant regions, more phosphosites were identified in the
   malignant regions (Figure S8 C-S8H).

   To illustrate the most active kinases in the malignant regions among
   the three glioma samples, we conducted the KSEA analysis. The results
   revealed that, in the sample with EGFR mutation, the most significantly
   activated kinase in malignant regions was MAPK7 (Figure S8I and S8 J).
   Previously, we demonstrated that MAPK7 could promote
   IDH1-WT/EGFR-mutant gliomas’ tumor growth through phosphorylating
   enzymes (PRPS1, PRPS2) that catalyzed the phosphoribosylation of ribose
   5-phosphate to 5-phosphoribosyl-1-pyrophosphate [[145]17]. We then
   evaluated the spatial distribution of phospho-substrates of MAPK7 and
   observed that phospho-substrates of MAPK7, including PRPS1_T255 and
   PRPS2_S41, which were involved in the nucleotide metabolism pathway,
   showed consistent expression patterns with the kinase activity of MAPK7
   (Figure S8 K and S8L).

   Intriguingly, apart from MAPK7, CDK2 was also observed to be activated
   in the malignant regions of the IDH1-WT/EGFR-mutant sample.
   Importantly, we found that although both MAPK7 and CDK2 were activated
   in the malignant regions, they exhibited spatial heterogeneity in
   kinase activity and protein expression. Specifically, CDK2 was
   activated in the center of malignant regions, while MAPK7 was activated
   in the peritumoral regions. Along with these findings, the
   phospho-substrates of CDK2, including XRCC6_T455 and CDK7_T170, showed
   enhanced expression in the malignant regions, specifically in the tumor
   centroid of IDH1-WT/EGFR-mutant sample (Figure S8M and S8 N).
   Furthermore, the pathway enrichment analysis indicated that the tumor
   center exhibited heightened activity in CDK-related G2M cell cycle
   pathways, while the boundary regions showed increased activation of
   EGFR-MAPK-related pathways (Figure S8O and S8P).

   For the sample with IDH1 mutation, the GSK3B kinase showed enhanced
   activity in the malignant regions compared to the non-malignant
   regions. Accordingly, the spatial kinase activity of GSK3B was also
   observed to show elevation in the malignant regions (Figure S9 A and
   S9B). To elucidate the phospho-signaling pathways activated by GSK3B,
   we investigated the spatial expression of its phospho-substrates. As a
   result, the phospho-substrate that presented the most significant
   consistent spatial expression patterns with GSK3B was HIF-1A_T554
   (Figure S9 C). Consistently, the hypoxia pathway was also observed to
   be activated in the malignant regions of the IDH1 mutant sample at the
   spatial proteomic level (Figure S9D). Thus, spatial phosphoproteomic
   data generated by PSERP revealed a mechanism through which the
   activation of the GSK3B-HIF-1 A phosphorylation signaling pathway could
   lead to the hypoxia features in the malignant regions of the
   IDH1-mutant sample.

   For the IDH1/EGFR-double-WT sample, we observed that the most
   significantly activated kinase in the malignant regions was PDGFRA
   (Figure S9E and S9 F). The activity of PDGFRA was also elevated in the
   malignant regions of the sample (Figure S9 F). Intriguingly, the
   phospho-substrates of PDGFRA, such as MUC1_Y1218 and ABL2_Y272,
   exhibited spatial expression patterns consistent with that of PDGFRA,
   and were enriched in the angiogenesis pathway (Figure S9G and S9H).
   Additionally, by calculating the micro-vessel density (MVD) score, we
   observed significantly higher MVD score in the malignant regions of the
   IDH1/EGFR-double-WT sample compared to the malignant regions of the
   IDH1-WT/EGFR-mutant and IDH1-mutant samples (Figure S9I), further
   confirming the angiogenesis characteristics of the malignant regions in
   the IDH1/EGFR-double-WT sample. In conclusion, our data revealed the
   concordance between spatial phosphoproteomic data and proteomic data in
   portraying the specific features of diverse tumor regions. These
   results underscored the ability of our PSERP method to systematically
   elucidate the spatial heterogeneity of tumors.

Profiling of the glioma ecosystem reveals region-specific cellular
composition and cell‒cell communications

   Previous single-cell transcriptomic studies of gliomas have provided
   insight into the intra-tumoral heterogeneity and dynamic cellular
   plasticity in gliomas [[146]34–[147]36]. Combining spatial omic data
   with single- nuclei transcriptomic data could offer a more
   comprehensive understanding of the glioma ecosystem. However, previous
   spatial proteomic studies were limited by insufficient protein coverage
   or detection throughput, which hindered their interactive analysis with
   single-cell transcriptomic data. Using our PSERP approach, which can
   provide in-depth proteomic profiling, we analyzed the multicellular
   ecosystems in different glioma samples. Specifically, we performed
   single-nuclei transcriptomic analysis on the same glioma samples (the
   IDH1-WT/EGFR-mutant sample, the IDH1-mutant sample, and the
   IDH1/EGFR-double-WT sample) used in our spatial proteomic analysis.
   After performing quality control (Methods), we profiled a total of
   31,872 cells from three glioma specimens with different mutational
   statuses (Fig. [148]6A-B), with an average of 6,781 genes detected per
   cell.

Fig. 6.

   [149]Fig. 6
   [150]Open in a new tab

   The integrative analysis of spatial tumor ecosystem. A UMAP plot of all
   cells from the IDH1-WT/EGFR-mutant G1 sample, IDH1-mutant G2 sample,
   and IDH1/EGFR-double-WT G3 sample, colored by their sample source. B
   UMAP plot of major cell clusters. These plot were colored based on
   deconvoluted cell types using the sc-RNAseq data. C The compositions of
   diverse malignant cell types (upper) and non-malignant immune cell
   types (lower) in the IDH1-WT/EGFR-mutant G1 sample, IDH1-mutant G2
   sample, and IDH1/EGFR-double-WT G3 sample, respectively. D, F, and H
   Heatmaps showing proteomic cluster maps of the IDH1-WT/EGFR-mutant G1
   sample (D), IDH1-mutant G2 sample (F), and EGFR & IDHI-WT G3 sample
   (H). Spots was colored by their clusters. E, G, and I Heatmaps showing
   spatial distribution maps of cell types in the IDH1-WT/EGFR-mutant G1
   sample (E), IDH1-mutant G2 sample (G), and IDH1/EGFR-double-WT G3
   sample (I), colored by their cell types. J Heatmaps showing the spatial
   distribution of the levels of EGFR (left) and CS1RF (right) in the
   IDH1-WT/EGFR-mutant G1 sample. K Heatmaps showing the spatial
   distribution of the levels of HIF-1 A (left) and NeuN (right) in the
   IDH1-mutant G2 sample. L Heatmaps showing the spatial distribution of
   the levels of PDGFRA (left) and MPO (right) in the IDH1/EGFR-double-WT
   G3 sample. M, O, and Q The cell-cell interactions among cell types
   indicated by the scRNA-seq data from the IDH1-WT/EGFR-mutant sample
   (M), IDH1-mutant sample (O) and EGFR & IDH1-mutant sample (Q),
   respectively. The line thickness represented the number of significant
   ligand-receptor pairs. N, P, and Q Dot plots indicating the
   correlations between ligands and receptors in the IDH1-WT/EGFR-mutant
   G1 sample (N), the IDH1-mutant G2 sample (P), and the
   IDH1/EGFR-double-WT G3 sample (Q). Cell types enriched in malignant
   regions are labeled in red text, while those enriched in non-malignant
   regions are labeled in blue text. Ligands and receptors upregulated in
   cell types enriched in malignant regions are labeled in red text, while
   those upregulated in cell types enriched in non-malignant regions are
   labeled in blue text

   To explore tumor cellular composition, we performed UMAP analysis and
   classified cells from each sample (Methods). This analysis revealed
   substantial intratumor heterogeneity, identifying 13 clusters in the
   IDH1-WT/EGFR-mutant sample, 7 clusters in the IDH1-mutant sample, and
   14 clusters in the IDH1/EGFR-double-WT sample. We also identified
   cluster-specific marker genes through differential gene expression
   analysis, which helped define the identity of each cell cluster
   (Fig. [151]6C, Methods). Comparative analysis revealed distinct
   malignant and immune cell compositions across samples. In the
   IDH1-WT/EGFR-mutant sample, the predominant malignant cells were
   AC-like (astrocytic-like), while the major non-malignant cells included
   tumor-associated macrophages (TAMs), cancer-associated fibroblasts
   (CAFs), CD8 + T cells, and microglia. In the IDH1-mutant sample,
   MES-like (mesenchymal-like) cells predominated, with smaller
   proportions of NPC-like and OPC-like (oligodendrocytic precursor)
   cells. The major non-malignant cells were astrocytes, oligodendrocytes,
   and neurons. In the IDH1/EGFR-double-WT sample, the major malignant
   cells were OPC-like, and the major non-malignant cells were
   neutrophils, endothelial cells, and CD4 + T cells.

   To characterize the spatial cellular composition of each sample, we
   integrated the spatial proteomic data with single-nuclei transcriptomic
   data using the Bulk2Space algorithm (Methods). The results revealed
   distinct cellular compositions across different spatial clusters within
   each sample. Specifically, in the IDH1-WT/EGFR-mutant sample, we
   observed significant enrichment of malignant cell types, such as
   AC-like cells, in the malignant regions, while non-malignant cell
   types, such as TAMs, were predominantly enriched in the non-malignant
   regions (Fig. [152]6D and E). In the IDH1-mutant sample, MES-like cells
   were enriched in the malignant regions, while neurons were enriched in
   the non-malignant regions (Fig. [153]6F and G). In the
   IDH1/EGFR-double-WT sample, OPC-like cells were enriched in the
   malignant regions, while neutrophils were more abundant in the
   non-malignant regions (Fig. [154]6H and I). These findings suggested
   different cell type distribution patterns between malignant and
   non-malignant regions were observed across IDH1-WT/EGFR-mutant,
   IDH1-mutant and IDH1/EGFR-double-WT samples.

   We next portrayed the spatial distribution of cell type markers within
   each sample. In the IDH1-WT/EGFR-mutant sample, we observed elevated
   expression of the AC-like marker EGFR in the malignant regions, while
   the non-malignant regions showed enhanced expression of CSF1R
   (Fig. [155]6J). In the IDH1-mutant sample, we observed dominant
   expression of the MES-like marker HIF-1A in the malignant regions,
   whereas the neuronal marker NeuN was enriched in the non-malignant
   regions (Fig. [156]6K). In the IDH1/EGFR-double-WT sample, the OPC-like
   marker PDGFRA was predominantly expressed in the malignant regions,
   while the neutrophil marker MPO showed significantly higher expression
   in the non-malignant regions (Fig. [157]6L).

   To explore the cellular relationships between malignant and
   non-malignant cell types, we analyzed cell–cell interactions by
   examining ligand-receptor pairs using CellPhoneDB [[158]37] based on
   single-cell transcriptomic data (Methods). Notably, further analysis
   revealed key ligand-receptor interactions between malignant and
   non-malignant cells across samples with different mutational statuses.
   Specifically, in the IDH1-WT/EGFR-mutant sample, AC-like tumor cells
   were observed to interact with non-tumor cells, including CAFs, through
   ligand-receptor pairs such as PLG-ITGB2 and PLAU-ITGB2 (Fig. [159]6M
   and 6 N). AC-like tumor cells also interacted with TAMs via
   ligand-receptor pairs C3-CD74, C3-CD81, and EGF-EGFR (Fig. [160]6O and
   P). In the IDH1-mutant sample, MES-like tumor cells were found to
   interact with non-tumor cells, including oligodendrocytes, through
   ligand-receptor pairs like ANXA1-FPR1 A and PON2-HTR2 A (Fig. [161]6Q
   and R). In the IDH1/EGFR-double-WT sample, OPC-like tumor cells
   interacted with non-tumor cells, including endothelial cells, through
   ligand-receptor pairs VEGFA-PDGFRA and PDGFA-PDGFRA (Fig. [162]6S and
   T).

Spatial integrative analysis revealed the impact of tumor microenvironment on
glioma tumorigenesis

   Building on these findings, we sought to elucidate whether the spatial
   distribution of distinct cell types in the tumor microenvironments
   impacted tumor growth. First, we performed survival analysis using our
   previously published glioma multi-omic dataset [[163]17]. The results
   revealed the enrichment of TAMs was significantly correlated to glioma
   patient survival (Figure S10 A), suggesting that higher TAM
   infiltration contributes to more aggressive gliomas. In the
   IDH1-WT/EGFR-mutant sample, combining PSERP-based spatial proteomic
   data with single-nuclei transcriptomic data, we found that TAMs were
   dominantly enriched in the tumor border regions of the
   IDH1-WT/EGFR-mutant sample (Figure S10B and S9 C). The TAM-associated
   protein, CSF1R, was observed to have elevated expression in the tumor
   border regions of the IDH1-WT/EGFR-mutant sample (Figure S10D). We also
   conducted multiplexed IHC (mIHC) staining for CSF1R + TAMs (TAM marker:
   CSF1R) and EGFR + tumor cells (AC-like tumor cell marker: EGFR). As a
   result, significant enrichments of CSF1R + TAMs were observed in
   non-malignant regions, specifically in the tumor border regions of the
   IDH1-WT/EGFR-mutant sample (Figure S10E). Further analysis focused on
   how TAM influenced the tumor cell growth. To this end, we screened for
   TAM-secreted factors potentially associated with glioma malignant
   growth using single-nuclei transcriptomic data. A total of 281
   transcripts were identified as specifically upregulated in TAMs, among
   which EGF was most significantly upregulated in TAMs (Figure S10 F-G).
   We then applied survival analysis using our glioma multi-omic data. The
   protein expression of EGF showed the most significantly negative
   correlation with patients’ overall survival (Figure S10H). These data
   suggested that EGF might be secreted by TAMs and that its expression
   was associated with poor survival in gliomas. In accordance with this
   finding, cell–cell interaction analysis also confirmed a positive
   bidirectional interaction between TAMs and AC-like glioma tumor cells
   through the EGF-EGFR ligand-receptor pair (Figure S10I). We then
   examined the spatial distributions of TAMs and AC-like glioma tumor
   cells using the Bulk2Space algorithm. The analysis revealed that the
   malignant tumor region near the tumor border were predominantly
   enriched with AC-like glioma tumor cells (Figure S10 J and S10 K).

   Furthermore, we investigated how EGF-EGFR signaling promoted glioma
   cell growth. In previous work, we demonstrated that EGFR could promote
   glioma tumor cell growth by activating MAPK7-mediated phosphorylation
   signaling pathway [[164]17]. We found that both EGFR and MAPK7 were
   enriched in the malignant tumor regions (Figure S10L and S10M).
   Correlation analysis revealed a positive correlation between the
   spatial expression levels of EGFR and MAPK7 (Figure S10 N). To further
   validate the impact of EGF-EGFR signaling on glioma tumor cell growth,
   we constructed EGFR overexpressed glioma tumor cells (EGFR-OE-U87).
   These cells were stimulated with EGF or left untreated as a control,
   and tumor growth was assessed using a CCK-8 assay. The results
   demonstrated that EGF stimulation significantly increased tumor cell
   growth in EGFR-OE-U87 cells (Figure S10O). These findings confirmed the
   role of EGF-EGFR signaling in promoting glioma tumor cell growth. In
   conclusion, the distribution of TAMs in tumor border regions of glioma
   may promote tumor cell growth through the EGF-EGFR signaling pathway.

   In addition, we also investigated the influence of the tumor
   microenvironment on glioma tumorigenesis in both the IDH1-mutant and
   IDH1/EGFR-double-WT samples. In the IDH1-mutant sample, the application
   of the Bulk2Space algorithm revealed a significant enrichment of
   neuronal cells in the non-malignant regions (Figure S11 A and S11B). To
   explore the impact of neuronal cell enrichment on the tumor
   microenvironment, we analyzed the spatial expression of neuron-specific
   proteins and found that the spatial expression of dopamine receptor
   DRD1 showed the most significantly positive correlation with the
   spatial distribution of neurons (Figure S11 C and S11D). Previous
   research on dopamine receptors has shown that their high expression can
   inhibit NLRP3, thereby reducing inflammation [[165]38]. We then
   examined the enrichment of inflammasomes in these regions and found
   that areas with higher neuronal enrichment had relatively lower
   inflammasome enrichment (Figure S11E). Correlation analysis revealed a
   negative correlation between inflammasome enrichment and DRD1
   expression, suggesting that the elevated neuronal enrichment in this
   region might reduce inflammation by modulating NLRP3 (Figure S11 F). To
   further confirm the elevated expression of DRD1 + neuron cells in the
   non-malignant regions and validate the impact of DRD1 on inflammasome,
   we employed multiplexed IHC (mIHC) staining for DRD1 + neuronal cells,
   IDH1^R132H + tumor cells, and NLRP3 + inflammasome activated cells. The
   results confirmed significant enrichment of DRD1 + neuron cells in
   non-malignant regions, accompanied by a reduction in
   inflammasome-positive cells (Figure S11G). These findings indicated
   that the enrichment of dopaminergic-like neurons in the non-malignant
   regions was associated with reduced inflammasome activation.

   In the IDH1/EGFR-double-WT sample, neutrophils were the most abundant
   cell type in the non-malignant regions based on the single-nuclei
   transcriptomic data (Figure S11H). Consistent with this observation,
   PSERP-based spatial proteomic data also revealed elevated expression of
   neutrophil marker MPO in the non-malignant regions of the
   IDH1/EGFR-double-WT sample (Figure S11I). To explore the biological
   functions associated with the enrichment of neutrophils in the
   non-malignant regions, we identified pathways correlated with the
   spatial distribution of neutrophils. The results revealed that the
   neutrophil degranulation pathway had the strongest positive correlation
   with neutrophil enrichment, suggesting substantial activation of this
   pathway in the regions with high neutrophil content (Figure S11 K and
   S11L). Building on this finding, we screened for neutrophil-specific
   molecules and found that the spatial distribution of proteinase 3 (PR3)
   was strongly associated with neutrophil degranulation (Figure S11M and
   S11 N). This is consistent with previous reports indicating that PR3, a
   key enzyme on the neutrophil surface, plays a crucial role in
   modulating chemotactic migration and culminating in a process known as
   NETosis [[166]39]. Furthermore, we examined the spatial enrichment of
   the NETosis pathway and found it to be significantly upregulated in the
   neutrophil-enriched regions (Figure S11O). The significant positive
   correlation between PR3 expression and NETosis further highlighted the
   role of PR3 in neutrophil degranulation and NETosis in the
   non-malignant regions of the IDH1/EGFR-double-WT sample (Figure S11P).

Prospective cohort analysis confirmed the effectiveness of PSERP

   To further demonstrate that our PSERP-based spatial proteomic methods
   could effectively decipher the tumor spatial heterogeneity, we
   collected additional FFPE samples of 20 glioma patients from Zhongshan
   Hospital, Fudan University (Methods). Samples were classified into two
   groups (10 IDH1-mutant, 10 IDH1-WT). No significant inter-group
   differences were observed in age distribution (median age: 50 years
   (IDH1-mutant group) vs. 52 years (IDH1-WT group); Wilcoxon rank-sum
   test, p = 0.75) and gender proportion (60% male in both groups). To
   uncover molecular heterogeneities of the tumor boundary between
   IDH1-mutant and IDH1-WT samples, we preformed both bulk proteomic
   analysis and spatial proteomic analysis on representative areas,
   including the tumor core, tumor boundary, non-malignant region near the
   tumor boundary and non-malignant region distal to the tumor. The
   detailed information of these samples were presented in Table S1.

   To be more specific, IDH1-mutant samples exhibited lower immune cell
   enrichment scores compared with IDH1-WT samples (Figure S12 A). The
   results also showed that NLRP3 and CXCR3 proteins exhibited reduced
   expression in IDH1-mutant samples (Figure S12 C). Conversely,
   IDH1-mutant samples showed higher neuron infiltration scores, with
   elevated expression of dopamine signaling pathway proteins (DRD1, DRD2,
   MAOA, MAOB) (Figure S12B-C). These results were consistent with our
   spatial proteomic observations, indicating that the immune
   microenvironment features identified via PSERP-based spatial proteomics
   may reflect a universal mechanism underlying IDH1 mutation-driven
   alterations in the tumor microenvironment.

   We then systematically characterized the spatial proteomic features of
   IDH1-mutant and IDH1-WT samples. In the malignant regions, IDH1-mutant
   protein was only detected in IDH1-mutant samples (Figure S12D-E). As
   for the non-malignant regions, consistent with the findings from the
   bulk proteomic analysis, the neuronal cell related protein, dopamine
   receptor DRD1, was significantly enriched in non-malignant regions of
   IDH1-mutant samples, especially near the tumor boundaries (Figure
   S12D-E). In contrast, the inflammasome-related protein, NLRP3, was
   predominantly localized in non-malignant regions adjacent to tumor
   boundaries in IDH1-WT samples (Figure S12D-E). Multiplexed IHC staining
   for DRD1 + neuron cells and IDH1^R132H + tumor cells confirmed
   pronounced enrichment of DRD1 + neurons in non-malignant regions
   (Figure S12 F-G). These results validated our hypothesis that the
   expanded population of DRD1 + neurons may establish an
   immunosuppressive niche in peritumoral regions through the suppression
   of NLRP3 inflammasome in IDH1-mutant gliomas. Specifically, dopamine
   signaling via DRD1 receptors inhibits NLRP3-mediated pyroptosis in
   myeloid cells.

   To investigate how neuronal cells at the tumor boundary impact tumor
   cell progression through neuron-tumor interactions, we performed
   further analysis and functional experiments. We analyzed
   neuron-secreted proteins and identified NLGN3 as the most significantly
   upregulated protein in non-malignant regions near tumor boundaries of
   IDH1-mutant gliomas (Figure S13 A-C). NLGN3 plays a important role in
   promoting glioma progression. In our study, correlation analysis
   revealed that NLGN3 expression was strongly associated with mTOR
   signaling activity (Figure S13D). Spatial mapping results further
   demonstrated the co-enrichment of NLGN3 and the mTORC1 complex at tumor
   boundaries (Figure S13E). Collectively, the spatial distribution of
   NLGN3 and its correlation with mTORC1 suggested its potential
   functional role in linking neuronal cells to glioma tumor cells, and
   implied that it might promote tumor cell progression via the mTOR
   signaling pathway.

   To validate this mechanism, IDH1-mutant PDCs were treated with NLGN3
   (100 ng/mL) and/or the mTOR inhibitor rapamycin (20 nM). CCK-8 assays
   showed that NLGN3 significantly promoted PDC proliferation (p < 0.01),
   while rapamycin (m) completely abolished this effect (p = 0.92 vs.
   untreated controls, Figure S13 F-G). These results supported our
   assumption that elevated NLGN3 in the tumor boundary might enhance
   glioma tumor cell growth by activating the mTOR signaling pathway.
   Furthermore, these findings provided actionable insights into
   niche-specific therapeutic targeting for glioma in the future.

   Through the prospective cohort validation and functional experiments,
   we have demonstrated that our PSERP platform reliably resolves
   spatially organized proteomic patterns in gliomas, particularly the
   IDH1 mutation-dependent enrichment of neuronal signaling components
   (e.g., DRD1, NLGN3) at tumor boundaries. The identification of the
   NLGN3-mTORC1 axis as a driver of proliferation in IDH1-mutant gliomas
   not only validates PSERP’s analytical power but also reveals a
   targetable niche-specific mechanism.

   In summary, our PSERP approach provides comprehensive spatial proteomic
   coverage across the entire slide, enabling integrative analysis with
   single-cell transcriptomic data to elucidate the spatial distribution
   and interactions of distinct cell types. Using the PSERP approach, we
   elucidated both the cell type heterogeneity in gliomas with different
   mutational statuses and the role of spatially specific recruitment of
   cell types in promoting tumor cell communication and fostering a
   supportive microenvironment for tumor growth. Moreover, we observed
   several receptor-ligand interactions that primarily occurred at the
   protein level, highlighting the critical role of proteomic analysis in
   understanding cellular communication. Our PSERP-based spatial proteomic
   approach directly identified these interactions, enabling a more
   comprehensive understanding of biological pathways within the tumor
   microenvironment. By pinpointing the high expression of DRD1 and its
   correlation with inflammasome enrichment in the IDH1-mutant sample, as
   well as the observation of NETosis in the IDH1/EGFR-double-WT sample,
   we provided clear examples of how spatial proteomic data uncovers
   intricate molecular interactions within tumor microenvironments, which
   may not be fully captured by spatial transcriptomic data alone.
   Therefore, while spatial transcriptomics is a powerful tool for gene
   expression analysis, the PSERP-based spatial proteomic method provides
   a unique perspective on the functional molecular landscape of tissues,
   which is particularly valuable in understanding complex diseases with
   spatial heterogeneity.

Spatial-resolved proteomics revealed the intra-tumor heterogeneity of
tumor-specific mutational peptides

   Gene mutations can generate neoantigen candidates for personalized
   immunotherapy [[167]12, [168]13]. By integrating tumor-specific genomic
   events, our PSERP-based spatial proteomic approach can portray spatial
   distributions of tumor specific peptides at proteomic level.
   Specifically, in our study, we integrated tumor-specific genomic events
   from our previously published glioma multi-omics dataset, which
   included the glioma specimens utilized in this dataset [[169]17]. The
   tumor-specific genomic events for each tumor were translated into 1, 3,
   or 6 frames and used to build a glioma-specific spectral library. Based
   on this library, we searched for tumor-specific peptides using the raw
   LC–MS files of spatial proteomics. A total of 479 tumor-specific
   mutational peptides were identified (FDR < 0.1 at peptide level and
   protein level). The MHC-I binding affinities for all the detected
   tumor-specific mutational peptides were predicted using NetMHCpan v4.0,
   and 75 tumor-specific mutational peptides with a predicted binding
   affinity < 500 nM were identified (Fig. [170]7A; Methods).

Fig. 7.

   [171]Fig. 7
   [172]Open in a new tab

   The heterogeneity of mutated peptide identified by PRSEP. A The
   strategy of mutated peptide identification and quantification. B Mutant
   peptides whose distribution showed significant spatial-specificity in
   specific cluster of each sample. C, D Heatmaps (upper) showing spatial
   distribution of specific mutant peptides in the IDH1-WT/EGFR-mutant G1
   sample (C) and IDH1-mutant G2 sample (D). Boxplots (bottom) showing
   differential expression levels among clusters in the
   IDH1-WT/EGFR-mutant G1 sample (C) and IDH1-mutant G2 sample (D),
   respectively. E, F IFNy ELISPOT output of PBMCs co-cultured with
   synthesized tumor-specific peptides identified in the
   IDH1-WT/EGFR-mutant G1 sample (E) and IDH1-mutant G2 sample (F).
   Scrambled peptides (control), tumor-specific peptides, or combined
   peptides were performed. G, H Line graphs showing progression rates of
   PBMCs co-cultured with synthesized tumor-specific peptides identified
   in the IDH1-WT/EGFR-mutant G1 sample (G) and IDH1-mutant G2 sample (H)
   by CCK-8 assay. Scrambled peptides (control), tumor-specific peptides,
   or combined peptides were performed. H Images of PDX tumors from mice
   treated with specific mutant peptide or control

   We then analyzed the spatial distribution of these 75 tumor-specific
   peptides, and focused on identifying those whose spatial distribution
   covered the maximal area of malignant regions in tumor samples. To this
   end, five tumor-specific peptides (CTNNA2_p. V401 A, EGFR_p. G598
   VPTEN_p. K125E, SLC25 A5_p. T254M, and CD36_p. G326 V) were screened
   out for the IDH1-WT/EGFR-mutant sample, and five tumor-specific
   peptides (TP53_p. R248 W, IDH1_p. R132H, PML_p. P410L, EIF3 A_p. R700
   K, and HDAC6_p. G518 V) were screened out for the IDH1-mutant sample
   (Fig. [173]7B-D). The combination of these selected peptides covered
   more than 95% of the malignant regions in their respective samples.

   T cells from healthy donors provide a rich source of T-cell repertoires
   that can specifically recognize neoantigens presented on tumors. To
   evaluate the immunogenicity of tumor-specific mutational peptides, we
   synthesized ten tumor-specific peptides (Methods), and evaluate T-cell
   responses by treating them with these candidate peptides, followed by
   assessment of the responses using an ELISpot assay (Methods). IFN-γ
   secretion was significantly increased when treated with the combination
   of peptides (Fig. [174]7E and F). Furthermore, an in vitro cytotoxicity
   assay (Methods) revealed that vaccines derived from the combination of
   peptides showed enhanced tumor cell growth inhibition efficacy
   (Fig. [175]7G and H).

   Based on these findings, we constructed a PDX mouse model via
   subcutaneous injection (Methods), to assess whether tumor-specific
   mutational peptides can activate T cells in vivo. Histological
   examination, proteomic profiling, and phosphoproteomic analysis
   demonstrated a high degree of biological similarity between the
   original tumors and the PDX tumors (Methods, Supplementary Data 1). A
   single type of tumor-specific peptide-derived vaccine or combination of
   peptide-derived vaccines was then intraperitoneally injected into the
   PDX mouse model (Methods). The results revealed that the combination of
   peptide-derived vaccines significantly delayed xenograft growth
   (Fig. [176]7I). In summary, our data highlighted that the PSERP
   approach could directly reveal the spatial expression profiles of
   tumor-specific mutational peptides (potential neoantigens) and their
   nearby microenvironments. This is important since current T-cell
   reactivity detection methods can only be applied to a limited set of
   neoantigens per tumor. Information on the spatial expression of
   neoantigens and tumor microenvironments can help us efficiently filter
   out low-expression or nonspecifically expressed peptides. Moreover,
   combinational targeting of tumor-specific mutational peptides is a
   trend in neoantigen-based treatments. Therefore, selecting appropriate
   neoantigen combinations and retaining highly confident and actionable
   candidates is crucial [[177]13]. Our approach can help identify
   tumor-specific mutational peptides with spatial expression specificity.
   Thus, we can select combinations of tumor-specific mutational peptides
   whose spatial expression covers the largest areas of tumor regions,
   providing insights into potential neoantigen combinations that could
   most efficiently target and kill tumor cells.

Discussion

   Spatial omics is a rapidly developing research field that aims to
   enhance our understanding of health and disease within the spatial
   context of biological structures, such as organs, tissues, or tumors.
   The use of spatial technologies in disease research aids in
   understanding drug action, drug delivery, and immune system activities,
   while also providing critical location information [[178]1, [179]4].
   Spatial transcriptomic studies have begun to derive diagnostic and
   prognostic value in several oncology contexts, demonstrating that
   spatial omic data can facilitate the pathological assessment of tumor
   tissue at the point of care [[180]40]. To date, spatial transcriptomic
   techniques have made significant progress, and several landmark
   methods, such as Stereo-seq [[181]41], Seq-Scope [[182]42], and ZipSeq
   [[183]43], have been employed to perform high-resolution spatial
   analysis across entire tissue sections. Despite the extraordinary
   achievements of spatial transcriptomics, it is important to note that
   proteins, rather than DNA or RNA, govern most cellular processes and
   execute biological functions. To accurately elucidate tissue
   heterogeneity and the underlying molecular mechanism, spatial
   proteomics is critical.

   To assess the advantages and limitations of various spatial proteomic
   methods, we compared PSERP with widely used techniques, including CODEX
   [[184]44], LCM [[185]45], and ProteomEX [[186]46]. CODEX, a multiplexed
   imaging technique, enables high-resolution analysis of tissue
   architecture by simultaneously detecting multiple proteins. However, it
   requires advanced equipment, is limited by available antibodies and
   markers, and has a restricted scanning range. LCM allows for the
   precise microdissection of specific tissue regions for proteomic
   analysis, but it is labor-intensive and requires significant
   histological expertise. ProteomEx combines tissue expansion with manual
   microdissection, but it cannot provide comprehensive spatial proteomic
   profiling across the entire tissue, potentially overlooking gradual
   molecular changes. Additionally, its manual approach may lead to
   reproducibility issues. In summary, although current spatial proteomic
   techniques have made significant progress, they still face limitations
   in terms of protein coverage and throughput (Supplementary Data 2).
   Under these circumstances, we developed PSERP, a matrix-based spatial
   proteomic technique. PSERP enables rapid, quantitative profiling of
   proteomic spatial variability across whole tissue sections with
   submillimeter resolution. Its key advantage lies in its ability to
   capture protein distributions without bias, offering flexible
   field-of-view selection and comprehensive protein coverage. PSERP also
   offers significant time savings by automating much of the process
   through robotic workstations. With standardized protocols, consistent
   quality control, and quantification metrics, PSERP ensures reproducible
   results across different laboratories and large-scale studies,
   enhancing its reliability for research applications. As a result, PSERP
   can effectively capture dynamic alterations in protein expression
   across entire tissue sections. Moreover, its in-depth protein coverage
   and flexible capturing capabilities facilitate subsequent integrative
   analyses with single-cell transcriptomic data or other spatial omic
   data.

   Gliomas exhibit complex intratumoral and intertumoral heterogeneity,
   making them an ideal model for a proof-of-principle study of PSERP. The
   selection of glioma samples was based on several scientific and
   clinical factors. The WHO 2021 classification integrates genetic
   mutations into glioma diagnosis, prompting us to focus on common
   mutations such as EGFR and IDH1, which are frequently observed in
   gliomas [[187]47–[188]49]. More importantly, EGFR and IDH1 mutations
   have been confirmed to be associated with patient’s prognosis, with
   EGFR mutations leading to poor clinical outcomes and IDH1 mutations
   associated with improved clinical outcomes [[189]17]. In addition to
   mutational frequencies, IDH1 and EGFR mutations have been shown to have
   distinct impacts on immune microenvironments. Specifically, compared to
   IDH1-WT/EGFR-mutant samples, IDH1-mutant samples displayed molecular
   characteristics indicative of an immunologically “desert-like”
   microenvironment [[190]17]. Additionally, patients with IDH1 mutations
   can be classified into those with and without1p/19q co-deletion,
   according to the WHO 2021 guidelines. Our preliminary research on IDH1
   mutations included both samples with and without 1p/19q co-deletion
   [[191]17]. Comparative analysis revealed that IDH1-mutant samples,
   regardless of 1p/19q co-deletion status, exhibited a “desert-like”
   microenvironment in contrast to IDH1-WT samples, in contrast to IDH1-WT
   samples [[192]17]. Based on these findings, for this study, we
   specifically selected IDH1-mutant gliomas without 1p/19q co-deletion
   for further analysis. For this proof-of-principle study, we analyzed
   IDH1-mutant and IDH1-WT (IDH1-WT/EGFR-mutant and IDH1/EGFR-double-WT)
   glioma samples to analyze tumor spatial heterogeneity using PSERP-based
   spatial proteomic data.

   In this study, we used the PSERP approach to profile the spatial
   proteome of gliomas with varying mutational statuses. In addition to
   characterizing the distinct features of malignant and non-malignant
   tumor regions, we also illustrated the dynamic expression patterns of
   proteins from the tumor center to the border and paratumoral regions, a
   feature predicted by computational methods [[193]50]. The diverse
   proteomic patterns across the core tumor, tumor border, and paratumoral
   regions provide valuable clinical insights for identifying potential
   therapeutic targets in gliomas. Notably, the strong immune infiltration
   observed at the tumor border enhances our understanding of the
   mechanisms driving glioma metastasis and may inform the development of
   more effective therapeutic strategies. While previous single-cell
   transcriptome studies have highlighted the cellular heterogeneity of
   gliomas [[194]51], the spatial resolution of cellular interactions
   remains poorly understood. By integrating analysis of single-cell
   transcriptome and spatial proteome data, we not only depicted the
   spatial diversity of cellular distribution in gliomas with different
   genetic mutations but also demonstrated how specific cell types reshape
   the tumor ecosystem through cell–cell communication at spatial
   resolution. For instance, OPC-like glioma cells, characterized by
   elevated PDGFRA expression, were found to recruit endothelial cells and
   promote angiogenesis. This observation aligned with our previous
   findings that gliomas with increased expression of PDGFRA were
   characterized by angiogenesis [[195]52]. Additionally, we applied our
   PSERP-based spatial proteomic analysis tool to profile both proteomic
   and phosphoproteomic data, revealing distinct kinase activation
   patterns across samples. These findings underscored the robustness of
   our spatial proteomic and phosphoproteomic data in capturing the
   molecular features of different glioma regions, highlighting the
   potential of PSERP to elucidate tumor spatial heterogeneity.

   In addition, the brain is a complex organ with distinct regional
   functional characteristics that may influence tumor behavior. To
   explore whether the spatial proteomic patterns can reflect the
   characteristic distribution of tumors in the brain, we reviewed the
   clinical records of patients whose tumor samples were included in this
   study. The three IDH1-WT/EGFR-mutant samples and three IDH1-mutant
   samples in our study were derived from distinct brain regions.
   Specifically, for the IDH1-WT/EGFR-mutant samples (G1, G4 and G5), G1
   sample originated from the temporal lobe, while G4 and G5 samples came
   from the frontal lobe. For the IDH1-mutant samples (G2, G6 and G7), G2
   and G6 samples originated from the frontal lobe, while G7 sample came
   from the parietal lobe. We compared tumor characteristics among samples
   with the same mutation features in different brain regions. For
   IDH1-mutant tumors, those in the frontal lobe showed significant
   activation of the Notch signaling pathway, while tumors in the parietal
   lobe exhibited significant activation of the HIF-1 signaling pathway
   (Figure S14 A). Consistent with this finding, proteins such as NOTCH2,
   NOTCH1, and SLC4 A4 were elevated in the frontal lobe, while HIF-1 A,
   TP53, and PDK1 were elevated in the parietal lobe (Figure S14B). For
   IDH1-WT/EGFR-mutant tumors, those in the frontal lobe exhibited
   significant activation of the TGF signaling pathway, while tumors in
   the temporal lobe showed significant activation of the IGF pathway
   (Figure S14 C). Proteins such as SMAD3, SMAD4, and TGFβ were elevated
   in the frontal lobe, while IGF1R, IGF1, and IGF2R were elevated in the
   temporal lobe (Figure S14D). The differences in molecular
   characteristics were also observed in our previous glioma cohort
   studies (Figure S14E-F). These results indicated that tumors in
   different brain regions may exhibit distinct molecular characteristics.
   Furthermore, our methods enhance the understanding of glioma
   infiltration in the brain, potentially aiding in the precise surgical
   resection of tumors in a clinical context. Future research should focus
   on expanding the range of brain tumor regions studied and increasing
   the sample size for more comprehensive comparisons.

   Importantly, our PSERP includes a spatially resolved workflow for
   tumor-specific peptidome identification. Advanced neoantigen prediction
   algorithms can identify millions of potential mutation-derived
   neoantigens. For example, the TSNAdb database compiles neoantigens
   predicted by NetMHCpan from somatic mutations in The Cancer Genome
   Atlas (TCGA) tumor samples [[196]53]. Meanwhile, the GNIFdb database
   covers glioma-specific neoantigens across various subtypes [[197]54].
   Despite these algorithmic advancements, most predicted neoantigens are
   still difficult to validate at the proteomic level, limiting their
   broader application in immunotherapy. To date, no neoantigen-directed
   therapy has been granted regulatory approval for use in cancer
   patients. Several reasons might explain this. One key reason is that
   early research on neoantigens primarily relied on genomic and
   computational approaches, which could not directly capture their
   expression patterns, tumor specificity, or spatial distributions. As a
   result, many predicted neoantigens are not stably expressed,
   compromising their therapeutic efficacy [[198]55]. Our PSERP workflow
   is a powerful tool with several advantages. First, it allows for the
   quantitative detection of tumor neoantigens. Second, it enables the
   spatial mapping of these neoantigens. This helps identify
   tumor-specific antigens expressed exclusively within the tumor region,
   thus minimizing potential damage to normal tissue during targeted
   therapy. Another reason is the diverse immune evasion mechanisms that
   tumors use to escape immune surveillance, which are closely linked to
   intratumoral heterogeneity [[199]56]. Thus, understanding the potential
   impact of neoantigens on the tumor microenvironment is crucial. Using
   our PSERP methods, we observed diverse spatial expression patterns of
   neoantigens. By combining this with the spatial distribution of immune
   cells, we were able to predict the impact of neoantigens on tumor
   immune environments, providing valuable information for therapeutic
   neoantigen selection. Moreover, this method can help select tumor
   neoantigens that cover the tumor most extensively, improving the
   targeted therapeutic effect. To our knowledge, our PSERP provides the
   first spatial proteomic landscape of neoantigens. Therefore, these
   findings could have a significant impact on the clinical translation
   and therapeutic application of neoantigens.

   Furthermore, our method can be integrated with the Tandem Mass Tag
   (TMT) multiplexing strategy, which enables the labeling of peptides
   extracted from each voxel. This integration has the potential to
   accelerate the detection process by allowing parallel detection of
   multiple samples. For example, when using 16-plex TMT (Tandem Mass Tag)
   for sample labeling, the detection time for 600 samples could be
   significantly shortened to 1–2 days. In future studies, by combining
   PSERP with TMT, we could enhance the throughput and efficiency of
   spatial proteomic profiling, reducing overall detection time without
   compromising data quality or resolution.

   Although we have highlighted many advantages of the PSERP method, there
   are still areas for improvement. While PSERP provides submillimeter
   spatial resolution, it may be less sensitive than methods like CODEX,
   particularly for detecting low-abundance proteins in very small or
   sparse tissue samples. This limitation could impact its ability to
   achieve single-cell or subcellular resolution in certain contexts, thus
   affecting its sensitivity in applications that demand ultra-high
   resolution. Furthermore, although PSERP minimizes the need for
   expensive specialized equipment compared to techniques like CODEX, it
   still requires access to mass spectrometry platforms and robotic
   workstations. These technical requirements may present a challenge for
   laboratories lacking such infrastructure, potentially limiting the
   broader adoption and application of PSERP in diverse research settings.

   In conclusion, we propose the PSERP approach, which integrates tissue
   expansion and matrix-based sample preparation with an optimized
   high-throughput mass spectrometry platform. The proposed framework
   offers an efficient and cost-effective approach for characterizing
   intra-tumoral cellular heterogeneity in both tumor and normal tissues.
   We demonstrated the value of our framework in exploring tumor
   ecosystems, making it possible to characterize the spatial
   heterogeneity present in cancer, which could inform the development of
   more precise and effective treatments for brain tumors and other solid
   tumors.

Methods

Specimens and clinical data

   Frozen primary tumor and germline blood samples from qualified patients
   were collected from Zhongshan Hospital. One 4 μm thick slide from each
   tissue sample was sectioned for hematoxylin and eosin (H&E) staining.
   Histopathologically defined adult glioma tumors were considered for
   analysis. Clinical data were obtained from Zhongshan Hospital and
   reviewed for correctness and completeness. The Research Ethics
   Committees of Zhongshan Hospital, Fudan University, approved this study
   (B2019-200R), and all patients provided written informed consent.
   Tissue expansion in PSERP method was performed following the previously
   published protocol [[200]15]. The 10 μm thick slides were
   deparaffinized with xylen and washed with ethanol for PSERP-based
   proteomic analysis. We used 10 μm thick slides from FFPE sample blocks
   for bulk proteomic analysis. We used 10 μm thick slides from FFPE
   sample blocks. The slides were macro-dissected, deparaffinized with
   xylene, and washed with ethanol. The extracted tissues were then lysed
   in a buffer comprising 0.1 M Tris–HCl (pH 8.0), 0.1 M DTT, and 4% SDS
   at 99 °C for 30 min. The crude extract was then clarified via
   centrifugation at 16,000 xg for 10 min, and the supernatant was loaded
   into a 10 kD Microcon filtration device (Millipore), centrifuged at
   12,000 × g for 20 min, and then washed twice with Urea lysis buffer (8
   M Urea, 100 mM Tris–HCl pH 8.0) and twice with 50 mM NH4HCO3. The
   samples were digested using trypsin at an enzyme to protein mass ratio
   of 1:25 overnight at 37 °C. Finally, the peptides were extracted and
   dried SpeedVac (Eppendorf).

Cohort sample acquisition

   For cohort validation, glioma FFPE samples were obtained from Zhongshan
   Hospital, Fudan University. A total of 20 participants underwent
   surgical resection from January 2018 to December 2024. Clinical data,
   including gender, age, tumor grade, and IDH1 mutational status were
   obtained from Zhongshan Hospital and are summarized in Table S1.

Sample processing for whole exome sequencing (WES) and RNA-seq

   Our study sampled a single site of the primary tumor from surgical
   resections due to the internal requirement to process a minimum of 125
   mg of tumor tissue and 50 mg of adjacent normal tissue. DNA and RNA
   were extracted from normal tumor and blood specimens via a coisolation
   protocol using Qiagen QIAsymphony DNA Mini Kit and QIAsymphony RNA Kit.
   Genomic DNA was also isolated from peripheral blood (3–5 mL) to serve
   as a matched normal reference material. A Qubit™ dsDNA BR Assay Kit was
   used with a Qubit® 2.0 fluorometer to determine the concentration of
   dsDNA in the aqueous solution. Any sample that passed quality control
   and produced enough DNA to be subjected to various genomic assays was
   subjected to genomic characterization. RNA quality was quantified using
   a NanoDrop 8000, and quality was assessed using an Agilent Bioanalyzer.
   A sample that passed RNA quality control and had a minimum RNA
   integrity number (RIN) of 7 was subjected to RNA sequencing. Identity
   matches for germline tissue, normal adjacent tissue, and tumor tissue
   were assayed at the BCR using the Illumina Infinium QC array.

Whole exome sequencing

Library construction

   Library construction was performed as described in [[201]57], with the
   following modifications: initial genomic DNA input into shearing was
   reduced from 3 μg to 20–250 ng in 50 μL of solution. For adapter
   ligation, Illumina paired-end adapters were replaced with palindromic
   forked adapters purchased from Integrated DNA Technologies with unique
   dual-indexed molecular barcode sequences to facilitate downstream
   pooling. Kapa HyperPrep reagents in 96-reaction kit format were used
   for end repair/A-tailing, adapter ligation, and library enrichment PCR.
   In addition, during postenrichment SPRI cleanup, the elution volume was
   reduced to 30 μL to maximize the library concentration, and a vortexing
   step was added to maximize the amount of template eluted.

In-solution hybrid selection

   After library construction, the libraries were pooled into groups of up
   to 96 samples. Hybridization and capture were performed using the
   relevant components of Illumina’s Nextera Exome Kit following the
   manufacturer’s suggested protocol, with the following exceptions.
   First, all libraries within a library construction plate were pooled
   prior to hybridization. Second, the Midi plate from Illumina’s Nextera
   Exome Kit was replaced with a skirted PCR plate to facilitate
   automation. All hybridization and capture steps were automated on an
   Agilent Bravo liquid handling system.

Preparation of libraries for cluster amplification and sequencing

   After postcapture enrichment, library pools were quantified using qPCR
   (automated assay on the Agilent Bravo) using a kit purchased from KAPA
   Biosystems with probes specific to the ends of the adapters. For qPCR
   quantification, the concentrations of the libraries were normalized to
   2 nM.

Cluster amplification and sequencing

   Cluster amplification of DNA libraries was performed according to the
   manufacturer’s protocol (Illumina) using exclusion amplification
   chemistry and flow cells. Flow cells were sequenced utilizing
   sequencing-by-synthesis chemistry. The flow cells were then analyzed
   using RTA v.2.7.3 or later. Each pool of whole exome libraries was
   sequenced on 76 paired cycle runs with two 8 cycle index reads across
   the number of lanes needed to meet coverage for all libraries in the
   pool. Pooled libraries were subjected to HiSeq 4000 paired-end runs to
   achieve a minimum of 150 × target coverage per sample library. The raw
   Illumina sequence data were demultiplexed and converted to fastq files;
   adapter and low-quality sequences were trimmed. The raw reads were
   mapped to the hg38 human reference genome, and the validated BAMs were
   used for downstream analysis and variant calling.

Genomic data analysis

Somatic variant calling

   Somatic variants were called from WES tumor and normal paired BAMs
   using SomaticWrapper v1.5, a pipeline designed for the detection of
   somatic variants from tumor and normal exome data. The pipeline merges
   and filters variant calls from MuTect v1.1.7 [[202]58]. SNV calls were
   obtained from Strelka, Varscan, and Mutect. Indel calls were obtained
   from Stralka2, Varscan, and Pindel. The following filters were applied
   to obtain variant calls of high confidence: normal VAF ≤ 0.02 and tumor
   VAF ≥ 0.05; read depth in tumor ≥ 14 and normal ≥ 8; and indel length
   < 100 bp. All variants must be called by 2 or more callers; all
   variants must be exonic; and all variants must exclude variants in
   dbSNP but not in COSMIC.

Single-nuclei RNA sequencing

Single-nuclei RNA library preparation and sequencing

   Approximately 20–30 mg of cryopulverized powder from GBM specimens was
   resuspended in lysis buffer (10 mM Tris–HCl (pH 7.4), 10 mM NaCl, 3 mM
   MgCl2, and 0.1% NP-40). This suspension was pipetted gently 6–8 times,
   incubated on ice for 30 s, and pipetted again 4–6 times. The lysate
   containing free nuclei was filtered through a 40 μm cell strainer. We
   washed the filter with 1 mL of wash and resuspension buffer (1 × PBS
   + 2% BSA + 0.2 U/μL RNase inhibitor) and combined the mixture with the
   original filtrate. After a 6-min centrifugation at 500 × g and 4 °C,
   the nuclear pellet was resuspended in 500 μL of wash and resuspension
   buffer. After DRAQ5 staining, the nuclei were further purified by
   fluorescence-activated cell sorting (FACS). The FACS-purified nuclei
   were centrifuged again and resuspended in a small volume (approximately
   30 μL). After counting and microscopic inspection of nuclear quality,
   the nuclei were diluted to approximately 1,000 nuclei/μL. Approximately
   20,000 nuclei were subjected to single-nuclei RNA sequencing
   (snRNA-seq) on the 10X Chromium platform. We loaded single nuclei onto
   a Chromium Chip B Single Cell Kit, 48 rxns (10 × Genomics, PN-1000073),
   and processed them through a Chromium Controller to generate GEMs (Gel
   Beads in Emulsion). We then prepared the sequencing libraries with the
   Chromium Single Cell 3’ GEM, Library & Gel Bead Kit v3, 16 rxns (10
   × Genomics, PN-1000075) following the manufacturer’s protocol.
   Sequencing was performed on an Illumina NovaSeq 6000 S4 flow cell. The
   libraries were pooled and sequenced using the XP workflow according to
   the manufacturer’s protocol with a 28 × 8x98 bp sequencing recipe. The
   resulting sequencing files were available as FASTQ files per sample
   after demultiplexing.

snRNA-seq data preprocessing

   For each sample, we obtained the unfiltered feature-barcode matrix per
   sample by passing the demultiplexed FASTQs to the Cell Ranger v3.1.0
   ‘count’ command using default parameters, and a customized pre-mRNA
   GRCh38 genome reference was built to capture both exonic and intronic
   reads. The customized genome reference modified the transcript
   annotation from the 10 × Genomics prebuilt human genome reference 3.0.0
   (GRCh38 and Ensembl 93).

   Seurat v3.1.2 was used for all subsequent analyses. We constructed a
   Seurat object using the unfiltered feature-barcode matrix for each
   sample. A series of quality filters were applied to the data to remove
   those cell barcodes that fell into any one of the following categories
   recommended by Seurat: too few total transcript counts (< 300);
   possible debris with too few genes expressed (< 200) and too few UMIs
   (< 1,000); possibly more than one cell with too many genes expressed (>
   10,000) and too many UMIs (> 10,000); possibly dead cells or a sign of
   cellular stress and apoptosis with too high of a proportion of
   mitochondrial gene expression relative to the total transcript count (>
   10%).

   Each sample was scaled and normalized using Seurat’s ‘SCTransform’
   function to correct for batch effects (with parameters: vars.to.regress
   = c("nCount_RNA","percent.mito"), variable.features.n = 3000). We then
   merged all the samples and repeated the same scaling and normalization
   method. All cells in the merged Seurat object were then clustered using
   the original Louvain algorithm [[203]59]: Fast unfolding of communities
   in large networks, J. Stat. Mech. Theor. Exp., 2008 (2008), p. P10008)
   and the top 30 PCA dimensions via Seurat’s ‘FindNeighbors’ and
   ‘FindClusters’ (with parameters: resolution = 0.5) functions. The
   resulting merged and normalized matrix was used for subsequent
   analysis.

snRNA-seq cell type annotation

   Cell types were assigned to each cluster by manually reviewing the
   expression of marker genes. The marker genes used were referred to
   CellMarker2.0, and previously published glioma single cell researches.

snRNA-seq analysis

   Differentially expressed genes within among different cell types were
   identified by the FindMarkers function, which compares cells belonging
   to one cluster to those belonging to the other clusters. The Wilcoxon
   statistical test was used. A log2 FC > 0.3 and an FDR < 0.05 were used
   to filter DEGs.

Analysis of cell–cell interactions

   To analyze the cellular cross-talk between cell types in different
   regions, CellPhoneDB (version 2.0.1) [[204]60, [205]61], was utilized
   to identify significant ligand-receptor interactions. The interaction
   score refers to the mean of the average expression values for all
   individual ligand-receptor partners in corresponding interacting pairs
   for the for the different cell types.

Expansion microscopy (ExM) for proteomics

Tissue expansion

   Tissue expansion for proteomics was performed following a previously
   published protocol [[206]52, [207]62]. Here, we further explain the
   specific steps in the protocol. Firstly, we used paraformaldehyde-fixed
   tissues, which preserve proteins and tissue structure to prevent
   diffusion. Secondly, we applied the succinimidyl ester of 6 ((acryloyl)
   amino) hexanoic acid (acryloyl-X (AcX), 0.1 mg/mL in PBS) to the
   tissue, which was then incubated overnight in a humid chamber at room
   temperature. AcX can form a linker that covalently binds to the primary
   amine groups on proteins. This interaction facilitates the
   incorporation of AcX into the swellable hydrogel polymer during
   gelation. Consequently, the bonds formed by AcX on the proteins are
   tethered to the hydrogel polymer chains, preserving the spatial
   organization of biomolecules relative to one another. Thirdly, the
   tissue was incubated in the expansion gel solution. Initially, tissues
   were incubated in the gel solution on ice for 1 h, which included 0.01%
   (w/v) 4-hydroxy-TEMPO, 2 M NaCl, 8.6% (w/v) sodium acrylate, 30% (v/v)
   acrylamide/bisacrylamide (30% solution; 37.5:1), and 0.2% (w/v) APS in
   1 × PBS. After removing the gel solution, tissues were then incubated
   with an activated gel solution containing the same components plus 0.2%
   (w/v) TEMED. The activated solution was spread into a gelation chamber
   and incubated at 37 °C for 2 h. The gelation chamber was constructed
   using a glass microscope slide and two coverslip spacers (Fisher
   Scientific, catalog no. 12-541B), forming a gel with a thickness of
   approximately 170 μm.

   After gelation, the gel was peeled off from the chamber and homogenized
   in denaturation buffer (200 mM SDS, 200 mM NaCl, 50 mM Tris in Milli-Q
   water) at 75 °C for 4 h. SDS served as a homogenization agent, aiding
   lipid removal while preserving protein integrity by maintaining their
   anchoring to the hydrogel. This step ensured that the sample maintained
   its structural integrity throughout the expansion process. The
   homogenized samples were then washed in on a shaker for three times (30
   min per wash) at room temperature.

   To facilitate the visualization of the samples after expansion, we
   applied Coomassie staining for gel samples. For Coomassie staining,
   gels were washed 3 times with 50% methanol for 15 min per wash and 1
   time with 100% methanol for 1 h at room temperature. Methanol was
   replaced with Coomassie stain solution and incubated for 1.5 h at 95 °C
   then replaced with fresh Coomassie stain solution and incubated for
   another 1.5 h at 95 °C. Although gels may shrink to the original size
   or smaller during Coomassie staining, the mechanical stability of the
   hydrogel in the expanded state was not effected. For gel destaining,
   samples were serially washed in 50 mM Tris, 25 mM Tris, 10 mM Tris and
   2.5 mM Tris (all buffers pH = 8.8) for 1 h each at 60 °C.

   Finally, the samples were washed in Milli-Q water to initiate
   expansion, with the water changed every one hour until the samples
   reached a needed expansion factor, which is a sixfold expansion factor
   in this study. The expanded tissue gel can retain the spatial
   distribution of proteins. Technically, the expansion factor can be
   adjusted according to user requirements, such as threefold, sixfold,
   ninefold and so on. In our proof-of-principle study, we used the
   sixfold expansion factor for all samples, which was corresponded to a
   216-fold (6 × 6 × 6) increase in volume of both the sample and the
   surrounding gel. The original sample thickness was 10 μm, which
   expanded to 60 μm after expansion. The initial gel thickness of
   approximately 170 μm increased to around 1,202 μm after the expansion
   process.

   In our study, we ensure that the gel is homogeneous during the hydrogel
   expansion procedures. Specifically, (1) as for the hydrogel
   composition, we referred to previously published articles [[208]15].
   Our hydrogel composition includes acryloyl-X (AcX), a linker that has
   been reported to covalently binds to the primary amine groups on
   proteins. This interaction facilitates the incorporation of AcX into
   the swellable hydrogel polymer during the gelation phase [[209]15].
   Consequently, the connections formed by AcX on the proteins are
   tethered to the hydrogel polymer chains, effectively preserving the
   spatial organization of biomolecules relative to one another. (2) In
   addition, we included a mechanical homogenization step, using sodium
   dodecyl sulfate (SDS) to maintain the structural integrity and
   organization of the sample throughout the expansion process. SDS serves
   as a homogenization agent that not only aids in lipid removal but also
   preserves protein integrity by affecting their conformation and charge
   without disrupting their anchoring to the hydrogel [[210]15]. As the
   hydrogel-tissue hybrid expands uniformly upon immersion in water, this
   physical magnification of the tissue ensures that the original
   localization of proteins and the overall organizational structure are
   retained. Taken together, the hydrogel composition and the hydrogel
   expansion procedures ensure the homogeneity of the gel and the tissue.

Details on tissues before and after expansion

   In this proof-of-principle study, we applied a sixfold expansion factor
   consistently across all samples, resulting in an isotropic increase in
   length, width, and height. This corresponded to a 216-fold (6 × 6 × 6)
   increase in volume of both the sample and the gel. The original
   thickness of samples used in our study is 10 μm, and the sample
   thickness after expansion is 60 μm. For gel thickness, the gels started
   with a thickness of approximately 170 μm, which expanded to around
   1,202 μm after the expansion process. Moreover, the images of samples
   before and after the expansion were presented to show the tissue sizes
   of the pre-expanded and post-expanded samples (Figure S2 A).

Sample segmentation

   The expanded tissue gel retained the spatial distribution of proteins.
   Subsequently, the expanded tissue was physically segmented into
   numerous volumetric units utilizing a custom-designed 3D-printed
   cutting toolkit. This toolkit facilitated the precise division of the
   expanded tissue into thousands of cubic voxels, with an adjustable
   cutting pitch. The utilization of the 3D-printed cutting toolkit
   ensured uniform and stable sample segmentation. By precisely
   controlling the cutting parameters, the toolkit enabled the division of
   tissue into voxels of stable dimensions. Here, we divided sixfold
   expended tissues into 2 mm × 2 mm voxels (~ 333 μm × 333 μm original
   regions).

The details of the cutting toolkit

   In our study, the cutting toolkit is used for cutting the expanded gel
   into voxels. To better explain the details of the cutting toolkit, we
   describe the preprocessing steps for the samples before cutting
   briefly. Firstly, the sample was embedded in the hydrogel, followed by
   homogenization of the tissue-hydrogel composite to ensure uniformity
   before expansion. Secondly, the tissue-hydrogel composite was expanded
   in water. Thirdly, the expanded tissue-hydrogel composite was placed in
   the cutting toolkit. Subsequently, the expanded tissue-hydrogel
   composite was cut into voxels of the specific size using the cutting
   toolkit.

   The cutting toolkit includes a specialized 3D-printed platform, which
   can support the gels and facilitate the cutting process (Figure S1
   A-B). The 3D-printed platform is a grid base with evenly spaced grooves
   that divide the surface into multiple small squares. The users can
   place the gel onto the platform, and then the gel can be cut into
   voxels along these grid lines. The distance between the centers of two
   adjacent grooves is designed to match the specific voxel size (e.g.,
   2 mm in this proof-of-principle study). We can accomplish a precise
   division of the expanded tissue into hundreds or thousands of uniform
   cubic voxels during the cutting process. The specific size of the
   voxels can be adjusted by using cutting toolkit with different sizes.
   This cutting toolkit enables precise cutting of the gel into uniform
   voxels based on the defined sizes.

   After placing the gel onto the cutting platform, we utilized a robotic
   arm (Tecan Freedom EVO) equipped with a fine, tensioned stainless steel
   cutting wire (30 µm diameter) to uniformly divide the gel into
   specified voxels along the grooves of the cutting toolkit (Figure S1B).
   The cutting size can be customized by the researchers by adjusting the
   cutting toolkit and the robotic arm, allowing for the generation of
   voxels with varying sizes.

   The design of the cutting toolkit aimed to provide flexibility in
   accommodating both the original size of the target samples and the
   specific requirements for the resolution by researchers. To this end,
   we have utilized this customizable and replaceable 3D-printed cutting
   toolkit that was able to align with the cutting procedure and the gel
   samples, according to the specific experimental needs.

The method for recording the position of voxels within the sample

   The expanded gel sample was first imaged using Biorad’s imaging system
   (ChemiDoc, Biorad) (Figure S15 A). The length and width of the image
   were recorded by image J (Biorad). The recorded image was then
   segmented according to the custom-determined voxel size. Each segmented
   voxel will be recorded with its corresponding coordinates (x,y), and
   marked with a barcode (Figure S15 A). The expanded gel sample was then
   segmented by the robotic arm according to the marked image. The
   segmented voxels were then transferred to a 96-well plate in the order
   of barcodes, and processed with further peptide extraction and LC–MS/MS
   analysis. The barcode of each voxel will then be utilized to locate the
   proteomic data of the voxel back to its location on the expanded gel
   sample as well as the pre-expanded sample.

The process of transferring the cut tissue samples

   After cutting, the voxels were transferred using a robotic arm. The
   robotic arm was equipped with a vacuum suction tip that allowed for the
   gentle attachment and movement of the hydrogel. Since the tissue pieces
   remained embedded within the hydrogel, the transfer of gel pieces also
   transferred the cut tissue pieces. This suction mechanism ensured that
   the hydrogels were moved without compromising their integrity,
   minimizing any potential damage during the transfer process. Then, the
   tissue pieces were transferred to the multi-well plates, such as the
   96-well plate, for further processing. The robotic arm can execute
   continuous transfers of the hydrogel in a predetermined sequence,
   enabling fast and high-throughput operations. During the transfer, the
   (x,y) coordinates of the voxels were assigned to the specific wells on
   the 96-well plate, ensuring that after transfer and subsequent
   analysis, the data could be mapped back to the original locations. The
   robotic arm was capable of executing continuous transfers in a
   predetermined sequence, enabling fast and high-throughput operations
   (Figure S1 C). We accurately recorded the positions of voxels within
   the sample to ensure that the spatial information of the original
   tissue was faithfully represented during analysis. After segmenting the
   expanded tissue image and the corresponding H&E image, each voxel was
   assigned a (x, y) coordinate. These coordinates were then mapped to the
   corresponding locations in the corresponding H&E image (Figure S15B-D).

Proteomic analysis

Protein extraction and tryptic digestion

   The excised tissue-gel samples were washed in ddH[2]O three times at 25
   °C and incubated with 50%/50% (v/v) acetonitrile (ACN)/ddH[2]O for 30
   min at 30 °C on a shaker. The samples were then washed in 50%/50% (v/v)
   ACN/100 mM ammonium bicarbonate (ABB) for 10 min and dried in a
   SpeedVac at 45 °C for 3 min. The dried samples were treated with 20 mM
   TCEP for 30 min in darkness at 32 °C, followed by alkylation by the
   addition of 55 mM iodoacetamide (IAA) and incubation for 30 min in the
   dark at 25 °C. The sample pieces were further washed with 100 mM ABB.
   The samples were dehydrated by washing with 50%/50% (v/v) ACN/100 mM
   ABB 2 times for 5 min each at 25 °C on a shaker and then dried in a
   SpeedVac. Protein digestion was performed with 10 ng/μL trypsin in 25
   mM ABB (pH = 8.0) by incubating twice at 37 °C for 4 h and 8 h. The
   digested peptide solutions were collected in the following steps and
   combined: 1) 30–40 μL of the supernatant was collected; 2) 100 μL of 25
   mM ABB was added and vortexed for 10 min at 25 °C, and the supernatant
   was collected; 3) 100 μL of 50% ACN/2.5% formic acid was added, and the
   mixture was vortexed for 10 min; the supernatant was collected, and the
   process was repeated three times; and 4) 100 μL of 100% ACN was added,
   and the mixture was vortexed until the gel pieces turned white and
   sticky. Peptide samples were placed under vacuum to reduce the volume
   to 20–30 μL. The peptides were then desalted using C18 spin columns
   (Pierce™ C18 Spin Tips, Thermo Fisher Scientific, US) and dried in a
   SpeedVac. The cleaned peptides were stored at –20 °C until further
   analysis.

LC‐MS/MS acquisition

   For FFPE derived bulk proteomic analysis, peptide samples were analyzed
   on an Easy-nLC 1200 liquid chromatography system (Thermo Fisher
   Scientific) coupled to a Q Exactive HFX via a nano-electrospray ion
   source (Thermo Fisher Scientific).

   The dried peptides were redissolved in 10 μL loading buffer (5%
   methanol and 0.2% formic acid), and 5 μL of the sample was loaded onto
   a trap column (100 μm × 2 cm, home-made; particle size, 3 μm; pore
   size, 120 Å; SunChrom) with a maximum pressure of 280 bar using solvent
   A, then separated on home-made 150 μm × 12 cm silica microcolumn
   (particle size, 1.9 μm; pore size, 120 Å; SunChrom) with a gradient of
   5–35% mobile phase B (acetonitrile and 0.1% formic acid) at a flow rate
   of 600 nL/min for 75 min. MS analysis was conducted with one full scan
   (300–1,400 m/z, R = 120,000 at 200 m/z) at an automatic gain control
   (AGC) target of 3e6 ions, followed by up to 20 data-dependent MS/MS
   scans with higher-energy collision dissociation (target 5e4 ions, max
   injection time 20 ms, isolation window 1.6 m/z, normalized collision
   energy of 27%). Detection was performed using Orbitrap (R = 7,500 at
   200 m/z) and data were acquired using Xcalibur software (Thermo Fischer
   Scientific).

   For ExM derived proteomics, samples were measured using LC‐MS
   instrumentation consisting of an EASY‐nLC 1200 ultra‐high‐pressure
   system (Thermo Fisher Scientific) coupled via a nano‐electrospray ion
   source (Thermo Fisher Scientific) to a Fusion Lumos Orbitrap (Thermo
   Fisher Scientific). The high-performance liquid chromatography (HPLC)
   gradient was set to 15 min, which enabled efficient peptide separation
   tailored for our samples. For each voxel analyzed, one run was
   performed. The total number of runs depended on the specific paraments
   during sampling, such as the original size of the tissue section, the
   cutting size, and the expansion factor. In this study, for example, we
   collected around 600 voxels on average per sample, which resulted in a
   total of 600 runs on average for each sample. Given that each
   instrument in our platform can analyze approximately 40 samples per
   day, the estimated total measurement time for each sample, involving
   around 600 runs, was approximately 15 days in one instrument. It is
   important to note that the LC–MS/MS analysis time can be adjusted
   according to the desired protein coverage. In our study, we utilized a
   15-min gradient to optimize the number of identified proteins within a
   relatively short time. If the objective is to increase the number of
   identified proteins, the gradient time can be extended. Conversely, for
   studies focusing on a larger number of samples with less stringent
   protein identification requirements, the gradient duration can be
   shortened to enhance throughput.

   MS analysis was performed in a data‐independent manner (DIA). The DIA
   method consisted of MS1 scans from 300–1,400 m/z at 60 k resolution
   (AGC target 4e5 or 50 ms). Then, 30 DIA segments were acquired at 15 k
   resolution with an AGC target of 5e4 or 22 ms for maximal injection
   time. The setting “inject ions for all available parallelizable time”
   was enabled. HCD fragmentation was set to a normalized collision energy
   of 30%. The spectra were recorded in profile mode. The default charge
   state for MS2 was set to 3. For quality control (QC) of the performance
   of the LC‒MS/MS proteomics system, standard HEK293 T QC samples were
   inserted into the sample test sequence at intervals of every 20 runs.

Peptide identification and protein quantification

   For FFPE derived bulk proteomic analysis, MS raw files generated by
   LC–MS/MS were processed with.

   MaxQuant, version 2.3.1.0 with a false discovery rate set at 0.01
   [[211]63]. Proteins were identified with at least 2 peptide of a
   minimum length of 6 amino acids by searching against the human Uniprot
   database (UniprotKB, downloaded in 2022). Protease was Trypsin/P. The
   maximum number of missed cleavages was set to two. A mass tolerance of
   ± 10 ppm for precursor was allowed. The fixed modification was
   Carbamidomethyl (C), and the variable modification was oxidation (M).
   LFQ quantification [[212]64] was used for protein quantification.

   For ExM derived proteomic analysis, the raw data were processed using
   DIA-NN 1.7.10 in the “robust LC (high precision)” mode with
   RT-dependent median-based cross-run normalization enabled Ref: DIA-NN:
   neural networks and interference correction enable deep proteome
   coverage in high throughput). The MS2 and MS1 mass accuracies were set
   to 20 and 12 ppm, respectively, and the scan window size was set to 6.
   Although DIA-NN can automatically optimize such parameters, we fixed
   them to these values to ensure comparability. To generate a
   project-specific spectral library, a 24-fraction high-pH reversed-phase
   fractionated precursor library was created from the same tissue
   specimen and other glioma tissue specimens from a previously published
   study and acquired in DDA mode. The raw files were analyzed with
   MSFragger [[213]65]: MSFragger: ultrafast and comprehensive peptide
   identification in mass spectrometry–based proteomics) under default
   settings (with the exception that cysteine carbamidomethylation was
   removed from fixed modifications) to generate the library file used in
   DIA-NN. The human UniProt (UniProt Consortium, 2022) isoform sequence
   database (3 AUP000005640) was used to annotate the library. The library
   was first automatically refined based on the dataset in question at a
   0.01 global q-value (using the “Generate spectral library” option in
   DIA-NN). DIA-NN performs such refinement by identifying the highest
   score for each library precursor across all runs in the experiment and
   then replacing the library data with the empirically observed spectrum
   and retention time. The resulting report was stringently filtered at a
   0.01 precursor-level q-value, 0.005 precursor-level library q-value and
   0.05 protein group-level q-value. Protein quantification was performed
   using the MaxLFQ algorithm [[214]66]. In our study, 11,700 proteins
   were identified in G1 sample (EGFR-mutant sample), 11,329 in G2 sample
   (IDH1-mutant sample), and 11,465 in G3 sample (IDH1 and EGFR-wild type
   sample).

Patient specific database construction and variant peptide identification

   The proteogenomic database tool CustomProDB [[215]67] was used to
   incorporate the somatic SNVs, indels from our previous published 213
   patients’ glioma cohort [[216]68], into a searchable protein specific
   database with default parameters for variant peptide identification.
   The identified MS2 spectra were removed from the original raw files
   were searched against customized databases using DIA-NN. The database
   searching parameters were almost identical to those described above
   except that Oxidized methionine and protein N-term acetylation were set
   as variable modifications and the data were filtered at 1% PSM FDR.

HLA typing and neoantigen prediction

   Neoantigen prediction was performed with NetMHCpan v4.0 (NetMHCpan-4.0:
   improved peptide-MHC class I interaction predictions integrating eluted
   ligand and peptide binding affinity data), testing all novel 9-mer
   mutant peptides containing the mutated amino acids for binding to the
   patient’s germline HLA alleles. A peptide was defined as a putative
   neoantigen with a predicted binding affinity < 500 nM.

The bulk multi-omics data utilized in the PSERP study

   We utilized the multi-omics data of 213 glioma patients from our
   previously published paper [[217]17] and built a glioma-specific
   spectral library. The detailed methods for library construction were as
   follows: We integrated tumor-specific genomic events (SNVs and InDels),
   from our previously published glioma multi-omics dataset, which also
   included the multi-omics data of the same glioma samples utilized for
   the PSERP-based spatial proteomic analysis in this study. The
   tumor-specific genomic events for each tumor were translated into 1, 3,
   or 6 frames and were constructed to the glioma-specific spectral
   library by CustomProDB [[218]67].

   The raw files generated by PSERP-based spatial proteomic analysis were
   then searched against this glioma-specific spectral library. In this
   way, we identified 479 tumor-specific mutational peptides in total.
   These results proved that by including multi-omics data, PSERP could
   not only reveal the spatial heterogeneity of proteins but also could
   show the impact of mutations on proteins at the spatial level.

   Based on phosphoproteomic data from our previously published paper
   [[219]17], we built a phospho-spectral library, utilizing the FragPipe
   computational platform [[220]44]. The raw DIA data of the three glioma
   samples (IDH1-WT/EGFR-mutant G1, IDH1-mutant G2, and
   IDH1/EGFR-double-WT G3) generated in our research by PSERP were then
   searched against the phospho-spectral library.

The criteria for selecting voxel sizes

   Based on the previous explanation, the voxel size used in our method
   can be adjusted according to the custom requirement. The selection of
   voxel size is primarily based on two factors: proteomic identification
   depth and detection throughput. 1. Regarding identification depth:
   during the experimental design, to evaluate proteomic identification
   depth at different voxel sizes, we performed expansion and sampling of
   a homogeneous glioma tumor sample. We selected samples with
   post-expansion sizes of 1 mm × 1 mm, 2 mm × 2 mm, and 5 mm × 5 mm
   (corresponding to actual sizes of approximately 166 μm, 333 μm, and 833
   μm before expansion) for peptide extraction and mass spectrometry
   analysis. Then, we compared the number of identified peptides and
   proteins in voxels with different sizes. The results showed that on
   average, 946 (ranging from 759 to 1110), 2350 (ranging from 2022 to
   2671), and 3079 (ranging from 2882 to 3311) proteins were identified at
   voxel sizes of 1 mm × 1 mm, 2 mm × 2 mm, and 5 mm × 5 mm (corresponding
   to 166 μm × 166 μm, 333 μm × 333 μm, and 833 μm × 833 μm before
   expansion), respectively. These results suggested that the number of
   identified proteins increased with tissue volume and plateaued at a
   voxel size of 2 mm × 2 mm after expansion (corresponding to 330 μm
   × 330 μm before expansion) (Figure S15E). Using larger voxel sizes may
   lead to a reduction in spatial resolution without a significant
   increase in the number of identified proteins. 2. Regarding detection
   throughput: in our method, the voxel size and the target sample area
   size together determine the number of voxels that need to be analyzed.
   In our study, for a target sample area of approximately 5 cm × 5 cm
   after expansion, selecting a 2 mm voxel size resulted in approximately
   600 voxels. Further reducing the voxel size would significantly
   increase the number of voxels, thereby greatly increasing the time and
   equipment costs for detection. Therefore, voxel size selection also
   depended on the actual detection throughput and costs. Here, we
   presented the corresponding results for different voxel sizes from the
   sample after the six-fold expansion to provide researchers with a
   reference for the suitable selection (Table S1).

   Taken together, based on the considerations of both identification
   depth and detection throughput, we selected a voxel size of 2 mm × 2 mm
   (corresponding to the actual size of 333 × 333 μm) for our
   proof-of-principle study. For each voxel, we also evaluated its peptide
   concentration using BCA analysis. The result revealed that average 104
   ng peptides could be extracted from each 2 mm × 2 mm voxel,
   corresponding to the actual size of 333 × 333 μm.

   Moreover, both the specific voxel size and the expansion factor can be
   adjusted according to the user’s specific needs for further studies.
   The flexible combinations of the voxel size and the expansion factor
   can result in different actual sizes of the voxel and different number
   of voxels (Table S1).

The deal of the boarder of tissue

   To deal with the boarder of the tissue where no full squared voxels can
   be obtained, we used the two quality control procedures for border
   voxels based on two methods used in spatial transcriptomic studies
   [[221]60, [222]69–[223]72]. 1.The alignment with H&E image: in multiple
   spatial transcriptomic studies [[224]60, [225]69, [226]71], spots were
   automatically aligned to the paired H&E images by software, and all
   spots under tissue detected by the software were included in downstream
   analysis. According to this method, we captured the H&E image of each
   sample, and determine the tissue area according to the image (Figure
   S16 A). We recorded the coordinates of all voxels containing tissue as
   mentioned above. After that, all the voxels containing tissue,
   including the border voxels which were not fully populated with tissue,
   were processed for peptide extraction, and mass spectrometry analysis.
   2. The cutoff of identification number for downstream analysis: several
   spatial transcriptomic studies also added a cutoff to filter out spots
   with low identification number. For instance, in Kuppe C et al.’s study
   [[227]71], spots with less than 300 measured genes and less than 500
   UMIs were filtered out. In Kaufmann M et al.’s study [[228]70] spots
   transcriptomes consisting of < 300 expressed genes were removed. And in
   Qi et al.’s study [[229]72], spots were filtered for minimum detected
   gene count of 200 genes. In our study, voxels with fewer than 500
   identified proteins were filtered out from the following analysis.

   Additionally, we compared the number of identified proteins among
   normal voxels which were fully populated with tissue, border voxels
   which were not fully populated with tissue, and blank gel voxels
   (Figure S16B). The result indicated that border voxels which were not
   fully populated with tissue had a relatively lower number of identified
   proteins compared to the fully populated voxels (Figure S16 C).
   However, the number of identified proteins in these boundary voxels was
   significantly higher than that of the control blank gel voxel located
   outside the tissue border that contained no tissue at all (Figure S16
   C). Moreover, during our unsupervised clustering analysis, the
   malignant region boundary voxels with relatively low protein
   identification did not segregate clearly from other malignant voxels;
   instead, they clustered together, indicating that these boundary voxels
   retained meaningful biological information (Figure S16D-E). Similarly,
   non-tumor boundary voxels were also grouped within the non-tumor
   cluster (Figure S16D-E), further suggesting that the boundary voxels
   can contribute valuable insights into the biological context of the
   samples. In summary, to ensure the quality of data, we integrated the
   protein identification results with the image recognition outcomes to
   filter the border voxels, which allowed us to utilize high quality data
   for subsequent analysis and enhance the reliability of our findings.

Bioinformatic analysis

   Proteomic data analysis was performed within the R environment
   ([230]https://www.r-project.org/). Proteins with more than 30% missing
   values (those displayed as 0 in the MaxQuant output) were removed.
   Missing values were subsequently imputed using a normal distribution
   approach with parameters set to a width of 0.3 and a downshift of 1.8
   prior to statistical analysis. To assess differences in protein
   distribution between two groups, we employed the Wilcoxon rank-sum test
   to compare differences between two groups of voxels. The paired t-test
   was used to compare the means of two different groups from the same
   samples. Statistical significance was determined using a Bonferroni
   correction to account for multiple comparisons. In addition, we applied
   a permutation-based false discovery rate (FDR) adjustment with a
   significance threshold of 5% for multiple hypothesis testing.

Functional enrichment analysis of multiomics data using GSVA/ssGSEA

   To further analyze the biological characteristics of the different
   samples, we performed single-sample gene set enrichment analysis
   (ssGSEA/GSVA). Gene expression data of proteome data across different
   samples were used to determine enrichment scores across gene sets
   (i.e., KEGG, Reactome, Gene Ontology and HALLMARK) downloaded from the
   MsigDB (v7.4,
   [231]https://data.broadinstitute.org/gseamsigdb/msigdb/release/7.4/)
   with at least 10 overlapping genes via the R/Bioconductor package GSVA.
   GSEA was performed with GSEA software
   ([232]https://www.gseamsigdb.org/gsea/index.jsp). Gene sets including
   KEGG, Reactome, Gene Ontology and HALLMARK downloaded from the MsigDB
   (v7.4,
   [233]https://data.broadinstitute.org/gsea-msigdb/msigdb/release/7.4/)
   were used as the background.

Estimation of stromal and immune scores

   xCell analysis ([234]https://xcell.ucsf.edu/) [[235]73] was used to
   infer stromal and immune scores based on the RNA-seq data. The
   abundance of 64 different cell types was also computed via xCell.

Application of Cottrazm

   The Cottrazm clustered spots of the G1 sample using the K-nearest
   neighbor (KNN) algorithm based on the normalized spatial protein
   expression matrix. The results revealed the six clusters (Figure S5 C).
   The CNV scores for each spot were calculated using the InferCNV
   algorithm. Two clusters with the greatest CNV scores were identified as
   the core spots of malignant spots (Figure S5D). Spatial spots were then
   arranged, and the Manhattan distance between spots was calculated. The
   neighboring spots of the tumor core where the Manhattan distance was
   less than a linear model-fitted radius were identified. Spots as
   malignant spots or tumor boundary spots were then determined based on
   the UMAP distance to the tumor centroid (Figure S5E and S5 F).

Fuzzy C-means Clustering

   Proteins from spatial proteomic data were grouped into different
   clusters using Mfuzz package in R with fuzzy c-means algorithm
   [[236]74].

Integrated analysis of single-cell and spatial proteomic data by Bulk2Space

   We generated spatially resolved single-cell expression profiles from
   bulk proteomes via Bulk2Space
   ([237]https://github.com/ZJUFanLab/bulk2space) [[238]75] using
   high-quality scRNA-seq data and spatial proteomics as references.