RNA-sequencing (RNA-Seq) expression profiles for LSCC were downloaded from the TCGA database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) and the full TCGA-LSCC dataset (501 patients with LSCC, phs000178) was included in the present study. The ferroptosis-related genes presented in the present study were obtained from Liu et al (16). All of the analysis methods were implemented in R (version 4.0.3; R Foundation for Statistical Computing), as previously described (17–19).

The method used was similar to that of previous studies (17–20). The R software's ‘limma’ package was used to study differentially expressed mRNA. To identify mRNAs with differential expression, P<0.05 and log₂ (fold change) >1 or log₂ (fold change) <-1 were used as thresholds. An enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) can be used to study gene function in conjunction with high-level genome functional information. KEGG pathways were enriched using the ‘clusterProfiler’ package (version 3.18.0) in R software to describe the carcinogenesis of mRNA. In addition, a box plot was created using R software with the ‘ggplot2’ package and a heat map with the ‘pheatmap’ package.

In the present study, Kaplan-Meier survival analysis was divided into two parts: The first part was a Kaplan-Meier survival analysis (R software version 4.0.3) of the TCGA data. To generate Kaplan-Meier curves, log-rank tests were used to calculate P-values, hazard ratios (HRs) with 95% confidence intervals (CIs) and survival time between distinct groups. Disease-free survival (DFS) time was the period between the start of randomization and the time of the last follow-up visit or the time of death (from any cause). The progression-free survival (PFS) time was defined as the time between the start of randomization and the onset of (any aspect of) tumor progression, death (from any cause) or the last visit after the last randomization. The overall survival (OS) time was calculated from the time of randomization until death (from any cause) or the last follow-up visit. All of the analysis methods and R packages were implemented in R (version 4.0.3). P<0.05 was considered to indicate a statistically significant difference (21). The second part of the analysis was a Kaplan-Meier survival analysis of the cohort data of the present study. PFS and OS time for metastatic LSCC were compared for high and low expression of HSPA5 in the cohort and calculated using Kaplan-Meier analyses. The final follow-up was performed on May 1st, 2022.

As previously described, multivariate and univariate Cox regression analyses were performed to identify the suitable terms required for building a nomogram (22). A forest plot was used to show the P-value, HR and 95% CI of each variable with the ‘forestplot’ package in R software (version 4.0.3), and a nomogram was developed according to the results of the multivariate Cox proportional hazards analysis to predict the 1-, 3- and 5-year overall recurrence. The nomogram was used to provide a graphic representation of the factors used to calculate the risk of recurrence for an individual patient using the points associated with each risk factor through the ‘rms’ package in R software (version 4.0.3).

The mRNA stemness index was calculated using the machine-learning method of one-class logistic regression (OCLR) developed by Malta et al (23). Based on mRNA expression signatures, 11,774 genes were detected in the gene expression profile. Spearman's correlation coefficient was utilized to analyze the relationship between variables, followed by the transformation of the dryness index to a standardized range of [0,1] through a linear conversion involving subtraction of the minimum value and division by the maximum value. The aforementioned analysis methods and R packages were implemented in R software (version 4.0.3) as previously described (23).

SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3 and PDCD1LG2 were selected as immune checkpoint-relevant transcripts (24–29), and all the expression values of these eight genes were extracted. The aforementioned analyses and R packages were implemented in R software (version 4.0.3) using the ‘ggplot2’ and ‘pheatmap’ packages (version 4.0.3) (30–32).

First, the potential ICB response was predicted using a tumor immune dysfunction and exclusion (TIDE) algorithm as previously described (33). The TIDE algorithm uses a set of gene expression markers to assess two different mechanisms of tumor immune escape: Dysfunction of tumor-infiltrating cytotoxic T lymphocytes (CTLs) and exclusion of CTLs by immunosuppressive factors. High TIDE scores were associated with poor efficacy of ICB therapy and short survival time after administration of ICB therapy.

The ‘GSVA’ package in R software (version 4.0.3) was used for the analysis, selecting the single-sample gene set enrichment analysis (‘ssGSEA’) method as a parameter. The correlation between genes and pathway scores was then analyzed using Spearman's correlation coefficient (34).

The chemotherapeutic response for each sample was predicted based on the largest publicly available pharmacogenomics database, the Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/). The prediction process was implemented with the R package ‘pRRophetic’ (R Foundation for Statistical Computing, version 4.0.3). The IC₅₀ of the samples was estimated by ridge regression. All parameters were set as the default values using the batch effect of combat and tissue type of all tissues, and the duplicate gene expression was summarized as the mean value.

Metastatic LSCC tissues and para-cancerous tissues (n=100) were obtained from the Fifth Affiliated Hospital of Zhengzhou University (Zhengzhou, China). A total of 100 patients were included in the study, which was performed from July 2015 to December 2018. All patients had undergone surgery and eventually developed lung metastases. These patients were treated according to the National Comprehensive Cancer Network guidelines for LSCC (35). Metastatic patients who have not undergone surgery are excluded. Formalin-fixed paraffin-embedding (FFPE) tissues were used in the present study and the detailed procedure was as described previously (36). In brief, the first step was to dewax and rehydrate formalin-fixed and paraffin-embedded sections of LSCC and para-cancerous tissues (5 µm thick). Antigen retrieval was performed by heating the slides in 10 mM Tris buffer with 1 mM EDTA (pH 9) in a streamer for 20 min. Samples were immersed in 3% H₂O₂ for 5 min to inhibit endogenous peroxidase activity. Subsequently, after Tris-buffered saline containing 0.1% Tween had been used for washing, endogenous biotin was inhibited by sequential incubation with 0.1% anti-biotin protein and 0.01% biotin (Dako; Agilent Technologies, Inc.) for 10 min at room temperature. Subsequently, 3% skimmed milk powder was applied for 30 min at room temperature to block unspecific binding. Thereafter, LSCC and para-cancerous tissue sections were incubated with HSPA5 monoclonal mouse anti-human antibodies (cat. no. ab21685; dilution, 1:1,000; Abcam) at 4°C overnight. The sections were then serially rinsed and incubated with secondary antibody (cat. no. ab98799; dilution, 1:3,000; Abcam) at room temperature 1 h. Finally, color was developed by incubation with 3,3′-diaminobenzidine (DAB) substrate [ImmPACT DAB Peroxidase (HRP) substrate; cat. no. ab64238; Abcam], followed by counterstaining with hematoxylin for 5 min and mounting with Vecmount (cat. no. H-5000; Vector Laboratories). The IHC staining was evaluated independently by two experienced pathologists blinded to patient characteristics and results via microscopy (BX53; Olympus Corp.). Based on a combination of the % positive stained cells and the intensity of the staining, a semi-quantitative scoring system (H-score) was calculated, as previously described (32–34). The H-score was calculated as follows: [H-score=∑ % (0–100%) × intensity (1–3)]=(% weak-intensity cells ×1) + (% moderate-intensity cells ×2) + (% strong-intensity cells ×3). The median H-score was selected as the cut-off value for high or low HSPA5 expression.

R software (version 4.0.3) was used for all statistical analyses of the TCGA data. The statistical details of all the experiments are reported in the Materials and methods section, including the statistical analyses performed and statistical significance. GraphPad Prism 9.0 (Dotmatics) and IBM SPSS Statistics 24.0 (IBM Corp.) were used to statistically analyze the cohort data. Quantification of the HSPA5 density was performed using an unpaired t-test, and PFS and OS were calculated using the Kaplan-Meier estimator. Univariate and multivariate Cox proportional hazards regression models were used to estimate the HRs, along with the associated CIs and P-values. The Student's unpaired t-test and the χ² test were used for inferential statistical analysis. For all data, P<0.05 was used to indicate a statistically significant difference.

An analysis of RNA-Seq data from the TCGA database was performed to determine the presence of ferroptosis-related genes in LSCC. Prior studies have identified 22 genes that have been reported to serve crucial roles in regulating ferroptosis (16). These ferroptosis regulator genes are CDKN1A, HSPA5, TTC35/EMC2, SLC7A11, NFE2L2, MT1G, HSPB1, GPX4, FANCD2, CISD1, FDFT1, SLC1A5, SAT1, TFRC, RPL8, NCOA4, LPCAT3, GLS2, DPP4, CS, ALOX15 and ACSL4. Following this, the expression patterns of ferroptosis-related genes in LSCC and adjacent tissues were compared. Of these 22 ferroptosis-related genes, the mRNA expression levels of HSPA5, HSPB1, GPX4, FANCD2, CISD1, FDFT1, NFE2L2, SLC1A5, RPL8, NCOA4, TFRC and SLC7A11 were significantly increased in LSCC compared with those in the adjacent tissues (Fig. 1). However, the expression levels of CDKN1A, EMC2, SAT1, LPCAT3, GLS2, DPP4, CS and ACSL4 were significantly decreased in LSCC compared with those in the adjacent healthy tissues (Fig. 1). However, there was no significant difference in the expression levels of MT1G and ALOX15 in LSCC compared with those of the adjacent tissues (Fig. 1).

Ferroptosis is involved in the metabolism of cells and the immune system, so it has been suggested that it and its regulators could be linked to cancer survival (37). In light of this, the prognostic value of increased ferroptosis-related gene expression was evaluated in LSCC. To determine whether ferroptosis-related gene expression is prognostic, the DFS, PFS and OS were evaluated. There were 12 genes related to ferroptosis with increased expression, and only the higher expression of HSPA5 was significantly associated with a shorter DFS (Fig. 2A) and PFS (Fig. 2B), but not OS (Fig. 2C). There was no significant difference in DFS, PFS or OS associated with high and low expression of the other 11 ferroptosis-related genes (Fig. S1). However, HSPA5 was significantly elevated in lung adenocarcinoma tissue compared with that of normal tissue (Fig. S2A). However, there was no significant difference in the DFS, PFS and OS between high and low expression of HSPA5 (Fig. S2B-D). Of the 22 ferroptosis-related genes, only HSPA5 predicted DFS, PFS and OS in patients with LSCC, suggesting that it may be a potential biomarker for predicting LSCC.

A multivariate and univariate Cox analysis was performed to evaluate the independent prognostic value of HSPA5 in terms of DFS, PFS and OS of patients with LSCC. According to the univariate analysis, the group with high HSPA5 expression had a significantly worse DFS time (Fig. 3A) and the multivariate analysis demonstrated that high HSPA5 expression was independently associated with a significantly shorter DFS time (Fig. 3B). This result indicates that HSPA5 expression may serve as an independent prognostic factor for LSCC. To develop a clinically applicable method to predict a patient's survival probability, a nomogram was generated to construct a predictive model that considered clinicopathological covariates. With the multivariate and univariate analysis of DFS rates, a nomogram was created to predict the 1-, 3- and 5-year DFS rates in the discovery group, which demonstrates the similar results (Fig. 3C). The probability of the death of patients was also predicted using generalized linear regression. Among the predictors under evaluation was when HSPA5^hi vs. HSPA5^low met the P<0.05 risk assessment criterion. According to an ideal model, DFS rates in the entire cohort were well predicted at 1-, 3- and 5-year intervals (Fig. 3D). The independent prognostic significance of HSPA5 was then measured in terms of the PFS in patients with LSCC and the same results were obtained (Fig. 3E-H). By contrast, in terms of the OS of patients with LSCC, there were no significant differences between the HSPA5 high and low expression groups (Fig. S3). In summary, in both the univariate and multivariate Cox analysis, as well as the nomograms, HSPA5 was predicted to be an unfavorable prognostic factor.

Immune checkpoint inhibitors have become a common type of immunotherapy for patients with lung cancer in recent years (38,39). In the present study, TIDE scores were used to evaluate the ICB responses of patients with LSCC with high (HSPA5^hi) and low (HSPA5^low) expression of HSPA5. HSPA5^hi patients had significantly higher TIDE scores than HSPA5^low patients, indicating lower responses to ICB (Fig. 4A). Due to the strong association between immune checkpoint molecules and ICB responses, the expression of immune checkpoint molecules in patients with LSCC with HSPA5^hi and HSPA5^low was then assessed. PDCD1LG1, CTLA4 and PDCD1LG2 were all significantly upregulated, but not TIM3, LAG3, PDCD1, TIGIT and SIGLEC15, in the HSPA5^hi group (Fig. 4B). Previous studies have indicated that cancer progresses through the gradual loss of a differentiated phenotype and the acquisition of stem-cell-like characteristics (23). Consequently, numerous treatment approaches, including immunotherapy, do not kill cancer stem cells effectively (40,41). Therefore, the stemness features of patients with LSCC with HSPA5^hi and HSPA5^low were tested using an OCLR machine-learning algorithm, as previously described (23). Patients with LSCC with HSPA5^hi were demonstrated to have similar OCLR scores to those with HSPA5^low (Fig. 4C). In summary, these data indicate that patients with LSCC with HSPA5^hi are not sensitive to ICB.

As HSPA5 can be used as a prognostic biomarker for lung adenocarcinomas, high and low expression levels of HSPA5 of patients with lung carcinomas were compared. The volcano plot in Fig. 5A depicts the 237 upregulated and the 27 downregulated genes in the HSPA5^hi vs. HSPA5^low groups (Table SI) and Fig. 5B presents a heat map of the differentially expressed genes. The KEGG signaling pathways enriched with differentially expressed genes were analyzed to predict their primary biological actions. Fig. 5C shows the KEGG pathways enriched by the upregulated genes in the HSPA5^hi vs. HSPA5^low group, including ‘glutathione metabolism’, ‘focal adhesion’, ‘extracellular matrix (ECM)-receptor interaction’ and numerous cancer pathways. Fig. 5D presents the KEGG pathways enriched by the downregulated genes in the HSPA5^hi vs. HSPA5^low group, including ‘tyrosine metabolism’, ‘fatty acid degradation’ and ‘ATP-binding cassette transporters’. Correlations between the expression of HSPA5 with enriched signaling pathways with the use of ssGSEA analyses were also assessed, as previously described (42). Notably, positive correlations between HSPA5 expression and ‘ferroptosis’, ‘cellular responses to hypoxia’, ‘tumor proliferation signature’, ‘G2M checkpoint’, ‘MYC targets’ and ‘TGFB’ were demonstrated (Fig. 6).

Cisplatin, gemcitabine, docetaxel, paclitaxel and vinorelbine are the five most commonly used chemotherapy agents in the treatment of LSCC (43,44). The present study therefore assessed whether high or low HSPA5 expression was associated with sensitivity to treatment with the five different chemotherapeutic agents commonly used for LSCC by IC₅₀ analysis (Fig. 7). The IC₅₀ scores for cisplatin, gemcitabine, docetaxel, paclitaxel and vinorelbine were tested according to previous reports (45). Patients with LSCC in the HSPA5^hi group exhibited significantly lower gemcitabine, docetaxel, paclitaxel and vinorelbine IC₅₀ scores than those of the HSPA5^low group (Fig. 7B-E); however, the IC₅₀ scores were similar between groups for cisplatin (Fig. 7A). According to these results, patients with LSCC with HSPA5^hi may be more sensitive to chemotherapy.

A total of 100 metastatic LSCC tissues and para-cancerous tissues were obtained from the Fifth Affiliated Hospital of Zhengzhou University (Zhengzhou, China) to confirm the expression pattern of HSPA5. A total of 100 patients, of whom 47 (47%) were male and 53 (53%) were female, with a median age of 63 years, were included in the study, which was conducted from July 2015 to December 2018. All patients had undergone surgery and eventually developed lung metastases. These patients were treated according to the National Comprehensive Cancer Network guidelines for LSCC (35). Table I lists the detailed clinicopathological characteristics of the patients. To detect HSPA5 expression in each patient, IHC was performed with an anti-HSPA5 antibody. In Fig. 8A, representative IHC images of HSPA5 in cancerous and paraneoplastic tissues are provided. H-scores were used to perform a semiquantitative analysis of HSPA5 expression, as previously described (42). According to the quantitative analysis, the H-scores of HSPA5 in cancerous tissues were significantly higher than those in paraneoplastic tissues (Fig. 8B). Furthermore, it was demonstrated that the HSPA5 expression levels were different in every LSCC sample (Fig. 8B). The prognostic value of HSPA5 expression was then evaluated. The median H-score was used as a cut-off value for determining high and low expression of HSPA5. On the basis of the cutoff value for LSCC, the LSCC group of patients was divided into high-expression (HSPA5^hi) and low-expression (HSPA5^low) subgroups. The PFS and OS of the patients with HSPA5^hi and HSPA5^low metastatic LSCC was then evaluated. In line with the results for the TCGA database cohort, higher expression of HSPA5 in metastatic LSCC was significantly associated with shorter PFS (Fig. 8C) and OS (Fig. 8D). Furthermore, univariate analysis revealed that an advanced tumor-node-metastasis (TNM) stage (46) was a risk factor for PFS and OS (Table II). These risk factors from the univariate analysis were adopted as covariates in a multivariate Cox proportional hazards model, and an advanced TNM stage and HSPA5 were determined to be independent prognostic indicators for PFS and OS (Table III).