A scoring system based on copper-related long RNAs to predict outcomes and immune features in pancreatic cancer

By searching the Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov↗), transcriptome expression data (including 4 normal and 179 tumor samples), clinical information data (n = 185, including survival, age, grade, and TNM stage), and gene mutation data (n = 169) were extracted from the TGCA-PAAD cohort. Progression-free survival (PFS) data were retrieved from the Xena Explore (https://xenabrowser.net/↗) database at the University of California, Santa Cruz. Ten cuproptosis-related genes (CRGs) selected for this study were derived from previously reported studies¹¹. A protein–protein interaction (PPI) network for those 10 CRGs was constructed using the String (https://cn.string-db.org/↗) database. Strawberry perl (version 5.32.1.1) was used to integrate matrix files for transcriptome (including mRNA and lncRNA) expression and mutation data. Autophagy related genes from MSigDB database (http://www.gsea-msigdb.org/gsea/msigdb/↗) of the Human Gene Set: GOBP_REGULATION_OF_AUTOPHAGY.

First, we extracted the expression of ten CRGs in the TCGA-PAAD cohort. The correlation coefficient between CRGs and lncRNA was calculated by Pearson correlation analysis. We used | correlation coefficient | of > 0.3 and a p-value of < 0.001 as the threshold to obtain cuproptosis-related lncRNA. Visualization was performed with Cytoscape (version 3.8.0) software. Additionally, we integrated expression and survival data of cuproptosis-related lncRNA in the TGCA-PAAD cohort for subsequent analyses. The above-mentioned analyses were performed by using R software “limma,” “dplyr,” “ggalluvial,” and “ggplot2” packages.

First, we divided the TCGA-PAAD integrated datasets randomly into a training cohort and a test cohort in a 1:1 ratio. Next, we screened for survival-associated lncRNA by performing Cox survival analysis on the training cohort. Then Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was performed on training cohort with survival-related lncRNA. Specifically, we constructed a penalty function to get a more refined model through cross-validation to find the minimum λ value mapping of lncRNAs²⁵, and these lncRNAs were used to build CRLss. The risk score was calculated as follows:

In the equation, the expression (i) and corresponding coefficient (i) represent the expression and Cox regression coefficient in CRLss. According to the median risk score in the training cohort, the training and testing cohorts were divided into high-risk and low-risk groups, and the testing cohort was used as a validation set to evaluate the predictive performance of the CRLss.

We first assessed the clinical baseline variability of the entire cohort, training cohort, and testing cohort to validate the predictive performance of the CRLss. With the R software “limma,” “scatterplot 3d” package, allGene, cuproptosis-related gene, cuproptosis-related lncRNA, and risk lncRNAs were used as the main characteristics to perform principal component analysis (PCA) of CRLss. Then, K-M survival analysis (including OS and PFS), risk curves, and risk heat maps for different cohorts were used to further verify the predictive performance of CRLss. Given the close association between cuproptosis and autophagy, we also explored the correlation between risk lncrnas and autophagy-related genes using Pearson correlation analysis ( | correlation coefficient |> 0.5 and P < 0.001 ).

In addition, we also used the R software “survival” and “survminer” to identify the survival correlation of different clinical characteristics, including age, gender, pathological stage, and grading, in the high- and low risk-groups of CRLss, which was used to evaluate whether the constructed CRLss was applicable to different clinical groups of PAAD patients. Univariate and multifactorial COX regression analyses were used to assess whether risk score and other clinical characteristics were independent prognostic factors. We calculated their concordance index (C-index) through the R package “dplyr,” “survival,” “rms,” and “pec,” which was used to evaluate their predictive ability in the model.

With the help of the R package “TimeROC,” “Survival,” and “Survminer,” the 1-, 3-, and 5-year ROC survival curves in the CRLss were plotted. Area under the curve (AUC) were used to evaluate the clinical prognostic value of the CRLss. Clinicopathological parameters were also stratified as subgroups for analysis. Additionally, based on the results of uni- and multi-factor Cox analysis, logistic model and Cox proportional risk model, we constructed a nomogram consisting of risk score, clinical features, and survival prognosis to predict 1-, 3-, and 5-year OS in PAAD patients. A calibration curve based on the Hosmer–Lemeshow goodness of fit test was used to assess the clinical credibility of the nomogram.

With | log2 fold change (FC) | of > 1 and p of < 0.05 as the threshold, differentially expressed genes (DEGs) in the high- and low-risk groups of CRLss were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)²⁶ enrichment analyses were performed to explore the biological functions of DEGs.

TIME is closely related to the occurrence and progression of cancer. Immune cell infiltration and stromal cell metabolism in the tumor microenvironment have a profound influence on the TIME²⁷. The ESTIMATE algorithm was used to calculate the abundance of immune cells and stromal cells in tumor tissue, as well as the purity of tumor tissue²⁸. Next, the immune function scores of the tumor samples were calculated using ssGSEA²⁹, and the differences between tumor microenvironment and immune function scores in different risk groups of CRLss were analyzed using the R software “reshape2” and “ggpubr” packages.

In addition, the degree of immune cell infiltration is one of the indicators to predict the immunotherapy response, which is closely related to the prognosis and survival of pancreatic cancer³⁰. CIBERSORT is the most frequently cited tool for estimating immune cell infiltration³¹. We used the R language “CIBERSORT” package for deconvolution analysis of the gene expression matrix of immune-related cell subtypes. We set the perm value to 1000 to ensure the accuracy of the results. Then, the infiltration abundance of different immune cells in tumor samples was calculated, and the correlation between CRLss risk score and immune cells was analyzed.

Immunotherapy response refers to immune checkpoint inhibitors binding with corresponding immune checkpoint genes on tumor cells to activate the immune recognition and immune response of T cells to tumor cells to kill tumor cells³². Based on this, we explored the differential expression levels of immune checkpoint genes in high- and low-risk groups of the CRLss.

Tumor Immune Dysfunction and Exclusion (TIDE) is a newly developed computational method for predicting an immunotherapy response³³. We obtained the scores in the TCGA-PAAD samples from the TIDE (http://tide.dfci.harvard.edu/↗) database and analyze the variability of scores in different risk groups of the CRLss for predicting the immunotherapy response in the model.

Tumor mutation burden (TMB) is defined as the total number of somatic gene coding errors, base substitutions, and insertion or deletion errors detected per million bases. There is growing evidence that TMB expression levels correlate with the efficacy and prognosis of PD-1/PD-L1 inhibitors in selected tumors^34,35. Therefore, we explored the difference in TMB expression in CRLss high- and low-risk groups. Next, the R “survival” and “survminer” packages were used to obtain the optimal cutoff of TMB, which was used to plot the K-M survival curves of TMB in different risk groups of the CRLss.

R “pRRophetic,” “ggpubr,” and “limma” packages were used to obtain the half-maximal concentration (IC50) of the drug in the high- and low-risk groups to identify the difference in drug sensitivity of different risk groups of the CRLss. Then we screened clinically commonly used drugs for presentation by drawing box plots.

Strawberry version of perl (version 5.32.1.1), R software (version 4.1.2), and related packages were used for statistical analysis of data and graphing. Cytoscape (version 3.8.0) was used to visually demonstrate the correlation between 10 cuproptosis-related genes and 34 cuproptosis-related lncRNAs. Wilcoxon’s and Kruskal–Wallis were used to compare differences between groups using Pearson’s correlation coefficient to assess correlations between variables, and Kaplan–Meier and Cox regression models were used for survival correlation analysis. All statistical P-values were bilateral, and a P of < 0.05 was considered statistically significant without special note.

A flow chart is shown in Fig. 1 to directly reflect the ideas and details of this research. The PPI network of cuproptosis-related genes showed that these genes are closely related to various biological processes (Fig. 2A). and then we integrated the TCGA-PAAD transcriptome data and extracted the expression levels of 10 cuproptosis-related genes. A total of 180 cuproptosis-related lncRNAs were obtained by Pearson correlation analysis (Fig. 2B). Next, we randomised the TCGA-PAAD patients into a training and testing cohort in a 1:1 ratio. In the training cohort we screened 34 lncRNAs associated with PAAD survival by Cox survival analysis. We found that except for CASC8 (hazard ratio = 1.644), the other 33 lncRNAs were low-risk lncRNAs [hazard ratio (HR) < 1] (Fig. 2C). Figure 2D shows the expression landscape of 34 survival-related lncRNAs in each TCGA-PAAD sample.

Table 1 demonstrates the baseline characteristics of the clinical features (including age, gender, grade, and stage) for the different subgroups. Then, LASSO-COX regression analysis was performed on 34 survival-related lncRNAs (Fig. 3A, B). Cross-validation yielded the minimum λ value and finally mapped six cuproptosis-related lncRNAs (AL117335.1↗, AC044849.1↗, AL358944.1↗, ZNF236-DT, Z97832.2↗, and CASC8) (Table 2). According to the risk scoring system established above, the prognostic model was constructed. Table 3 shows the association of risk lncRNAs and Curoptosis-Related Genes. The PCA showed that lncRNAs involved in model construction were more obvious than cuproptosis gene, cuproptosis lncRNA, and allGene in the high-low risk group of model differentiation (Fig. 3C–F). Moreover, significant differences were also observed in the expression of six lncRNAs in the high-low risk group (Fig. 3G–L). Given the close association between Copper metabolism and autophagy³⁶, we demonstrated the correlation between risk lncRNAs and autophagy-related genes using Sankey plots. (Supplementary Fig. 1).

Figure 4A shows the heat map of different clinical features in the high- and low-risk group. K-M survival analysis of clinical feature subgroups showed that the low-risk group was superior to the high-risk group in age, gender, and grade (Fig. 4B–G). Although the same phenomenon as other clinical features was observed in stages I–II (Fig. 4H), no statistically significant P values were observed in stages III–IV (Fig. 4I). The reason might be the small sample size of stage III–IV cases in the TCGA-PAAD cohort. Nevertheless, we observed a trend toward longer survival in the low-risk group than in the high-risk group. Overall, the model was applicable to PAAD patients with different clinical characteristics.

We assessed heatmaps of expression, risk score, and survival status for the entire cohort, the training cohort, and the testing cohort using the scoring system to further verify the performance of the model. The results showed that the three cohorts were consistent (Fig. 5A–I). Survival analysis showed that the low-risk group had better OS and PFS than the high-risk group (Fig. 5J–O). Age, grade, and risk score were independent factors affecting the prognosis of PAAD patients (Fig. 6A, B). The C-index curve indicated that the risk score was superior to other clinical features in predictive performance (Fig. 6C). In addition, the AUC values at 1, 3, and 5 years were 0.707, 0.762, and 0.880 for the ROC survival curves, respectively (Fig. 6D). Consistent with the C-index curve, the AUC of risk score at 1, 3, and 5 years was also significantly higher than that of other clinical features. Taken together, these results confirmed the reliable clinical predictive accuracy of this model.

Based on the results of Cox analysis, we integrated the risk score, clinicopathological parameters (age, grade), and survival data of TCGA-PAAD patients and constructed the nomogram by concretized Cox regression model. The predicted OS of PAAD patients at 1, 3, and 5 years is shown in Fig. 6H. Calibration curves showed that nomogram-predicted values were reliably consistent with actual values.

Through differential gene analysis of the high- and low-risk groups, we finally obtained 1318 DEGs, including 99 upregulated and 1219 down-regulated genes. Then, we explored the biological function of these DEGs by GO and KEGG analyses. GO analysis showed that DEGs were responsible for such activities as T cell activation, T cell receptor complex, metal ion transmembrane transporter activity, and channel activity (Fig. 7A, B). KEGG analysis showed that DEGs were enriched in multiple signaling pathways, such as cytokine-cytokine receptor interaction, cell adhesion molecules, chemokine signaling pathway, and T cell receptor signaling pathway (Fig. 7C, D). In general, DEGs were closely related to ion transport, tumor proliferation and metastasis, and immune regulation.

First, we performed a correlation analysis between risk score and immune cell infiltration abundance. The results showed that a total of 65 immune cells were significantly associated with risk scores (Supplementary Fig. 2). Only a few immune cells were positively correlated with the risk score (Fig. 8A). Then, ssGSEA analysis showed that the content of B cells, CD8⁺ T cells, immature dendritic cells, mast cells, neutrophils, NK cells, plasmacytoid DCs, T helper cells, Th1 cells, and tumor-infiltrating lymphocytes in the low-risk group was higher than that in the high-risk group (Fig. 8B). Immune function analysis also showed that the low-risk group was superior to the high-risk group in CCR (a chemokine receptor), checkpoint, cytolytic activity, promoting inflammation, T cell, T cell costimulation, and type II IFN response (Fig. 8C, D). Based on the ESTIMATE algorithm, the StromalScore and ImmuneScore analysis in different risk groups also obtained consistent results (Fig. 8E). Figure 8F shows a scatter plot of immune cells with the top 10 correlation coefficients. In addition, the expression of 30 immune checkpoint-related genes in the low-risk group was significantly higher than that in the high-risk group (Fig. 9A). The same trend was observed in the TIDE analysis (Fig. 9F). Collectively, these results suggested that patients in the low-risk group had a higher degree of immune cell infiltration and were a potential population to benefit from immunotherapy.

There is growing evidence of a close relationship between TMB and immunotherapy response. Therefore, the mutation landscape of PAAD was also under our attention. The waterfall plot showed a lower mutation frequency in the low-risk group (68.54%) than in the high-risk group (98.63%). Additionally, the top three mutated genes in both groups were KRAS (H/L: 82%/44%), TP53 (H/L: 73%/43%), and SMAD4 (H/L: 26%/18%) (Fig. 9B, C). This is consistent with previous complete exome sequencing of pancreatic cancer³⁷. Figure 9D shows that the TMB of the high-risk group was higher than that of the low-risk group. K-M survival analysis showed that patients with low-risk scores and low levels of TMB had better outcomes (Fig. 9E).

We predicted clinical drug response by using the R software package “pRRophetic.” The results showed that there were 56 chemotherapeutic and targeted drugs with differential IC50 values in high- and low-risk groups (Supplementary Fig. 2). We screened clinically common drugs such as chemotherapy and molecular targeting drugs to demonstrate. Compared to low risk groups, We found that low-risk group was better suited to Lenalidomide, Metformin, Nilotinib, Pazopanib, Temsirolimus. While high-risk group was more suitable for Bicalutamide, Epothilone.B, Lapatinib, Paclitaxel, Sorafenib (Fig. 10A–J).

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.