We downloaded RNA sequencing (RNA-seq) and expression files and mutation files from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository). The data included tumor tissues of 343 STAD patients and 30 matched normal tissues. Data were downloaded and handled according to TCGA guidelines.

According to the study by Tsvetkov et al. (2022), 19 cuproptosis-associated genes were evaluated (Supplementary Table S1). Correlation between cuproptosis-related genes and differentially expressed lncRNAs was evaluated. Pearson’s correlation coefficients (R) of gene expression patterns were used as a measure of gene coexpression. The PCC threshold to retrieve cuproptosis-related lncRNAs was 0.4 (|R| > 0.4), with a p value < 0.001.

The downloaded clinical and demographic data of STAD patients were analyzed with univariate Cox regression analysis to identify lncRNAs associated with patient overall survival (OS) and those associated with cuproptosis were further identified as candidate lncRNAs for the construction of prognostic signature. Lasso regression was performed to screen lncRNAs that were truly correlated with a patient’s survival on the basis of 10-fold cross-validation. Based on the nine optimal lncRNAs identified, the risk scores of patients were calculated according to the following formula: r i s k   s c o r e = ∑ i n X i * Y i Where X was regression coefficient and Y was expression level of cuproptosis-related lncRNAs.

A total of 343 STAD patients were allocated to either the training cohort or the test cohort randomly in a 1:1 ratio for constructing and validating the cuproptosis-related lncRNAs signature. Patients in each cohort were classified into either low--risk group or high-risk group according to the cut-off value, which was the median risk score (Meng et al., 2019; Hong et al., 2020). The Chi-square test and the receiver operating characteristics (ROC) curves were used to help determine if observed OS was in line with expected OS, and the 1-year, 3-years, and 5-years OS rates were compared between the low-risk group and the high-risk group by Kaplan–Meier analysis. We further constructed a nomogram with cuproptosis-related lncRNA risk score and established clinical risk factors to calculated patient survival time. Then concordance index (C-index) and calibration curves were used to evaluate the prediction power of the nomogram. Finally, stratified analysis was used to assess whether the signature retained its predictive ability in subgroups of patients (stages I–II and stages III–IV). The “survival”, “rms”, “survminer” and “timeROC” R packages were used.

We used principal component analysis (PCA) to characterize cuproptosis-related lncRNAs expression patterns. PCA is a common unsupervised method for the analysis of gene expression data. 3D scatter plots were used to visualize the relationship between the three variables of samples. The analysis of differentially expressed genes (DEGs) was performed with the glm method of the “edgeR” R package. We set the threshold value of log fold change (log2FC) at |log2FC| ≥ 1, with a false discovery rate (FDR) < 0.05, to identify important DEGs. Gene Ontology (GO) was used to interpret DEGs for the relevant cellular components, biological processes, and molecular functions. Differential Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between the high-risk group and the low-risk group were screened using Gene Set Enrichment Analysis (GSEA), with a FDR <0.25.

Single-sample GSEA (ssGSEA), an extension of GSEA, was used to calculate separate enrichment scores for immunological pathways by the normalized enrichment score (NES) (Subramanian et al., 2005). Each ssGSEA enrichment score represents the degree to which the genes are coordinately upregulated or downregulated within a sample.

We downloaded the somatic mutation file and calculated each patient’s tumor mutation burden (TMB) score. The influence of TMB on patient OS was evaluated by Kaplan–Meier analysis and compared between the high- and low-risk groups by t-test. Maftools R package was used.

To predict treatment response of immune checkpoint blockades (ICBs), tumor immune dysfunction and exclusion (TIDE) algorithm was used to identify signatures of T cell dysfunction and signatures that exclude T cell infiltration into tumors (Jiang et al., 2018). To predict treatment response of the most important groups of drugs again STAD, the half-maximal inhibitory concentrations (IC₅₀) were calculated using pRRophetic as described in Genomics of Drug Sensitivity in Cancer (GDSC) (Geeleher et al., 2014).

Figure 1 illustrates the results of the search and the process of screening. A total of 16,773 lncRNAs that may be associated with 19 cuproptosis-associated genes were found. Among these lncRNAs, 298 lncRNAs met the pre-defined criteria ((|R| > 0.4). All 298 lncRNAs upregulated the expression of cuproptosis genes in the Sankey diagram (Figure 2A. Univariate Cox regression analysis found that 13 lncRNAs were prognostic factors of patient survival (Figure 2B).

To construct a risk model with cuproptosis-related lncRNAs in STAD, we randomly allocated 343 STAD cases into the training set and the test set at 1:1 ratio. The chi-square test showed that the two groups were comparable in terms of both clinicopathologic and demographic parameters (Table 1).

To avoid overfitting, nine lncRNAs were further identified by LASSO regression method (Figures 2C,D. A formula was established with the expression levels of nine lncRNAs:

Risk score = LINC01094 × (0.5250) + AC022182.1 × (2.02146) + AC011747.1 × (0.1655) + LINC02476 × (0.1295) + AC005014.2 × (−0.6903) + AC090809.1 × (0.2959) + AC084781.2 × (0.3942) + SENCR × (0.6958) + AC010422.4 × (−0.8166) (Meng et al., 2019).

As expected, the high-risk group had worse survival in each sample set (Figure 3).

The areas under the l-, 3- and 5-years ROC curves (AUC) were 0.719, 0.773, and 0.755 respectively (Figure 4A). The AUC of risk score was 0.719 and the C-index in the risk model was 0.726, indicting a perfect predictive ability (Figures 4B,C). In the uni-Cox regression, the hazard ratios (HR) of the risk score was 1.0726 (p < 0.001), and in the multi-Cox regression, HR of the risk score was 1.092 (p < 0.001) (Figures 4D,E).

A nomogram model was drawn to predict OS of patients (Figure 5A). The calibration plots showed the predicted l-, 3- and 5-years OS was consistent with the actual OS (Figure 5B). Thus the nomogram was well calibrated, with good prediction of patient survival. The high value of C index (0.726) indicated that the nomogram has excellent discriminative ability.

The results of decision curve analyses to compare the performance of the nomogram are shown in Figure 5C. The nomogram has greater net benefit than other clinical parameters in all patients.

The 3D scatter diagram showed the low-risk group and the high-risk group had distinct aggregation features of PCA (Figures 6A–C). GO analysis indicated related biological processes included B cell activation signaling pathway, antigen receptor−mediated signaling pathway, and immune response−regulating signaling pathway; related cellular components included immunological synapse, endocytic vesicle membrane, endocytic vesicle, T cell receptor complex, and immunoglobulin complex, and related molecular functions included immune receptor activity, heparin binding, glycosaminoglycan binding, sulfur compound binding, immunoglobulin receptor binding, and antigen binding (Figures 6D,E). GSEA identified genes involved in PI3K−Akt signaling pathway, cell adhesion, cytokine−cytokine receptor interaction and chemokine signaling pathway were differentially expressed between the low--risk group and high-risk group (Figures 6F,G).

Somatic mutations between the two groups were compared. The ten most mutated genes were TP53, TTN, PCLO, ZFHX4, CSMD3, SYNE1, ARID1A, LRP18, MUC16, and ACVR2A. The high-risk group had more frequent TP53 mutation (Figures 7A,B) but overall lower TMB (Figure 7C). Patients with higher scores and lower TMB had the worst prognosis among the four groups (Figures 7D,E).

The TIDE scores were significantly higher in the high-risk group compared to the low-risk group. This indicated that TIDE could be used to evaluate sensitivity to ICB therapy for STAD patients (Figure 8A). Indeed, several immune-related pathways had different activities between the two groups. Patients in the high-risk group had higher activities in terms of T cell co−inhibition and check−point (Figure 8B). Drug sensitivity comparison showed most drugs have similar IC₅₀ between the two groups, and there were eight drugs that had lower IC₅₀ in the high-risk group: PD−173,074, AZD8055, BEZ235, CGP-60474, Dasatinib, Pazopanib, TGX221, and HG-6-64-1 (Figure 8C).