Prognostic value and immunotherapy analysis of immune cell-related genes in laryngeal cancer

Data sources

RNA-seq data of LC were extracted from The Cancer Genome Atlas-Head and Neck Squamous Cell (TCGA-HNSC) cohort (https://portal.gdc.cancer.gov/), comprising 112 LC samples and 11 normal tissue samples. After excluding samples with incomplete survival or clinical data, 83 LC samples were retained for subsequent analyses. Additionally, the GSE27020 dataset, serving as an external validation set, was sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo). This dataset includes expression profiles from 59 LC and 50 normal samples, along with transcriptome and disease-free survival (DFS) data and clinical features. Moreover, the GSE51985 dataset, utilized for model gene expression validation, contains the expression profiles of 10 LSCC tissue samples, also obtained from GEO.

Identification of DEGs and DE-ICs

Differential expression analysis was first conducted to identify differentially expressed genes (DEGs) between the 11 normal and 112 LC samples using the limma R package (version 3.44.3) (Ritchie et al., 2015), with a statistical significance threshold of P-value <0.05 and |Log2FC| >1. Single-sample Gene Set Enrichment Analysis (ssGSEA, version 1.03) was then performed to compare the proportions of 29 immune cell types between LC and normal samples, calculating the immune score for each cell type. Differentially expressed immune cells (DE-ICs) were defined based on these immune scores. DE-ICs were categorized into high and low expression groups according to the median immune score, and survival probabilities for each DE-IC group were compared.

Weighted gene co-expression network analysis

For further analysis, the goodSamplesGenes function in the weighted gene co-expression network analysis (WGCNA) R package was used to filter samples with a TRUE status, and immune cell infiltration levels with significant survival differences were incorporated as WGCNA trait data to identify co-expressed immune-related genes in LC. In the network construction process, a soft thresholding power β was selected as the lowest power. Modules were defined by the dynamic cutting method, with a minimum of 80 genes per module. Similar gene expression profiles were clustered, and modules were visualized in a tree diagram with a MEDissThres setting of 0.2. To identify modules associated with LC immune characteristics, a heatmap of module-feature relationships was constructed, displaying correlation coefficients and P-values. Modules with P < 0.05 and |cor|>0.3, associated with immune-related phenotypes (e.g., parainflammation, Tfh, NK cells), were selected as modules of interest, and genes within these modules were considered module genes.

Identification and functional enrichment of key genes

Overlap analysis was performed between DEGs and module genes, with the intersecting genes regarded as key genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were subsequently conducted on the key genes using the clusterProfiler R package (version 4.0.2; R Core Team (2020)), with a significance threshold of P-value <0.05.

Construction and validation of risk model

The 83 TCGA-LC samples were randomly divided into a training set (59 samples) and an internal validation set (24 samples) at a 7:3 ratio, utilizing the Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) for key gene expression levels. The risk model was constructed through univariate Cox analysis followed by Least Absolute Shrinkage and Selection Operator (LASSO) regression, with 10-fold cross-validation. Specifically, univariate Cox analysis was conducted to identify survival-related genes (P < 0.05) using the survival package (version 3.1-12) (Friedman, Hastie & Tibshirani, 2010). These survival-related genes were then input into LASSO logistic regression analysis using the glmnet package (version 4.1-1) (Borgan, 2001), with family set to binomial and type.measure set to class, for the selection of prognostic biomarkers.

To assess the prognostic value of the risk model, risk scores for each patient with LC were calculated by multiplying the risk coefficients from the LASSO model with the expression levels of the respective genes, following the formula: RiskScore = (−0.69348) × RENBP + 0.390337× OLR1. Risk curves, scatter plots, and heat maps of model gene expression were generated based on the risk scores. Overall survival curves were plotted for patients with LC in both high- and low-risk groups, determined by whether their risk scores exceeded the median, using the survminer package (version 0.4.6). The model’s efficacy was further assessed by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curves. ROC curves for 1, 3, and 5-year survival time nodes were generated using the survival ROC package (version 3.1-12). Additionally, the same evaluation procedures were applied to both internal and external validation sets (GSE27020) to further validate the model’s effectiveness.

Establishment of a nomogram

To assess the prognostic value of the risk model in relation to clinical factors, a correlation analysis was performed between the clinicopathological characteristics (Age, T, N, M, Gender, and Stage) of the 83 TCGA-LC samples and the model genes, with the results visualized in a heatmap. Survival differences across subgroups of clinicopathological characteristics (M: M0, MX; Age: >60 and ≤60; T: T1, T2, T3, T4; Gender: Male, Female; lymph node metastasis status: N, N0 indicating no lymph node metastasis, N1, N2, and N3 indicating varying degrees of metastasis, and NX indicating undefined. Stage: Stage I, II, III, IV) were examined between the two risk groups. Subsequently, univariate Cox regression analysis was conducted to evaluate the prognostic impact of clinicopathological factors and the risk score. Clinicopathological variables with P < 0.05 were selected for multivariate Cox regression analysis to identify independent prognostic factors.

A nomogram combining the risk model with independent prognostic factors was then developed using the rms package (version 6.1-0). The discriminatory power of the nomogram was assessed using the concordance index (C-index), and its predictive accuracy was evaluated through calibration curves and ROC curve AUCs.

Expression analysis of model genes

The expression levels of the model genes were extracted from the TCGA-LC dataset and compared with those in normal samples. To validate these expression levels, the gene expression levels were further validated in the GSE51985 dataset.

Functional enrichment analyses of high- and low-risk groups

To analyze the impact of high- and low-risk groups on LC progression, hallmark and KEGG (C2) gene sets were downloaded from MSigDB (Broad Institute), and Gene Set Variation Analysis (GSVA) was performed for all genes in the high- and low-risk groups of TCGA-LC. The aim was to identify common functions and related pathways associated with the gene sets. For hallmark enrichment, the criteria |NES|>1 and FDR <0.05 were applied, while for KEGG pathways, the thresholds of |NES|>1 and FDR <0.25 were used to define significant enrichment.

TME analysis

The Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm was applied to compare stromal scores, immune scores, and ESTIMATE composite scores between the high- and low-risk groups.

Additionally, the ssGSEA algorithm was utilized to calculate the proportions of immune cells in each sample, and the ssGSEA score comparison between the risk groups was visualized using violin plots and heatmaps.

To further evaluate the potential immunotherapy effects in LC, the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was applied to predict the response to immune checkpoint inhibitors (ICIs). The differences in immunotherapy sensitivity between the high- and low-risk groups were assessed. TIDE scores, including TIDE, IFNG, Merck18, CD8, Dysfunction, and Exclusion scores, were compared between the groups using rank-sum tests.

Additionally, the correlation between TIDE scores and risk score was analyzed. Finally, the expression levels of PD-1 were extracted and the differences between high- and low-risk groups were investigated (Roh et al., 2017).

Clinical specimen collection

The clinical specimens used in this experiment were all from the tissues excised from patients with LC in the First Hospital of Shanxi Medical University. Written informed consent had been obtained from all patients. This study was conducted in accordance with the Declaration of Helsinki. The protocol was approved by the Ethics Committee of the First Hospital of Shanxi Medical University on November 8, 2022 with the ethical approval number (2022) Ethical review (K184). We collected 29 clinical specimens for RT-qPCR and 12 specimens for immunohistochemical staining. The cancer group was obtained from laryngeal cancer tissue, and the normal group was taken from normal tissue five centimeters away from the edge of laryngeal cancer. Adopting the macroscopic dissection method, the minimum volume of the sample was 0.5 * 0.5 * 0.5 mm³, which may vary slightly depending on the size of the primary tumor. As for RT-qPCR, the specimen was immediately placed in an enzyme-free EP 2ml tube containing RNA later protective solution and stored at 4 °C. After 24 h, the protective solution was removed, and the specimen was then frozen in a −80 °C ultra-low temperature refrigerator. As for IHC, the specimen was immediately soaked in polyformaldehyde and stored in a refrigerator at 4 °C.

RT-qPCR

In this experiment, one mL of Trizol reagent (Ambion) was added to 50–100 mg of LC or adjacent normal tissue for total RNA extraction. RNA concentration was measured using a GENESYS 150 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was synthesized via reverse transcription using the EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TranSgen, Beijing, China) according to the manufacturer’s instructions. One microgram of RNA was used for cDNA synthesis in a final volume of 20 µL, which was stored at −20 °C. A 20 µL qPCR reaction system was prepared using PerfectStart Green qPCR Supermix, and qPCR was conducted in triplicate on a LightCycler 96 real-time PCR system (Roche, Basel, Switzerland). The two-step PCR protocol consisted of an initial denaturation at 94 °C for 30 s, followed by 40 cycles of 94 °C for 5 s and 60 °C for 30 s. GAPDH served as the internal reference gene for the prognostic biomarkers, and relative expression levels were calculated using the 2^−ΔΔCt method. RT-qPCR primers are listed in Table 1.

Immunohistochemistry (IHC)

Paraffin-embedded tissue samples were sectioned to four µm, dewaxed using xylene and ethanol, and subjected to antigen retrieval in sodium citrate solution. After three 5-minute PBS washes, endogenous peroxidase activity was blocked using a 3% H2O2 methanol solution, followed by blocking with goat serum for 20 min. Sections were incubated overnight at 4 °C with primary antibodies (1:400 dilution, Bioss, BS-2044R and BS-8497R, Beijing). Afterward, secondary antibodies (ZSGB-BIO, PV-6000, Beijing) were applied for 1 h at 37 °C. DAB (CST 8059S, USA) was used for staining, followed by counterstaining with hematoxylin (Nanchang Yulu). The sections were subsequently dehydrated, mounted, and scanned.

Statistical analysis

Data are expressed as mean ± standard deviation. Statistical analyses were performed using GraphPad Prism 9 (GraphPad Software, San Diego, CA, USA). For paired samples, normality was assessed using the Shapiro–Wilk test, and a T-test was applied for normally distributed data, while the Wilcoxon test was used for non-normally distributed data. Fisher’s exact test was employed for analyzing clinicopathological parameters. A P-value of <0.05 was considered statistically significant.

Identification of DEGs and DE-ICs

A total of 1,509 DEGs were identified between LC and normal samples in the TCGA-LC dataset, comprising 914 up-regulated and 595 down-regulated DEGs (Figs. 1A, 1B). The ssGSEA algorithm revealed the proportions of immune cells (Fig. 1C). Differential analysis of immune cell distributions indicated significant differences in the proportions of 21 immune cell types, including aDCs, APC.co.stimulation, CCR, Check.point, Cytolytic.activity, HLA, iDCs, Inflammation.promoting, Macrophages, MHC.class.I, NK.cells, parainflammation, pDCs, T.cell.co.inhibition, T.cell.co.stimulation, Tfh, Th1.cells, Th2.cells, TIL, Treg, and Type.I.IFN.Response, categorized as DE-ICs (Fig. 1D). These 21 DE-ICs were classified into high and low expression groups based on the median value. Survival analysis revealed significant survival differences associated with three immune cell types, where high expression of Parainflammation correlated with poorer prognosis, while elevated levels of Tfh and NK cells were associated with better prognosis (Fig. 1E).

Identification of module genes

In the WGCNA, the soft threshold was set to seven (R² = 0.85) to construct a scale-free network (Fig. 2A), which helped to identify gene modules that might have had an important role in the development of laryngeal cancer. Fifteen modules except for the grey module were identified based on the setting of MEDissThres = 0.2 and dynamic tree cutting (Fig. 2B). These modules represented collections of genes with similar expression patterns that might have been associated with different biological processes or features of laryngeal cancer. The cyan and pink modules were the most relevant modules to Parainflammation, Tfh, NK.cells, and case (|cor|>0.3 and P < 0.05; Fig. 2C), which suggested that the genes in these modules might have played important roles in the regulation of the immune microenvironment, inflammatory response, and disease progression in laryngeal cancer. Therefore, a total of 1,051 genes (224 in cyan, 827 in pink) in the two modules were obtained as module genes (Fig. 2D, Tables S2A, S2B).

Identification and enrichment analyses of key genes

The overlap analysis showed that 124 intersection genes were acquired between 1,051 module genes and 1,509 DEGs (Table S3A, Fig. 3A), which were defined as key genes for subsequent analyses. The enrichment analysis of the key genes illustrated that 866 GO terms were enriched, such as defense response to virus, defense response to symbiont, response to virus, regulation of response to biotic stimulus, etc. in BP; secretory granule membrane, endocytic vesicle, ficolin−1 −rich granule, etc. in CC; and peptide binding, antigen binding, integrin binding, peptide antigen binding, etc. in MF (Fig. 3B). Moreover, 47 KEGG pathways were enriched, which mainly contained phagosome, Epstein-Barr virus infection, viral myocarditis, natural killer cell mediated cytotoxicity, etc. (Fig. 3C). This suggested that key genes might have influenced immune escape from laryngeal cancer, cancer cell clearance, and disease progression by participating in these pathways.

RENBP and OLR1 were prognostic biomarkers

Based on the FPKM of the 124 key genes, the univariate Cox analysis results were visualized in a forest map. As can be seen, two key genes relevant to survival were screened out from the training set, including RENBP and OLR1 (Fig. 4A). This suggested that the expression levels of these two genes might have been strongly associated with the survival prognosis of laryngeal cancer patients. The LASSO regression analysis of the two key genes suggested that these two genes were screened out as prognostic biomarkers when the cross-validation error was lowest (lambda.min = 0.03) (Fig. 4B). The risk score of the two prognostic biomarkers = − 0.69348 ×RENBP + 0.390337×OLR1. It could be found that patients in the high-risk group had worse survival, and OLR1 and RENBP were highly expressed in the high- and low-risk groups respectively (Fig. 4C). The Kaplan–Meier (K-M) curve illustrated that the low-risk patients had higher survival probability (P = 0.0098) (Fig. 4D). Additionally, ROC curve analysis showed AUCs greater than 0.7 for 1-, 3-, and 5-year survival in the training set, indicating the robust predictive capacity of the risk model (Fig. 4E). Moreover, the validation results of both internal and external validation (GSE27020) sets showed the consistent results with those of the training set (Figs. 5A–5F), and the AUCs of its 1-, 3- and 5-year survival rates were all greater than 0.65. The stability and reliability of the risk model, which could effectively predict the prognosis of laryngeal cancer patients on different sample sets, had been fully demonstrated, and it had certain potential for clinical application.

Prediction of nomogram was accurate

The correlations between clinicopathological characteristics (Age, T, N, M, Gender, and Stage) of 83 TCGA-LC samples and the prognostic biomarkers were shown in the heat map (Fig. 6A), and it could be found that there were significant survival differences in N1-N3, M0, T3-T4, Stage III–IV, male, and age≤60 subgroups between different risk groups (Fig. 6B). This suggested that these clinicopathological features were closely related to the prognosis of patients with laryngeal cancer and could be used as an important reference index for assessing the prognostic risk of patients. The univariate Cox independent prognostic analysis result illustrated that N, Age and risk score could be considered independent prognostic factors (P < 0.05) (Fig. 6C). Next, the Multivariate Cox independent prognostic analysis showed that the N (P = 0.04) and risk score (P = 0.001) could be regarded as independent prognostic factors for LC (Fig. 6D).

Moreover, it could be seen that the C-index of the nomogram based on N and risk score was 0.6886, and the calibration curves revealed that the prediction of this nomogram was accurate (Figs. 6E, 6F). Moreover, in the 1-, 3-, and 5-year DFS survival ROC curves, the AUCs of the nomogram were optimal, reaching 0.827, 0.775 and 0.763, respectively (Fig. 6G).

Expression validation of model genes

It could be found that, in the training set, the expressions of RENBP and OLR1 were both significantly higher in LC samples (Fig. 7A). And the validation results from the GSE51985 showed consistency with those in the training set (Fig. 7B).

GSVA

To analyze the effects of grouping on LC progression, GSVA was employed. The hallmark enrichment results revealed that 19 and 13 pathways were enriched in high- and low-risk groups, respectively. Enrichment of these pathways suggested that laryngeal cancer cells might have been more invasive and metastatic in the high-risk groups. It could be observed that in the high-risk group, the significantly enriched pathways mainly included epithelial-mesenchymal transition (EMT), coagulation, myogenesis, TGF-beta-signaling, angiogenesis, and hypoxia. In the low-risk group, MYC-targets-VI, reactive-oxygen-species-pathway, E2F-target, and DNA-repair were the major enriched pathways (Fig. 8A). In terms of KEGG pathways, 23 and 26 pathways were significantly enriched in the high- and low-risk groups respectively. For instance, in the high-risk group, ECM-receptor-interaction, focal adhesion, cardiomyopathy-HCM, pathways in cancer, etc. were the major enriched KEGG pathways. In addition, autoimmune thyroid disease, ribosome, and oxidative phosphorylation were the main enriched KEGG pathways in the low-risk group (Fig. 8B).

Immune infiltration and therapy

ESTIMATE analysis revealed a significant difference in stromal scores between the two risk groups, with the high-risk group exhibiting a higher stromal score (Fig. 9A). This suggested that the stromal component of the tumor microenvironment might have been more abundant in the high-risk group, which might have been associated with tumor characteristics such as aggressiveness, growth and metastasis. Additionally, the violin plot of ssGSEA showed that the differences in the proportions of the majority of immune cells were insignificant between the two groups. It was noteworthy that the percentages of CD8.T.cells, Cytolytic.activity, Inflammation.promoting, and pDCs were significantly higher in the low-risk group, while the proportions of T.cell.co.stimulation and Tfh were higher in the high-risk group (Fig. 9B), and the heat map showed the similar results (Fig. 9C).

The TIDE results demonstrated that there were significant differences in the TIDE score, CD8 score and Exclusion score between the two risk groups. Specifically, TIDE and Exclusion scores were higher in the high-risk group, but the CD8 score showed the opposite result (Fig. 9D). Additionally, except that positive correlations were found between the TIDE score and risk score, and between Exclusion score and risk score, other scores showed negative correlations with risk score. The strongest correlation was found between Exclusion score and risk score, which was 0.47 (Fig. 9E). Finally, the expression level of CD274 (PD-L1) was significantly higher in the low-risk group (Fig. 9F). Higher PD-L1 expression in the low-risk group might have meant that immune cells in this group were in a relatively active but suppressed state.

The expression of OLR1 and RENBP in tissues

The results of RT-qPCR showed that the expression of OLR1 (n = 29, P = 0.0003) in LC tissues was higher than that in adjacent normal tissues (Fig. 10A); and the expression of RENBP (n = 29, P = 0.0431) in LC tissues was slightly higher than that in adjacent normal tissues (Fig. 10B). Capture three fields of view randomly under 200× magnification for positive expression area analysis in each group. The results showed that OLR1 (n = 10, P = 0.0068) and RENBP (n = 9, P = 0.0078) were highly expressed in LC tissues (Fig. 10C). Subsequently, we used IHC to detect the expression of PD-L1 in cancer and adjacent tissues. The expression of PD-L1 (n = 12, P < 0.0001) was increased in cancer tissues compared with adjacent tissues (Fig. 10D). Finally, we conducted the analysis of clinical pathological parameters, and the result was presented in Tables S4A, S4B. This further confirmed the up-regulation of OLR1 and RENBP expression in laryngeal cancer identified in the bioinformatics analysis, and provided strong evidence at the experimental level for the use of these two genes as laryngeal cancer-associated biomarkers, suggesting that their aberrant expression in laryngeal cancer tissues might have been closely related to the developmental process of laryngeal cancer.