Systematic analysis of glutamine metabolism family genes and exploration of the biological role of GPT in colorectal cancer

Data download and processing

The RNA-seq expression data and survival information from the CRC cohort (TCGA-COAD/READ) were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Expression data and prognostic information of 31 cancers for pancancer analysis were downloaded from UCSC database (https://xena.ucsc.edu/). Single-cell RNA sequencing (scRNA-seq) data of 12 CRC patients were downloaded from GSE166555. The scRNA-seq data in GSE166555 were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. The glutamine metabolism-related genes involved in this study were obtained from previous studies published by He. et al. [23]. Analysis of the differential expression of GMRGs in the TCGA-CRC cohort using the R package “Limma”, and the screening condition was set to “|Foldchange > 2, P < 0.05|”.

Machine learning screening for differentially expressed GMRGs

For the differentially expressed GMRGs, we further used SVM-REF and RandomForest algorithm for separation to screen out the characteristic genes. SVM-REF is a powerful classification algorithm [24], that can be applied to characterize cancer, allowing the identification of biomarkers for diagnosis. Similarly, RandomForest is also a powerful classification algorithm [25] that ranks gene importance, and in this study, we picked genes with an importance score greater than 2. Finally, we use the common genes identified by these two machine learning algorithms for further research and consider them as characteristic genes of CRC.

scRNA-seq analysis of characteristic genes

The scRNA-seq data were preprocessed using the “Seurat” package. “PercentageFeatureSet” is used to calculate the percentage of mitochondrial genes. Cells with a gene count below 50 and a mitochondrial gene percentage above 5 were eliminated. The 1500 genes with the largest variation between cells were analyzed as cell clusters, the principal component was set to 20, and the resolution was 0.5.

Evaluation of diagnostic performance of characteristic genes

We further used the receiver operating characteristic curve (ROC) based on specificity and sensitivity to evaluate the diagnostic performance of the characteristic gene. The AUC indicated the diagnostic performance of the characteristic gene. Different AUC values mean different diagnostic values [26]. AUC = 1.0 means perfect diagnostic value; AUC = 0.9-1.0 means high diagnostic value; AUC = 0.7-0.9 means relative diagnostic value; AUC = 0.5-0.7 means low diagnostic value; AUC = 0.5 means no diagnostic value.

Evaluation of the prognostic value of characteristic genes

We performed a further analysis on the prognostic value of the characteristic genes. First, we observed the survival curves of different expression levels of the characteristic genes using the Kaplan-Meier method. Second, we further evaluated the ability of characteristic genes to affect the prognosis of CRC using univariate Cox regression analysis and displayed it using a forest plot. In the prognostic analysis, we excluded the samples with missing survival data and survival time ≤30.

External validation and experimental validation of core gene GPT expression

The GEPIA2 database analyzed the differential expression of the key gene GPT at the mRNA level of colon cancer and rectal cancer online (http://gepia2.cancer-pku.cn/#index) [27]. The UCLCAN database analyzed the expression difference of GPT at the protein level in the CPTAC database online (http://ualcan.path.uab.edu/index.html) [28]. The HPA database (https://www.proteinatlas.org/) obtained the immunohistochemical images of GPT in normal colorectal and CRC tissues [29], and obtained the basic information of the corresponding tissue samples.

CRC specimen collection, RNA extraction and quantitative real-time PCR (qRT-PCR) reaction

The 20 pairs of paracancerous normal tissues and matched CRC tissues were all from Shanghai Shidong Hospital. All patients provided written informed consent to participate in this study.

The qRT-PCR experiment was performed according to the previous procedure [30]. Briefly, total RNA from normal and CRC tissues was extracted using RNA-easy isolation reagent (Vazyme, China), and then the RNA was reversed to cDNA using PrimeScript™ RT Master Mix (Takara Bio, Japan). Finally, we performed qRT-PCR reactions using SYBR-Green qPCR Master Mix (Vazyme, China) and normalised using GAPDH. The experimental results were quantified relative to the 2-ΔΔCt method. Primer sequences were designed in this study as follows, The GAPDH primer sequences, forward: GTCTCCTCTGACTTCAACAGCG, reverse: ACCACCCTGTTGCTGTAGCCA. The GPT primer sequences, forward: CAAGAAGGTGCTCATGGAGATGG, reverse: TGTTCACCACCTCCACATAGCC.

Co-expression analysis and enrichment analysis

The co-expression network of GPT was analyzed using the LinkedOmics database (http://www.linkedomics.org/). GO analysis and KEGG analysis explored the biological processes involved in the co-expression network.

Assessment of the level of immune cell infiltration

Based on the gene expression profile of TCGA-CRC, we used CIBERSORT algorithm to evaluate the proportion of twenty-two types of immune cell [31], and samples with p < 0.05 were used for further analysis. Then, this research evaluated differences in the infiltration levels of 22 types of immune cells in CRC patients with different GPT expression.

Somatic mutation spectrum analysis

We used the R package “maftools” [32] to analyze the somatic gene mutation profile of TCGA-CRC patient samples. First, we analyzed the mutation status of characteristic genes in TCGA-CRC samples. Second, we further analyzed the 10 most mutated genes in samples with different levels of GPT expression.

Drug sensitivity analysis

We used R package “oncoPredict” [33] to analyze the half maximal inhibitory concentration (IC50) of TCGA-CRC samples based on gene expression data. The lower the IC50 value, the higher the drug sensitivity, and the larger the IC50 value, the lower the drug sensitivity.

Pancancer analysis of GPT

As previously described, we used the R package “Limma” to explore the expression differences of GPT in pancancer. Likewise, the diagnostic performance of GPT in pancancer was evaluated using ROC curves. The prognostic value of GPT in pancancer was evaluated using univariate Cox analysis and Kaplan-Meier method.

Statistical analysis

All bioinformatics research in this study was done by R language (version 4.1.2). Differences between groups were analyzed using the Wilcoxon test. The Kaplan-Meier method was used to draw survival curves. qRTPCR results were completed using Graphpad prism 9 software. p < 0.05 means that the results are statistically significantly different.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Characteristic gene screening of GMRGs

First, we analyzed the differential expression of 118 glutamine metabolism family genes in the TCGA-CRC cohort, and the results showed that 31 GMRGs were differentially expressed genes (DEG), of which 4 GMRGs were down-regulated and 27 GMRGs were up-regulated. Figure 1A shows the volcano map of GMRG as well as the heat map of DEG (Figure 1B). Then, we performed two machine learning algorithms on DEG to screen out characteristic gene: RandomForest (Figure 1C, 1D) and SVM-RFE (Figure 1E). Detailed results are detailed in Supplementary Table 1. According to the results of the above two machine learning algorithms, we used the Venn diagram to obtain 9 overlapping genes and defined them as characteristic gene (Figure 1F). Finally, we explored the mutation status of CRC signature genes and found that CAD (6%), SLC3A2 (5%) and MET (4%) had the highest mutation rates (Figure 1G).

Figure 1. Characteristic gene screening. (A) Volcano map of differential expression of GMRGs. (B) Heatmap of DEG. (C, D) RandomForest evaluates the relative importance of DEG. (E) SVM–RFE algorithm. (F) Venn diagram of two machine learning algorithms. (G) Somatic mutation profiles of characteristic genes.

scRNA-seq analysis of characteristic genes

First, we processed single-cell sequencing data of CRC cancer tissues to classify cells into 18 clusters (Figure 2A). We then displayed the 10 genes that were significantly overexpressed in each cluster in Figure 2B. In addition, we automatically annotated cell types through the “SingleR” package, and finally determined that there are 9 main types of cells in CRC tissues, including T cells, B cells, Epithelial cells, Monocyte, Tissue stem cells, Smooth muscle cells, Endothelial cells, NK cells, Fibroblasts (Figure 2C).

Figure 2. scRNA-seq analysis of CRC tissue. (A) tSNE analysis of CRC tissue scRNA-seq data to classify cell clusters. (B) Heatmap showing the top 10 genes highly expressed in each cell cluster. (C) The “Sin-gleR” package annotates the cell clusters into 9 cell types. (D) The distribution of the characteristic genes in the 9 cell types. (E) The bubble plot shows the expression levels of the characteristic genes in the 9 cell types.

Second, we explored the connections between cells using cell communication analysis. We visualized the number (Figure 3A) and weight (Figure 3B) of interactions between the 9 cell types. And centering on each cell type, the communication relationship between it and 9 cell types is shown separately (Figure 3C–3K). Furthermore, we use heatmaps to visualize the ligand-receptor molecules that mediate interactions between cells (Figure 3L).

Figure 3. Cell communication analysis of scRNA-seq data. (A) Network diagram showing the number of connections between different cell types. (B) Network diagram showing the weight of connections between different cell types. (C–K) Diagram of the communication network between different cells and other cells. (L) Bubble plot showing genes involved in cell communication.

Finally, we explored the distribution of characteristic genes in the TME (Figure 2D). Bubble plots show the expression levels of signature genes in each cell type (Figure 2E). Combining the above results, we found that SLC7A5 was highly expressed in monocyte, B cells and T cells. MET is mainly expressed in endothelial cells and epithelial cells. SHMT2 has similar expression level in other cell types except T cells. CAD and PPAT were expressed at low levels in all cell types. PYCR1 is highly expressed in tissue stem cells, epithelial cells and fibroblasts. GPT is almost exclusively expressed in epithelial cells. SLC3A2 was highly expressed in monocytes, but very low in fibroblasts. MYC is mainly expressed in endothelial cells, smooth muscle cells and epithelial cells.

Diagnostic and prognostic value of characteristic genes

First, we analyzed the diagnostic performance of the characteristic genes using ROC curves. We found that the 9 characteristic genes all had high AUC values (Figure 2A–2I), specifically: CAD (AUC = 0.976), GPT (AUC = 0.973), MET (AUC = 0.983), MYC (AUC = 0.974), PPAT (AUC = 0.962), PYCR (AUC = 0.972), SHMT2 (AUC = 0.981), SLC3A2 (AUC = 0.968) and SLC7A5 (AUC = 0.992). These results imply that all characteristic genes are a potentially powerful class of diagnostic markers.

At the same time, we analyzed the survival curves of the 9 characteristic genes with different expression levels using the Kaplan-Meier method (Figure 2J–2U). We found that high expression of CAD and SLC7A5 had a worse prognosis, while high expression of GPT, MET, MYC, and PPAT had a better prognosis. The results showed that high expression of CAD and SLC7A5 had worse prognosis, while high expression of GPT, MET, MYC and PPAT had better prognosis. However, the prognosis of patients with different expression levels of PYCR, SHMT2, and SLC3A2 was not significantly different.

Finally, we used univariate Cox analysis and found only GPT as a protective factor affecting the prognosis, and none of the remaining 8 characteristic genes were factors affecting the prognosis in CRC (Figure 4A). Combining the above results, we regard GPT as the most important GMRGs in CRC and conduct further analysis.

Figure 4. Evaluation of diagnostic performance and prognostic value of characteristic genes. (A–I) ROC curves of 9 characteristic genes. (J–U) KM survival curves of 9 characteristic genes.

Co-expression network and enrichment analysis of GPT in CRC

In order to fully understand the biological function of GPT in CRC, we analyzed the co-expression network of GPT in the LinkedOmics database online. The results showed that GPT expression was positively correlated with 3462 genes and negatively correlated with 5207 genes (Figure 5B). The heatmap shows the top 50 most significantly positively and negatively correlated genes with GPT expression, respectively (Figure 5C, 5D). We further performed GO analysis and KEGG analysis on the top 200 significantly positively correlated genes. The bubble charts show the top 15 GO analysis terms (including 5 BP, CC and MF terms respectively) and the 10 KEGG analysis terms (Figure 5E, 5F). GO analysis results showed that, in terms of BP, mainly enriched in hormone metabolic process; one−carbon metabolic process; cellular modified amino acid metabolic process; primary alcohol metabolic process; steroid metabolic process. in terms of CC, mainly enriched in brush border; brush border membrane; cluster of actin−based cell projections; apical part of cell; microvillus. in terms of MF, mainly enriched in carbonate dehydratase activity; oxidoreductase activity, acting on the CH−OH group of donors, NAD or NADP as acceptor; oxidoreductase activity, acting on CH−OH group of donors; steroid dehydrogenase activity, acting on the CH−OH group of donors, NAD or NADP as acceptor; carboxylic ester hydrolase activity. KEGG analysis mainly enriched in Pancreatic secretion; Bile secretion; Steroid hormone biosynthesis; Retinol metabolism; Glycerolipid metabolism. In short, the results suggest that GPT co-expressed genes are mainly involved in various metabolic pathways.

Figure 5. Co-expression analysis and enrichment analysis. (A) Univariate Cox regression analysis of characteristic genes. (B) Correlation of GPT ex pression with other genes. (C) Heatmap showing the 50 genes with the strongest positive correlation with GPT expression. (D) Heatmap showing the 50 genes most negatively correlated with GPT expression. (E) GO analysis. (F) KEGG analysis.

External validation and experimental validation of GPT expression

We first explored the differential expression of GPT in COAD and READ using the GEPIA2 database. We found that GPT mRNA expression levels were low in COAD and READ (Figure 6A). In addition, we detected GPT mRNA levels of GPT using qRT-PCR. The results also showed that compared with CRC tissues, GPT was highly expressed in normal tissues (Figure 6B). In addition, we also evaluated the expression difference of GPT at the protein level. First, we used the CAPTC database on UALCAN to analyze the differences in GPT protein expression levels, and found that GPT protein levels were significantly decreased in COAD tissues compared with normal tissues (Figure 6C). In addition, we obtained immunohistochemical images of normal colon and rectal tissues and CRC tissues from the HPA database. We found that the intensity of GPT in CRC tissues was significantly reduced compared to normal tissues (Figure 6D–6G).

Figure 6. External data and experimental data verify GPT expression of GPT in CRC. (A) GEPIA2 database showing differential expression of GPT in COAD and READ (mRNA level). (B) qRT-PCR showed that GPT was significantly lower in CRC tissues than in adjacent normal tissues. (C) The UALCAN database showed that the low expression of GPT in COAD was correlated with cancer tissues (protein level). (D, E) HPA database showing IHC images of GPT in COAD and colon normal tissues. (F, G) HPA database showing IHC images of GPT in READ and rectal normal tissues.

Association analysis between GPT expression and TME

Our research assessed the infiltration levels of 22 immune cells in TCGA-CRC tissue samples using the CIBERSORT algorithm. We found that compared to the GPT low expression group, the infiltration levels increased significantly increased in the GPT high expression group, while the infiltration levels of T cells CD4, gamma delta of T cells gamma delta, Macrophages M0 and Macrophages M2 decreased significantly decreased (Figure 7A).

Figure 7. Correlation between GPT expression and the tumor microenvironment. (A) Differential analysis of immune cell infiltration levels with different levels of GPT expression. (B) Correlation between GPT expression and immune checkpoint molecules. (C, D) Somatic mutation profiles at different GPT expression levels.

We also analyzed the correlation between GPT expression and 48 immune checkpoint molecules (ICM), in which GPT expression was correlated with 26 ICM, including positive correlations with 17 ICM and negative correlations with 9 ICM (Figure 7B).

Furthermore, we explored the somatic mutation profile in patients with CRC and identified the 10 most mutated genes in different GPT expression level groups. The results showed that APC, TP53, KRAS, TTN, PIK3CA, MUC16, and SYNE1 were the 8 most mutated genes in different groups of expression level groups, and APC had the highest mutation rate (Figure 7C, 7D).

Prediction of potentially sensitive drugs

In order to explore the sensitivity of patients with different GPT expression groups to common chemotherapy drugs, we used the R package “oncoPredict” to analyze the IC50 values of CRC patients in the TCGA-CRC cohort to 5−Fluorouracil, Oxaliplatin, Gefitinib, Tamoxifen, Sorafenib and Dabrafenib. The results showed that the IC50 values of 5−Fluorouracil (Figure 8A), Oxaliplatin (Figure 8B), Gefitinib (Figure 8C), Tamoxifen (Figure 8D) and Sorafenib (Figure 8E) were higher in the GPT low expression group, whereas those of Dabrafenib (Figure 8F) were lower.

Figure 8. Drug sensitivity analysis. Sensitivity analysis of 5-Fluorouracil (A), Oxaliplatin (B), Gefitinib (C), Tamoxifen (D), Sorafenib (E), Dabrafenib (F) in different GPT expression levels.

Diagnostic performance of GPT in pancancer

As mentioned above, the diagnostic performance of GPT in pancancer is evaluated by the ROC curve. We found that GPT had excellent diagnostic performance in CHOL (AUC: 0.985). In BRCA (AUC: 0.711), CESC (AUC: 0.790), ESCA (AUC: 0.785), GBM (AUC: 0.780), HNSC (AUC: 0.894), KICH (AUC: 0.775), KIRP (AUC: 0.827), LUAD (AUC: 0.767), LUSC (AUC: 0.748), PCPG (AUC: 0.887), THCA (AUC: 0.716) and THYM (AUC: 0.786) had higher diagnostic performance (Figure 9).

Figure 9. RUC curve to evaluate the diagnostic performance of GPT expression in 22 cancer types.

Expression of GPT in pancancer

GPT was differentially expressed in 16 cancers compared to normal tissues. Specifically, GPT was low expressed in BRCA, CHOL, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PCPG, STAD, UCS, and highly expressed in BLCA, LUAD, LUSC, THCA (Figure 10A). Furthermore, we sorted the expression of the tumor tissues, we found that GPT had the lowest expression level in TGCT and the highest expression level in LIHC (Figure 10B).

Figure 10. GPT expression patterns in pancancer. (A) Different expression of GPT in pancancer. (B) Ranking map of GPT expression in tumor tissues (from high to low).

Prognostic value of GPT in pancancer

As previously mentioned, we explored the role of GPT expression in the prognosis of pancancer using univariate Cox analysis and drew a forest plot. We found that GPT expression can affect the prognosis of ACC, BRCA, KIRC, KIRP, LGG, LIHC, and THYM. Specifically, GPT acts as a protective factor affecting the prognosis of ACC (HR: 0.373), BRCA (HR: 0.832), KIRC (HR: 0.869), KIRP (HR: 0.720), LGG (HR: 0.739), LIHC (HR: 0.840) and acts as a risk factor affecting the prognosis of THYM (HR: 2.495) (Figure 11A). Finally, for these 7 cancer types, we used the Kaplan-Meier method to analyze the survival curves of different expression levels of GPT, and the results showed that the GPT high expression group had a better prognosis in ACC (Figure 11B), BRCA (Figure 11C), KIRC (Figure 11D), KIRP (Figure 11E), LGG (Figure 11F), LIHC (Figure 11G), and in THYM (Figure 11H) The prognosis is poor.

Figure 11. Prognostic value of GPT in pancancer. (A) Univariate Cox analysis to evaluate the correlation between GPT expression and OS in the patient. (B–H) Survival curves of GPT expression in 7 cancers.