The cancer names and abbreviations involved in this study include Adrenocortical Carcinoma (ACC), Bladder Urothelial Carcinoma (BLCA), Breast Invasive Carcinoma (BRCA), Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma (CESC), Cholangiocarcinoma (CHOL), Colon Adenocarcinoma (COAD), Lymphoid Neoplasm Diffuse Large B-cell Lymphoma (DLBC), Esophageal Carcinoma (ESCA), Glioblastoma Multiforme (GBM), Head and Neck Squamous Cell Carcinoma (HNSC), Kidney Chromophobe (KICH), Kidney Renal Clear Cell Carcinoma (KIRC), Kidney Renal Papillary Cell Carcinoma (KIRP), Acute Myeloid Leukemia (LAML), Brain Lower Grade Glioma (LGG), Liver Hepatocellular Carcinoma (LIHC), Lung Adenocarcinoma (LUAD), Lung Squamous Cell Carcinoma (LUSC), Mesothelioma (MESO), Non-Small Cell Lung Cancer (NSCLC), Ovarian Serous Cystadenocarcinoma (OV), Pancreatic Adenocarcinoma (PAAD), Pheochromocytoma and Paraganglioma (PCPG), Prostate Adenocarcinoma (PRAD), Rectum Adenocarcinoma (READ), Sarcoma (SARC), Skin Cutaneous Melanoma (SKCM), Stomach Adenocarcinoma (STAD), Testicular Germ Cell Tumors (TGCT), Thyroid Carcinoma (THCA), Thymoma (THYM), Uterine Corpus Endometrial Carcinoma (UCEC), Uterine Carcinosarcoma (UCS), and Uveal Melanoma (UVM).

The foundational data for this study were sourced from the TCGA and GTEx databases, comprising 15,901 normal or adjacent tissue samples and 10,201 tumor tissue samples. These include gene expression data for SLC38A2 and clinical data from patients. Additionally, single-cell sequencing data from the National Institutes of Health (NIH) were utilized, covering datasets GSE176078, GSE168652, GSE138709, GSE166555, GSE162631, GSE171306, GSE166635, GSE162025, GSE162708, and GSE141445. Moreover, data from the omicsdi database (EMTAB6149) and the NGDC database (CRA001160) were included. The study also employed immunohistochemistry slides and protein 3D structure data from the Human Protein Atlas (HPA).

For TCGA data, we downloaded and organized RNAseq data from 33 cancer projects processed by the STAR pipeline and extracted the data in TPM format. We discarded data without corresponding clinical information and converted it to log₂(value + 1). For expression matrices downloaded from the GEO database, we used the ComBat_Seq function from the “sva” package to remove batch effects, preserving grouping information using the mod option.

From the TCGA expression database, we used the R packages “ggplot2”, “stats”, and “car” for data visualization, and apply the Wilcoxon rank sum test for statistical analysis.

The scRNA-seq data for single-cell analysis was derived from 12 datasets, and the expression matrix transformation and batch effect correction were completed as described earlier. Subsequently, the “Seurat” package was used for quality control. Specifically, genes expressed in at least three cells and cells with at least 200 genes expressed were extracted for data initialization. During the filtering stage, cells with gene counts between 200 and 7,500 were retained, and cells with mitochondrial gene proportions exceeding 20% were removed (Zhang et al., 2024). The “FindVariableFeatures” function was used to select 3,000 highly variable genes for analysis, the “FindNeighbors” function was used to construct a KNN graph based on the top 30 PCs, and UMAP was employed for visualization. For cell annotation, the “FindAllMarkers” function was used to identify differentially expressed genes for each cluster and perform automatic annotation. Additionally, the “FeaturePlot” function was used to label and visualize the expression of SLC38A2.

For the forest plot of SLC38A2 across various cancers, select OS (Overall Survival), DSS (Disease-Specific Survival), and PFI (Progression-Free Interval) as outcome variables. Use the “survival” package to test the proportional hazards assumption and perform univariate Cox regression analysis.

For multivariate Cox regression analysis, we use Overall Survival (OS) as the outcome variable. The default clinical indicators included are pathologic M stage, pathologic T stage, pathologic N stage, gender, and age. For cancer types without TNM staging, we use the WHO classification standards. The proportional hazards assumption is tested using the survival package, followed by multivariate Cox regression analysis. Variables meeting the set p-value threshold in the univariate analysis are included in the multivariate Cox model construction. Visualize the forest plot using the “ggplot2” package. For survival curves, use the “survminer” and “ggplot2” packages for visualization.

To analyze the correlation between SLC38A2 expression and immune infiltration matrix data across 33 types of cancer in the TCGA database, the ssGSEA algorithm provided by the R package “GSVA” was utilized. This involved calculating immune infiltration using markers of 24 immune cell types. Additionally, the “estimate” package was used to compute stromal and immune scores. Visualization was performed using the “ggplot2” package.

Based on previous cell annotation results, divide the single-cell data in BRCA and PAAD into three major cell groups: malignant cells, immune cells, and stromal cells for independent trajectory analysis. Use the R package “Monocle3” for pseudotemporal analysis, apply “learn graph” to learn cell development trajectories, and use “order cells” to set the root node based on biological markers. Use “graph test” to identify differentially expressed genes along the trajectory.

We obtained 580 samples from the LUAD cohort (58 normal tissues and 522 tumor tissues) and 503 samples from the COAD cohort (41 normal tissues and 462 tumor tissues) for differential analysis based on the expression profiles of 947 metabolism-related genes. Significant differentially expressed genes were identified using the DESeq2 algorithm. In COAD, 518 out of 947 metabolism-related genes (54.7%) were differentially expressed between normal (41 samples) and colorectal cancer (462 samples) tissues, with 253 genes upregulated and 265 genes downregulated in tumors. In LUAD, 244 out of 947 metabolism-related genes (25.7%) were differentially expressed between normal (58 samples) and lung adenocarcinoma (522 samples) tissues, with 160 genes upregulated and 84 genes downregulated in tumors.

Subsequently, WGCNA was performed on the 244 and 518 differentially expressed genes in the two cancers, respectively. A scale-free topology network was constructed based on the soft threshold power, and co-expression modules were identified using the dynamic tree cut algorithm. After calculating module eigengenes, their Pearson correlation with clinicopathological features (T stage, N stage, pathological stage, and overall survival) was assessed. The clusterProfiler was used for KEGG pathway enrichment analysis of key modules, analyzing the blue module in LUAD and the yellow module in COAD.

Based on the tumor samples from both, we first performed molecular subtyping using consensus clustering (ConsensusClusterPlus R package), determining the optimal number of clusters k = 3 through cumulative distribution function (CDF) analysis. Differential expression analysis (DESeq2 R package) was then conducted for the three subtypes (C1/C2/C3), and Kaplan-Meier survival analysis (survival R package) was used to evaluate prognostic differences among subtypes.

We collected the half-maximal inhibitory concentration (IC₅₀) values of 265 small molecules and their corresponding mRNA gene expression data from 860 cell lines in the Genomics of Drug Sensitivity in Cancer (GDSC). Similarly, we gathered IC₅₀ values of 481 small molecules and their corresponding mRNA gene expression data from 1,001 cell lines in the Cancer Therapeutics Response Portal (CTRP). We then integrated the mRNA expression and drug sensitivity data. Subsequently, we performed a Pearson correlation analysis to determine the correlation between SLC38A2 expression and IC₅₀ values. Finally, a volcano plot was created using the “ggplot2” package, with the horizontal axis representing the correlation coefficients and the vertical axis representing the adjusted p-values.

The Wilcoxon test was used to calculate expression differences between groups. Log-rank and COX regression were employed to examine the significance of survival curve differences between high and low SLC38A2 expression groups. Spearman correlation analysis was used to assess the correlation between SLC38A2 expression and the TME (including immune score, stromal score, and immune cell infiltration). All data were log-transformed. All p-values were adjusted for false discovery using the Benjamini–Hochberg procedure. Significant associations in analyses required an FDR-adjusted p-value (q-value) < 0.05, and a p-value of less than 0.05 was considered statistically significant.

The results of the unpaired sample differential analysis indicate that SLC38A2 exhibits differential expressions in 14 types of cancer tissues compared to adjacent non-cancerous tissues. Specifically, SLC38A2 is overexpressed in eight types of cancer which include BLCA (p < 0.05), COAD (p < 0.001), ESCA (p < 0.05), HNSC (p < 0.001), KIRC (p < 0.001), KIRP (p < 0.05), LUSC (p < 0.001), and STAD (p < 0.001). Conversely, it is underexpressed in six types of cancer that include BRCA (p < 0.001), LIHC (p < 0.001), LUAD (p < 0.01), PRAD (p < 0.001), THCA (p < 0.001) and UCEC (p < 0.001) (Figure 1A). In paired sample analysis, SLC38A2 shows differential expressions in 10 types of cancer, including BLCA (p < 0.01), COAD (p < 0.01), HNSC (p < 0.001), KIRC (p < 0.001), KIRP (p < 0.05), LUSC (p < 0.001), BRCA (p < 0.05), LIHC (p < 0.001), THCA (p < 0.001), and UCEC (p < 0.001). The first six cancers correspond to overexpression as observed in the unpaired sample analysis, while the last four correspond to underexpression (Figure 1B).

Immunohistochemistry results indicate that SLC38A2 expression is significantly lower in LIHC and PAAD compared to normal liver and pancreatic tissues. In contrast, its expression is higher in COAD compared to normal colonic epithelium. Additionally, the expression of SLC38A2 in lung, breast, and oral epithelial tissues, as well as their corresponding cancers (LUAD, BRCA and HNSC) were also examined (Figures 1E–J).

The mutation analysis results indicate that the main types of mutations for SLC38A2 in cancer are amplification, mutation, and deep deletion. The top three cancers with the highest probability of these mutations are BLCA, DLBC, and UCEC (Figure 1D).

Using overall survival (OS) as the evaluation metric, SLC38A2 serves as a risk factor for BRCA (HR = 1.597, p < 0.05), LUAD (HR = 1.650, p < 0.01), MESO (HR = 2.007, p < 0.05), and PAAD (HR = 1.761, p < 0.05), while acting as a protective factor for KIRC (HR = 0.625, p < 0.05) (Figure 2A). For disease-specific survival (DSS), SLC38A2 is a risk factor for LUAD (HR = 1.762, p < 0.05) and PAAD (HR = 1.999, p < 0.05), and a protective factor for KIRC (HR = 0.524, p < 0.01) (Figure 2B). In terms of progression-free interval (PFI), SLC38A2 is a risk factor for LUAD (HR = 1.666, p < 0.01) and PAAD (HR = 1.766, p < 0.05) (Figure 2C). Survival curve results indicate that high expression of SLC38A2 is associated with shorter overall survival in BRCA, LUAD, MESO, and PAAD, suggesting a correlation with poor prognosis (Figure 2D). Additionally, the results of the multivariate Cox regression further confirmed that SLC38A2 acts as an independent risk factor in these cancers (Supplementary Material 1).

The analysis results from the Estimate algorithm indicate that the expression of SLC38A2 is positively correlated with stromal scores in most cancer types, including BRCA, COAD, ESCA, KIRC, LIHC, LUAD, OV, PAAD, PCPG, PRAD, READ, STAD, TGCT, THCA, and UCS. In contrast, for immune scores, SLC38A2 exhibits a negative correlation pattern in cancers such as BLCA, BRCA, CESC, HNSC, KIRC, KIRP, LGG, LUAD, LUSC, SARC, SKCM, TGCT, THYM, and UCEC (Figure 3A).

The analysis results from the ssGSEA algorithm indicate that the expression of SLC38A2 is generally negatively correlated with immune cells in most cancers. This includes aDCs, B cells, CD8⁺ T cells, cytotoxic cells, DCs, NK CD56bright cells, NK CD56dim cells, NK cells, pDC, T cells, TFH, Th17 cells, and Tregs. Additionally, SLC38A2 shows a generally positive correlation with certain immune cells, such as T helper cells, Tcm, Tem, Tgd and Th2 cells (Figure 3B). Further analysis was conducted on the correlation between SLC38A2 and infiltrating immune cells in four types of cancer which include BRCA, LUAD, MESO, and PAAD. The results indicate that SLC38A2 shows a moderate positive correlation with infiltrating Tcm (R = 0.435, p < 0.001) and T helper cells (R = 0.361, p < 0.001) in BRCA, while it exhibits a moderate negative correlation with pDCs (R = −0.325, p < 0.001) in BRCA. Additionally, SLC38A2 is moderately negatively correlated with infiltrating Th17 cells (R = −0.349, p < 0.001) in LUAD. In MESO, SLC38A2 demonstrates a moderate positive correlation with infiltrating T helper cells (R = 0.384, p < 0.001). In PAAD, SLC38A2 shows a strong positive correlation with infiltrating T helper cells (R = 0.588, p < 0.001) and a moderate positive correlation with infiltrating Tcm (R = 0.405, p < 0.001), Neutrophils (R = 0.402, p < 0.001), Th2 cells (R = 0.380, p < 0.001), Th1 cells (R = 0.374, p < 0.001), and Macrophages (R = 0.314, p < 0.001) (Figure 3D).

Additionally, we investigated the relationship between high and low expression of SLC38A2 and immune cell infiltration. In BRCA, the subgroup with high SLC38A2 expression corresponds to higher infiltration of T helper cells (p < 0.001) and Tcm (p < 0.001), and lower infiltration of pDCs (p < 0.001). In LUAD, the high SLC38A2 expression subgroup corresponds to higher infiltration of Th2 cells (p < 0.001) and lower infiltration of Th17 cells (p < 0.001) and pDCs (p < 0.001). In MESO, the high SLC38A2 expression subgroup is associated with increased infiltration of T helper cells (p < 0.01), and Th2 cells (p < 0.05). Finally, in PAAD, the subgroup with high SLC38A2 expression corresponds to greater infiltration of T helper cells (p < 0.001), Tcm (p < 0.001), and Neutrophils (p < 0.001) (Figure 3C).

By performing dimensionality reduction and clustering on high-throughput sequencing datasets from 12 types of cancer, we categorized the cells within the samples into the following subpopulations: Alveolar, B cells, Cholangiocytes, CD4⁺ Tconv, CD8⁺ T cells, CD8⁺ Tex, Dendritic Cells (DCs), Endothelial, Epithelial, Endometrial stromal cells, Fibroblasts, Hepatocytes, Monocytes/Macrophages (Mono/macro), Mast cells, Myofibroblasts, Microglia, Mural cells, Natural Killer (NK) cells, Neutrophils, Plasma cells, Pericytes, Smooth Muscle Cells (SMC), Trpolif and so on (Figures 4A,C).

We found that SLC38A2 is expressed across various cell lineages. Notably, there is a certain level of enrichment of SLC38A2 in tumor-associated stromal cells, such as epithelial cells, endothelial cells, and fibroblasts. In contrast, the expression of SLC38A2 is relatively low in immune-related cells (Figures 4D,E). Additionally, we discovered that SLC38A2 expression is associated with tumor immune subtypes. In BRCA, SLC38A2 is expressed at higher levels in tumors with C1 and C4 immune subtypes, while lower levels are observed in C2, C3, and C6 subtypes. In LUAD, SLC38A2 shows higher expression in C1 and C6 immune subtypes, with lower expression in C2, C3, and C4 subtypes. In MESO, higher expression is seen in C1, C3, and C6 immune subtypes, and lower expression in C2 and C4 subtypes. In PAAD, SLC38A2 is more highly expressed in C1, C3, and C6 subtypes, with reduced expression in the C2 subtype (Figure 4B).

In BRCA, the expression of SLC38A2 is positively correlated with GSTT4, GGT1, GGLC, HPGDS, GSR, GSTT2, GSTK1, MGST2, RRM2B, and GSTT2B, while it is negatively correlated with G6PD, GPX1, IDH2, CHAC2, ODC1, PGD1, RRM2, SMS, CHAC1, GSTO1, and PRDX6 (Figure 5A). In LUAD, SLC38A2 shows positive correlations with GGLC, GPX8, CHAC2, RRM2B, LAP3, TXNDC12, RRM2, and SMS, and negative correlations with GGT1, OPLAH, HPGDS, GPX1, GSTK1, and MGST2 (Figure 5B). In MESO, SLC38A2 is positively correlated with G6PD, GSTT4, GCLC, IDH1, GPX8, ODC1, RRM2B, TXNDC12, RRM2, and SMS (Figure 5C). In PAAD, SLC38A2 exhibits positive correlations with G6PD, GCLC, HPGDS, GSR, IDH1, GPX8, CHAC2, ODC1, RRM2B, LAP3, PGD, RRM2, GSTT2B, SMS, and PRDX6 (Figure 5D).

In BRCA and PAAD, we categorized cells into three main types that include malignant cells, stromal cells, and immune cells. We analyzed the expression patterns of SLC38A2 and genes involved in glutathione metabolism within these cell types. The results indicate that SLC38A2 is highly expressed in the later stages across most cell subtypes, including malignant and stromal cells in BRCA, as well as malignant, stromal, and immune cells in PAAD. It is worth mentioning that SLC38A2 shows early enrichment in the immune cell subtype of BRCA (Figure 6).

To further elucidate glutamine metabolism in pan-cancer and analyze its association with SLC38A2, we conducted further validation in COAD and LUAD. First, in COAD, we identified 518 differentially expressed genes out of 947 metabolic genes in the TCGA patient cohort, with 253 genes upregulated and 265 downregulated in tumors (Figure 7A). WGCNA results indicated that these genes could be well clustered into six distinct gene modules (Figure 7B). The yellow gene module was significantly associated with clinical indicators (pathologic_N = −0.15, pathologic_stage = −0.17) (Figure 7D). Functional enrichment analysis of the yellow module revealed enrichment in the glutamine metabolism pathway (Figure 7C). Next, we performed clustering analysis based on the expression levels of the 518 differentially expressed genes in the TCGA COAD cohort, grouping patients into three clusters (C1, C2, C3) (Figures 7E,F). We analyzed the expression levels of SLC38A2 in these patient subgroups and found higher expression in clusters C1 and C3 compared to normal samples (Figure 7G).

Similarly, in LUAD, we identified 244 differentially expressed metabolic genes and used WGCNA to cluster these genes into four gene modules (Figures 8A,B). The blue module was significantly associated with cancer prognosis (pathologic_N = 0.12, pathologic_T = 0.16, OS.time = −0.11, pathologic_stage = 0.14) (Figure 8D). Functional enrichment analysis of the blue module indicated that these genes are related to the glutamine metabolism pathway (Figure 8C). Further clustering results showed significant survival differences among the C1, C2, and C3 patient cohorts, with patients in the C3 group exhibiting normal survival rates (Figures 8E–G). Additionally, the expression levels of SLC38A2 in all three patient groups were lower than in normal samples (Figure 8H).

In BRCA, SLC38A2 expression is positively correlated with ADORA2A, CD160, CD274, CSF1R, HAVCR2, IL10, IL10RB, KDR, PDCD1LG2, TGFBR1, and VTCN1, and negatively correlated with IDO1, KIR2DL1, LAG3, LGALS9, and PDCD1. In LUAD, SLC38A2 expression is positively correlated with CD160, CD274, IL10, KDR, PDCD1LG2, and TGFBR1, and negatively correlated with LGALS9 and PDCD1. In MESO, the expression of SLC38A2 is positively correlated with CD160, IL10, KDR, TGFB1, and TGFBR1, and negatively correlated with LAG3 and LGALS9. And in PAAD, the expression of SLC38A2 is positively correlated with ADORA2A, BTLA, CD160, CD244, CD274, CD96, CSF1R, CTLA4, HAVCR2, IDO1, IL10, KDR, KIR2DL3, LAG3, PDCD1, PDCD1LG2, PVRL2, TGFB1, TGFBR1 and TIGIT (Figure 5E).

The GO enrichment analysis indicates that SLC38A2 is primarily associated with the following pathways in cancer. These include cell surface receptor protein serine/threonine kinase signaling pathway, transforming growth factor beta receptor superfamily signaling pathway, regulation of cell-matrix adhesion, regulation of bone formation and development, regulation of monocyte chemotaxis, and transmembrane transport of amino acids (Figure 9A).

First, we evaluated the correlation between SLC38A2 expression and patient age. We found that SLC38A2 expression was higher in younger patients, especially in six types of cancer: ESCA, KIRC, PAAD, READ, TGCT, and THCA (Figure 9B). Time-dependent AUC results indicate that SLC38A2 has a higher predictive value in three types of cancer that include DLBC (1-year AUC = 0.792, 5-year AUC = 0.757), KICH (1-year AUC = 0.898, 3-year AUC = 0.766, 5-year AUC = 0.790), and PAAD (5-year AUC = 0.739) (Figure 9C). Diagnostic ROC results show that SLC38A2 has diagnostic value in COAD, GBM, HNSC, LAML, LIHC, LUAD, LUSC, PAAD, PCPG, PRAD, SKCM, THCA, THYM, UCEC, and UCS (Figure 9D).

Results from the CTRP database indicate that the expression of SLC38A2 is significantly correlated with the sensitivity of 62 small-molecule compounds (Figure 10A), as well as the GDSC database results show a significant correlation between SLC38A2 expression and the sensitivity of 87 small-molecule compounds (p < 0.05) (Figure 10B). Additionally, SLC38A2 expression exhibits a moderate negative correlation with the sensitivity to Dabrafenib and Vemurafenib, a weak positive correlation with the sensitivity to Pyrazoloacridine, and a weak negative correlation with the sensitivity to Selumetinib (Figure 10C). Besides, we analyzed the relationship between the sensitivity of six chemotherapy drugs and the grouping based on SLC38A2 expression (Figure 10D).