Comprehensive analysis of tertiary lymphoid structures-related genes for prognostic prediction, molecular subtypes and immune infiltration in gastric cancer

Download of GC dataset and acquisition of TLS-RGs

The RNA expression data, clinical characteristics information, and somatic mutation data about GC were downloaded from the The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The GC-related dataset GSE84437 with complete clinical information and sufficient number of samples was located and downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Those samples without complete information were excluded, and the gene ID-transformed datasets were batch corrected and merged by the “sva” package and the “Combat” package, and the merged datasets were used for subsequent analysis. 39 TLS-RGs were obtained from previous studies [23]. The clinical information and data involved in this study were obtained from public databases. Therefore, written informed consent from patients and approval from the ethics committee were not required for this study.

Differential expression analysis and mutation analysis

To identify DEGs between normal and GC samples or between different subgroups, differential analysis was performed by the “limma” package based on |foldchange| ≥ 1 and P < 0.05. Mutation data were analyzed to obtain sample tumor mutational burden (TMB) scores. Copy number variations (CNV) frequencies of the samples were obtained from UCSC Xena (https://xena.ucsc.edu/) and the gain and loss of CNV frequencies were analyzed.

Consensus clustering analysis based on correlated genes

Consensus clustering analysis is often used to target subgroups of one trait, which in turn allows understanding the differences between different subgroups. In our study, we performed consensus clustering analysis on TLS-RGs and co-DEGs by the “ConsensusClusterPlus” package. Kaplan-Meier (KM) curves can describe the differences in survival between different subgroups and can be performed with the “survival” package.

Enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis provides an understanding of the signaling pathways involved at the molecular level. Gene ontology (GO) enrichment analysis enables the qualification and description of gene and protein functions, including biological process (BP), cellular component (CC) and molecular function (MF). Gene set variation analysis (GSVA) enables differential analysis at the level of signaling pathways. In our study, they can be implemented by the “clusterProfiler” package and the “GSVA” package.

Immune landscape assessment

In order to understand the immune landscape of each patient, we apply multiple algorithms to assess the immune cell infiltration and immune score. Single-sample Gene Set Enrichment Analysis (ssGSEA) is able to analyze the infiltration of 23 types of immune cells in each sample, which can be implemented through the “ssGSEA” package. The ESTIMATE package provides ESTIMATE scores, immune scores and tumor purity for each sample.

Principal component analysis (PCA)

PCA analysis enables the assessment of the classification of different subgroups, and to a certain extent can reflect the differences between subgroups. In addition, TLS scores can be obtained by PCA analysis of TLS-RGs.

Construction and verification of nomograms

Nomograms are widely used for prognostic prediction of tumors. In our study, nomograms based on clinical characteristics and TLS scores were constructed to predict 1-, 3-, and 5-OS of GC patients. Concordance index (C-index) and decision curves were able to assess the predictive accuracy of nomograms. Univariate and multifactorial regression analysis were used for the identification of prognostic marker genes.

Immunotherapy response and drug sensitivity analysis

The tumor immune dysfunction and exclusion (TIDE) score allows assessment of the likelihood of tumor immune escape based on gene expression data. We obtained TIDE scores for each sample through the TIDE database (http://tide.dfci.harvard.edu/login/). The Cancer Imaging Archive (TCIA) results can reflect the outcome of patients treated with anti-PD-1 or anti-CTAL4. Microsatellite instability (MSI) is considered one of the important indicators of tumor immune response, so we further analyzed the TCIA results to obtain MSI results for each sample. Drug sensitivity analysis is performed through the “pRRophetic” package, where the half maximal inhibitory concentration (IC50) is used as an indicator of drug sensitivity.

Statistical analysis

The data processing, analysis, and visualization involved in this study were performed by R software (version 4.1.2). Wilcoxon rank sum test and analysis of variance (ANOVA) were used for the analysis of differences between two and three groups, respectively. p < 0.05 was considered a statistical difference between the comparison groups.

Data availability statement

The datasets in this study were obtained from the TCGA database (https://portal.gdc.cancer.gov/) or the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Raw data and original codes are included in the article or in the supplementary material and further inquiries can be addressed to the respective authors.

Expression and mutation analysis of TLS-RGs in GC

According to the TCGA dataset, 39 TLS-RGs were obtained and the expression differences of these 39 TLS-RGs in normal tissue and GC tumor tissue were analyzed (Supplementary Table 1). The results are shown in Figure 1A, a total of 19 TLS-RGs were considered to have significantly altered expression in GC, 14 of which were up-regulated and 5 were down-regulated. The PPI network results indicated a significant association and interplay of TLS-RGs in GC (Figure 1B). TLS-RGs were also mutated in GC (11.55%), notably, IL1R2 (3%), STAT5A (2%), and CD4 (2%) had the highest mutation frequency in GC (Figure 1C). CNV frequency results showed that there were significant copy number changes in TLS-RGs in GC and more genes lost their copy number (Figure 1D). These results confirmed that the expression and mutation frequency of TLS-RGs were altered in GC indicating that TLS-RGs were associated with GC.

Figure 1. Analysis result of TLS-RGs in expression and mutational landscape. (A) Result of DE-TLS-RGs expression in normal and GC tissues. (B) PPI network analysis of TLS-RGs. (C) Somatic mutational landscape of TLS-RGs. (D) CNV frequencies of TLS-RGs in GC. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Identification and immune infiltration of TLS subgroups in GC

To clarify the role of TLS-RGs in GC, we performed unsupervised clustering analysis of GC patients based on the expression of 39 TLS-RGs. Good results were obtained when the number of subgroups was 3 (Figure 2A). PCA results showed the separation of the three clusters from each other (Figure 2B). The KM curves showed that the survival curves were different among three clusters, with cluster C being the best and cluster A or cluster B being worse (Figure 2C). The expression heat map of TLS-RGs in three clusters is shown in Figure 2D, and it should be noted that the expression of TLS-RGs in cluster B is the most specific. To understand the reasons for the survival differences between these three clusters, differential analysis was performed for two of these three clusters and the biological signals involved in the differentially expressed genes were identified by enrichment analysis (Figure 2D–2G). The results showed that immune-related pathways were more enriched in cluster A and cluster C, and cluster C was higher than cluster A (Figure 2D–2G). This suggests that immunity is involved in the regulatory of TLS-RGs in GC.

Figure 2. Identification of TLS subgroups in GC. (A) Results of unsupervised clustering analysis. (B) PCA results of three clusters. (C) KM curves of three clusters. (D) Heat map of TLS-RGs expression in three clusters. (E) Enrichment analysis of differential genes between cluster A and cluster B. (F) Enrichment analysis of differential genes between cluster A and cluster B. (G) Enrichment analysis of differential genes between cluster A and cluster B.

To investigate the relationship between the altered immune environment and TLS subgroups, we applied multiple immune algorithms to assess the immune landscape of the samples. ssGSEA results showed significant differences in the abundance of 23 immune cell in the 3 clusters (Figure 3A). It is worth noting that the infiltration of Neutrophils and Type 17 T helper cells in TLS cluster A or TLS cluster C is different from other immune cells. And cluster C had a higher ESTIMATE score and higher immune score but lower tumor purity compared to cluster B (Figure 3B–3E). TIDE results also pointed to a better immune response in cluster C (Figure 3F). In addition, we performed immunotherapy analysis on these 3 clusters and found that cluster C showed a greater response to anti-PD-1/anti-PD-1+ anti-CTLA4 treatment (Figure 3G–3I). In conclusion, the better survival of cluster C is closely related to its own immune environment.

Figure 3. Immune landscape analysis results of TLS subgroups in GC. (A) Results of 23 immune cell abundance analyses by ssGSEA. (B–E) ESTIMATE algorithm results. (F) TIDE score results. (G–I). Results of immunotherapy analysis. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Subgroup identification based on co-differential expressed genes

To garner subgroups marked by more conspicuous distinctions in prognostic outcomes, a reiteration of consensus clustering analysis was executed upon a cohort of GC patients. This secondary consensus clustering analysis was undertaken with the intention of amplifying the repertoire of DEGs, with a specific focus on identifying those that bear relevance to prognosis. The overarching objective herein was to refine the construction of a more robust prognostic risk model. Subsequently, a compendium of common DEGs was discerned within the context of TLS clusters through the adept utilization of the “limma” computational package. Thereafter, a meticulous enrichment analysis was conducted, yielding pertinent insights into the functional implications of these genes (Supplementary Table 2). These co-DEGs were found to be mainly enriched in the biological functions of the immune process (Figure 4A, 4B), suggesting that TLS-RGs-mediated alterations in the immune process are critical for the regulation of GC. Further, the role of co-DEGs in GC was further discussed by dividing patients into different genetic subgroups based on the expression of co-DEGs using a consensus clustering approach. The results showed that the optimal grouping results could be obtained when the number of subgroups was 2. The KM curves indicated that there was a significant difference in prognosis between the two gene clusters, with gene cluster B having a better prognosis than gene cluster A (Figure 4C), and the difference was greater than that of the TLS cluster (P_{gene cluster} = 0.006 < P_{TLS cluster} = 0.018). Expression heat map of co-DEGs is shown in Figure 4D, and it is clear that their expression is altered in these two gene clusters. We also analyzed the expression of TLS-RGs between these two gene clusters and the results showed that the expression of most TLS-RGs differed and that the expression of TLS-RGs was higher in gene cluster A than in gene cluster B (Figure 4E). This indicates that TLS-RGs in gene clusters appear significantly different.

Figure 4. Results of subgroup identification based on co-differentially expressed genes. (A) KEGG results of co-DEGs. (B) GO results of co-DEGs. (C) KM curves of gene clusters. (D) Expression heat map in gene clusters of co-DEGs. (E) Expression differences of TLS-RGs in gene clusters. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Assessment and classification of TLS scores

In pursuit of enhanced precision in forecasting both OS and treatment responsiveness among GC patients, a TLS scoring framework was meticulously crafted, leveraging the outcomes of Principal Component Analysis (PCA) applied to gene clusters. This strategic development seeks to furnish a robust predictive tool capable of furnishing more refined prognostic and treatment response insights for the benefit of GC patient care. As shown in Figure 5A, 5B, TLS scores differed in either gene clusters or TLS clusters. Note that TLS scores are highest in gene cluster A and TLS cluster A, which indicates that TLS scores are associated with immune regulation. Based on the TLS score, the optimal threshold was selected and the patients were divided into a low TLS score group and a high TLS score group. The distribution of GC patients among TLS clusters, gene clusters and TLS score clusters is shown in Figure 5C. The scatter plot showed that GC patients in the low TLS score group had a longer survival time (Figure 5D). GC patients in the low TLS score group tended to have better OS (Figure 5E). PCA results showed that patients with different TLS scores were separated from each other (Figure 5F).

Figure 5. Assessment and classification results of TLS scores. (A) Differences in TLS scores among gene clusters. (B) Differences in TLS scores among TLS clusters. (C) Distribution of TLS scores among gene clusters, TLS clusters and TLS score subgroups. (D) Scatter plot of GC patients with different TLS scores. (E) KM curve results of GC patients with different TLS scores. (F) PCA results of GC patients with different TLS scores. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Construction of nomogram based on clinical characteristics and TLS scores

In order to accurately predict the prognosis of GC patients, we constructed a prognostic nomogram based on clinical characteristics and TLS scores (Figure 6A). Concordance index (C-index) results showed that the nomogram had good prognostic ability (Figure 6B). The decision curves showed that the nomogram was more accurate in predicting OS of GC patients than other clinical characteristics or TLS scores (Figure 6C). We further assessed the predictive independence of clinical characteristics and TLS scores by univariate and multifactorial Cox regression. The results showed that age (hazard ratio (HR) = 1.026, p < 0.001), T (HR = 1.255, p = 0.001), N (HR = 1.549, p < 0.001) and TLS score (HR = 1.180, p < 0.001) were associated with OS in GC patients, whereas age (HR = 1.031, p < 0.001), N (HR = 1.496, p < 0.001) and TLS score (HR = 1.201, p < 0.001) were able to serve as independent prognostic factors for GC patients (Figure 6D, 6E).

Figure 6. Results of nomogram based on clinical characteristics and TLS scores. (A) Prognostic nomogram based on clinical characteristics and TLS scores. (B) C-index results of the prognostic nomogram. Result is close to 45° indicating that the nomogram prediction is accurate. (C) Decision curve results of the prognostic nomogram. (D) Results of univariate Cox regression analysis of clinical characteristics and TLS scores. (E) Results of multifactorial Cox regression analysis of clinical characteristics and TLS scores. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Immune landscape and immunotherapy in different TLS scoring groups

The abundance of 23 immune cells in the TLS score subgroups was analyzed by ssGSEA, and the results showed that there were significant differences in the distribution of the other 20 immune cells in different TLS scoring groups, except for CD56dim natural killer cell, Type 17 T helper cell and Type 2 T helper cell (Figure 7A). As shown in Figure 7B, 7C, SD/PD group had higher TLS scores. ESTIMATE results indicated that the high TLS score group had higher ESTIMATE score, immune score as well as lower tumor purity (Figure 7D–7G). Notably, there was more potential for immune escape in the high TLS score group due to the higher TIDE score in the high TLS score group compared to the low TLS score group (Figure 7H). In the analysis of ICP, we found that the expression of ICP was higher in the high TLS score group than in the low TLS score group, with the exception of HHLA2 (Figure 7I). In conclusion, we found significant immunological differences between subgroups of TLS score groups, with patients in the low TLS score group receiving immunotherapy with better outcomes than the high TLS score group.

Figure 7. Results of immune landscape and immunotherapy in different TLS score groups. (A) Abundance results of 23 immune cells by ssGSEA in different TLS score groups. (B) Results of anti-PD-L1 immunotherapy responses in different TLS score groups; Abbreviations: CR: complete remission; PR: partial remission; SD: stable disease; PD: progressive disease. (C) Distribution of TLS scores in different immunotherapy responses. (D–G) ESTIMATE results in different TLS score groups. (H) TIDE results in different TLS score groups. (I) Results of ICP in different TLS score groups. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.

Analysis of tumor mutation landscape and drug sensitivity in different TLS score groups

Microsatellite Instability (MSI) is considered as one of the important indicators of tumor immune response. In our study, patients in the low TLS score group were associated with MSI-L and MSI-H, while patients in the high TLS score group were associated with MISS (Figure 8A, 8B). Tumor mutation burden (TMB) can help predict patient response to immunotherapy. In our study, TMB was negatively correlated with TLS score (Figure 8C). Notably, TMB was higher in the low TLS score group than in the high TLS score group (Figure 8D), indicating that patients in the low TLS score group had a better immunotherapy response. Subsequently, we performed survival analysis on patients with different TMB, and the results showed that patients with high TMB tended to have better OS (Figure 8E). Further combined with the TLS score, patients with low TLS score accompanied by high TMB had the best prognosis (Figure 8F). The results of somatic mutation analysis for TLS score showed that the incidence of mutations in TTM, TP53, MUC16 and ARID1A were all higher than 20%, and it should be noted that the incidence of mutations in these genes was higher in low TLS score than in high TLS score (Figure 8G, 8H). These results suggest that TLS-RGs can influence the response to tumor immunotherapy. Drug sensitivity can reflect the patient’s response to drug therapy, and drug sensitivity analysis for TLS score was performed in our experiments. The results showed that Rapamycin, Dasatinib, Sunitinib and Saracatinib had a higher IC50 in the low TLS score group of patients, in contrast to AKT inhibitor VII (Figure 8I–8M). In conclusion, TLS-RGs was also associated with drug sensitivity in GC patients.

Figure 8. Results of tumor mutation landscape and drug sensitivity analysis in different TLS score groups. (A) MSI results in different TLS score groups. (B) Correlation results of MSI and TLS score. (C) Correlation results of TMB and TLS score. (D) Distribution of TMB in different TLS score groups. (E) KM curves of patients with different TMB. (F) KM curves of patients with different TMB and different TLS scores. (G) Somatic mutation analysis results in low TLS scores. (H) Somatic mutation analysis results in high TLS scores. (I–M) Drug sensitivity analysis results in different TLS score groups. ^*p < 0.05; ^**p < 0.01; ^***p < 0.001.