Integrated analysis of TCGA data identifies endoplasmic reticulum stress-related lncRNA signature in stomach adenocarcinoma

Rationale for utilizing TCGA data analysis methods

This study utilized data from TCGA due to its comprehensive and well-validated nature, rendering it an ideal resource for cancer research. TCGA offers extensive genomic data, including transcriptomic, mutational, and clinical datasets, essential for identifying potential biomarkers and therapeutic targets in cancer research. Specifically, our focus on lncRNAs in the context of ERS in STAD necessitated a database with broad and detailed genomic information. The integrative bioinformatics approach adopted in this study, which combined gene expression profiles with clinical data, facilitated a comprehensive understanding of the molecular mechanisms underlying STAD. This method also aided in identifying potential prognostic lncRNAs and therapeutic targets. This approach also facilitated the development and validation of a prognostic model based on these lncRNAs, emphasizing the utility of TCGA data in clinical applications for STAD.

Data preprocessing

The study acquired 448 transcriptome expression datasets, 434 tumor mutation datasets, and 443 clinical datasets from the TCGA database, which is accessible at https://tcga-data.nci.nih.gov/tcga/. This clinical data included gender, age, and survival details. Using the Perl programming language, this data was consolidated into a matrix file. The research identified 295 ERS-related genes (Table S1), by querying the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp ) using ERS-specific keywords (Table S2).

Evaluation of differentially expressed ERS-related lncRNAs

The “limma” R package was utilized to identify 1364 ERS-associated lncRNAs, applying cutoff values of |r|>0.6 and p<0.001 using Pearson’s correlation coefficient. Differentially expressed ERS-related lncRNAs were further selected by applying a threshold of adjusted p<0.05 and |log₂ (fold change, FC)|>1, in a comparison between STAD samples and non-carcinoma samples. The ‘limma’ R package was employed to identify 1364 ERS-associated lncRNAs, using cut-off values of |r|>0.6 and p<0.001 as determined by Pearson’s correlation coefficient. The selection of differentially expressed ERS-related lncRNAs was further refined by applying a threshold of adjusted p<0.05 and |log₂ (fold change, FC)|>1, in comparisons between STAD samples and non-carcinoma samples.

Development of a prognostic model for ERS-related lncRNAs

A prognostic model for ERS-related lncRNAs was constructed in the following manner: initially, survival-related lncRNAs were identified from the pool of differentially expressed ERS-related lncRNAs, employing univariate Cox regression analysis and taking into account clinical data from STAD patients. Subsequently, the identified prognostic lncRNAs were subjected to LASSO-Cox regression analysis to mitigate overfitting. This process involved one thousand cycles and was conducted with a significance level set at p<0.05. The risk scores were then calculated using the following formula:

where Σ represents the sum of all included lncRNAs, coef _i represents coefficient index, and expr _i represents expression level.

Verification of the risk model

Initially, the entire TCGA cohort underwent random division into training and testing groups at a 1:1 ratio. Observations were allocated to either the low- or high-risk group based on their median risk scores, and the predictive efficacy of the risk model was assessed using the Chi-square test. Survival curves were generated according to the number of events in each group, utilizing the “survival” and “survminer” packages. The reliability of the risk scores was evaluated using the testing set. Furthermore, univariate Cox regressions, LASSO, and multivariate Cox regressions were conducted using STAD patients’ clinical data, including ages, gender, and stage. The prognostic effectiveness was assessed using time-dependent Receiver Operating Characteristic (ROC) curve analysis and Harrell’s concordance index (C-index) with the “survival”, “survminer”, and “timeROC” R packages.

Nomogram establishment

A nomogram was constructed to evaluate the reliability of the prediction outcomes, incorporating factors such as gender, age, grade, risk score, and tumor stage. The “rms” R package was utilized for this purpose, and Hosmer-Lemeshow goodness-of-fit test curves were created as line graphs for 1-, 2-, and 3-year overall survival (OS). The predictive performance of the nomogram was assessed using the C-index and time-dependent ROC curve, computed with the “survival” and “timeROC” R packages.

Gene ontology and Kyoto encyclopedia of genes and genomes enrichment analysis

In this study, an enrichment analysis was conducted to elucidate various signaling pathways and their biological effects, based on the findings obtained. This analysis was facilitated using the ‘clusterProfiler’ package, an R-based tool that offers a streamlined way to perform statistical analysis and visualization of functional profiles for genes and gene clusters. The package can be accessed via the Comprehensive R Archive Network at CRAN URL for clusterProfiler. For the analysis, both the p-value and q-value thresholds were set at 0.05. Statistically significant enrichment was defined as cases where both the False Discovery Rate (FDR) and p-value were less than 0.05. For reference to the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, please visit:

GO: http://geneontology.org/; KEGG: https://www.genome.jp/kegg/.

Gene set enrichment analysis

Several R packages were employed to further explore changes in gene expression among different risk groups in the study. These include: ‘org.Hs.e.g. db’: An R package providing mappings between Entrez Gene IDs and various gene identifiers. This package is crucial for gene annotation in our analysis. Available at Bioconductor URL for org.Hs.e.g. db. ‘limma’: utilized for analyzing gene expression data from microarray or RNA-seq technologies. It uses linear models for differential expression in complex experiments. Available at Bioconductor URL for limma. ‘enrichplot’: this package assists in visualizing functional enrichment results, making it easier to interpret gene set enrichment analysis. Available at Bioconductor URL for enrichplot. ‘clusterProfiler’: as previously mentioned, this package was used for both statistical analysis and visualization of functional profiles. For Gene Set Enrichment Analysis (GSEA), refer to the official GSEA website: http://www.gsea-msigdb.org/gsea/index.jsp.

Tumor mutation burden analysis

The R package “maftools” was utilized to investigate the Mutation Annotation Format (MAF) data on somatic mutations, and the tumor mutation burden (TMB) score for each STAD sample was calculated using a specified methodology:

Immune cell infiltration analysis and immune checkpoints analysis

The R package ‘ESTIMATE’ assesses stromal cell and immune cell scores, enabling the evaluation of tumor purity. Immune cell infiltration in tumor samples was analyzed using single-sample GSEA for 29 different immune cell types. Furthermore, the ‘ggpubr’ R package was used to evaluate immune checkpoint molecules and visualize them through boxplots in R.

Drug sensitivity analysis

Chemosensitivity analysis was conducted to explore potential differences in the sensitivity to chemotherapeutic drugs between low- and high-risk groups among STAD patients. The analysis of drug sensitivity was performed using the ‘limma’ package with a p-value threshold of 0.001, and diagrams illustrating the results were created using ‘ggplot2,’ with the ‘ggpubr’ package aiding in the visualization.

Consensus clustering

Consensus clustering was conducted on patients from TCGA database using the Consensus Cluster Plus R package, with a focus on the five selected ERS-related lncRNAs. The optimal number of clusters was determined through the analysis of the cumulative distribution function and consensus matrix.

Statistical analysis

Pearson correlation analysis was determined to investigate the association between ERS-related genes and lncRNAs using R version 4.2.2. The Wilcoxon signed-rank test was applied to assess immune infiltrating cells and p-value less than 0.05 was considered significant.

Identification of LncRNAs co-expressed with RES in STAD

As depicted in Figure 1, the analysis was initiated by acquiring transcriptome data from 412 STAD tissues and 36 normal tissues, followed by conversion of Ensembl ID into official gene names and categorization of transcripts as either mRNAs or lncRNAs. Subsequently, 295 ERS-related genes were identified based on their mRNA expression. Potential ERS-related lncRNAs were identified by filtering out 958 lncRNA genes that exhibited a correlation coefficient greater than 0.6 and a significance level of p<0.001.

Figure 1:

The workflow of this study.

Construction and the validation the prognostic risk signature

The WilcoxTest was employed to analyze the differential expression of lncRNAs co-expressed with ERS genes between tumor and normal groups in TCGA. In the tumor group, 523 genes exhibited upregulation, while 58 genes showed downregulation when compared to the normal group (Figure 2A and B). To construct a prognostic risk score model, clinical data of TCGA-STAD was downloaded and merged with the above differential expression data. The samples were divided into two groups: a training group for constructing the model and a testing group for validating the model. Univariate Cox regression analysis identified 32 lncRNAs (AL352979.4, AL137798.1, LINC02533, AC079467.1, LINC01354, AL022326.1, SNHG14, AC084880.3, AC100849.1, LINC01990, AC110491.1, AL161457.1, ADAMTS9-AS2, LINC01537, AC012055.1, AP006248.3, AC027243.3, RBMS3-AS3, AL033543.1, LINC01235, AC008808.1, GAS1RR, MRGPRF-AS1, LINC00571, AC073548.2, RMDN2-AS1, LINC02073, DIRC3, ARHGEF26-AS1, AL583829.1, AC109927.1 and CFAP61-AS1) that were associated with prognosis (Figure 2C and D). Using LASSO regression analysis, 12 lncRNAs (AL352979.4, LINC02533, AC110491.1, AL161457.1, AC012055.1, AP006248.3, AC027243.3, LINC01235, AC008808.1, LINC00571, LINC02073 and CFAP61-AS1) were selected (Figure 2E and F). Finally, five prognostic lncRNAs were identified using multivariate Cox regression analysis. The regression coefficients of the five lncRNAs were obtained using the R-package “Survival”. The risk score was calculated based on the screened significant lncRNAs using the following formula: risk score = (1.2059 × AC012055.1) + (0.32184 × LINC01235) + (−1.3766 × LINC00571) + (1.0984 × LINC02073) + (0.3024 × CFAP61-AS1).

Figure 2:

Identifying differentially expressed ERS-related lncRNAs in STAD. (A) Volcano plot depicting differential expression of ERS-related lncRNAs in STAD tissues compared to those in normal tissues; (B) heatmap illustrating the expression levels of ERS-associated differentially expressed lncRNAs; (C) heatmap displaying the expression levels of the 32 prognostic lncRNAs in all patients; (D) univariate Cox regression analysis of the 32 lncRNAs associated with prognosis; (E) distribution of LASSO coefficients for the ERS-related lncRNAs; (F) cross-validation procedure for optimizing LASSO regression parameters.

The prognostic power of the model was evaluated through several analyses, including Kaplan-Meier (KM) survival curve, independent prognostic, receiver operating characteristic (ROC) curve, and nomogram. The KM analysis showed that low-risk group patients had a better prognosis compared to the high-risk group (p<0.001), with an increasing death rate as the risk score increased (Figure 3A). Furthermore, a heatmap showed that LINC00571 was a protective lncRNA, while AC012055.1, LINC01235, LINC02073, and CFAP61-AS1 were risk lncRNAs (Figure 3B). ROC analysis revealed that the risk score exhibited superior predictive performance compared to other clinical indicators (Figure 3C), achieving an area under the ROC curve (AUC) of 0.656 for 1-year survival, 0.697 for 3-year survival, and 0.678 for 5-year survival (Figure 3D). Univariate and multivariate analyses of independent prognostic factors revealed that the risk score obtained from the prognostic prediction model was an independent prognostic factor (Figure 3E and F). A nomogram was developed to enhance surgical prediction for STAD patients. The calibration curves for 1-, 2-, and 3-year OS demonstrated strong agreement between the nomogram’s predictions and actual observations (Figure 3G and H). The C-index for this risk model was 0.695 (95 % CI: 0.671–0.719). These findings indicate that the newly constructed risk model is a more accurate predictor of prognosis (Figure S1A–F).

Figure 3:

Prognostic value of the five ERS-related lncRNAs model. (A) KM survival curves comparing the OS of patients in these two groups in the training cohort; (B) the reference risk curve and heat map depicting the expression levels of the five ERS-related lncRNAs in the training cohort; (C) comparison of the predictive accuracy of the clinical characteristics and the risk model using ROC curves; (D) time-dependent ROC curve indicating the overall survival rates at 1, 2, and 3 years; (E) univariate Cox regression analysis of clinical features and the risk model for prognosis; (F) multivariate Cox regression analysis of clinical features and the risk model for prognosis; (G) nomogram integrating the risk score and other clinical factors to predict overall survival; (H) calibration curves for 1-, 2-, and 3-year survival predictions demonstrating the high agreement between predictions and observations.

Identify the differentially expressed genes (DEGs) and functional enrichment analysis between high- and low-risk groups

Differential expression analysis was conducted to investigate the biological functions and pathways associated with the risk score. This analysis compared the high- and low-risk groups, applying a cutoff value of |logFC|>1 and an FDR<0.05. Our analysis revealed 696 up-regulated genes and 54 down-regulated genes (Figure 4A and B). Subsequently, GO, KEGG enrichment analysis, and GSEA were conducted to identify enriched biological functions and pathways. Our GO enrichment analysis revealed that the DEGs were significantly enriched in collagen-containing extracellular matrix and receptor ligand activity (Figure 4C and D). KEGG analysis showed that the DEGs were significantly enriched in neuroactive ligand-receptor interaction and dilated cardiomyopathy (Figure 4E and F). Furthermore, GSEA analysis revealed that amino, nucleotide, and lipid metabolism were enriched in the low-risk group. In contrast, the high-risk group was enriched in vascular smooth muscle contraction, dilated cardiomyopathy, hypertrophic cardiomyopathy, focal adhesion, and neuroactive ligand receptor interaction (Figure 4G and H).

Figure 4:

Identification of DEGs and functional enrichment analysis. (A) Volcano plot displaying the DEGs between the high- and low-risk groups; (B) heatmap of the top 50 DEGs; (C) GO analysis of DEGs enriched in biological processes; (D) GO analysis of DEGs enriched in molecular functions; (E) KEGG analysis of DEGs; (F) KEGG pathway analysis of DEGs; (G) GSEA showing enriched pathways in low-risk group; (H) GSEA showing enriched pathways in high-risk group.

Analysis of mutation pattern between low- and high-risk groups

The relationship between the risk score and mutation patterns was investigated by conducting TMB analysis on 434 STAD samples obtained from TCGA. Significantly mutated genes associated with ERS were identified using Fisher’s exact test with a p-value threshold of 0.01. It was observed that there were diverse somatic mutation profiles in both the low- and high-risk groups. The genes with the most mutations were TTN, TP53, MUC16, LRP1B, ARID1A, CSMD3, SYNE1, FAT4, FLG and PCLO, with TTN ranked first (Figure 5A and B). A significant negative correlation was observed between the risk score and TMB (R=−0.12, p=0.013) (Figure 5C and D). Additionally, patients with higher TMB scores displayed improved overall survival in the STAD cohort. Within the high-risk category, individuals with low TMB scores exhibited the poorest overall survival (Figure 5E). Moreover, patients in the low-risk group with higher TMB scores (Figure 5F).

Figure 5:

Illustrates the TMB analysis between the low- and high-risk groups. (A, B) The somatic mutations in each group; (C) the relationship between the risk score and TMB; (D) TMB distribution in STAD cases categorized by risk groups; (E, F) the correlation between TMB, which includes different risk scores, and survival probability.

Drug sensitivity analysis and immune cell infiltration analysis

Our analysis of the tumor microenvironment has revealed significant differences in immune and stromal scores between high- and low-risk groups (Figure 6A–C; p<0.001). Correlations between immune cell populations and risk scores were also investigated, with T cells and NK cells demonstrating significant negative correlations with risk scores, while other immune cell populations showed significant positive correlations (Figure 6D). Additionally, survival analyses performed on samples divided by immune cell content showed that cancer-associated fibroblast_MCPCOUNTER, Endothelial cell_XCELL, and Macrophage_TIMER had a higher prognosis in the low-risk group (p<0.05), while T cell follicular helper_CIBERSORT and uncharacterized cell_EPIC had a higher prognosis in the high-risk group (p<0.05) (Figure 6E–I). Our analysis of immune cell and function scores using ssGSEA further demonstrated significant differences between low- and high-risk groups (Figure 6J and H). Immune function scores were also significantly higher in the high-risk group (Figure 6L; p<0.05). Specifically, HLA (-DRA, -DOA, -DPB1, -DPA1, -DPB2, -DQA1, and -DRB6) expressions were significantly higher in the high-risk group than in the low-risk group (Figure 6L; p<0.05). These results highlight the important role of the tumor microenvironment and immune cell regulation in STAD and provide valuable insights for future research in this area.

Figure 6:

Presents the results of the immune cell infiltration analysis. (A–C) Box plots reveal a statistically significant difference in immune score, stromal score, and estimat score; (D) correlation analysis indicates a significant relationship between immune cell infiltration and risk score; (E–I) the correlation between immune cell infiltration and survival probability is demonstrated through KM curves; (J) ssGSEA enrichment analysis suggests the involvement of various immune cells in the tumor microenvironment; (K) functional enrichment analysis identifies immune-related pathways that are significantly enriched; (L) box plots show differential expression of immune checkpoint genes in the high-risk group compared to the low-risk group. *p<0.05, **p<0.01, ***p<0.001 vs. low-risk group.

Our tumor microenvironment analysis revealed that individuals in the high-risk group displayed a pronounced immune-hot phenotype, whereas those in the low-risk group exhibited immune-cold characteristics. These findings suggest that the high-risk group may benefit more from immune checkpoint inhibitor (ICI) treatment. Drug sensitivity results demonstrate that AZD8055, AZD8186, BMS-75807, CZC24832, Dactolisisb, AMG-319, AT13148, AZ960, AZDBB2, AZD2014, Dasatinib, Entospletinib, Foretinib, GNE-317, Luminespib, GSK2606414, IGFIR-3801, JAK-8517, JAKI_8709, Mirin, Nilotinib, NU7441, PCI-34051, PD173074, Pictilisib, PLX-4720, PRIMA-IMET, RVX-208, Staurosporine, Uprosertib, WZ4003 and XAV939 were effective in the low-risk group (Figures 7, S2, and S3). Conversely, Afatinib, Dabrafenib, Dihydrorotenone, Erlotinib, Gefitinib, Lapatinib, Osimertinib, PF-4708671, Ribociclib, Sapitinib and TAFI-5496 were more effective in the high-risk group (Figures 7, S2, and S3).

Figure 7:

Analysis of drug sensitivity across high- and low-risk groups. (A) Sensitivity of Afatinib differentiated between high- and low-risk cohorts; (B) comparative sensitivity of AMG-319 in high-vs. low-risk groups; (C) differential sensitivity to AT13148 between the risk categories; (D) assessment of AZ960 sensitivity across the risk strata. (E) Sensitivity variation of AZD1332 between high- and low-risk groups; (F) AZD2014 sensitivity contrasted between the risk groups; (G) sensitivity profile of AZD8055 in the two risk classifications; (H) disparity in AZD8186 sensitivity between high- and low-risk patients; (I) evaluation of BMS-754807 sensitivity across high- and low-risk groups; (J) sensitivity differentiation of CZC24832 among the risk cohorts; (K) dabrafenib sensitivity comparison between high- and low-risk classifications; (L) examination of dactolisib sensitivity between the identified risk groups.

Using an unsupervised consensus clustering algorithm to validate the prognostic risk signature

This study employed a combination of statistical analyses and clustering algorithms to identify and characterize prognostic risk signatures in a cohort of cancer patients. Specifically, the approach involved the utilization of univariate Cox regression analysis, LASSO, and multivariate Cox regression analysis to develop the risk signature. This signature was then validated using unsupervised consensus clustering. This approach led us to identify three distinct gene subtypes, referred to as gene C1-3 (Figure 8A), with gene C3 showing the best survival prognosis (Figure 8B). To gain additional insights into the relationship between the risk score and genotypes, a Sankey diagram was employed. This diagram revealed that gene C3 was closely associated with the low-risk group, while the other two clusters were more closely associated with the high-risk group (Figure 8C). Dimensionality reduction analyses using PCA and t-SNE further confirmed these findings, with clear clustering observed between high and low-risk groups based on risk score, as well as between the three gene subtypes (Figure 8D and E). The study also investigated the relationship between the risk signature and the tumor microenvironment, including the evaluation of immune cell expression and scores for immune infiltration, stromal cells, and ESTIMATE scores. Notably, gene C3 exhibited significantly lower immunescores, stromalscores, and ESTIMATEscores compared to C1 and C2, suggesting a less immunogenic and stromal-rich tumor microenvironment in this group (Figure F-H; p<0.001). In line with this, the study observed that stroma score_XCELL, endothelial cell_XCELL, cancer-associated fibroblast_MCPCOUNTER, cancer-associated fibroblast_XCELL, and endothelial cell_MCPCOUNTER were all lower in C3 when compared to the other two gene subtypes (Figure 8J). The study’s analysis of immune checkpoints revealed significant differences in the expression of several key immune checkpoint genes. These included ADORA2A, BTNL2, BTLA, CTLA4, C10orf54, CD (27, 28, 40, 40LG, 44, 48, 70, 80, 86, 160, 200, 200R1, 244, 274, 276), IDO (1, 2), ICOS, ICOSLG, LAG3, LAIR1, TNFRSF (4, 8, 9, 18), TNFSF (4, 9, 14, 15, 18), TIGIT, TMIGD2, PDCD1LG2, PDCD1, KIR3DL1, HAVCR2, VTCN1, and NRP1 (Figure 8I). Notably, gene C3 exhibited lower expression of several immune checkpoint genes, including CTLA4, LAG3, IDO1, and PDCD1LG2, which could potentially make it more responsive to immune checkpoint inhibitor therapy (Figure 8I). Overall, our findings highlight the potential clinical utility of our prognostic risk signature in guiding personalized treatment decisions for cancer patients, particularly those in the low-risk and gene C3 subgroups.

Figure 8:

The relationship between gene subtypes and various aspects of the TME, including immune cell infiltration, immune checkpoint expression, TMB, and drug sensitivity. (A) The consensus clustering of STAD patients with K=3; (B) overall survival curves for three different gene subtypes; (C) the Sankey diagram illustrates the distribution of different risk groups and gene subtypes; (D) PCA and tSNE analysis of the high- and low-risk groups; (E) PCA and tSNE analysis of the different gene subtypes; (F–H) boxplots depicting the statistical differences in immune score, stromal score, and estimat score; (I) boxplots showing the differential expression of immune checkpoint genes; (J) the expression levels of immune cells in different gene subtypes.

Our analysis on drug sensitivity underscored that distinct compounds manifested efficacy tailored to specific gene subtypes, either C3 or C2. For the C3 group, therapeutic potential was evident with compounds including AZ960, AZD1332, AZD2014, AZD8055, Dasatinib, BMS-754807, Cediranib, Doramapimod, Entospletinib, IGFIR_3801, JAK_8517, JQ1, NU7441, Nutlin-3a, PCI-34051, RO-3306, RVX-208, SB216763, Staurosporine, WZ4003, and XAV939 (Figures 9, S4, and S5). On the other hand, the C2 group exhibited responsiveness to Acetalax, Afatinib, BMS-345541, Dabrafebib, Dihydrorotenone, Erlotinib, Gefitinib, Lapatinib, Leflunomide, ML323, OSI-027, Oxaliplatin, Pevonedistat, PF-4708671, Sapitinib, and TAFFI_5496 (Figures 9, S4, and S5). These findings suggest that personalized therapeutic interventions based on gene subtypes may lead to better clinical outcomes for patients with different prognostic risk signatures.

Figure 9:

Analysis of drug sensitivity in different gene subtypes. (A) Sensitivity of Acetalax differentiated among C1, C2, and C3 gene subtypes; (B) comparative sensitivity of Afatinib in these groups; (C) differential sensitivity to AZ960 among the gene subtypes; (D) assessment of AZD1332 sensitivity across the gene strata. (E) Sensitivity variation of AZD2014; (F) AZD8055 sensitivity contrasted among groups; (G) sensitivity profile of BMS-345541 in the three gene classifications; (H) disparity in BMS-754807 sensitivity among C1, C2, and C3 patients; (I) evaluation of Cediranib sensitivity; (J) sensitivity differentiation of dabrafebib among the gene cohorts; (K) dasatinib sensitivity; (L) examination of dihydrorotenone sensitivity.