Identification of a novel necroptosis-related LncRNA signature for prognostic prediction and immune response in oral squamous cell carcinoma

2.1The data extraction

The clinical information and RNA sequences of OSCC patients’ were retrieved from the TCGA database (as of March 1, 2022, https://portal.gdc.cancer.gov/ repository). Later, RNA sequence (RNA-seq) data from the TCGA database were normalized to million fragments per thousand bases (FPKM). This study only included samples with complete clinicopathological data, which were then analyzed and preprocessed using the Strawberry Perl programming language (version -5.30.0.1; https://www.perl.org) to increase the study’s accuracy. Following that, patient data were randomized to a training or testing set using the R software “Caret” package (version of 4.1.0; https://www.r-project). GSEA and previous reports yielded a total of genes associated with necroptosis.

2.2Identification of necroptosis-related signatures (NRLs) in OSCC

Using the GENCODE website (https://www.gencodegenes.org/, as of September 3, 2021) to retrieve signature GTF annotation files, 14,056 signatures were identified from the TCGA–seq HNSC RNA data to distinguish the mRNA and signature. Pearson’s correlation analysis which was considered an accepted method was used to investigate the correlation between coding genes and lncRNAs to examine the co-expression of necroptosis-related genes (NRGs) in each signature among OSCC samples. These Signatures are thought to be significantly associated with necroptosis mRNAs (for the cutoff, | Cor |> 0.4, P< 0.001). To compare the differences in necroptosis-related gene expression in normal and tumor tissue, the “Edger” and “DEseq” R packages were used to generate the p-value for each gene, and a fold change (Fc), with the P< 0.05 and |log2 Fc |> 1 and the genes were identified as NRLs.

2.3Construction and verification of NRLs signature for OSCC

For the purpose of screening lncRNAs in relationship to prognosis, Univariate Cox regression analysis was adopted. In addition, LASSO regression analysis was used to screen for lncRNAs that were significantly associated with overall survival (OS) in OSCC patients. This study conducted a multifactorial Cox regression analysis to generate the best model. Finally, five lncRNAs associated with necroptosis were considered as prognostic factors. RS was determined using the following formula: Risk score=∑i=1nCoef⁢(i)×x⁢(i), (n : number of genes incorporated into the signature; Coef(i): the NRL’s coefficient, and x(i): the gene expression level). Subsequently, all the OSCC samples were divided into low-or high-risk subgroups according to the median RS value. Then, the R software “survival” package was used for constructing the survival curve to compare overall survival (OS) between two risk subgroups from both datasets. Images were used to depict the risk score distribution, survival status, and prognostic signature expression profile of OSCC patients.

2.4Establishment of nomogram and correlation between prognostic markers and clinicopathological signatures

To ensure that the NRLs signal was independent, univariate and multivariate Cox regression was used to assess the relation between clinical characteristics and NRLs markers with OS. Specifically, the role of NRLs markers as risk factors independent of additional clinicopathological characteristics like age, sex, stage, grade, and TNM classification was investigated. Survival analysis was conducted in the present study based on different clinicopathological signatures to confirm the applicability of the NRLs-based set. The “RMS” package used NOMO plots for independent prognostic factors employed to predict the OS rates of OSCC patients at 1, 3, and 5 years. Finally, receiver operating characteristic (ROC) curves were plotted to ascertain if the above factors were precise and specific in predicting the OSCC prognosis.

2.5Levels of tumor-infiltrating immune cells (TIICs)

The most recent technology was used to calculate the immune penetration status in TCGA samples, including QUANTISEQ (http://icbi.at/quantiseq), XCELL (http://xcell.ucsf.edu/), MCPCOUNTER, TIMER (version 2.0; http://timer.cistrome. org/), EPIC (http://epic. gfellerlab.org), CIBERSORT (http://ciber sort.stanford. edu/), and CIBERSORT-abs to determine the correlation between TIICs levels and a risk score. Moreover, this work also utilized the Wilcoxon test for assessing the difference in tumor invasive immune cell content of both risk groups. In addition, this study conducted Spearman’s correlation analysis to find a relationship between RS and TIICs levels. The R software limma, scales, GGplot2, and GGText packages are used for data visualization. In addition, a single-sample gene set concentration analysis (SsGSEA) score was determined to assess the concentration levels of immune-related functions in both groups. This work also detected differences in immune checkpoint (ICP) molecules between the two populations for evaluating the gene set in predicting immunotherapy for OSCC.

2.6Pathway enrichment

To clarify different pathways enriched in two risk groups, this work utilized Gene Set Enrichment Analysis (GSEA) software (version 4.10) for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment.

2.7Significance of NRLs-based signatures in antitumor drugs

This study used a “pRophetic” program package to evaluate the median inhibitory concentration (IC₅₀) of 138 common antitumor agents, including paclitaxel, cisplatin, and imatinib, in order to predict the response to chemotherapy agents in patients with OSCC in two different risk subgroups. (AJCC) The IC50 values between the two risk groups were analyzed by the Wilcoxon signed-rank test.

2.8Consensus clustering

For exploring potential molecular subtypes in MIBC patients,the “ConsensusClusterPlus” R package was used to sort out the optimal cluster value based on determined prognostic NLRs. After that, we used the “Rtsne” R package to accomplish principal component analysis (PCA).

2.9Cell culture

Human oral squamous cell lines CAL-27, SCC-25 and normal oral keratinocytes cell lines HOK were used in this study. The cell lines were purchased from the United States Type Culture Repository (ATCC, Manassas, USA). All cells were subcultured and stored in Shanxi Province Key Laboratory of Oral Diseases Prevention and New Materials, and tested regularly to ensure mycoplasma negative. All cells were cultured in high glucose DMEM medium (Gibco, CA, USA) containing 10% fetal bovine serum (Gibco, CA, USA) and 1% penicillin/streptomycin solution at 37^∘C and 5% CO₂.

2.10Real-time PCR

Totol RNA was extracted from cells using the M5 Universal RNA Mini Kit(Mei5 Biotechnology,China) according to the instructions. Totol RNA(Mei5 Biotechnology,China) was then reverse-transcribed into cDNA using the M5 Hiper lncRNA cDNA Synthesis Kit with gDNA remover lncRNA cDNA. Real-time quantitative PCR amplification was performed using the M5 Hiper lncRNA Fluorescence quantitative detection kit(Mei5 Biotechnology,China). Using GAPDH as reference, the mRNA relative expression was measured by 2 ^{- ΔΔCT}. Primer sequence is shown in Supplementary Table S1.

2.11Statistical analysis

Before conducting statistical analyses, the data were checked for normality to ensure they follow the assumptions required for parametric tests. The normality of the data was assessed using the Shapiro-Wilk test and by visual inspection of Q-Q plots. If the data did not conform to a normal distribution, appropriate non-parametric tests were used. This step was critical to ensure the validity of the statistical analyses applied in the study. Statistical analysis was performed using R software (version 4.1.0). Student’s t-test and one-way ANOVA were used to calculate differences between two groups or more. The association between necroptosis-related signature genes, risk scores, and clinical signatures was calculated by the chi-square test. A log-rank test with the best cutoff value was used to plot a Kaplan-Meier (KM) survival curve to analyze OS. Further, a Cox regression model (univariate or multivariate) was used to establish the relationships between various variables with clinical outcomes. The statistical significance level was P< 0.05.

3.1Identification of necroptosis-related signatures (NRLs) within OSCC

Figure 1.

Flowchart. General flow is shown in the figure.

Figure 2.

Identification of different expressions of signatures associated with necroptosis in OSCC. (A) mRNA co-expression network diagram of Signatures and necroptosis. (B–C) Volcano map and heat map of long non-coding RNA (signature) with differential expression in OSCC compared with para-carcinoma samples. (|logFC|> 1, P< 0.05).

Figure 1 depicts the flow chart created in the present work. This study retrieved transcriptome RNA-Seq data from 330 TCGA-OSCC cases, including 330 OSCC tissues and 30 nearby normal tissues, with pertinent clinical information. The present study included samples with adequate clinical data for further analysis. There were 14,056 signatures detected from the TCGA-HNSC gene expression file based on the GTF annotation file for human signatures. GSEA and previous reports yielded 67 genes (nrmRNA) associated with necroptosis (Supplementary Table S2). Pearson’s correlation analysis yielded 666 signature necroptosis-related genes (nrmRNA) with significant correlation (|cor|> 0.5, P< 0.05). Based on these results, we created an mRNA-signature co-expression network (Fig. 2A) to determine the potential impact of necroptosis-related signatures. Further, 307 signatures with abnormal expression in OSCC were identified as necroptosis-related signatures (NRLs) by differential analysis (|logFC|> 1, P< 0.05), and gene expression heat map and volcano map were drawn (Fig. 2B,C).

Figure 3.

NRLs prognostic gene set establishment. (A) Forest map showing prognostic signatures. 8 OS-related signatures among OSCC cases were screened through univariate Cox proportional hazards regression (P< 0.05). (B) 10-fold cross-validation. (C) Changes in different gene coefficients. (D) Calorigram of NRLs expression in 330 OSCC tissues and normal tissues. (E) Sankey diagram showing the regulatory relationship and connection degree of NRLs with NRGs.

To identify lncRNAs with prognostic significance, single-factor Cox analysis was adopted for detecting 8 NRLs (SLC16A1-AS1, GAS1RR, AC099850.3, STARD4-AS1, AC011978.1, LINC01503, CDKN2A-DT, LINC00973) (P< 0.05, Fig. 3A) (|cor|> 0.4 and P< 0.001) (Fig. 3A). Subsequently, LASSO regression analysis detected 8 NRLs (Fig. 3B,C). The heatmap depicting the expression levels of the eight signatures across tumor and non-tumor tissues was generated (Fig. 3D). The Sankey diagram demonstrated the interactions between NRLs and NRGs and showed a positive correlation between OS in OSCC patients and the prognostic gene sets of the 8 NRLs (Fig. 3E). Multifactor Cox analysis further determined 5 NRLs (AC099850.3, StarD4-AS1, AC011978.1, LINC01503, and CDKN2A-DT) with prognostic significance.

3.2Construction of NRLs signatures

To assess the prognostic risk of OSCC patients, a risk model was established with eight cuproptosis-related lncRNAs. Each OSCC patient in the TCGA database was set a risk score based on the formula: RS = (0.406300575347132 × AC099850.3) + (-1.49039102223428 × StarD4-AS1) + (1.39186509 330249 × AC011978.1) + (0.288732320378307 × LINC01503) + (2.26941578067647 × CDKN2A-DT). The cases were classified as a high- or low-risk group based on the median RS.

3.3Analysis and validation of risk score based on 5-NRLs signature

Figure 4.

Analysis and verification of necrotic apoptosis-related signature risk score in OSCC patients. (A–C) KM analysis on low-and high-risk OSCC cases from all samples, Train group, and Test group. (D–F) Distribution of risk scores. (G–I) Scatter plot showing the relation of survival time with RS. (J–L) Heat map of 5NRLs expression profile in OSCC patients between low-and high-risk groups.

Figure 4.

continued. (M) 1-, 3-and 5-year t-ROC curves. (N) Analysis of receiver operating characteristic curve of risk score. Relative to additional clinicopathological signatures, our signature-based RS showed increased AUC.

Kaplan-Meier survival analysis was used to evaluate the independent prognostic performance of the 5-NRLs prognostic gene set among OSCC cases. The high-risk group exhibited a significantly lower overall survival (OS) compared to the low-risk group (P< 0.001, Fig. 4A). Additionally, the high-risk group had lower survival rates, higher mortality rates (Fig. 4D), a greater number of cases and death events as the risk score (RS) increased (Fig. 4G), and abnormal expression in tumor tissues (Fig. 4J). To validate the accuracy of RS, all samples were randomized into a training group (n= 161) and a testing group (n= 161), and further analysis yielded consistent results with the initial findings (Fig. 4A–L).

Time-dependent ROC (t-ROC) curves for 1-, 3-, and 5-year survival were plotted to evaluate RS’s predictive performance, with the area under the curve (AUC) measuring 0.719, 0.646, and 0.664, respectively. The first-year AUC value was higher than those for the third and fifth years, indicating that the signature was particularly accurate for predicting OS in OSCC (Fig. 4M). Moreover, compared to age (0.579), gender (0.506), tumor pathological grade (0.561), and tumor stage (0.555), the AUC score for OS (0.719) was considerably higher (Fig. 4N), underscoring the utility of the constructed signature in predicting OSCC prognosis.

3.4Correlation of 5-NRLs prognostic markers with clinicopathological parameters and their independent prognostic value

Figure 5.

An Independent assessment of the association of NRLs markers with patient outcomes. (A) Univariate Cox regression showing the significant relationship between the survival and signature-based RS (P< 0.001), age (P= 0.009), stage (P< 0.001), and grade (P= 0.043). (B) Multivariate Cox regression identified signature-based RS (P< 0.001), age (P= 0.008), and stage (P< 0.001) as factors independently associated with patient survival. The risk score independently predicted OSCC survival in the TCGA database. ^**** P< 0.0001, ^*** P< 0.001, ^** P< 0.01, ^*P< 0.05. (C–D) The KaplanMeier curve showed that patients in different stage with high risk displayed a shorter overall survival than those with low risk. (E) A Nomogram illustrating the relationships between RS and clinicopathological signatures.

To further evaluate the role of signatures associated with necrotic apoptosis in the development of OSCC, univariate Cox regression analysis demonstrated that the risk score (RS) independently predicted OSCC prognosis (hazard ratio [HR] = 1.577, 95% confidence interval [CI] = 1.367–1.820, P< 0.001). Additionally, multivariate Cox regression indicated that even when other confounding factors were adjusted, RS remained a significant predictor of OSCC prognosis (HR = 1.502, 95% CI: 1.303–1.732, P< 0.001) (Fig. 5A, B). Furthermore, the low-risk group showed longer overall survival (OS) compared to the high-risk group in both stages I and II (P= 0.002, Fig. 5C) and in stages III and IV (P< 0.001, Fig. 5D).

Additionally, by adopting a stepwise Cox regression model, we combined the necrotic apoptosis-related RS with clinicopathological characteristics to create a clinically adaptive nomogram to estimate and quantify the 1-, 3-, and 5-year survival probabilities for OSCC patients (Fig. 5E). These findings suggest that this signature can be used for prognostic prediction.

3.5Functional analysis based on the risk model

Figure 6.

5-GSEA pathway enrichment analysis of NRLs. (A) The top 5 enriched pathways for both risk groups (B) Enriched pathways for the high-risk group. (C) Enriched pathways for the low-risk group.

This study used GSEA to perform KEGG pathway enrichment to elucidate different pathways enriched into high-and low-risk groups. As a result, the five most significantly enriched pathways in the high-risk group included RNA polymerase, pyrimidine metabolism, spliceosome, ribosome, and RNA degradation (Fig. 6A), as well as multiple oncogenic pathways, such as oxidative phosphorylation, P53 signaling pathway, mutual conversion of pentose and glycoside, starch and sucrose metabolism (Fig. 6B). Forty-four pathways, were markedly related to the low-risk subset. Among them, the top five were arrhythmogenic right ventricular cardiomyopathy, aldosterone-regulated sodium reabsorption, dilated cardiomyopathy, vascular smooth muscle contraction, and transendothelial migration of white blood cells (Fig. 6A). B-cell receptor signaling, T-cell receptor signaling, antigen processing and presentation, and natural killer cell-mediated cytotoxicity were associated with the low-risk group (Fig. 6C). Supplementary Table S3 and Supplementary Table S4 shows all the enrichment pathways for both risk groups.

3.6TME features and TIICs for both risk groups based on 5-NRLs

Figure 7.

TME features and TIICs between two risk subpopulations based on 5-NRLS. (A) Study on the correlation of TIICs with RS. (B) TIICs in high and low-risk subsets. (C) ICP in high and low-risk subsets. (D–F) The difference in ESTIMATE score,immune score and stomal score between low-risk and high-risk groups. (G) Comparison analysis between two risk groups according to the ICP gene expression profiles (^***P< 0.001; ^**P< 0.01; ^*P< 0.05).

Figure 7.

continued.

The TME comprises stromal cells, immune cells, and tumor cells. The present study used the “ESTIMATE” package to measure the significant difference in TME between two risk groups and calculate expression levels of 60 common TIICs. As a result, the infiltration levels of 10 common TIICs showed positive relations to RS, like non-regulatory T cell CD4 + (COR = 0.31, P= 8.59E-09), helper T cell CD4 + (COR = 0.28, P= 3.23E-07), eosinophils (COR = 0.15, P= 0.004), etc. The invasion level of 50 infiltrating immune cells, such as aDC, was statistically significant and negatively correlated with RS (Fig. 7A). The 14 immune cell infiltrates (aDC, B cells, CD8+ T cells, mast cells, DC, iDC, NK cells, neutrophils, pDC, helper T cells, Tfh, Th1, Th2, Til, and Tregs) showed significant differences between two risk subsets (Fig. 7B). In addition, this study performed ssGSEA to assess the 47 immune-related function enrichment degrees of both risk groups. As a result, nine ICPs (including APC costimulation, CD8+ T-cell checkpoint, CCR, cytolytic activity, type II IFN response, proinflammatory response, T-cell costimulation, and T-cell coinhibition) of both subpopulations varied significantly (Fig. 7C). Based on the TME score, the low-risk subgroup showed a higher TME score than the high-risk subgroup (Fig. 7D–F), which indicates immunosuppressive TME may be partially responsible for the low survival rate in the high-risk subgroup. In addition, box plots were used to highlight the differences in the levels of 26 common TIICs between the two risk groups (Fig. 7G).

3.7The importance of NRLs in antitumor drugs

Figure 8.

Sensitivity analysis of 24 anticancer drugs among high-risk OSCC cases.

Figure 8.

continued.

OSCC exhibits relatively low sensitivity to various antitumor drugs, limiting their widespread use in clinical practice. We used the “pRRophetic” algorithm to evaluate the half-maximal inhibitory concentration (IC50) for different chemotherapy agents in two risk groups to tailor chemotherapy accordingly. The results indicated significant differences in response to 24 anticancer drugs between the groups (Fig. 8). The low-risk subgroup showed higher sensitivity to 13 anticancer drugs, including ABT.263, AKT inhibitors, AZD6482, CCT007093, DMOG, GDC0941, Imatinib, LFM.A13, MK.2206, PAC.1, Temsirolimus, WO2009093972, and Z.llnle.CHO, suggesting that low-risk patients might benefit from these chemotherapy agents (Table 1).

Conversely, the IC50 values for drugs such as Bi.d1870, BIRB.0796, Doxorubicin, Epothilone.b, Erlotinib, GW.441756, Obatoclax.mesylate, PD.0325901, SL.0101.1, Sorafenib, and VX.680 were higher among high-risk patients. This might indicate that these drugs are more suitable for cases with an elevated risk score based on the 5-NRLs characteristics (Table 2).

3.85-NRLs cluster analysis and the relationship between clusters and overall survival and tumor microenvironment

Figure 9.

5-NRLs cluster analysis and the relationship between clusters and overall survival and tumor microenvironment. (A) Clustering sample distribution. (B) Four clusters and the relationship between the overall survival rate. (C) The sankey diagram, each subgroup and the relationship between the high and low risk distribution. (D–J) PCA and T-SNE distribution in each subgroup and the relationship between the tumor microenvironment score. (K) Each subgroup immune cells expressing quantity heat map. (L) Immune checkpoint expression in each subgroup analysis of the differences. (M–N) Sensitivity analysis of PD.0325901 and PD.0332991 among clusters.

Figure 9.

continued.

To explore the relationship between genes associated with necrotic apoptosis and OSCC subtypes, we performed an unsupervised and consistent cluster analysis based on the expressions of necrotic lncRNAs by changing the cluster variable (k) from 2 to 9 through Consensus ClusterPlus R package (Fig. 9A). The distribution diagram of cluster samples with k= 2–9 is shown in Supplementary Fig. S1, and the classification table of each sample is shown in Supplementary Table S5. The results showed that when k= 4, the intra-group correlation was the highest, while the inter-group correlation was low. Samples in the TCGA data set could be separated stably, indicating that OSCC patients could be well divided into four clusters. We also show the classification of each sample under different cluster numbers (k). Therefore, we divided the data sets for OSCC patients into clusters C1, C2, C3, and C4. Through survival analysis, it can be seen that the overall survival of patients in four different clusters is different (P= 0.005) (Fig. 9B). The Sankey chart shows the relationship between four clusters and high and low risk of cancer (Fig. 9C), with most of the C2 subgroup belonging to the low risk group and most of the C3 subgroup belonging to the high risk group. Through PCA and T-SNE analysis, it can be observed that C1, C2, C3 and C4 samples can be separated according to the expression level of 5-NRLs, that is to say, samples of different clusters can be distinguished according to the expression level of 5-NRLs (Fig. 9D–G). Next, to determine the relationship between 5-NRLs and tumor microenvironment in OSCC, we observed differences in tumor microenvironment scores between C1 and C2, and between C3 and C1, C2, and C4 (Fig. 9H), and partial differences in immune cells among different subtypes(Fig. 9I)., with the highest number of immune cells in C3. There were some differences in stromal cells among different types, and stromal cells had the highest content in C3 (Fig. 9J). Type to intuitive see different points between the expression of immune cells, through seven different prediction software for drawing heat maps (Fig. 9K). In addition, the expression of 30 immune checkpoint genes differed among different genotypes (Fig. 9L). Next, the semi-inhibitory concentration levels of anti-tumor drugs on different subgroups were analyzed, and it was found that 93 antitomor drugs had significant differences among different subgroups (Supplementary Table S6). Among them, PD.0332991 and PD.0325901 showed significant differences among multiple subgroups (Fig. 9M, N).

3.9Verification of 5- NRLs expression levels in OSCC

Figure 10.

The expression of 5-NRLs in oral normal keratinoid epithelial normal cell lines and OSCC cell lines was determined by RT-qPCR. (A) AC011978.1; (B) AC099850.3; (C) LINC01503; (D) CDKN2A-DT; (E) STARD4-AS1; (^*P< 0.05, ^**P< 0.01 and ^***P< 0.001, ns: non-significant).

In order to further explore the expression of 5-NRLs, we used RT-qPCR analysis. The data were analyzed using one-way ANOVA to compare the expression levels across the human normal oral keratinocyte line (HOK), SCC25, and Cal27 cells. The results showed that compared with HOK, the expression levels of AC099850.3, AC011978.1, and LINC01503 were up-regulated in SCC25 and Cal27 cells (Fig. 10A–C). However, the expression levels of STARD4-AS1 and CDKN2A-DT were down-regulated in OSCC cell lines (Fig. 10D–E). Overall, these experimental results further validate the reliability of our established risk model and suggest that 5-NRLs may play an important role in the occurrence and progression of OSCC.