Ferroptosis-related lncRNAs signature to predict the survival and immune evasion for lung squamous cell carcinoma

Introduction: the investigation on the interactions between ferroptosis and lncRNAs for lung squamous cell carcinoma (LUSC) has been scare, and its impact on tumor immune microenvironment remained unknown. We aim to not only identify a ferroptosis-related lncRNAs signature for LUSC prognosis, but also evaluate its correlation to tumor immune evasion.

Methods: RNA sequencing data and survival information were obtained from The Cancer Genome Atlas database. A ferroptosis-related lncRNAs signature (FerRLSig) was developed and validated by univariate Cox regression, Least Absolute Shrinkage and Selection Operator regression and multivariate Cox regression. The tumor immune microenvironment and immune evasion were subsequently evaluated based on the FerRLSig stratification.

Results: the FerRLSig consisted of 10 ferroptosis-related lncRNAs and significantly associated with overall survival with satisfactory area under curve (HR = 2.240, 95% CI: 1.845–2.720, p < 0.001, 5-years AUC: 0.756). Based on the FerRLSig stratification, the high-risk group demonstrated not only significantly higher immune infiltration, but also more profound T cell dysfunction and immune evasion, which might ultimately lead to the resistance to current immune checkpoint inhibitors.

Conclusion: a robust prognostic FerRLSig for LUSC has been developed and validated, demonstrating a close association not only with tumor immune cell infiltration, but also with T cell dysfunction and immune evasion. Further investigation is warranted to better improve the survival of LUSC patients based on the FerRLSig stratification.

Introduction

Lung cancer has been the top one cause of cancer death with one of the highest incidence rates second only to breast cancer worldwide (Sung et al., 2021). Lung squamous cell carcinoma (LUSC) accounts for approximately 20% of all lung cancer cases and constitutes the bulk of non-small cell lung cancer with lung adenocarcinoma (Barta et al., 2019). There has been numerous effective targeted therapies for lung adenocarcinoma, which has significantly prolonged the survival of lung adenocarcinoma patients with certain mutations (Sordella et al., 2004; Solomon et al., 2014; Ramalingam et al., 2020; Shaw et al., 2020; Wu et al., 2020). Moreover, a large number of researches have utilized the transcriptome data of lung adenocarcinoma to build various prognostic and predictive tools, to derive useful risk stratification, and to provide valuable insights on the development of sensitive drugs (Guo et al., 2021; Lu et al., 2021; Zheng et al., 2021). However, compared to lung adenocarcinoma, LUSC lacks effective targeted therapy and is generally less well-defined in terms of gene expression profile.

On the other hand, tumor immune evasion has been identified as one of the hallmarks of cancer and closely related to the tumor immune microenvironment (Hanahan and Weinberg, 2011; Gajewski et al., 2013). And fortunately, immune checkpoint inhibitors tackling tumor immune evasion have made remarkable breakthroughs in LUSC and significantly improved the LUSC patients’ survival (Brahmer et al., 2015; Paz-Ares et al., 2018). But still some LUSC patients were resistant to the current treatments including immunotherapy, leading to intractable progression or relapse, and ultimately cancer death. Multiple studies have identified increased CD8⁺ T cell infiltration as a favorable prognostic factor (Lee and Ruppin, 2019; Morad et al., 2021). However, T cell dysfunction has also been recognized as an important mechanism of immunotherapy resistance (Thommen and Schumacher, 2018). Therefore, further evaluation of the LUSC gene expression pattern’s impact on tumor immune microenvironment and immune evasion are still needed for LUSC patients.

Ferroptosis is an iron-dependent, oxidatively regulated cell death activated by extrinsic blockade of the cystine/glutamate transporter or intrinsic blockade of intracellular antioxidants. Recent studies have demonstrated that ferroptosis plays a significant part in tumorigenesis and various treatment sensitivity, thus might be a useful tool in cancer prognosis and patient stratification (Dixon et al., 2012; Chen et al., 2021). Moreover, long non-coding RNA (lncRNA), with more than 200 nucleotide and without functional protein translation, has been found to be closely related to tumorigenesis and tumor progression via ferroptosis in recent studies (Mao et al., 2018; Wang et al., 2019; Zhang et al., 2020; Statello et al., 2021). Therefore, interactions between ferroptosis and lncRNAs are likely to be critical to overcome cancer progression. However, investigation on the ferroptosis-related lncRNAs signature on LUSC has been scare and its impact on tumor immune microenvironment remained unknown, thus warranting further investigation.

In this study, we utilized The Cancer Genome Atlas database on lung squamous cell carcinoma (TCGA-LUSC) (Cancer Genome Atlas Research Network, 2012), and we aimed to develop a ferroptosis-related lncRNAs signature for LUSC prognosis. Furthermore, explorations on tumor immune microenvironment, tumor immune evasion and T cell dysfunction were also performed based on the FerRLSig stratification.

Materials and methods

Data acquisition

RNA sequencing data of TCGA-LUSC patients and all available clinical data were downloaded from the Genomic Data Commons portal (https://portal.gdc.cancer.gov/) (Grossman et al., 2016). Overall survival (OS) was calculated from the time of lung cancer diagnosis to the time of death or the last follow-up. Patients with incomplete survival information or with OS less than 30 days were excluded.

This is a retrospective study based on publicly available TCGA database. The Ethics Committee of our hospital has confirmed that no additional ethical approval or informed consent is required.

A list of 108 validated ferroptosis genes was obtained from FerrDb database (http://www.zhounan.org/ferrdb) (Zhou and FerrDb, 2020). LncRNAs was identified in TCGA-LUSC RNA sequencing data through GENCODE annotation (https://www.gencodegenes.org/) (Frankish et al., 2021). The correlations between 108 ferroptosis genes and the lncRNAs expression in the entire set of TCGA-LUSC were analyzed with Pearson’s correlation, and 4,259 ferroptosis-related lncRNAs were identified by the selection criterion of |Pearson R|) > 0.3 and p < 0.001 with any one of the 108 ferroptosis genes.

Development and validation of the ferroptosis-related LncRNAs signature

The TCGA-LUSC patients were randomized into a training set and a testing set by the “caret” R package at the ratio of 7:3 (Kuhn, 2008). The FerRLSig was established with the training set while the validation was performed with the testing set and the entire set. Univariate Cox regression, Least Absolute Shrinkage and Selection Operator (LASSO) regression and multivariate Cox regression were applied in order to establish the final FerRLSig while avoiding overfitting with cross-validation (Friedman et al., 2010). The correlation between FerRLSig lncRNAs and ferroptosis genes was further visualized in both heatmap and Sankey diagram. The risk score of every patient was calculated as ${\displaystyle \sum }_{i=1}^{n}Coef\left(i\right)\times Expr\left(i\right)$ , with Coef(i) and Expr(i) representing the regression coefficient and expression level for each FerRLSig lncRNA respectively. The entire set of TCGA-LUSC patients were stratified into low-risk and high-risk groups by the median of the calculated risk scores.

Further evaluation of the FerRLSig and establishment of the prognostic nomogram

The prognosis effect of the FerRLSig was further evaluated with Kaplan–Meier survival plot, subgroup analysis by universal clinical characteristics, correlation with universal clinical variables, univariate and multivariate Cox regression, and the area under the receiver operating characteristic curve (AUC). In addition, the discriminative ability of the FerRLSig was further evaluated through the comparison to the Principal component analysis (PCA) and the t-distributed stochastic neighbor embedding (t-SNE).

A nomogram incorporating the FerRLSig, age, gender and TNM stage, was developed to visualize the prognostic model and facilitate its clinical application for the OS of LUSC patients. The prognostic power of the nomogram was evaluated with AUC curves at 1, 3, and 5 years’ OS, calibration curves were also plotted and visually assessed. Furthermore, decision curve analysis (DCA) was also conducted to demonstrate the actual net benefit gain of the nomogram (Vickers and Elkin, 2006).

Tumor immune microenvironment exploration and drug sensitivity screening based on the FerRLSig stratification

To analyze and compare the tumor immune microenvironment based on the FerRLSig stratification, xCell and ESTIMATE were both applied to inferred the immune and stromal cell infiltration (Yoshihara et al., 2013; Aran et al., 2017). In addition, to further characterize the potential underlying molecular pathways, gene set enrichment analysis (GSEA) was performed to identify the significantly enriched pathways in low-risk and high-risk groups respectively. The hallmark gene sets and C5 gene sets from the Gene Ontology (GO) were downloaded from the Molecular Signatures Database as the reference files (Subramanian et al., 2005). A nominal p value < 0.05 and a false discovery rate (FDR) q value < 0.25 were set as the statistically significant thresholds for GSEA GO analysis. Moreover, the expression level of several immune checkpoints and immune inhibitory factors, including CCL2, CD274 (PD-L1), CTLA4, CXCR4, IL6, LAG3, PDCD1 (PD-1), and TGFB1 were compared between low-risk and high-risk group. To further evaluate the predictive application for immune checkpoint inhibitors of the FerRLSig, TIDE score were compared between low-risk and high-risk groups (Fu et al., 2020). Furthermore, drug sensitivity screening was also performed with the 198 compounds available from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). And the half-maximal inhibitory concentration (IC50) of the available compounds on low-risk and high-risk groups was extracted with the Oncopredict R package (Maeser et al., 2021).

Statistical analysis

Discrete variables were described as counts and percentages, their differences between groups were statistically evaluated with Pearson chi-square test or Fisher’s exact test (any expected values less than 5). On the other hand, continuous variables were described as median, mean and/or interquartile range, and their differences between groups were compared with Mann–Whitney–Wilcoxon test.

All statistical analyses and visualizations were performed with R version 4.1.0 (http://www.R-project.org) and corresponding packages. The Kaplan-Meier method was utilized in survival analysis and survival curves were compared with log-rank test. Two-sided p < 0.05 was considered statistically significant.

Results

Development and validation of the FerRLSig

In total, 493 patients with RNA sequencing data were downloaded from the TCGA-LUSC database, and 473 patients with necessary survival information were included and randomly assigned into training set and testing set by the ratio of 7:3 as demonstrated in the study workflow (Figure 1). The mean age of the patients in the entire set was 67 years old, 74.2% were male, 71.0% were Caucasian and 67.7% were early stage (I-II). No statistically significant difference was found between the training set and the testing set (Table 1). In the entire set, 4,259 lncRNAs were significantly correlated with any one of the 108 ferroptosis genes (|Pearson R|>0.3 and p < 0.001, Supplementary Figure S1A), among which 43 were significantly associated with OS by univariate Cox regression in the training set (p < 0.01, Figure 2A). Within the training set, 21 ferroptosis-related lncRNAs were further identified as significant prognostic factors via LASSO regression with the λ set at lambda. min that gave the minimum mean squared error (Figures 2B,C). The final FerRLSig was established with multivariate Cox regression, consisting of 10 lncRNAs (Figure 2D). The correlation between FerRLSig lncRNAs and ferroptosis genes were visualized in both heatmap and Sankey diagram for the entire set (Supplementary Figures S1B,C). The training set was divided into high-risk group and low-risk group by the median of the FerRLSig. Compared to low-risk group, high-risk group had evidently more death event and shorter overall survival time (Figure 3A). Within the FerRLSig, RP11-65J21.3, ST3GAL5-AS1, ADAMTS9-AS2, RP5-940J5.8, RP11-535M15.1, and RP1-32I10.10 were over-expressed in the high-risk group with positive coefficients and classified as risk promoters. On the other hand, RP11-1085N6.3, KB-1836B5.4, LCMT1-AS1, and LINC01426 were over-expressed in the low-risk group with negative coefficients and classified as risk inhibitors (Figure 3A; Supplementary Figure S1C). And the final FerRLSig Formula equal to 0.305*(RP11-65J21.3) + 0.934*(ST3GAL5-AS1) + (−1.366)*(RP11-1085N6.3) + 1.942*(ADAMTS9-AS2) + 1.674*(RP5-940J5.8) + (−0.331)*(KB-1836B5.4) + (−1.229)*(LCMT1-AS1) + 0.404*(RP11-535M15.1) + (−0.572)*(LINC01426) + 1.365*(RP1-32I10.10). The high-risk group was significantly associated with worse OS compared to low-risk group in the training set (HR = 3.345, p < 0.001, Figure 3C). For validation, the same model was applied to the testing set and the entire set, all FerRLSig lncRNAs demonstrated similar expression profiles in the high-risk and the low-risk groups. In addition, the high-risk group in both testing and entire set was also significantly associated with worse OS compared to low-risk group (testing set HR = 2.290, p < 0.001, Figures 3B,D; entire set HR = 2.606, p < 0.001, Supplementary Figure S2).

FIGURE 1

FIGURE 1. Flow chart of this study; *FerRLSig: the ferroptosis-related lncRNAs signature.

TABLE 1

TABLE 1. The baseline characteristics of LUSC patients in TCGA database.

FIGURE 2

FIGURE 2. Development of prognostic ferroptosis-related lncRNAs signature for LUSC patients; (A) Forest plot of the 43 selected lncRNAs significantly associated with overall survival by univariate Cox regression analysis; (B) The coefficient profile of 21 OS-related lncRNAs chosen by LASSO regression; (C) The mean-squared error curve with different tuning parameters (logλ) and perpendicular dotted lines were drawn at the logλ corresponding to the minimum mean squared error (lambda.min) and the most regularized model within one standard error of the minimum mean squared error (lambda.1-se); (D) Forest plot of the 10 final selected lncRNAs by multivariate Cox regression analysis.

FIGURE 3

FIGURE 3. Overall survival analysis and validation of the ferroptosis-related lncRNAs signature; (A–B) Distribution of risk score, OS time, OS status and heatmap of the 10 prognostic ferroptosis-related lncRNAs signature in the TCGA-LUSC training set (A) and TCGA-LUSC testing set (B). (C–D) Kaplan-Meier survival curves of the OS of the patients in the high- and low-risk groups for the TCGA-LUSC training set (C) and TCGA-LUSC testing set (D).

Evaluation of the FerRLSig and the establishment of prognostic nomogram with the FerRLSig

To compare the whole expression profile and the FerRLSig, both PCA and t-SNE were applied to the RNA sequencing of the entire set and annotated with the FerRLSig risk stratification. The high-risk group and low-risk group demonstrated distinctly different distribution in both PCA and t-SNE, indicating that the FerRLSig risk stratification recapitulated the major variability of the TCGA-LUSC RNA sequencing (Figures 4A,B). The correlations of the FerRLSig risk score with clinical characteristics including age, gender, and TNM stage were explored and no statistically significant correlation was found (Supplementary Figure S3). Univariate and multivariate Cox OS analysis were further performed with FerRLSig and other universal clinical characteristics. Both TNM stage and FerRLSig (univariate: HR = 2.281, p < 0.001; multivariate: HR = 2.240, p < 0.001) demonstrated significant prognostic effect in both univariate and multivariate Cox regression (Figures 4C,D). In addition, time-dependent AUC of the FerRLSig and other clinical characteristics were plotted, and the FerRLSig demonstrated consistently higher AUC compared to other clinical characteristics, including TNM stage (5-years AUC: FerRLSig 0.756, TNM stage 0.607, Figures 4E,F).

FIGURE 4

FIGURE 4. Evaluation of the ferroptosis-related lncRNAs signature (FerRLSig) in the entire set of TCGA-LUSC; (A) Principal component analysis of TCGA-LUSC RNA sequencing annotated with the FerRLSig stratification; (B) t-distributed stochastic neighbor embedding (t-SNE) analysis annotated with the FerRLSig stratification; (C) Univariate Cox regression overall survival analysis of the FerRLSig score and universal clinical characteristics; (D) Multivariate Cox regression overall survival analysis of the FerRLSig score and universal clinical characteristics; (E) Time-dependent area under curve plot of the risk score and clinical characteristics. (F) Receiver operating characteristic (ROC) curves of the universal clinical characteristics and risk score of the 5-years overall survival.

To further demonstrate the model applicability in different population, subgroup OS analysis was performed in different age, gender, ethnicity and stage groups. The high-risk group was consistently associated with worse OS in all subgroups, which not only validated the model’s wide applicability in different population, but also demonstrated its independent prognostic role (Figure 5).

FIGURE 5

FIGURE 5. Subgroup overall survival analysis of the FerRLSig stratification by different age, gender, ethnicities and TNM stage subgroups in the TCGA-LUSC entire set.

Moreover, the FerRLSig and other universal clinical characteristics were incorporated into a prognostic nomogram to better predict the one, three and 5 years’ OS probabilities (Figure 6A). Overall, the nomogram demonstrated satisfactory AUC on 1, 3, and 5 years’ OS and calibration (AUC: 1-year 0.668, 3-years 0.761, 5-years 0.779, Figures 6B,C). In addition, the nomogram also demonstrated net benefit gain in decision curve analysis compared to both “intervention to none” and “intervention to all” in both 3 and 5 years’ OS prediction (Figure 6D).

FIGURE 6

FIGURE 6. Development and evaluation of the prognostic nomogram; (A) A clinical prognostic nomogram was developed to predict the 1-, three- and 5-years overall survival (OS) probability. A vertical line between each variable and point scale can be drawn to determine the points for each variable, then all the points are summed up as the total points, and the predicted overall survival rate of the 1-, three- and 5-years were calculated by drawing a vertical line from the total points scale to the 1-, three- and 5-years survival scales; (B) Area under curve plot of the nomogram for the 1-, three- and 5-years OS; (C) Calibration curves of the nomogram for 1-, three- and 5-years overall survival: nomogram-predicted overall survival is plotted on the x-axis, actual overall survival is plotted on the y-axis, a plot along the 45-degree line indicates a satisfactory model in which the predicted probabilities are identical to the actual outcomes; (D) Decision curve analysis demonstrating the clinical benefit gain of the nomogram for the three- and 5-years OS: the y-axis measures the net benefit, which is calculated by summing the benefit (true positives) and subtracting the harms (false positives). The solid line indicates the prognostic model, and the two other lines indicate the “intervention for all” (dotted line) and “intervention for none” (black line). A model is considered of clinical value if it has a higher net benefit than other models at any given threshold.

The correlation between the FerRLSig and tumor immune microenvironment

To explore the correlation between the FerRLSig and the tumor immune microenvironment, immune cell infiltration was inferred and compared between the high-risk and low-risk group with xCell analysis. The two groups exhibited apparently different tumor immune microenvironment. The high-risk group demonstrated significantly higher dendritic cells, B cells, class-switched memory B cells, CD8⁺ T cells, and multiple myeloid cells infiltration while the low-risk group had significantly higher pro B cells, Th1 cells and Th2 cells infiltration (Figures 7A,B). Moreover, both immune score and stromal score were significantly higher while tumor purity was significantly lower in high-risk group, confirming higher immune cell infiltration in high-risk group (Figures 7C–E). Furthermore, hallmark gene set enrichment analysis was performed to further characterize the different molecular pathways activated by different risk groups. Multiple immune-related hallmarks were significantly enriched in high-risk group, including complement, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, inflammatory response, interferon-alpha response, interferon-gamma response, TGF-BETA signaling and TNFA signaling via NFKB were enriched in high-risk group. On the other hand, low-risk group demonstrated multiple proliferation-related hallmarks and DNA damage repair hallmarks (Figures 7F, G). In addition, Gene Ontology pathway enrichment analysis also corroborated the above finding that high-risk group demonstrated multiple significantly enriched immune-related biological process and molecular function, while low-risk group was correlated with multiple proliferation and transcription ones(Figures 7H, I).

FIGURE 7

FIGURE 7. Evaluation of the tumor immune microenvironment based on the ferroptosis-related lncRNAs signature stratification; (A) The volcano plot depicted the different cellular landscapes of the tumor based on the ferroptosis-related lncRNAs signature (FerRLSig) stratification by the xCell analysis, the infiltration of green dot cells were significantly higher while red dot cells were significantly lower in low-risk group compared to high-risk group. (B) The boxplots depicted the significantly different infiltration of immune cell based on the FerRLSig stratification. (C–E) The boxplots compared the immune scores (C), stromal scores (D) and tumor purity (E) based on the FerRLSig stratification by the ESTIMATE analysis. (F–G) GSEA results demonstrated the differential gene set enrichment in Hallmark with high-risk (F) and low-risk (G) group. (H–I) GSEA results demonstrated the differential gene set enrichment in C5 of biological process (BP), cellular component (CC), and molecular function (MF) in high-risk (H) and low-risk (I) groups based on the FerRLSig stratification. (J) Comparison of the immune checkpoints and immune inhibitory factors, including CCL2, CD274, CTLA4, CXCR4, IL6, LAG3, and PDCD1 and TGFB1based on the FerRLSig stratification in violin plots and boxplots. (K) Comparison of TIDE score based on the FerRLSig stratification in boxplots. *: p < 0.05, **: p < 0.01, ***: p < 0.001.

Apparently, the FerRLSig was strongly correlated with tumor immune microenvironment and high-risk group demonstrated increased immune activities compared to low-risk group. We continued to compare the expression level of immune checkpoints and immune inhibitory factors between the two groups to evaluate the potential predictive application of the FerRLSig on the current immune checkpoint blockade therapy. We found that the high-risk group was significantly associated with higher expression level of immune checkpoints and immune inhibitory factors, including CCL2, CTLA4, CXCR4, IL6, LAG3, PDCD1, and TGFB1, indicating tumor immune evasion (Figure 7J). We utilized the TIDE to further evaluate the tumor immune microenvironment of different FerRLSig groups. And surprisingly, although high-risk group consistently demonstrated higher effector T cell signatures including IFNG (interferon gamma) and CD8⁺ T cell infiltration, the T cell dysfunction signature were significantly higher while microsatellite instability score was significantly lower in high-risk group. This ultimately led to significantly higher overall TIDE score in high-risk group, indicating resistance to immune checkpoint inhibitors (Figure 7K).

Drug sensitivity screening based on the FerRLSig

Besides immune checkpoint inhibitors, we also aimed to evaluate the predictive application of the FerRLSig on other drugs available from GDSC database. Therefore, drug sensitivity screening with GDSC was performed based on the FerRLSig. We identified 119 compounds from GDSC to have statistically significant different half-maximal inhibitory concentration (IC50) based on the FerRLSig stratification (Supplementary Table S1). Notably, low-risk group was significantly more sensitive to platinum and taxane compared to high-risk group (Figures 8A–C), which are the backbones for LUSC chemotherapy, therefore might account for the superior OS of the low-risk group. On the other hand, high-risk group was seemingly more intractable with fewer clinically available systemic therapies compared to low-risk group. Three representative drugs with significantly lower IC50 in high-risk group compared to low-risk group were identified, targeting WNT signaling, MAPK signaling and PI3K signaling pathway (Figures 8D,E).

FIGURE 8

FIGURE 8. Drug sensitivity screening based on the ferroptosis-related lncRNAs signature (FerRLSig) stratification with the Genomics of Drug Sensitivity in Cancer (GDSC) database; (A) Boxplot of the cisplatin half-maximal inhibitory concentration (IC50) based on the FerRLSig stratification; (B) Boxplot of the paclitaxel IC50 based on the FerRLSig stratification; (C) Boxplot of the docetaxel IC50 based on the FerRLSig stratification; (D) Boxplot of the SB216763 (targeting WNT signaling pathway) IC50 based on the FerRLSig stratification; (E) Boxplot of the Selumetinib (targeting MAPK signaling pathway) IC50 based on the FerRLSig stratification; (F) Boxplot of the AZD8186 (targeting PI3K signaling pathway) IC50 based on the FerRLSig stratification.

Discussion

Previous studies have successfully developed prognostic ferroptosis-related lncRNAs signatures in lung adenocarcinoma (Guo et al., 2021; Lu et al., 2021; Zheng et al., 2021). However, the investigation of the ferroptosis-related lncRNAs signature in lung squamous cell carcinoma (LUSC) has been scare and has not evaluated the signature’s impact on tumor immune microenvironment yet. In this study, we developed and validated a ferroptosis-related lncRNAs signature (FerRLSig) for the prognosis stratification of lung squamous cell carcinoma (LUSC). High-risk group had significantly worse OS compared to low-risk group (HR = 2.240, 95%CI: 1.845–2.720, p < 0.001), which was further corroborated in different age, gender, ethnicities and TNM stages subgroups, indicating the wide applicability and independent prognostic effect of the FerRLSig. And notably, compared to TNM stage, the FerRLSig demonstrated consistently improved AUCs (5-years AUC: FerRLSig 0.756, TNM stage 0.607, Figures 4E,F) on OS. Thus through this study, we have developed and validated a robust prognostic ferroptosis-related lncRNAs signature for LUSC.

A previous retrospective study utilizing TCGA database identified 29 ferroptosis-related lncRNAs with univariate Cox regression and constructed a prognostic ferroptosis-related lncRNAs. The 1-, 2-, and 3-years area under curve (AUC) of their signature were 0.658, 0.693, and 0.687 respectively (Yao et al., 2022). In our study, the TCGA-LUSC patients were randomized into a training set and a testing set by the “caret” R package at the ratio of 7:3. The FerRLSig was established with the training set while the validation was performed with the testing set and the entire set. Univariate Cox regression, Least Absolute Shrinkage and Selection Operator (LASSO) regression and multivariate Cox regression were applied in order to establish the final FerRLSig while avoiding overfitting with cross-validation. A more concise ferroptosis-related lncRNAs signature comprising 10 ferroptosis-related lncRNAs with an AUC of 0.756 for 5-years OS was established and validated. Compared to the previous study, we applied more stringent statistical methods to identify ferroptosis-related lncRNAs to avoid overfitting. Moreover, we underwent additional internal validation, which was absent in the previous study. Therefore, we believed our FerRLSig to be more statistically stringent with better prognostic effect compared to the previous study.

More importantly, we have also evaluated the correlation between the FerRLSig and the tumor immune microenvironment. Generally, high-risk group demonstrated significantly higher immune cell infiltration in the xCell analysis, notably by dendritic cells, CD8⁺ T cells, M1 macrophages and M2 macrophages. On the other hand, low-risk group demonstrated significantly higher pro B cells, Th1 cells and Th2 cells infiltration. In addition, the ESTIMATE and GSEA analysis also corroborated previous results that high-risk group had significantly higher immune and tumor scores with multiple immune-related gene sets enrichment compared to low-risk group, including complement, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, interferon-alpha response, interferon-gamma response, TGF-BETA signaling and TNFA signaling via NFKB. CD8⁺ cytotoxic T cells and its secreted interferon gamma are central to the tumor immune elimination, and increased immune cell infiltration indicates prominent immune response, but these do not necessarily lead to better tumor control or survival (Gajewski et al., 2013; Mojic et al., 2017; Morad et al., 2021). In addition, multiple immune inhibitory factors seen in high-risk group, including M2 macrophages and IL6, might render the infiltrating CD8⁺ T cell dysfunction and leading to immune evasion (Xue et al., 2014; Kumari et al., 2016). To further evaluate the mechanism of the immune evasion based on the FerRLSig stratification, several immune checkpoints and immune inhibitory factors including CCL2, CD274 (PD-L1), CTLA4, CXCR4, IL6, LAG3, PDCD1 (PD-1), and TGFB1 were compared between low-risk and high-risk group, and all were significantly higher in high-risk group except PD-L1, strongly suggesting the immune evasion and inhibitory microenvironment in high-risk group.

Immune checkpoint inhibitors have made remarkable breakthroughs in LUSC and significantly improved the LUSC patients’ survival (Brahmer et al., 2015; Paz-Ares et al., 2018). Considering the immune evasion and inhibitory microenvironment in high risk group based on the FerRLSig stratification, TIDE score was utilized to estimate the immune checkpoint inhibitor sensitivity. CD8⁺ T cell infiltration and interferon gamma signature were significantly higher in high-risk group while no statistically significant difference was found on CD274 (PD-L1) signature, which corroborated previous results. And high-risk group was significantly less sensitive to immune checkpoint inhibitor with significantly higher TIDE score and dysfunction score compared to low-risk group. Besides the inherent limitations of the TIDE analysis, one possible explanation would likely be that the T cell dysfunction with multiple alternative immune checkpoints including CTLA-4 and LAG3 within the tumor microenvironment is beyond the salvage of the single-target immune checkpoint inhibitor (Thommen and Schumacher, 2018). Therefore, the trials of combination immunotherapy targeting multiple immune checkpoints and further innovation are needed for the future improvement of LUSC patients. Drug sensitivity screening was also performed based on the FerRLSig with drugs available from GDSC database. Low-risk group was significantly more sensitive to platinum and taxane compared to high-risk group, which might partially account for its better OS. On the other hand, high-risk group were significantly more sensitive to three representative drugs that targeting WNT signaling, MAPK signaling and PI3K signaling pathways. And intriguingly, all three pathways are known to cancer immune evasion (Sumimoto et al., 2006; Dituri et al., 2011; Martin-Orozco et al., 2019) and these kinds of drugs might be combined with immunotherapy to further improve the survival for LUSC patients. However, further investigation is needed to verify these possibilities.

Several limitations are worth mentioning in this study. Firstly, this is a retrospective study from a single database, thus external validation and further prospective study are required. Secondly, several important clinical variables including the extent of resection, resection margin, comorbidities are currently unavailable, which warrant further investigation in future studies. Thirdly, the distance to actual clinical application remains long as whole transcriptome RNA sequencing for lncRNAs identification has not been easily accessible in clinical practice yet. The last but not the least, both in vitro and in vivo experiments are required to further explore the molecular mechanism underlying the ferroptosis-related lncRNAs signature.

In conclusion, a robust prognostic FerRLSig for LUSC has been developed and validated, demonstrating a close association not only with tumor immune cell infiltration, but also with T cell dysfunction and immune evasion. Further investigation and innovation are required to validate the results from our study and better improve the survival of LUSC patients based on the FerRLSig stratification.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/.

Author contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by RZ, XZ, HY, YL, YW, DZ, and LC. The first draft of the manuscript was written by RZ and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Funding

This study was funded by the National Key Research and Development Program of China (Grant number 2021YFC2500905).

Acknowledgments

The results shown here are in whole or part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.968601/full#supplementary-material

References

Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18, 220. doi:10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Barta, J. A., Powell, C. A., and Wisnivesky, J. P. (2019). Global epidemiology of lung cancer. Ann. Glob. Health 85, 8. doi:10.5334/aogh.2419

PubMed Abstract | CrossRef Full Text | Google Scholar

Brahmer, J., Reckamp, K. L., Baas, P., Crinò, L., Eberhardt, W. E. E., Poddubskaya, E., et al. (2015). Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer. N. Engl. J. Med. 373, 123–135. doi:10.1056/NEJMoa1504627

PubMed Abstract | CrossRef Full Text | Google Scholar

Cancer Genome Atlas Research Network (2012). Comprehensive genomic characterization of squamous cell lung cancers. Nature 489, 519–525. doi:10.1038/nature11404

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Kang, R., Kroemer, G., and Tang, D. (2021). Broadening horizons: the role of ferroptosis in cancer. Nat. Rev. Clin. Oncol. 18, 280–296. doi:10.1038/s41571-020-00462-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Dituri, F., Mazzocca, A., Giannelli, G., and Antonaci, S. (2011). PI3K functions in cancer progression, anticancer immunity and immune evasion by tumors. Clin. Dev. Immunol. 2011, 947858. doi:10.1155/2011/947858

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, S. J., Lemberg, K. M., Lamprecht, M. R., Skouta, R., Zaitsev, E. M., Gleason, C. E., et al. (2012). Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell 149, 1060–1072. doi:10.1016/j.cell.2012.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Frankish, A., Diekhans, M., Jungreis, I., Lagarde, J., Loveland, J. E., Mudge, J. M., et al. (2021). GENCODE 2021. Nucleic Acids Res. 49, D916–D923. doi:10.1093/nar/gkaa1087

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J. H., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22. doi:10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 12, 21. doi:10.1186/s13073-020-0721-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajewski, T. F., Schreiber, H., and Fu, Y-X. (2013). Innate and adaptive immune cells in the tumor microenvironment. Nat. Immunol. 14, 1014–1022. doi:10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossman, R. L., Heath, A. P., Ferretti, V., Varmus, H. E., Lowy, D. R., Kibbe, W. A., et al. (2016). Toward a shared vision for cancer genomic data. N. Engl. J. Med. 375, 1109–1112. doi:10.1056/NEJMp1607591

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y., Qu, Z., Li, D., Bai, F., Xing, J., Ding, Q., et al. (2021). Identification of a prognostic ferroptosis-related lncRNA signature in the tumor microenvironment of lung adenocarcinoma. Cell Death Discov. 7, 190. doi:10.1038/s41420-021-00576-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of cancer: The next generation. Cell 144, 646–674. doi:10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhn, M. (2008). Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26. doi:10.18637/jss.v028.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumari, N., Dwarakanath, B. S., Das, A., and Bhatt, A. N. (2016). Role of interleukin-6 in cancer progression and therapeutic resistance. Tumour Biol. 37, 11553–11572. doi:10.1007/s13277-016-5098-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. S., and Ruppin, E. (2019). Multiomics prediction of response rates to therapies to inhibit programmed cell death 1 and programmed cell death 1 ligand 1. JAMA Oncol. 5, 1614–1618. doi:10.1001/jamaoncol.2019.2311

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, L., Liu, L-P., Zhao, Q-Q., Gui, R., and Zhao, Q-Y. (2021). Identification of a ferroptosis-related LncRNA signature as a novel prognosis model for lung adenocarcinoma. Front. Oncol. 11, 675545. doi:10.3389/fonc.2021.675545

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeser, D., Gruener, R. F., and Huang, R. S. (2021). oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief. Bioinform. 22, bbab260. doi:10.1093/bib/bbab260

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, C., Wang, X., Liu, Y., Wang, M., Yan, B., Jiang, Y., et al. (2018). A G3BP1-interacting lncRNA promotes ferroptosis and apoptosis in cancer via nuclear sequestration of p53. Cancer Res. 78, 3484–3496. doi:10.1158/0008-5472.CAN-17-3454

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin-Orozco, E., Sanchez-Fernandez, A., Ortiz-Parra, I., and Ayala-San Nicolas, M. (2019). WNT signaling in tumors: The way to evade drugs and immunity. Front. Immunol. 10, 2854. doi:10.3389/fimmu.2019.02854

PubMed Abstract | CrossRef Full Text | Google Scholar

Mojic, M., Takeda, K., and Hayakawa, Y. (2017). The dark side of IFN-γ: Its role in promoting cancer immunoevasion. Int. J. Mol. Sci. 19, E89. doi:10.3390/ijms19010089

PubMed Abstract | CrossRef Full Text | Google Scholar

Morad, G., Helmink, B. A., Sharma, P., and Wargo, J. A. (2021). Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell 184, 5309–5337. doi:10.1016/j.cell.2021.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Paz-Ares, L., Luft, A., Vicente, D., Tafreshi, A., Gümüş, M., Mazières, J., et al. (2018). Pembrolizumab plus chemotherapy for squamous non–small-cell lung cancer. N. Engl. J. Med. 379, 2040–2051. doi:10.1056/NEJMoa1810865

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramalingam, S. S., Vansteenkiste, J., Planchard, D., Cho, B. C., Gray, J. E., Ohe, Y., et al. (2020). Overall survival with osimertinib in untreated, EGFR-mutated advanced NSCLC. N. Engl. J. Med. 382, 41–50. doi:10.1056/NEJMoa1913662

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaw, A. T., Bauer, T. M., de Marinis, F., Felip, E., Goto, Y., Liu, G., et al. (2020). First-line lorlatinib or crizotinib in advanced ALK-positive lung cancer. N. Engl. J. Med. 383, 2018–2029. doi:10.1056/NEJMoa2027187

PubMed Abstract | CrossRef Full Text | Google Scholar

Solomon, B. J., Mok, T., Kim, D-W., Wu, Y-L., Nakagawa, K., Mekhail, T., et al. (2014). First-line crizotinib versus chemotherapy in ALK-positive lung cancer. N. Engl. J. Med. 371, 2167–2177. doi:10.1056/NEJMoa1408440

PubMed Abstract | CrossRef Full Text | Google Scholar

Sordella, R., Bell, D. W., Haber, D. A., and Settleman, J. (2004). Gefitinib-sensitizing EGFR mutations in lung cancer activate anti-apoptotic pathways. Science 305, 1163–1167. doi:10.1126/science.1101637

PubMed Abstract | CrossRef Full Text | Google Scholar

Statello, L., Guo, C-J., Chen, L-L., and Huarte, M. (2021). Gene regulation by long non-coding RNAs and its biological functions. Nat. Rev. Mol. Cell Biol. 22, 96–118. doi:10.1038/s41580-020-00315-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Sumimoto, H., Imabayashi, F., Iwata, T., and Kawakami, Y. (2006). The BRAF-MAPK signaling pathway is essential for cancer-immune evasion in human melanoma cells. J. Exp. Med. 203, 1651–1656. doi:10.1084/jem.20051848

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin. 71, 209–249. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Thommen, D. S., and Schumacher, T. N. T. (2018). T cell dysfunction in cancer. Cancer Cell 33, 547–562. doi:10.1016/j.ccell.2018.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Vickers, A. J., and Elkin, E. B. (2006). Decision curve analysis: a novel method for evaluating prediction models. Med. Decis. Mak. 26, 565–574. doi:10.1177/0272989X06295361

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Mao, C., Ouyang, L., Liu, Y., Lai, W., Liu, N., et al. (2019). Long noncoding RNA LINC00336 inhibits ferroptosis in lung cancer by functioning as a competing endogenous RNA. Cell Death Differ. 26, 2329–2343. doi:10.1038/s41418-019-0304-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y-L., Tsuboi, M., He, J., John, T., Grohe, C., Majem, M., et al. (2020). Osimertinib in resected EGFR-mutated non-small-cell lung cancer. N. Engl. J. Med. 383, 1711–1723. doi:10.1056/NEJMoa2027071

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, J., Schmidt, S. V., Sander, J., Draffehn, A., Krebs, W., Quester, I., et al. (2014). Transcriptome-based Network analysis reveals a spectrum model of human macrophage activation. Immunity 40, 274–288. doi:10.1016/j.immuni.2014.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, N., Zuo, L., Yan, X., Qian, J., Sun, J., Xu, H., et al. (2022). Systematic analysis of ferroptosis-related long non-coding RNA predicting prognosis in patients with lung squamous cell carcinoma. Transl. Lung Cancer Res. 11, 632–646. doi:10.21037/tlcr-22-224

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Wang, L., Li, H., Zhang, L., Zheng, X., and Cheng, W. (2020). Crosstalk between noncoding RNAs and ferroptosis: New dawn for overcoming cancer progression. Cell Death Dis. 11, 580. doi:10.1038/s41419-020-02772-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Z., Zhang, Q., Wu, W., Xue, Y., Liu, S., Chen, Q., et al. (2021). Identification and validation of a ferroptosis-related long non-coding RNA signature for predicting the outcome of lung adenocarcinoma. Front. Genet. 12, 690509. doi:10.3389/fgene.2021.690509

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, N., and FerrDb, Bao J. (2020). FerrDb: a manually curated resource for regulators and markers of ferroptosis and ferroptosis-disease associations. Database. 2020, baaa021. doi:10.1093/database/baaa021

PubMed Abstract | CrossRef Full Text | Google Scholar