The RNA-sequencing and clinicopathological data enrolled in this study were acquired from three public databases and an own cohort. Patients with primary gliomas included in this study. Here, the notion of gliomas is restricted to astrocytomas, oligodendrogliomas, and glioblastomas. Those with recurrent gliomas or age < 18 were excluded from this study. In total, 662 primary gliomas from the Cancer Genome Atlas (TCGA, of which 655 had survival data) were included in our study, and the fragments per kilobase million (FPKM) and survival data of them were downloaded from the TCGA website (https://portal.gdc.cancer.gov/). Another 415 primary gliomas were from the Chinese Glioma Genome Atlas (CGGA) 693 cohort, and the REMBRANDT cohort which consists of 369 primary gliomas. FPKM data of the CGGA cohort and array data of REMBRANDT cohort were obtained from the CGGA website (http://www.cgga.org.cn/). For data preprocessing, the genes with too low expression levels (FPKM maximum < 0.1 or standard deviation < 0.01) were excluded from subsequent analyses.

Our own cohort included 77 primary glioma patients from West China Hospital (WCH). We obtained their tumor tissues during craniotomy and sequenced mRNA of these tumor tissues. Subsequently, we used STAR to quantify the mRNA sequencing data and normalized them to FPKM. Survival data of these patients were acquired by regular follow-up, and the overall survival (OS) was defined as the period from surgery to death or the time of last follow-up (censored value). Additionally, the patients younger than 18 years old were excluded from analyses in all four cohorts. The detailed clinicopathological information of these patients is given in Table 1.

By searching the Molecular Signature Database (MSigDB) with the keyword “tryptophan metabolism” or “tryptophan metabolic process,” we identified 56 tryptophan metabolism-related genes (TrMGs), and 44 of which were kept after excluding the genes with low expression level. Subsequently, we performed unsupervised K-means clustering analysis to illuminate the distinctive tryptophan metabolism patterns in gliomas based on the expression patterns of tryptophan metabolism-related genes. We utilized the R package “factoextra” to determine the optimal number of clusters, which corresponded to the maximum average silhouette width (average distances of points to the centroids of clusters that it does not belong to minus the distance of points to the centroid of the cluster that it belongs to). To visualize the different expression patterns of TrMGs in each cluster, we conducted the t-Distributed Stochastic Neighbor Embedding (tSNE) analysis. Additionally, we utilized the TrMGs expression and cluster labels based on the TCGA cohort to construct a naïve Bayes classifier, and then stratified the patients of the other cohort into different clusters using this classifer.

Based on the expression of TrMGs, we constructed a risk signature system to elucidate the relationship between TrMGs and gliomas. First, we split the patients of TCGA cohort into training and validation tests with a ratio of 6:4. All the other three cohorts were utilized as validation sets. In the training set, we selected the TrMGs using the Least Absolute Shrinkage And Selection Operator (LASSO) Cox regression analysis. The TrMGs were determined as essential TrMGs in glioma if their coefficient was not zero at the lambdas corresponding to maximum C-index in over 80 random repetitions of LASSO Cox regression out of 100. Furthermore, a final multivariate Cox regression model was fitted to the training set with essential TrMGs. The tryptophan metabolism-related genes risk signature (TrMRS) was calculated using the following formula: T r M R S   R i s k   S i g n a t u r e = ∑ i = 1 ( β i ∗ Exp i )

In this calculating formula, β stand for the coefficient of each essential TrMG as fitted by the final multivariate Cox regression model, and Exp represented for the expression level of each essential TrMG. Subsequently, we determined the optimal TrMRS cut-off value by “surv_cutpoint” in the R package “survminer” with group proportion ≥ 0.3 for each dataset. Based on the optimal cut-off value, we allocated all patients into TrMRS low-risk or high-risk group. Additionally, we depicted the receiver operating characteristic (ROC) curve in validation sets of 1, 2, and 3-years survival and computed the area under the ROC curve (AUC) using the R package “timeROC.”

For the analyses of gene alternations and copy number variations (CNVs), we obtain the data of gene alterations and CNVs from the cBioPortal database (https://www.cbioportal.org/) for the TCGA cohort. The R package “maftools” was utilized to illustrate the different patterns of gene alterations and tumor mutation burdens (TMB) between different K-means clusters and TrMRS risk groups. Furthermore, the Genomic Identification of Significant Targets in Cancer (GISTIC) score was used to assess the different CNV levels in different clusters and risk groups.

To interpret the biological functions of the differential transcriptomes between clusters and TrMRS risk groups, we utilized the over-representation and gene set enrichment analysis (GSEA) to evaluate the differentially expressed genes (DEGs) with the R package “clusterProfiler”. R package “limma” was used to identify the DEGs between different clusters and risk groups. In the process of DEGs identification for GSEA, we stratified the patients into TrMRS high- and low-risk groups based on cut-off values of each cohort. Those genes with adjusted p-value < 0.05 and |log₂FC| > 0.5 were identified as DEGs. Furthermore, we used the R package “GSVA” to convert the logFPKM matrix of genes to pathway expression matrix. We identified the differentially expressed pathways between different clusters and risk groups with the “limma” package. To calculate the infiltration fraction of immune cells in glioma, we utilized the CIBERSORTx (https://cibersortx.stanford.edu/) and TIMER 2.0 [http://timer.comp-genomics.org/, which is based on signature genes correlated with estimated tumor purity and immune cell fractions (Li et al., 2016)]. Additionally, the Estimation of Stromal and Immune Cells in Malignant Tumor tissues using Expression data (ESTIMATE) was used to evaluate the infiltration of stromal and immune cells in tumor microenvironment and calculate the stromal, immune and ESTIMATE scores (Yoshihara et al., 2013). In this algorithm, the stromal-related genes were selected from the non-hematopoiesis-related genes that were differentially expressed between tumor cell fraction and match stromal cells fraction separated by laser capture microdissection in multiple cancers. Furthermore, the tumor purity data published by Aran et al. (2015) which included ESTIMATE score-based tumor purity and consensus purity estimation (CPE), were used to represent the tumor purity in gliomas. Another previously published algorithm was utilized to compute the tumor immunological phenotype (TIP) gene signature (Wang et al., 2021). Based on TIP signature, we identified immunological phenotypes of gliomas and distinguish relatively “hot” tumors from “cold” tumors. Furthermore, we also utilized the Tumor Immune Dysfunction and Exclusion (TIDE) suite (http://tide.dfci.harvard.edu/) to predict response to immune checkpoint inhibitors therapy in gliomas.

To determine prognostic factors in glioma, we used univariate and multivariate Cox regression analyses. We firstly enrolled TrMRS and other potential prognostic factors, including tumor grade, age, chemotherapy, radiotherapy, KPS, gender, 1p/19q codeletion, and IDH mutation, into the univariate Cox regression analysis. Subsequently, those factors with a p-value < 0.05 in univariate analysis were enrolled into multivariate Cox regression analysis to identify independent prognostic factors. Those factors with a p-value < 0.1 in multivariate Cox regression analysis were identified as independent factors and enrolled into construction of nomograms. The R package “rms” was utilized to construct the nomograms. The calibration curves were used to evaluate the efficacy of nomograms for prognosis prediction in glioma patients.

In all proceedings of bioinformatic analyses, we used the R software (version 4.2.1). To assess the differences between clusters or risk groups for continuous variables, we used the Wilcoxon rank sum test. To assess the differences for categorical variables, we used the chi-square test. The R package “survminer” was used to deliver survival analysis and generate Kaplan-Meier (K-M) curves, which were tested for differences by log-rank test. The “coxph” function in the R package “survival” was applied to conduct Cox regression analyses. The R package “glmnet” was utilized to perform the LASSO Cox regression analysis. T Iterative Grubbs tests were applied to expel the outliers in liner regression analysis.

Tumor samples and clinical data collection and use were performed strictly with ethics regulations and approved by the institutional review board of West China Hospital (No. 2018.569) based on local ethics regulations and the 1964 Helsinki declaration and its later amendments. In addition, the patients signed written consent for tumor tissue collection and processing. The sequencing data of West China Hospital generated in this study are available at the Genome Sequence Archive for Humans: accession code HRA002839 (access link: https://ngdc.cncb.ac.cn/gsa-human/s/JQssVoV1).

Based on 44 identified tryptophan metabolism-related genes, we conducted an unsupervised K-means clustering analysis in the TCGA cohort. Then the gliomas in the TCGA cohort were stratified into two clusters according to the average silhouette widths as described in the Material and Method section (Supplementary Figure S1A). Distinctions in TrMGs expression profile between two clusters were illustrated by tSNE analysis (Figure 1A). The expression levels of TrMGs in each glioma was given out in heatmap ordered by clusters (Supplementary Figure S1B). The critical metabolites and enzymes in two branches of tryptophan metabolism, kynurenine pathway and serotonin pathway, were interpreted in a schematic diagram (Figure 1B). The differences in the expression levels of ten important TrMGs involved in the tryptophan metabolism were also illustrated (Figure 1C). The cluster 2 was demonstrated with higher expression of IDO1, TDO2, KYNU, KMO, and HAAO, which were critical enzymes in the kynurenine pathway for tryptophan metabolism. Besides, KAT2, HADHA, GCDH, DDC, and ALDH2 were highly expressed in cluster 1, suggesting more consumption of metabolites of kynurenine pathway and more active serotonin pathway. The differences in the expression levels of each TrMG were shown in Supplementary Figure S1C.

Survival analysis revealed that the prognosis of cluster 1 was significantly better than cluster 2 (median OS, 63.8 vs. 17.7 months) (Figure 1D). The patients in the CGGA, REMBRANDT, and WCH cohorts were also stratified into two clusters based on the naïve Bayes clustering classifier train with the TCGA data. Survival analyses in these three cohorts also revealed significantly better prognosis in cluster 1 (Figures 1E–G), suggesting that the TrMG expression patterns were robust among different glioma cohorts.

The functional enrichment analyses based on two clusters depicted distinctive pathway alterations. The ECM receptor interaction, cell adhesion molecules, other pathways were identified as key pathway associated with differentially expressed genes (DEGs) between clusters in the KEGG gene sets (Figure 2A). In REACTOME gene sets, extracellular matrix organization and interferon gamma signaling pathways were found to be over-represented in the cluster DEGs (Figure 2B).

The gene alteration analyses revealed a different gene alteration landscape between two clusters. Cluster 1 gliomas were characterized by IDH1 mutations with co-occurring TP53 and ATRX mutations or CIC, FUBP1 and NOTCH1 mutations (Figures 2C,D). Higher proportion of EGFR could be found in cluster2, which were mutually exclusive with IDH1 mutations but co-occurrent with PTEN mutations (Figures 2C,E). The analysis of CNVs demonstrated that 1p/19q co-deletion, which was recognized as an essential diagnostic marker for oligodendroglioma, mainly occurred in the cluster 1 (Figure 2F). The gain of chromosome 7 and loss of chromosome 10 (+7/-10), which was a novel diagnostic marker for glioblastoma and indicated worse prognosis, occurred more frequently in cluster 2 compared to cluster 1, in line with the worse prognosis of cluster 2 in survival analysis. Analysis of clinicopathological features of these two clusters also demonstrated that cluster 2 had higher tumor grade, lower incidence of IDH mutation, and lower incidence of 1p/19q codeletion (Figure 2G). The differences in other clinicopathological features between two clusters were given in Supplementary Figure S2.

Based on these two K-means clusters, we conducted multiple analyses of immunological features to elucidate the differences in tumor immune landscape. The results of the TIMER score of TCGA cohort revealed that the tumor microenvironment of cluster 2 carried more macrophages, neutrophils, and CD8⁺ T-cells (Figure 3A). This conclusion could also be verified in the CGGA cohort (Figure 3B), suggesting a more complex tumor microenvironment in cluster 2. Furthermore, the results of the ESTIMATE revealed that the stromal score, immune score, and ESTIMATE score of the cluster 2 were significantly higher than cluster 1 in all four cohorts (Figure 3C). The tumor purity of cluster 2 was significantly lower than cluster 1 in all four cohorts (Figure 3C), indicating purer tumor cell microenvironment in gliomas of cluster 1. Through the calculation of TIP score, we were able to demonstrate that most gliomas of cluster 2 were more likely to be “hot” tumor immunological phenotype and expressed higher levels of gene markers for ‘hot’ tumors compared to gliomas in cluster 2 (Figure 3D). The analysis of TIP score in the CGGA cohort reached similar results (Figure 3E). Additionally, the TIP scores for the gliomas of cluster 2 were significantly higher than cluster 1 in all four cohorts (Figure 3F). All these results supported that the gliomas of cluster 2 had more complex tumor microenvironment, more immune cell infiltration, and presented with features of immunologically ‘hotter’ tumors.

To determine essential genes for the construction of TrMRS, we filtered the 44 TrMGs using the LASSO Cox regression in the training dataset. Subsequently, eight genes, including AOX1, CYP2E1, KMO, KYNU, ALDH2, OGDH, HSD17B10, and MAOB, were identified as essential TrMGs for the construction of TrMRS (Figure 4A). The formula to compute the TrMRS was as follows: 1.216*KYNU+0.254*KMO+0.231*AOX1+0.034*OGDH+0.006*HSD17B10 + 0.001*MAOB-0.066*ALDH2-1.064*CYP2E1.

Six genes of these essential TrMGs, including MAOB, HSD17B10, KYNU, AOX1, and OGDH, were determined as hazardous factors for glioma patients, and other two genes (CYP2E1 and ALDH2) were identified as protective factors (Figure 4B). Besides, representative immunohistochemical (IHC) staining for KYNU and ALDH2 from the Human Protein Atlas (Ponten et al., 2008) (https://www.proteinatlas.org/) was utilized to validate these results. The results of IHC demonstrated that the expression level of KYNU was higher in high-grade glioma compared to low-grade glioma (Supplementary Figure S3A), and the expression level of ALDH2 was higher in low-grade glioma (Supplementary Figure S3B), confirming the results from sequencing that KYNU was hazardous factor and ALDH2 was protective factor. Furthermore, we utilized the “surv_cutpoint” algorithm to calculate the optimal TrMRS cut-off for every cohort and allocated all the patients into TrMRS low- or high-risk groups according to this cut-off (Figure 4C). Survival analyses revealed that the overall survival of the patients in TrMRS high-risk groups was significantly poorer than TrMRS low-risk group in all four validation cohorts (Figures 4D–G), suggesting that TrMRS potentially functioned as factor for prognosis prediction. To test the efficacy of TrMRS in predicting prognosis of glioma patients, we performed ROC analyses to evaluate the performance of TrMRS alone in survival rate prediction at 1, 2, and 3 years. In the TCGA test set, the AUCs of TrMRS at 1, 2, and 3 years were 0.820, 0.818, and 0.819, respectively (Figure 4H). In other three validation cohorts, similar performances were also achieved (Figure 4H).

The expression pattern of eight essential TrMGs was exhibited with a heatmap ordered by TrMRS (Figure 5A). The relationship of other clinicopathological features, including tumor grade, histological diagnosis, IDH mutation status, 1p/19q codeletion, TERT promoter status, ATRX status, and MGMT promoter status, with the TrMRS were also given out (Figure 5A). Further analysis of gene mutations discovered that IDH1, TP53, and ATRX were top 3 most frequently mutated genes in TrMRS low-risk group (Figure 5B). TP53, PTEN, and EGFR were top 3 most frequently mutated genes in TrMRS high-risk group (Figure 5C). Additionally, analyses of tumor mutation burden (TMB) revealed that TrMRS high-risk group harbored higher TMB than low-risk group, and the correlation analysis demonstrated that TMB was positively correlated with the TrMRS (Figure 5D). Further analysis manifested the incidence of gene amplification, including EGFR, SEC61G, and LANCL2, was significantly higher in high-risk group than low-risk group (Figure 5E). Besides, the incidence of gene homozygous deletion, including CDKN2A and CDKN2B, was also higher in high-risk group (Figure 5F).

To explore the impact of TrMRS in terms of pathway, we utilized KEGG and REACTOME pathway gene sets in enrichment analyses of the DEGs between the TrMRS risk groups. The patient samples were stratified into TrMRS high- and low-risk groups. The analyses of GSEA were based on the differentially expressed genes (DEGs) between TrMRS high- and low-risk groups. The complement and coagulation cascades pathway (normalized enrichment score (NES) = 3.085, adjusted p-value < 0.001) and graft versus host disease pathway (NES = 2.855, adjusted p-value < 0.001) were ranked in the top five gene sets of the KEGG in the comparison between two risk groups (Figure 6A) using GSEA. Besides, the cytokine signaling in immune system pathway (NES = 2.950, adjusted p-value <0.001) and innate immune system pathway (NES = 3.248, adjusted p-value < 0.001) of REACTOME gene sets were ranked in the top five (Figure 6B). Furthermore, over-representation of the extracellular matrix receptor interaction of KEGG gene sets and the extracellular matrix organization of REACTOME gene sets were identified in the biological functions of DEGs between the two TrMRS risk groups (Figures 6C,D). Finally, the GSVA result demonstrating the top differentially expressed pathways in KEGG and REACTOME gene sets were also illustrated through heatmaps (Figures 6E,F).

To construct nomograms for the prediction of glioma patients’ prognosis, we first conducted univariate followed by multivariate Cox analyses to identify potential independent prognostic factors. Result demonstrated that the tumor grade, age, radiotherapy, TrMRS, chemotherapy, 1p/19q codeletion, and IDH mutation were significant univariate prognostic factor (Figure 7A). Subsequently, these factors were enrolled in multivariate Cox regression analysis and the result revealed that TrMRS, tumor grade, 1p/19q codeletion, and IDH mutation were independent prognostic factors in glioma (Figure 7B). Eventually, these factors were combined in the construction of a nomogram for personalized survival prediction (Figure 7C). The four factors were also used to construct a nomogram for the CGGA cohort (Figure 7D). The corrected C-indexes of the nomograms based on TCGA and CGGA cohort were 0.851 and 0.779, respectively. Additionally, the 1-, 2-, and 3-years calibration curves of the nomograms also validated their efficacy in predicting survival time of glioma patients (Figures 7E,F).

To elucidate the relationship between the TrMRS and tumor immune microenvironment in glioma, we conducted bundles of analyses on immunological features based on TrMRS. Firstly, the CIBERSORTx estimation of immune cell fractions depicted distinctive patterns of immune cell infiltration between two TrMRS risk groups. The TrMRS high-risk group manifested with more infiltration of resting NK cells, Macrophages (M0, M1, and M2), and neutrophils in its tumor microenvironment (Figure 8A). The low-risk group, on the other hand, had higher infiltration of activated NK cells and plasma cells. Further correlation analyses demonstrated that the infiltration of plasma cells and activated NK cells was negatively correlated with the TrMRS, and the infiltration of M2 macrophage and neutrophils was positively correlated with the TrMRS (Figure 8B). Besides, analyses of immune scores demonstrated that the TrMRS high-risk group presented with significantly higher stromal score, immune score, and ESTIMATE score (Figure 8C). The tumor purity of TrMRS high-risk group was remarkably lower than low-risk group (Figure 8D), suggesting more complicated tumor microenvironment in gliomas of TrMRS high-risk group with more immune cell infiltration. Correlation analyses revealed that the stromal score, immune score, ESTIMATE score, and tumor purity were strongly correlated with the TrMRS (R = 0.82, 0.762, 0.805, and 0.808) (Figure 8E).

To predict potential response to immunotherapy response, we also performed analysis on the expression of immunity-related genes. Gliomas of TrMRS high-risk group presented with significantly higher expression level of CD274 (PD-L1), CD276 (B7H3), HAVCR2 (TIM3), PD1 (CD279, PDCD1), and CD44 (Figure 9A). Correlation analysis also confirmed that the expression levels of these immunotherapy-related markers were positively correlated with the value of TrMRS (Figure 9B), suggesting potential ability of TrMRS to guide the choice of immunotherapy. Furthermore, immune phenotype analysis of high- and low-risk gliomas revealed that most gliomas of TrMRS high-risk group were identified as more likely to be immunological “hot” tumors, and most tumors of low-risk group were presumed “cold” tumors (Figure 9C). A strong positive correlation between TrMRS and TIP score was also confirmed by the correlation analysis (Figure 9D). The analyses of TIP score based on the CGGA cohort also manifested with similar results (Figures 9E,F). Additionally, we utilized TIDE algorithm to predict response to immune checkpoint inhibitors. The result revealed that the glioma patients of TrMRS high-risk group were more likely to benefit from therapy of immune checkpoint inhibitors in the TCGA and CGGA cohorts (Figures 9G,H). Most above findings can be validated in the other cohorts (Supplementary Figure S4).