The data related to mutation, copy number variation, and RNA-Seq profile of TCGA-STAD via TCGA GDC API were downloaded first (http://cancergenome.nih.gov/). GSE66229 expression data from the GEO database were obtained (https://www.ncbi.nlm.nih.gov/geo/). Samples need to be processed as follows (Brenner et al., 2009): excision of samples of primary tumors (Lambert et al., 2002); removal of incomplete samples with clinical characteristic information to ensure that samples have complete clinical prognostic information and transcriptome expression data. TCGA-STAD was the training set and the GSE66229 dataset was the validation set. We excluded samples without survival time or status. Finally, a sum of 350 primary tumor samples together with 32 normal matches were screened. We kept the protein-encoding genes for TCGA RNA-seq data analysis. Meanwhile, all data were log2 transformed, and RNA expression data were normalized. For the GEO dataset, 300 GC samples were finally screened. Specifically, the annotation information for each chip platform was acquired and subsequently utilized to map probes with all the detected genes. Then, we removed the probes that matched more than one gene. When one gene can be matched by more than one probes, the mean value of the gene expression was calculated and set as the value for that specific gene.

The tryptophan metabolism-related gene information is derived from a specific pathway named “KEGG_TRYPTOPHAN_METABOLISM”, which can be found in the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/).

We constructed a consensus matrix and clustered the samples through consensus clustering (Wilkerson and Hayes, 2010). The transcriptional expressions of prognostic genes, which are involved in the tryptophan metabolism score, were subsequently evaluated to determine the molecular subtypes. We executed 500 bootstraps employing the “hc” algorithm and “pearson” as the metric distance. Each bootstrap process involved around 80% of the training set patients. The cluster number was set within a range of 2–10, and we determined the optimal classification via cumulative distribution function (CDF). Notably, the consistency of CDF was carefully evaluated when constructing various molecular subtypes for GC samples.

The distinctly expressed genes were identified among the molecular subtypes, and then the distinctly expressed genes associated with statistical significance for prognosis were sorted out (|Log2FC|>1 and FDR<0.05). The gene numbers were further condensed by Lasso regression, and the prognostic genes potentially contributing to the tryptophan metabolism phenotype were filtered. Then, the risk score for each patient was calculated by the formula:

RiskScore = Σβi×Expi, where Expi is the expression value of each prognostic gene that determines tryptophan metabolism phenotype. The Cox regression coefficient of corresponding prognostic gene is referred to as β. Based on the calculated numbers, samples were then distributed into two subgroups, that is, RiskScore-high and RiskScore-low groups, with the threshold set as “0". The commonly used Kaplan-Meier method was utilized to analyze patient survival, and the patients’ statistical significance was calculated via a log-rank test.

In different molecular subtypes, pathways of different biological processes were investigated by performing GSEA for signaling pathway analysis based on the candidate gene sets from the KEGG/hallmark (https://www.gsea-msigdb.org/gsea/index.jsp).

Relative abundance of 22 immune cells in tumor tissues and percentage of immune cells was determined using CIBERSORT algorithm (https://cibersort.stanford.edu/) and ESTIMATE software (Wilkerson and Hayes, 2010), respectively. A total of 28 immune cells were scored using ssGSEA function (Charoentong et al., 2017).

The TIDE, as a widely used algorithm for immune checkpoint blockade (ICB) responsiveness prediction (Thorsson et al., 2018), was applied to verify the prediction of clinical responsiveness to ICB, which evaluated various cell types including tumor-related fibroblasts, which are responsible for excessive extracellular matrix deposition, immunosuppressive cells such as the M2 subtype of tumor-associated macrophages, and myeloid-derived suppressor cells that suppressed the T cell infiltration inside the tumor microenvironment, and two distinct mechanisms involved in escaping immune surveillance, including the score determining the dysfunctionality of tumor-infiltrating cytotoxic T lymphocytes (CTLs) and the score showing the rejection of CTLs by immunosuppressive factors.

A total of 40 genes were involved in tryptophan metabolism. To determine the genetic alteration of tryptophan metabolism in GC, mutation frequency of cells was analyzed among the tryptophan metabolism-associated genes. Among the 437 tumor samples, 121 (27.69%) samples had tryptophan metabolism mutations. As shown in Figure 1A, AOX1 and OGDHL genes had the highest mutation frequencies, and no mutation was found in the WARS1 gene. To understand the effect of mutations on tryptophan metabolism-related genes, we analyzed the biological signaling pathways in wild-type (WT) and mutant (Mut) groups through GSEA enrichment analysis. It shows that tumor-associated pathways, including TNFA_SIGNALING_VIA_NFKB, P53_PATHWAY, and MYC_TARGET, were enriched in the mutant group (Sanchez-Vega et al., 2018) (Figure 1B). We then examined the somatic copy number variations of these tryptophan metabolism-associated genes in GC tumor samples and discovered a lower frequency of copy number variation (CNV) deletion or amplification (Figure 1C). In order to explore mRNA expression of CNV value in tumor tissue, the samples were distributed into different groups relying on CNV value, including increase and loss of CNV, as well as no obvious variation in CNV. Comparison of the expression of genes correlating with tryptophan metabolism between these groups showed that patients with CNV gain had a higher mRNA expression level compared to those with CNV loss (Figure 1D). To determine the expression of the tryptophan metabolism-associated genes between tumor tissue samples and adjacent normal tissues. As indicated in Figure 1E, most tryptophan metabolism genes were significantly differentially expressed, such as AANAT, AFMID, HADH, IDO1, IDO2, IL4I1, KMO, KYNU, MAOA, MAOB, TDO2, and WARS2.

The tryptophan metabolism score in the TCGA dataset was calculated by the ssGSEA, and then Pearson was used to estimate the relationship between the protein-encoding genes and the tryptophan metabolism score. A total, 30 prognostic genes associated with tryptophan metabolism score were screened. Figure 2A showed the correlation between the 30 genes and tryptophan metabolism scores. We classified patients based on the consensus clustering on 30 prognosis-correlated gene expression profiles and selected the optimal cluster number based on the CDF. With relatively stable clustering results shown in Figures 2B, C, we finally chose k = 2 to acquire two molecular subtypes (Figure 1D). We further performed the prognostic analysis of these two molecular subtypes. As shown in Figure 2E, we found that the overall survival of C1 was significantly better than that of C2. In addition, a sample of GC patients from the GSE66229 dataset showed similar results. This suggests that GC patients in C1 would have a better prognosis relative to the C2 subtype (Figure 2F). Meanwhile, we determined the tryptophan metabolism scores of each sample in the TCGA and GSE66229 datasets, which showed the score of C2 was higher compared to that of the C1 subtype with a good survival benefit (Figures 2G, H). Furthermore, differences in clinicopathological characteristics of TCGA molecular subtypes were analyzed. As indicated in Figure 3, we found significant differences between the two molecular subtypes in terms of age, disease stages, grade classification, and patient survival status.

As shown in Figure 4A, some immune cell types were significantly different among the subtypes. Compared to the C2 subtype, activated dendritic cells (DCs) and CD4 positive memory T cells were significantly upregulated while the M2-phenotype macrophages were suppressed in the C1 subtype. Meanwhile, ESTIMATE was applied to evaluate the level of immune cells inside tumor tissues. Figure 4B showed that the ImmuneScore of the C2 subtype was significantly higher compared to that of other groups, indicating that C2 has a high level of immune cell infiltration. In addition, ssGSEA demonstrated significant differences in most immune cell scores between different subtypes (Figure 4C). We further analyzed the responsiveness to immunotherapy between different TCGA cohort molecular subtypes. As shown in Figure 4D, compared to C1 subtype, the gene expression of immune checkpoints such as IDO1, IDO2, and CD274 were increased dramatically in the C2 subtype.

We did this by employing the TIDE algorithm in order to assess the potential response of tumor samples to immune checkpoint inhibitors. Higher TIDE scores represent a higher likelihood of immune escape in response to immune checkpoint therapy (Jiang et al., 2018). As shown in Figure 4E, in the TCGA cohort, the C2 subtype showed a much higher TIDE score than the C1 subtype with best prognosis, indicating that the C2 subtype had greater potential of immune escape. In addition, analysis on the response of different molecular subtypes to conventional chemotherapeutic drugs showed that C2 patients were more sensitive to these drugs including sunitinib, MG-132, saracatinib, dasatinib, and CGP-60474 (Figure 4F).

Further, differences in genomic alterations between two different molecular subtypes in the TCGA cohort were analyzed. Here, we obtained molecular information of TCGA-STAD collected from a pan-cancer study (Thorsson et al., 2018). The C1 subtype showed a higher aneuploidy score, fraction altered, homologous recombination defects, non-silent mutation rate, tumor mutation burden (TMB), and number of segments (Figure 5A). Gene mutation differences between different molecular subtypes were also studied. The top 10 genes were exhibited in Figure 5B, which showed significant differences of TTN, TP53, and MUC16 genes in mutation frequency between the two molecular subtypes. Similarly to our results, it has been reported that TTN, TP53 and MUC16 are the most significantly mutated driver genes in GC, which is closely related to the prognosis of cancer patients (Dong et al., 2022).

We then performed GSEA analysis which showed that DNA replication, spliceosome and base excision repair signaling pathways were highly enriched in the C1 subtype, while the C2 subtype had high enrichment of the phagocytosis, chemokine, and leukocyte transendothelial migration signaling pathways (Figure 6A). We also evaluated the 10 oncogenic pathways from the previous study (Sanchez-Vega et al., 2018). The differences showed statistical significance in the rest of the pathways except for the TP53 and NRF1 signaling pathways (Figure 6B). Compared to the C2 subtype, the oncogenic pathways, such as the Wnt, PI3K, and RAS pathways (Zhang et al., 2001; Zhang et al., 2019; BY et al., 2020), were significantly downregulated in the C1 subtype, which indicated an association between activated oncogenic pathways and tryptophan metabolism that may result in poor prognosis in C2.

As shown above, we identified two different molecular subtypes with distinct prognostic, immune, mutational, and pathway signatures. Next, we screened genes which are differentially expressed in C1 subtypes compared with C2 subtypes (|Log2FC| > 1; FDR <0.05), and a sum of 618 distinctly expressed genes were obtained as shown in the volcano plot of differential analysis (Figure 7A) in which 36 genes were significantly increased and 582 genes were decreased. We next utilized univariate cox analysis on the 618 differentially expressed genes and identified a total of 218 genes that showed the greater impact on prognosis, including 217 Risk genes and 1 Protective gene (Figure 7B). Subsequently, we compressed these 218 significant differentially expressed genes through lasso regression. As shown in Figures 7C, D, 10-fold cross-validation was used to construct the model and eight genes at lambda = 0.0356 were screened for further analysis: EFNA3, GPX3, RGS2, CXCR4, SGCE, ADH4, CST2, GPC3. The final 8-gene signature formula was as follows:

RiskScore=(0.024*EFNA3)+0.079*GPX3+0.102*RGS2+0.051*CXCR4+0.014*SGCE+0.076*ADH4+0.067*CST2+0.066*GPC3.

The risk score for each TCGA sample was separately calculated, followed by conducting a receiver operating characteristic (ROC) analysis on the prognostic classification of the RiskScore. The prediction classification efficiency from 1 year to 4 years was calculated to have an area under the time-dependent ROC curves (AUC) of 0.7, which validated the prediction capability of the model (Figure 8A). We then performed the zscore on the RiskScore. When the RiskScore was lower than zero, the samples were separated into the low-risk group, while those with a RiskScore higher than zero were in the high-risk group. As shown in Figure 8B, the low-risk group showed prolonged survival time, indicating the good performance of this prognostic model. Additionally, GSE66229 dataset was utilized to test the robustness and validate the risk model constructed by these eight genes. As shown in Figures 8C, D, similar results were observed, which indicated the excellent predictive ability of this model.

We analyzed the differences in RiskScore between Gender, Age, TNM grades, Stage clinical grades, and Grade grades in the TCGA-STAD dataset to explore the relationship between RiskScore and clinical characteristics of GC. With the increase in the clinical grade, the RiskScore also increased (Figure 9A). Moreover, we analyzed the difference in RiskScores among different molecular subtypes. The RiskScore of the C2 subtype with worse prognosis was obviously higher compared to that of the C1 molecular subtype with the best prognosis. Clinicopathological characteristics between the RiskScore groups in the TCGA-STAD cohort were analyzed, and the high-risk group was found to have a higher clinical grade (Figure 9B), which was consistent with previous results.

Next, the enrichment of 22 immune cells in the high and low RiskScore groups was analyzed between RiskScore groups. As shown in Figure 10A, compared to the high-risk group, resting NK cells, activated CD4 T cells, and activated DCs were increased dramatically while M2 macrophages were substantially decreased in the low-risk group. ESTIMATE was also utilized to evaluate the level of immune cells in the tumor tissues. As shown in Figure 10B, the ImmuneScore in the high-risk group is much higher. Although it showed a higher level of immune cells in the high-risk group, the immunosuppressive cells such as MDSC and M2-phenotype macrophages also increased significantly, which might result in a poor prognosis. Further, we also used the ssGSEA function to analyze the scores of 28 types of immune cells, and most of the immune cell scores were significantly different between high- and low-risk groups (Figure 10C). The correlation between immune cells and RiskScore was further evaluated. As shown in Figure 10D, the RiskScore positively correlated with most immune cells, especially CD4 T cells, DCs, NK cells, and MDSCs, which could support the prediction of prognoses. This suggests that the infiltration of these immune cells increases as the risk score rises.

In Figures 10E–G, the most signaling pathways were positively correlated with the RiskScore of the samples and a significant positive correlation could be found between RiskScore and tryptophan metabolism ssGSEA scores, indicating that the activated tryptophan metabolism may induce poor prognosis.

We further explored the expression level of immune checkpoints between RiskScore groups. As shown in Figure 11A, some immune checkpoint genes, including NRP1, CD200, and CTLA4, were significantly downregulated in low-risk groups. We further analyzed the difference in immunotherapy among different RiskScore groups. TIDE software was applied to evaluate the immunotherapy response in the high and low RiskScore groups. We can find that in the TCGA cohort, the high-risk group showed a higher TIDE score (Figure 11B), indicating that the high-risk group showed greater potential of immune escape and less sensitive to immunotherapy. We further analyzed the relationship between RiskScore and TIDE score. Figure 11C showed a positive correlation between RiskScore and TIDE, IFNG, Exclusion, and Exclusion scores and a significant negative correlation with MDSC. In addition, we also analyzed the response of RiskScore groups to chemotherapeutic drugs in the TCGA cohort and found that the high RiskScore group was more sensitive to chemotherapeutic drugs including MG-132, dasatinib, CGP-60474, WH-4-023, and CMK (Figure 11D).

Combining RiskScore and clinicopathological features to optimize prognostic model and survival prediction.

Here, patient age, gender, TNM stage, Stage clinical grade, Grade, and RiskScore in the TCGA-STAD cohort were used to develop a decision tree. Only age, T stage, and RiskScore remained in the decision tree, where RiskScore was the most effective parameter (Figure 12A). Figure 12B showed a significant overall survival differences among the four risk subgroups. C2, C3, and C4 were all high-risk patients (Figure 12C). A significantly decreased survival benefit was found in the C2, C3, and C4 subgroups (Figure 12D). RiskScore, age, and Stage were significant prognostic factors, as confirmed by Univariate and multivariate Cox regression analysis on the clinical characteristics and RiskScore (Figures 12E, F). A nomogram combining the RiskScore with other clinicopathological characteristics was developed for the risk assessment and survival probability evaluation for patients. On survival rate prediction, the RiskScore showed the greatest impact. The accuracy of the prediction model was evaluated by the calibration curve (Figure 12G). As displayed in Figure 12H, the predicted calibration curve at 1 and 3 years(s) almost overlapped with the standard curve, indicating a strong prediction of the nomogram. Further decision curve analysis (DCA) demonstrated that the Riskscore and nomogram showed noticeably greater benefits than extreme curves. Additionally, both nomogram and RiskScore demonstrated the most powerful survival predictors compared to other clinicopathological features (Figures 12I, J).