The TCGA-GBM dataset was obtained from the UCSC Xena website (https://xenabrowser.net/datapages/). The Chinese Glioma Genomic Atlas (CGGA) dataset was acquired from the CGGA website (http://www.cgga.org.cn/).

The TCGA-GBM dataset which have 169 GBM patients and five healthy controls underwent preprocessing using the Affy package in R for normalization and RMA correction. To address the batch effect between the TCGA-GBM dataset and CGGA datasets, the Limma and sva packages were utilized. From the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/cards/GOMF_G_PROTEIN_COUPLED_RECEPTOR_ACTIVITY), a total of 870 GPCR-related genes were obtained, and 18 genes associated with GPCR exhibited correlation or survival status. Differentially expressed genes (DEGs) were identified based on log2 fold change (log2FC) > 1 and adjusted p-value <0.05. Ethical approval for the study was obtained from the local ethical committee at Shanghai Pudong New area People’s Hospital.

GBM patients with effective survival data ranging from 2 days to 3,881 days were included in the study. Univariate and multiple Cox models were employed to evaluate the relationship between GPCR genes and the survival status of GBM patients. The risk score was calculated using the formula: Risk score = Coefficient of mRNA × risk genes. The Least Absolute Shrinkage and Selection Operator (LASSO) method was applied to determine the precise coefficient of each gene based on survival status. The R packages “survival” and “survminer” were utilized to compare the survival status between high- and low-risk groups. The association between survival status, GPCR genes, and the TME was assessed using KM plots, which were further validated with in-house single-cell RNA sequencing data. The TME score, calculated based on infiltrated immune cells using CIBERSORT, determined high TME (indicating a hot tumor with increased infiltrated and functional immune cells) and low TME (indicating a cold tumor with reduced immune cell infiltration). The high GPCR and low GPCR groups were classified based on the LASSO risk score. The exact TME and GPCR scores for each sample can be found in Supplementary Table S2, S3.

Both Gene Set Enrichment Analysis (GSEA) and the fast gene set enrichment analysis (fGSEA) were employed to quantify the score of specific gene sets in each sample. fGSEA was used to calculate gene enrichment in GBM patients with different GPCR and TME status.

The CGGA dataset served as an external validation cohort for the risk model. GPCR scores for each patient were calculated using the formula from the TCGA model, and patients were subsequently categorized into high- and low-GPCR groups.

WGCNA analysis is a computational method for identifying co-expressed genes among different groups and identify marker genes based on the non-orientation analysis between gene sets and phenotypes. Here, WGCNA was utilized to identify gene modules related to different clusters in glioblastoma.

The single-cell data were obtained from our surgical samples: one is glioblastoma and the other is meningioma (serving as controls). The meningioma is a benign brain tumor without aggressive cell proliferation and invasiveness. The ethics was approved form the local ethical committee (2022-K18) and the patients’ approval has been obtained before the surgery for the willing on the surgery and taking part in this research. The consent form was also obtained from the patients and their relatives.

The exact procedure for scRNA-seq for brain tissues could be obtained from the 10X website (https://www.10xgenomics.com/cn/resources/document-library/4a0968). Fresh brain surgical samples from two brain tumor patients were collected and single-cell suspensions were isolated with the assay within 12 h after the surgery. Single-cell 3′libraries were constructed following the 10X Chromium protocol, and each single-cell 3’ library was sequenced using the Illumina Novaseq 6,000. The sequencing depth ranged from 35,000 to 50,000 mean reads per cell. And the estimated number of cells are 10,000–13,000. The median genes per cell is 2000–2,500. The database construction, sequencing, and data analysis were completed by Shanghai Ouyi Biomedical Technology Co., Ltd. Produced during high-throughput sequencing. The raw data (in fastq format), was processed using the CellRanger software from 10x Genomics. According to the data quality statistics and comparison with the reference genome, the software can distinguish cells by Recognition sequence. The high-throughput single cell Transcriptome was quantified by sequence markers and UMI markers of different mRNA molecules in each cell to obtain quality control statistical information such as high-quality cell count, gene median, and sequencing saturation. We processed the unique molecular identifier (UMI) count matrix using the R package Seurat (version 4.1.1). We first normalized the data with sctransform (SCT) in order to account for variance in sequencing depth across data points, detecting high-variance features, and stores the data in the SCT assay. Next, we carried out data quality control. We captured cells with less than 10% mitochondrial genes, and a total number of gene count ranging from 200 to 10,000, expressed in at least 3 cells were selected. Highly variable genes was set at 2000 (Zhang et al., 2021). The two samples (glioblastoma and meningioma) were integrated through single-cell transform (SCT) correction, and the uniform Manifold Approximation and Projection (uMAP) method was employed for data dimension reduction. Cellchat R package was applied to explore the cell-cell interaction between the cell clusters. The strength of cell-cell interaction is based on the total weight of ligand-receptor between different cell clusters, which includes the interaction number and communication probability (Jin et al., 2021). The monocle method was used for pseudo-time analysis to identify the stage and gene changes of T cells. GPCR extent was compared between glioblastoma and meningioma using the addmodulescore method.

Pseudo-time analysis was performed using the R package monocle v2.16.0 to obtain genes required for calculating pseudo-time. Differential gene expression analysis was conducted among T cells from different pathological stages. After calculating pseudo-time, differential gene expression analysis was repeated to determine genes that changed as a function of pseudo-time. The cell state containing the greatest number of S0-stage cells was considered the root state. The threshold of the q-value for multiple testing involved in the selection of differentially expressed genes (DEGs) was set at 0.01.

The maftools package (Mayakonda et al., 2018) was utilized for visualizing somatic mutations and determining the TMB of each patient. TMB status was compared between patients with high and low GPCR scores. Survival status, the expression levels of immune checkpoint molecules, and human leukocyte antigen (HLA) markers were also compared.

All data analyses were performed in R v4.0.3, and a p-value <0.05 was considered statistically significant.

Figure 1A shows a heatmap demonstrating the distribution of top GPCR genes between two groups. Lasso regression analysis was conducted and identified ten GPCR genes related to survival status (Figure 1B). Using the bootstrap method, a GPCR score was constructed with the following coefficients: 1.076*ADRB3+0.180*GPCR68 + 0.468*TPRA1+0.193*GPCR82 + 0.111*HTR7+0.583*CRCP+0.402*GPBAR1+(-0.226)*GLP1R+0.046*F2RL2+0.024*FZD1 (Figure 1C). Multiple cox analysis indicated that only CRCP was an independent risk factor in GBM (Figure 1D). However, we found that the GPCR_TME classifier was also an independent risk factor with a high odds ratio of 1.51, as illustrated in Figures 1E,F. Therefore, the classifier based on high GPCR/low GPCR or high TME/low TME can be used to predict the survival status in GBM, with the worst outcome observed in patients with high GPCR/low TME and a better survival status in those with low GPCR/high TME (p < 0.001). This classifier was also validated in the CGGA dataset, demonstrating that in patients with high GPCR and shorter survival status, the TME did not differentiate the outcome. However, in patients with low GPCR status, the TME could further classify patients based on the survival status (Figures 1G,H).

The co-expression network was identified using WGCNA, resulting in the identification of twelve modules (Figure 2A). The correlation between these twelve modules and the two groups was assessed by measuring the correlation between ME values and clinical features (Figure 2B). The results indicate that the brown and green modules positively correlate with the GPCR_high_TME_low phenotype (r = 0.51, p = 2e-12; r = 0.48, p = 7e-11) and negatively correlated with the GPCR_low_TME_high phenotype (r = -0.38, p = 5e-07; r = -0.27, p = 4e-04, Figure 2B). Genes in the green and brown modules were listed in Supplementary Table S4. Additionally, fGSEA indicate that the GPCR_high_TME_low group exhibit a positive correlation with positive regulation of T cell proliferation, while the GPCR_low_TME_high group show a negative association with T cell mediated immunity (Figure 2C). In the GPCR high and TME high group, there is also an association with the GPCR signaling pathway and a negative correlation with T cell-mediated immunity. Therefore, both GPCR and TME status might be involved in the T cell activity.

The TCGA samples were classified into three clusters: tumor, immune, and stromal parts (Figure 3A). Further analysis using the AddModuleScore method revealed significant enrichment of GPCR genes in immune cells (Figure 3B). This finding was consistent with the violin graph presented in Figure 3C. As GPCR-related genes were primarily expressed in immune cells, we investigated cell-cell interactions at a single-cell level using the CellChat R package. The analysis revealed strong interactions between several cell types (Figure 3D). The bubble plot of the network indicated that T cells with high GPCR levels were connected to tissue stem cells via PTN/NCL, while this connection was absent in T cells with low GPCR levels. Additionally, T cells with high GPCR levels were connected to macrophages via Tnf/Tnfrsf1b, which was not observed in T cells with low GPCR levels. These results suggest that immune cells with high GPCR levels may be involved in immune escape (Figure 3E). Comparison of GPCR levels between the GBM group and the meningioma group, based on in-house sc-RNA seq data, revealed significantly higher GPCR levels in the GBM group (Figure 3F).

Subsequently, T cells were reclustered, leading to the identification of 11 T sub-clusters (Figure 3G). Cluster eight and Cluster 10 exhibited significantly higher GPCR scores (Figure 3H). Pseudo-time analysis demonstrated that T cells had three stages, and there was a transition of T cell clusters from low GPCR to high GPCR levels (Figure 3I).

Considering the differences in GPCR scores between tumor and immune cells, we further compared the tumor mutation burden (TMB) level between different groups. Although not statistically significant, the TMB level was higher in the GPCR_low + TME_high group compared to the GPCR_high + TME_low group (Figure 4A). Consequently, the classifier based on GPCR level and TME level was unable to differentiate the survival status among these groups, although a trend was observed between the GPCR_low + TME_high group and the GPCR_high + TME_low group (Figure 4B). We then analyzed gene mutations between the two groups. In the GPCR_high + TME_low group (Figure 4C), the mutation level was lower compared to the GPCR_low + TME_high group (Figure 4D). The mutation percentage of each gene is also listed in Figures 4C,D. In the GPCR_low + TME_high group, the top five mutated genes were TP53 (39%), PTEN (34%), TTN (26%), EGFR (24%), and NF1 (21%), while in the GPCR_high + TME_low group, the top five mutated genes were EGFR (34%), TP53 (34%), TTN (29%), MUC16 (22%), and PTEN (22%). Moreover, expression of most immune checkpoints and HLA molecules was significantly higher in the GPCR_high + TME_low group, including BTLA, CD209, CD274, CD80, CD86, CTLA4, and PDCD1. Regarding the expression of HLA molecules, their levels were significantly higher in the GPCR_high + TME_low group.