We used multiple datasets of scRNA-seq and spatial transcriptome as well as bulk RNA-seq data of TCGA. Datasets used in this study is summarized in Supplementary Table 1 .

Using the ‘Recount2’ R package (16), we downloaded the RNA sequence data of 32 solid cancers obtained with the Illumina HiSeq 2000 mRNA platform (Illumina, San Diego, CA, USA). The coverage counts were scaled to estimate read counts using the scale_counts function of the recount2 package. For normalization, the read counts of all genes were converted into counts-per-million. We merged the transcriptome data of 32 solid cancers of TCGA projects into one large-scale expression matrix for pan-cancer analysis.

To evaluate overall immune cell enrichment in the TME, cell type enrichment scores were evaluated. A gene-signature based method, the xCell tool (http://xcell.ucsf.edu/), for inferring cell types from tissue transcriptome profiles was used (17). It infers 64 immune and stromal cell types of the TME. The composite score of immune cells (i.e., the ImmuneScore) was obtained for the analysis. More specifically, the immune score was defined as the sum of cell enrichment scores of B cells, CD4+ T-cells, CD8+ T-cells, dendritic cells, eosinophils, macrophages, monocytes, mast cells, neutrophils, and NK cells estimated by xCell. To estimate hypoxia score from RNA-seq data, single sample gene set enrichment analysis (ssGSEA) using hypoxia gene signature of biocarta pathway was applied (18).

The scRNA-seq data were scaled to log-normalization after the read counts were divided by the total number of transcripts and multiplied by 10,000. Two thousand highly variable genes were selected using the FindVariableFeatures function of Seurat (version 3.0) (19) based on a variance stabilizing transformation. The data were then scaled to z-scores. Principal component analysis was run on variable genes, and 10 principal components were selected for clustering analyses. The graph-based clustering approach was implemented using the FindClusters function. A key parameter to determine the number of clusters is the conservative resolution, which was set according to each dataset: 1.0 for the scRNA-seq data of breast cancer; 0.3 for the scRNA-seq data of HNSC; 0.2 for the scRNA-seq data of GBM; 0.5 for the scRNA-seq data of human melanoma treated with ICIs.

The scRNA-seq data were embedded by two-dimensional projection using t-SNE. To identify the marker genes of the clusters, the FindAllMarkers function of Seurat was used, and the first five high-ranked marker genes were identified according to the fold-change. The marker genes that were relatively highly expressed in a given cluster were extracted using the nonparametric Mann–Whitney U test.

Cancer cells and immune cells were identified by known marker genes of specific cell types as well as the extracted markers. For each dataset, at least four cell types were identified using specific well-known markers: cancer cells (EPCAM for epithelial cell origin tumors), T-cells (CD3D, CD8A, and CD4), B-cells (CD79A and IGHM), and myeloid cells (CD68 and LYZ). To identify cancer cells of GBM, we utilized chromosomal copy number alterations as described in previous studies (20, 21). Specifically, expression levels of genes located on chromosome 7 were extracted for all cells, and the scores for copy number alterations were calculated using a module score. The module score was evaluated using the AddModuleScore in the Seurat package. The cluster with a high copy number alteration score was regarded as a cancer cell cluster.

Spatial transcriptome data obtained from the 10X Genomics Visium platform were downloaded from publicly available datasets of 10X Genomics (https://www.10xgenomics.com/resources/datasets/). A human breast cancer dataset with spatial information, and H&E staining was used. The data were preprocessed with ‘SCTransfrom’ function of Seurat version 3. Immune cell enrichment scores were evaluated by xCell as described in the aforementioned section.

As a single gene feature count is sparse, we extracted module scores of GLUT1 and GLUT3 using each gene feature’s correlated genes. The GLUT1 and GLUT3 module scores were calculated using top-k positively correlated genes sorted by the Pearson’s correlation (k = 50). Gene ontology (GO) of GLUT1 (and GLUT3) correlated genes was evaluated using clusterProfiler (22). The module scores for each spot were calculated by AddModuleScore in the Seurat package.

To compare a specific cluster of scRNA-seq data with other clusters, we used the FindMarkers function of the Seurat package. Differentially expressed genes of the cluster were selected according to the following criteria: absolute log fold change > 0.25 and false discovery rate (FDR) corrected p-value < 0.05. GO functional pathway analysis was performed using the clusterProfiler (22). The GO terms were filtered with a cut-off of p-value < 0.05 and FDR of less than 0.2.

To examine the overall activities of glycolysis and OXPHOS of the TME, we used Reactome to select genes of certain metabolic pathways (23). The curated gene sets of the Reactome glycolysis and OXPHOS pathways to define metabolic profiles were obtained from MSigDB (Broad Institute, version 6.0). We used the AddModuleScore function of Seurat to calculate the scores of glucose metabolic pathways.

The correlations between variables were evaluated using Pearson’s correlation analysis. Comparison of the expression levels of gene features or module scores of two different clusters was performed using the Mann–Whitney U test. Further, comparison of the expression levels of multiple clusters was performed using the Kruskal–Wallis test. A p-value of less than 0.05 was considered statistically significant. All statistical analyses were performed using the R software package, version 4.0.2. (http://www.R-project.org).

The expression levels of GLUT1 and GLUT3, which are the main transporters for glucose uptake in cancer cells (24), were compared across different types of cancers using TCGA data. Most cancers showed higher levels of GLUT1 than GLUT3 expression, while testicular germ cell tumors, mesothelioma, sarcoma, diffuse large B-cell lymphoma, thyroid carcinoma, pheochromocytoma/paraganglioma, and liver hepatocellular carcinoma showed higher levels of GLUT3 than GLUT1 expression ( Figure 1A ). According to a previous study on lung cancer, GLUT3 is relatively highly expressed in immune cells within the TME, while GLUT1 is highly expressed in most cancer cells (13). Thus, we suggested that the GLUT3-to-GLUT1 ratio (GLUTratio) could reflect immune cell metabolism compared with cancer cell metabolism. We performed a correlation analysis between the GLUTratio and the composite immune enrichment score (ImmuneScore) across all cancer types. As a result, most cancer types showed a significant positive correlation between the GLUTratio and ImmuneScore ( Figure 1B ). Among the 32 cancer types, 24 presented with a significant moderate and positive correlation. Notably, an insignificant negative correlation was found in diffuse large B-cell lymphoma, which originated from immune cells. This positive association between the GLUTratio and immune profiles was also found between the GLUTratio and cytolytic score, which was calculated through granzyme B and perforin A expression levels to reflect effector T-cell functionality (25) ( Supplementary Figure 1 ). As it is well known that GLUT1 is upregulated in hypoxia, we evaluated the association between GLUTratio and the hypoxia signature. While the hypoxia signature was positively correlated with GLUT1 expression in most cancer types, GLUTratio was not ( Supplementary Figure 2 ).

To confirm whether the different expression patterns of GLUTs represent immune profiles of the TME, we employed multiple scRNA-seq datasets from various human solid tumors. The scRNA-seq data of human head and neck squamous cell carcinoma (HNSC) (26) were clustered, and immune cells of the TME and cancer cells were identified ( Figure 2A ) using marker genes of each cluster ( Supplementary Figure 2 ). Their expression patterns were then visualized using t-distributed stochastic neighborhood embedding (t-SNE) plots. GLUT1 and GLUT3 were differentially enriched according to clusters ( Figures 2B, C ). When clusters were categorized into two cell types (i.e., cancer and immune cells), GLUT3 and GLUT1 expression levels were significantly higher in immune cells and cancer cells, respectively ( Figure 2D ). Similar results were also found for other cancer types. The scRNA-seq data of human glioblastoma multiforme (GBM) were clustered ( Figure 2E ), and their marker genes are presented in Supplementary Figure 3 . Cancer cell clusters were identified using the chromosomal copy number alteration score calculated by the gene expression of chromosome 7 (21) ( Supplementary Figure 3 ). GLUT1 and GLUT3 also showed different patterns according to clusters ( Figures 2F, G ). As demonstrated in HNSC, GLUT1 and GLUT 3 were enriched in cancer cells and immune cells, respectively ( Figure 2H ). This different GLUT expression pattern in cancer and immune cells was also found in breast cancer (27) ( Supplementary Figure 4 ). In addition, the hypoxia score was evaluated for scRNA-seq data of HNSC and GBM ( Supplementary Figure 5 ). It showed that the hypoxia score was not specifically increased in a certain cell type such as cancer cell.

We additionally analyzed a spatial transcriptome dataset obtained from human breast cancer tissue. Spatial gene features were represented using spatially resolved transcriptome data of 3,813 spots ( Figure 3A ). The module scores of GLUT1 and GLUT3 were assessed using each gene feature’s associated genes ( Supplementary Figure 6A ). The GO results of GLUT1-associated and GLUT3-associated genes showed GLUT1 was spatially associated with metabolic process such as NAD metabolic process, while GLUT3 was spatially associated with extracellular matrix and regulation of inflammation response ( Supplementary Figure 6B ). The GLUT3-associated module score was spatially correlated with ImmuneScore, CD8-T-cells, and macrophages ( Figures 3A, B ), while the GLUT1-associated module score was spatially positively correlated with EPCAM expression and negatively correlated with ImmuneScore ( Figures 3A, B ). As the positive correlation between the GLUTratio and ImmuneScore was revealed by the pan-cancer analysis with TCGA data, results of the scRNA-seq and spatial transcriptome data showed that GLUT3 was highly expressed among immune cells in the TME compared with that in cancer cells.

We further analyzed whether glucose metabolic profiles (including GLUTs) were dynamically changed according to ICI treatment since reciprocal glucose uptake between cancer and immune cells reflects tumor immune function (3, 13, 28). A scRNA-seq dataset of human melanoma cells from pre- and post-treatment with ICI was analyzed (29). CD45+ immune cells were clustered into 11 cell types ( Figure 4A ). Consistent with our previous findings, GLUT3 was upregulated in all immune cells compared with GLUT1 ( Supplementary Figure 7A ). To characterize the subsets of immune cell clusters, the marker genes of subsets were depicted in Supplementary Figure 7B . ‘CD8Tcell_1’ and ‘CD8Tcell_2’ were related to exhausted T-cells. ‘CD8Tcell_3’ represented TCF7+ memory/effector cells (29). ‘Myeloid_1’ and ‘Myeloid_2’ types represented macrophages and ‘Myeloid_3’ represented markers related to plasmacytoid dendritic cells ( Supplementary Figure 7B ). Immune cell populations were differently distributed according to the response to ICI and pre- and post-treatment ( Figure 4B ). First, we compared GLUT3 expression levels according to ICI treatment. The GLUT3 level was increased in specific immune cell types after ICI treatment, and the degree of increment was different according to the response to ICI ( Figure 4C ). In particular, GLUT3 of the CD8+ T-cell cluster, which is the key player in ICI (‘CD8Tcell_1’), was increased after ICI treatment in both responders and non-responders ( Figure 4D ). Notably, the change in GLUT3 level after ICI treatment was different between responders and non-responders in myeloid cells, instead of CD8 T-cells ( Figure 4D ). In particular, GLUT3 expression in a myeloid cell cluster, ‘Myeloid_2’, was decreased among non-responders only. Furthermore, changes in GLUT3 expression in myeloid cells after ICI treatment were largely heterogeneous between responders and non-responders.

Considering that effector immune cells depend on glycolysis for energy consumption, we additionally analyzed enrichment scores of glucose metabolism of immune cells. Myeloid cells showed higher glycolysis and oxidative phosphorylation (OXPHOS) scores than other immune cells ( Supplementary Figure 7C ). The glycolysis score of ‘CD8 T-cell_1’ was increased after ICI in both responders and non-responders. The glycolysis score was also remarkably different according to the response to ICI in myeloid cells as the GLUT3 expression level changes. The difference in glycolysis according to the ICI response was also prominent in myeloid clusters. A cluster, ‘Myeloid_1’, showed increased glycolysis and ‘Myeloid_2’ showed decreased glycolysis after ICI treatment in responders, while the glycolysis of both clusters was not changed in non-responders ( Figure 4E ).

As the glucose metabolic profiles of myeloid cells were significantly modified after ICI treatment, we further analyzed myeloid cell subsets for a more detailed examination of metabolic changes. The myeloid subsets were re-clustered ( Figure 5A ) and the landscape of scRNA-seq data, as visualized with t-SNE plots, varied a lot according to ICI treatment ( Figure 5B ). In particular, non-responders showed a remarkable increase in a subtype of the myeloid cluster (‘subtype 2’) after ICI, while this cluster was decreased in responders after ICI ( Figure 5C ). According to the differential gene expression analysis, the marker genes of ‘subtype 2’ were identified as PLTP, MT1G, APOC1, APOE, and MT1H ( Figure 5D ). These differentially expressed genes of ‘subtype 2’ myeloid cells were associated with the lipid catabolic process and mitochondrial functions in GO analysis ( Figure 5E ). According to the KEGG pathway analysis, these genes were associated with ‘phagosome’ as the most significantly associated pathway ( Figure 5E ). GLUTs were differently expressed in the ‘subtype2’ cluster among myeloid cell clusters where GLUT3 was highly expressed in the TME according to the previous results. The ‘subtype 2’ cluster showed relatively low GLUT3 and high GLUT1 expression levels compared with ‘subtype 1’, which is the most abundant myeloid cell subtype in the TME ( Supplementary Figures 7D, E ). The OXPHOS score of ‘subtype 2’ was the highest among the five myeloid clusters ( Supplementary Figures 7D, E ). The ‘subtype 2’ cluster also showed high expression of classical M2 macrophage markers (MRC1 and CD163), and it also expressed some M1 macrophage markers, including CD86, ITGAX, HLA-DRA, and STAT1. Furthermore, it included PD-L1-positive cells ( Supplementary Figure 7F ). Another type, ‘subtype 4’, was relatively increased in responders after ICI. The ‘subtype 4’ was associated with IDO1, CD1C, and NRDG2 expression as well as M1 markers ( Supplementary Figure 7G ). According to the subset analysis of myeloid cells, non-responders showed increased myeloid cells with high OXPHOS, relatively high GLUT1 expression, low GLUT3 expression, and PD-L1-positive types. Thus, this non-responder associated myeloid cell cluster was different from the immune cells in TME which relatively showed high GLUT3 and low GLUT1 according to our previous results.