https://orcid.org/0000-0002-6602-6583
Abstract
Laryngeal squamous cell cancer (LSCC) is the second most prevalent malignancy occurring in the head and neck with a high incidence and mortality rate. Immunotherapy has recently become an emerging treatment for cancer. It is therefore essential to explore the role of tumour immunity in laryngeal cancer. Our study first delineated and evaluated the comprehensive immune infiltration landscapes of the tumour microenvironment in LSCC. A hierarchical clustering method was applied to classify the LSCC samples into two groups (high- and low-infiltration groups). We found that individuals with low immune infiltration characteristics had significantly better survival than those in the high-infiltration group, possibly because of the elevated infiltration of immune suppressive cells, such as regulatory T cells and myeloid-derived suppressor cells, in the high-infiltration group. Differentially expressed genes between two groups were involved in some immune-related terms, such as antigen processing and presentation. A univariate Cox analysis and least absolute shrinkage and selection operator analysis were performed to identify an immune gene-set-based prognostic signature (IBPS) to assess the risk of LSCC. The prognostic model comprising six IBPSs was successfully verified to be robust in different cohorts. The expression of the six IBPSs was detected by immunohistochemistry in 110 cases of LSCC. In addition, different inflammatory profiles and immune checkpoint landscape of LSCC were found between two groups. Hence, our model could serve as a candidate immunotherapeutic biomarker and potential therapeutic target for laryngeal cancer.
Introduction
Laryngeal cancer (LC) is among the most common malignant head and neck tumours. Laryngeal squamous cell carcinoma (LSCC) is the main pathological type and accounts for nearly 90% of LCs. According to the latest data released by the WHO, 184 615 new diagnoses and 99 840 deaths from LC will occur in 2020 (1). Data from the Global Cancer Observator (http://gco.iarc.fr) predict that the number of annual new cases of LC worldwide will increase to 28 million by 2040, suggesting that LC is a serious danger to human health. In China, the annual incidence of LC increased from 2.8 million in 2018 to 4.3 million, and ~40% of LC patients are diagnosed at an advanced stage. Considering many factors, such as lymph node metastasis, the survival rate of patients with LC has not been improved under modern treatment, including surgery, radiotherapy, chemotherapy and comprehensive treatment (2). Thus, it is urgent to find novel and specific treatments for LSCC.
Researchers have found that the immune microenvironment plays an important role in the occurrence and development of tumours (3,4). The human immune system protects the body from external diseases. Tumour cells target immune checkpoint molecules to prevent the body from producing an effective antitumour immune response, which activates the immune escape mechanism (5). In recent years, different degrees of tumour immune infiltration have been identified as prognostic factors for tumours, such as lung cancer, breast cancer and head and neck cancer (6–8). Although immunotherapy has achieved a breakthrough in the clinic, many patients still exhibit different degrees of toxicity and side effects under immune interventions. Individual differences across patients are partially due to the immune heterogeneity of tumours. Therefore, a comprehensive analysis of immune infiltration in tumours is essential for promoting the development of tumour immunotherapy.
In this study, we used 109 LSCC samples from the GEO database to delineate the comprehensive landscapes of the immune microenvironment, including inflammatory profiles and immune checkpoints in LSCC. Furthermore, we identified an immune gene-set-based prognostic signature (IBPS) comprising candidate biomarkers associated with immune infiltration patterns in LC. In total, 112 LSCC samples from the TCGA database and 110 LC tissue samples were applied to further verify the reliability of this gene set. The initial construction of the IBRS of patients with LSCC could help clarify the complex mechanisms underlying the relationship between immune molecules and the prognosis of LSCC and may help optimize immunotherapy strategies for treating LSCC.
Results
Unsupervised cluster analysis identified two immunogenomic subtypes in the LSCC samples
In total, 109 LSCC samples from GSE27020 were included in the current analysis. Based on 28 previously reported immune-associated gene sets representing diverse immune cell types, functions and pathways, we used the gene set enrichment analysis (ssGSEA) scores to quantify the activity or enrichment levels of the immune cells, functions and pathways in the LSCC samples. Then, an unsupervised hierarchical cluster analysis was performed to clearly divide the 109 LSCC samples into two clusters. Based on the heat map of the expression levels of the 28 immune-related gene sets (IRGs), we defined two clusters as two immunogenic subtypes of LSCC, namely, the high-infiltration group and the low-infiltration group (Fig. 1A). We used a principal component analysis (PCA) to show that the degree of infiltration of tumour-infiltrating immune cells (TIICs) had obvious clustering and individual differences (Fig. 1B). Simultaneously, activated CD4 T cells, activated CD8 T cells, myeloid-derived suppressor cells (MDSCs), type 1 T helper cells, type 17 T helper cells and regulatory T cells were considered as the representative immune cells. The degree of infiltration of the above immune cells in the high-infiltration group was significantly higher than that in the low-infiltration group (Fig. 1C–H).
Landscape of microenvironment TIIC composition in LSCC. (A) Enrichment of 28 immune gene sets in the two clusters yielded by hierarchical clustering. We named the two clusters high infiltration and low infiltration. (B) The degrees of TIIC composition in the high-infiltration and low-infiltration groups displayed distinct group-bias clustering and individual differences. (C–H) The degree of infiltration of activated CD4 T cells, activated CD8 T cells, myeloid-derived suppressor cells (MDSCs), type 1 T helper cells, type 17 T helper cells and regulatory T cells in the different immune clusters.
Prognostic significance of immune cells
First, the difference in survival between the two clusters was explored. We found that individuals with low immune infiltration characteristics had a significantly better disease-free survival (DFS) (log-rank P = 0.011) than those in the high-infiltration group (Fig. 2A). The multivariate Cox proportional hazard model showed that the patients in the low-infiltration group had better DFS (HR = 0.32; 95% CI = 0.11–0.97; P = 0.04; Fig. 2B). We performed a univariate Cox regression to determine whether the 28 human immune cell phenotypes were associated with patient prognosis (Fig. 2C). We found that regulatory T cells were significantly correlated with DFS (P = 0.03, HR = 78.353).
Prognostic significance of immune cells. (A) Disease-free survival (DFS) of all LSCC patients based on tumour microenvironment infiltration classes is shown as Kaplan–Meier curves, with the log-rank test showing an overall P = 0.011. (B) According to the age, tumour grade and different immune clusters of the LSCC patients, univariate and multivariate Cox proportional hazards models of DFS were established. (C) Forest plot of the ssGSEA result. A univariate Cox model was used to determine the statistical significance.
Inflammation and tumour immune features
ESTIMATE was used to determine the abundance of immune and stromal cell infiltration. The stromal and immune cell scores in the high-infiltration group were higher than those in the low-infiltration group (Fig. 3A and B). However, the low-infiltration group showed higher tumour purity (Fig. 3C). Then, we investigated whether there was a difference in cytotoxic activity (CYT) between the two immune infiltration groups. The results showed that the CYT score in the high-infiltration group was significantly higher than that in the low-infiltration group (P = 1.9e−10; Fig. 3D). Similarly, the tumour inflammation signature (TIS) and tumour-infiltrating lymphocyte score (TILs) in the high-infiltration group were significantly higher than those in the low-infiltration group (P = 1.5e−10; Fig. 3E and F).
Inflammation and tumour immune features. (A–C) The ESTIMATE model was used to determine the stromal and immune scores and tumour purity of the immune clusters. (D) Comparison of the relative cytotoxic activity scores (CYT) between the immune clusters. (E) Tumour inflammation signature (TIS) in the immune clusters. (F) Pathological evaluation of the percentage of tumour-infiltrating lymphocytes (TILs) in the immune clusters.
Identification of IRGs
Gene expression was compared between the samples of the high- and low-infiltration groups, and in total, 221 immune-related differentially expressed genes (DEGs) were obtained. The heat map shows the expression of the 221 DEGs between the high- and low-infiltration groups (Fig. 4A). The bubble chart shows the 20 pathways in which the DEGs were significantly enriched (Fig. 4B). The enrichment results show that the DEGs are mainly involved in allograft rejection, antigen processing and presentation and cell adhesion molecules. Subsequently, using a univariate Cox proportional regression model, it was found that 45 IRGs were statistically significantly correlated with DFS (Fig. 4C).
Analysis of immune-related genes (IRGs). (A) Heat map showing all 221 DEGs between the different immune clusters. (B) KEGG pathway analysis. (C) Forest plot of some DEGs. A univariate Cox model was used to determine the statistical significance, and 45 DEGs closely related to prognosis (P < 0.05) are displayed.
Verification of the prognostic classifier in TCGA cohort
We used a LASSO Cox regression model to calculate the most valuable prognostic genes and selected the minimum value of lambda as the standard to obtain a model with the following 6 genes: A2M, ANO1, C10orf10, IGFBP3, SERPINE1 and SPARCL1 (Fig. 5A and B). Subsequently, a risk score was established as follows: Risk Score = (0.09903569 × A2M expression) + (0.09367177 × ANO1 expression) + (0.01848474 × C10orf10 expression) + (0.06802683 × IGFBP3 expression) + (0.19349919 × SERPINE1 expression) + (0.12656874 × SPARCL1 expression).
Construction of a prognostic classifier of IRGs in TCGA validation set. (A) Five-fold cross-validation for tuning the parameter selection in the LASSO model. (B) LASSO coefficient profiles of the most useful prognostic genes. (C) The top shows the distribution of RS, the middle shows the patients’ survival time and status, and the bottom shows a heat map of the IRGs in the classifier. (D) Survival curve of the two groups. (E) ROC curve of TCGA cohort. (F, G) Heat map showing the hierarchical clustering results of the ssGSEA of the TCGA cohort. K–M analysis was performed according to the clustering results and was combined with RS.
The risk score of each patient was calculated using this formula. In TCGA validation cohort, the patients were divided into high-risk groups and low-risk groups based on the median risk score (Fig. 5C). The Kaplan–Meier survival analysis showed that the mortality rate in the high-risk group was significantly higher than that in the low-risk group (Fig. 5D, P = 0.015). The ROC curve showed that this model had a strong predictive ability in TCGA validation cohort, and the area under the curve (AUC) was 0.609 (Fig. 5E). Then, we used ssGSEA to calculate the abundance of 28 immune cells in TCGA validation set and performed unsupervised hierarchical clustering of the calculated results (Fig. 5F). The samples were divided into two groups with high- and low-infiltration levels after hierarchical clustering and were combined with the high- and low-risk scores in the K–M analysis (Fig. 5G). The results showed that the survival rate of patients with low-risk scores and low infiltration was better than that in the other three groups. Then, we explored the relationship between the risk score and clinicopathological features (Supplementary Material, Table S1) and found that the risk score was not correlated with sex or age. There was no significant correlation with the condition of the primary tumour, the involvement of the regional lymph nodes, or the stage of the tumour. This finding suggests that these clinicopathological indicators may serve as independent prognostic factors.
Inflammatory profiles and the immune checkpoint landscape in LSCC
First, we collected 9 metagene sets (9) representing different inflammation and immune responses. The details of metagene expression in GSE27020 are shown in Figure 6A. ssGSEA was used to assess the enrichment of the metagene sets. Then, according to the Spearman correlation coefficient between the risk score and the ssGSEA score of the 9 metagene sets, a correlation diagram was used to show the relationships between the variables (Fig. 6B). The analysis showed that a significant positive correlation exists between the risk score and chemokines and receptors (Spearman correlation coefficient = 0.330, P = 0.0005), interleukins and receptors (Spearman correlation coefficient = 0.363, P = 0.0001), other cytokines (Spearman correlation coefficient = 0.655, P < 0.0001), MHC II (Spearman correlation coefficient = 0.252, P = 0.0083) and other MHCs (Spearman correlation coefficient = 0.206, P = 0.0321) in the GEO cohort, with the highest correlation observed with other cytokines. To verify these findings from the GEO cohort, we also repeated the same procedure in TCGA cohort. As shown in Figure 6C and D, the risk score was only significantly positively correlated with other cytokines among the 9 metagenes (Spearman correlation coefficient = 0.489, P < 0.0001). As the risk scores are closely related to inflammatory active substances, we analysed the correlation between different levels of specific immune features and the risk scores (Supplementary Material, Fig. S1A–L). Among the features analysed in TCGA cohort, the risk score was significantly positively correlated with type 1 T helper cells (Spearman correlation coefficient = 0.22, P = 0.022) and regulatory T cells (Spearman correlation coefficient = 0.42, P = 6.4e−06) and significantly negatively correlated with activated CD8 T cells (Spearman correlation coefficient = −0.19, P = 0.049). In ESTIMATE, the risk score was significantly positively correlated with the stromal score (Spearman correlation coefficient = 0.38, P = 4.4e−05) and significantly negatively correlated with tumour purity (Spearman correlation coefficient = −0.22, P = 0.043). The risk score was significantly positively correlated with TIS (Spearman correlation coefficient = 0.21, P = 0.024) and significantly negatively correlated with TILs (Spearman correlation coefficient = −0.33, P = 0.00036).
Relationship between the risk score and inflammatory activity in LSCC. Risk score, clinical feature distribution and metagene expression status across the samples in the (A) GSE27020 and (C) TCGA cohorts. Metagene sets include MHC molecules (MHC I, MHC II and other MHCs), cytokines, immune costimulatory and coinhibitory molecules, chemokines and receptors, interferons and receptors and interleukins and receptors. Corrplots were plotted based on the Spearman correlation coefficient value between the risk score and nine metagenes in (B) GSE27020 and (D) TCGA. The score of each metagene set was calculated by the ssGSEA algorithm.
Survival analysis of the risk score and clinical indicators
We investigated the relationship between the clinical indicators (including sex, race, lymphovascular invasion, perineural invasion, tumour grade and primary tumour status) and the risk scores. As shown in Supplementary Material, Figure S2A, the survival rate of the men in the low-risk group was higher than that in the other three groups (P = 0.0042). Regarding race, the survival rate of white patients in the low-risk group was higher than that in the other three groups (P = 0.047; Supplementary Material, Fig. S2B). The patients in the low-risk group without lymphovascular invasion had a significantly higher survival rate than those in the other three groups (P = 0.0018; Supplementary Material, Fig. S2C), and the patients in the low-risk group with nonperineural invasion had a significantly lower mortality rate than those in the other three groups (P = 0.0038; Supplementary Material, Fig. S2D). Then, we placed the low-grade tumours (G1 and G2) in the favourable risk group and high-grade tumours (G3 and G4) into the unfavourable risk group. Supplementary Material, Figure S2E shows that the survival rate of the patients with grades G3 and G4 in the low-risk group was higher than that in the other three groups (P = 0.005). We placed primary tumours T1 and T2 in the early group and primary tumours T3 and T4 in the advanced group. We found that the survival rate of the patients with advanced tumours in the low-risk group was higher than that in the other three groups (P = 0.047; Supplementary Material, Fig. S2F).
Subsequently, the expression of A2M, ANO1, IGFBP3, SERPINE1, C10orf10 and SPARCL1 was detected by immunohistochemistry (IHC) in 110 LSCC samples. The details of the clinicopathological characteristics of the 110 LSCC samples are shown in Supplementary Material, Table S2. The data showed that 68 LSCC cases (61.82%) had a high level of ANO1 expression, 56 LSCC cases (50.91%) had a high level of A2M expression, 42 LSCC cases (38.18%) had a high level of IGFBP3 expression, 70 LSCC cases (63.64%) had a high level of SERPINE1 expression, 54 LSCC cases (49.09%) had a high level of C10orf10 expression, and 44 LSCC cases (40%) had a high level of SPARCL1 expression (Supplementary Material, Fig. S3A). We also found that the survival rate of the patients in the low-expression groups of A2M, ANO1 and IGFBP3 was better than that in the high-expression groups based on a K–M analysis (P < 0.05; Supplementary Material, Fig. S3B–G).
Discussion
Numerous studies have shown that the occurrence and development of tumours primarily depends on the tumour immune microenvironment (10–12). The tumour immune subtypes based on each individual patient’s immune heterogeneity have become another independent prognostic factor differing from classical tumour clinical parameters, such as the nuclear grade, stage and metastasis status. It has been proven that immune-related risk assessment models can provide a promising novel marker for tumour immunotherapy (13). Considering that the progression of LC is closely related to the immune microenvironment, it is vital to identify prognostic immune biomarkers of LC.
In the present study, 109 cases of LSCC from the GEO database were divided into the high-infiltration group and the low-infiltration group according to the expression levels of IRG sets. The high-infiltration group had worse outcomes than the low-infiltration group. In the high-infiltration group, we found that immune suppressive cells had an elevated infiltration level compared to those in the low-infiltration group. The concentration of immune suppressive cells including regulatory T cells and MDSCs can cause T cell dysfunction and may be a responsible for the poor prognosis in this group (14,15). Our univariate Cox regression analysis showed that the number of regulatory T cells was the only independent prognostic risk factor for immune cells in LSCC. Thus, a high-infiltration level of regulatory T cells contributes to this subtype’s poor prognosis. Regulatory T cells are immunosuppressive cells that play an immunosuppressive role by reducing the costimulatory activity and commission capacity of APC, releasing perforin to kill DCs, and secreting anti-inflammatory cytokines. In the tumour microenvironment, regulatory T cells are reprogrammed to obtain a higher immunosuppressive function (16). It has been confirmed that high regulatory T cell infiltration is associated with low survival rate in various tumours.
Then, we constructed the IBPS with the following 6 immune-related co-expressed genes: A2M, ANO1, C10orf10, IGFBP3, SERPINE1 and SPARCL1. The IBPS is a new prognostic evaluation method independent of the traditional evaluation method based on clinicopathological characteristics. Most genes in this model were found to be involved in tumour progression and tumour immune microenvironment remodelling (17–21). The protein encoded by A2M, as a protease inhibitor and cytokine transporter, can inhibit a series of proteases, including trypsin, thrombin and collagenase. This protein can also inhibit inflammatory cytokines and disrupt the inflammatory cascade. A2M has been reported to accelerate cancer cell proliferation, invasion and colony formation (20). Researchers have found that A2M is crucial for maintaining the stemness of CD9 leukaemia stem cells (22). Receptor-related protein-1 bound A2M was found to be associated with the inhibition of tumour cell proliferation, migration, invasion, spheroid formation and anchorage-independent growth, which might be tractable to therapeutic exploitation (23,24). ANO1, as a calcium-activated chloride channel, plays various roles in regulating epithelial fluid transport, gastrointestinal tract motility and saliva production (25). ANO1 is overexpressed in malignant tumours and promotes the proliferation and metastasis of tumour cells, such as breast cancer, oesophageal cancer, and head and neck cancer (21). Researchers found ANO1 involvement in head and neck cancer carcinogenesis, and ANO1 might be an actionable drug target (26). According to Godse, N. R et al.’s in vitro experiment and cohort analysis, ANO1 is also associated with chemoradiation failure and recurrence (27). The insulin-like growth factor (IGF) system plays a critical role in cell growth, survival and differentiation by paracrine or autocrine methods. IGFBP3 is a major carrier of IGFs in circulation which can regulate IGF activity in different cellular environments by its post-translational modifications, such as glycosylation, proteolysis and phosphorylation. Recent studies have shown that IGFBP-3/TMEM219 axis damage can have a significant impact on the survival outcomes of patients with specific cancers (28). Papadimitrakopoulou, V et al.’s clinical cohort analysis demonstrated that IGFBP3 is correlated with a significantly shorter disease-specific and DFS (29). Thus, IGFBP-3 can be used as a specific tumour diagnostic and prognostic marker and could provide a new target for cancer targeted therapy. C10orf10 has been closely associated with prognosis in cancer. C10orf10 expression is regulated by hypoxic stress and progesterone or androgen (30–32). The C10orf10 product DEPP, which is regulated by oxidative stress, is a major hypoxia-inducible gene involved in the activation of autophagy (33). SPARCL1 is a secreted glycoprotein of the SPARC family of matricellular proteins and plays different roles in tumour biology as an oncogene or tumour suppressor (34). SPARCL1 is also an important regulator promoting the angiogenesis function of cancer stem cells (35). SERPINE1 has a critical effect on sustaining proliferative signals, resisting cell death, angiogenesis, invasion and metastasis, and promoting the epithelial-to-mesenchymal transition and tumour-promoting inflammation (36,37). SERPINE1 affects cancer cell growth by modulating anti-fibrinolytic activity and cell adhesion. By inhibiting cell adhesion to vitronectin, SERPINE1 can also have a pro- and anti-apoptotic effect. Then, 112 LSCC patients from TCGA were chosen as a validation set to test the reliability and stability of the IBPS we constructed. The 6 genes also contribute to several regulatory pathways that affect tumour biology related features. For example, A2M and SERPINE1 are regulators of the body’s complement and coagulation cascade pathways and contribute to the body’s inflammation response (38–40). IGFBP3 and SERPINE1 play important regulatory roles in the cellular senescence pathway (41,42). To preliminarily understand the role of these 6 genes in LC, we detected the expression of the 6 genes in the prognostic model by IHC of 110 frozen LSCC tissues and found that a high expression of ANO1, A2M and IGFBP3 was closely associated with a poor prognosis in LSCC patients. Since inflammatory substances are an important constituent of the immune microenvironment, we further analysed the relationship between the risk score and inflammatory substances. The data indicated that the risk score was significantly positively correlated with TIS and significantly negatively correlated with TILs. Then, we found that high-risk patients were characterized by a high proportion of type 1 T helper cells and regulatory T cells and a low proportion of activated CD8 T cells. This finding suggests that the risk score calculated by the 6 co-expressed genes can serve as a biomarker of TIL, tumour inflammatory stimulus and the prognosis of LC. Previously, most studies focused on the biofunction of a single gene. However, a single-gene analysis may lack stability. In our study, we constructed a co-expressed IRG model, which was validated to be stable in the prognosis of LSCC.
In conclusion, we comprehensively explored the relationship between the composition of the immune microenvironment and LC and constructed a novel IBPS with 6 IRGs. The prognostic model was validated to be stable and reliable in different cohorts. This model may help clinicians make a precise prediction of the outcomes and the response to immunotherapy of LC patients. However, the study cohort we used was retrospective. Sufficient experimental verification and the potential mechanisms of IBPS in LC remain to be further explored.
Materials and Methods
Data collection and preprocessing
The expression profile of GSE27020 based on the GPL96 Affymetrix Human Genome U133A Array platform was downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/). A total of fresh frozen tumour tissue samples from 109 LSCC patients were included in this dataset. The patients’ clinical characteristics, including age, tumour grade, survival time and DFS status, were also obtained. Patients without clinical data were removed from this study. The UCSC Xena browser (GDC hub: https://gdc. xenahubs.net) was used to extract the gene expression data from The Cancer Genome Atlas (TCGA). We screened the data of 112 LSCC patients for further analysis, including all clinical information of these patients. In this study, we used the GSE27020 dataset in the GEO database as the training cohort and the LSCC sample cohort in TCGA-HNSC as the validation cohort.
In addition, the original data of the dataset obtained from the Affymetrix platform were adjusted in the background by the RMA algorithm in the Affy software package and normalized by the z-score. Regarding TCGA dataset, the RNA-sequencing data (FPKM values) were transformed into transcripts per kilobase million values and then normalized by the z-score. Finally, the different probe IDs were converted to their respective gene symbols to annotate them, and the repeated gene expression values were averaged.
Analysis of immune cell characteristics
The gene expression profile (GEP) data were analysed using a single sample ssGSEA, ESTIMATE, CYT, TIS and TILs. ssGSEA (43) in the R package ‘GSVA’ was applied to quantitatively analyse the infiltration levels of the immune cell types and predict the abundance of 28 TIICs in individual tissue samples based on gene sets of 782 genes (44). ESTIMATE (45) was used to calculate the immune and stromal scores. This algorithm can derive composite scores based on the level of immune cell infiltration in tumour tissues and the number of stromal cells present. Immune CYT (46) was determined based on a previously reported formula. The molecules perforin-1 (PRF1) and granzyme A (GZMA) were used to calculate CYT, and these molecules are known to be closely related to CD8+ T cell activation. The tumour inflammation signature (TIS) (47) is an 18-gene signature constructed for investigational purposes only that is able to quantify background adaptive immune responses that are suppressed in tumours. The TIS score is estimated as the average of gene expression of the marker genes. The density of T cells is known to be positively correlated with a better prognosis in many human cancers. To investigate tumour-infiltrating T cells, the proportion of TILs (48,49) was calculated.
Hierarchical clustering and identification of differentially expressed genes
Based on the ssGSEA results, we clustered the immune microenvironment of LSCC. Then, according to the results of the sample classification in hierarchical clustering, the R package limma was used to extract the DEGs to identify genes related to tumour immune microenvironment infiltration. The DEGs were determined by previously implemented significance criteria (| logFC | > 2/3, P < 0.05).
Functional enrichment and pathway analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) (50) and DAVID (51) analyses (http://david.abcc.ncifcrf.gov/home.jsp) were used to explore the functions of the DEGs.
Construction of the risk assessment model
After the IRGs were identified in the GEP of LSCC in GSE27020, a risk assessment model was constructed. First, a univariate Cox regression analysis was used to explore the effects of the DEGs on DFS. DEGs with P < 0.05 were considered significantly associated with survival and were used in the subsequent analyses. Then, the least absolute shrinkage and selection operator (LASSO) method was used to select variables in the Cox regression model by decreasing the dimensionality of the features to determine the best prognostic genes, and the minimum value of lambda was used (52,53). Finally, the following risk score formula was calculated by considering the expression of optimized genes and the correlation estimated Cox regression coefficients: Risk score = (exp Gene1 × coef Gene1) + (exp Gene2 × coef Gene2) + … + (exp Gene6 × coef Gene6). By sorting the calculated risk scores, the LSCC patients were classified into high-risk or low-risk groups based on the median risk score. The difference in survival between the high-risk and low-risk groups was assessed using the Kaplan–Meier method and a two-tailed log-rank test.
Tissue samples
In this study, 110 paraffin-embedded specimens of LSCC were collected from the Department of Pathology of the Second Affiliated Hospital of Harbin Medical University. From February 2008 to October 2012, partial or total laryngectomy was performed in the Department of Otolaryngology of the Second Affiliated Hospital of Harbin Medical University. All 110 patients were followed up for at least 5 years. All patients provided written informed consent in accordance with the Declaration of Helsinki. This study was approved by the Ethics Committee of the Second Affiliated Hospital of Harbin Medical University.
Immunohistochemistry
Paraffin-embedded specimens from 110 cases of LSCC were cut into 4 mm thick sections. The tissue sections were dewaxed in xylene and rehydrated in alcohol. The antigen was recovered by microwaving under different powers. Endogenous peroxidase was removed by treatment with 3% hydrogen peroxide for 10 min at room temperature. The immunohistochemical analysis was performed with an A2M antibody (1/500, ab109422, Abcam Biochemicals, UK), ANO1 antibody (1/500, ab190803, Abcam Biochemicals, UK), IGFBP3 antibody (1/100, ab217205, Abcam Biochemicals, UK), SERPINE1 antibody (1/50, ab66705, Abcam Biochemicals, UK), C10orf10 antibody (1/50, ab204420, Abcam Biochemicals, UK) and SPARCL1 antibody (1/1000，ab255597, Abcam Biochemicals, UK). The two pathologists adopted a unified standard and single blinded method. At least 1000 tumour cells were counted on each slide and evaluated by IHC. We evaluated the IHC results as previously described (54).
Statistical analysis
The Cox regression model and K–M analysis were based on the ‘survival’ package. The ‘forestplot’ package was used to summarize the survival results. The limma package was applied in the differential analysis. A Wilcoxon rank sum test with nonparametric statistical hypotheses was used to analyse data from two groups. A Spearman correlation analysis was performed to determine the correlation between two variables. The statistical analysis was implemented in R software (version 3.6.1). In addition, when we used P to evaluate the results, P < 0.05 indicated that the results were statistically significant.
Acknowledgements
We are grateful to the contributors of the data to public databases, including TCGA, GEO, and the Affymetrix platform.
Conflict of Interest statement. The authors declare no conflicts of interest.
Funding
The China Postdoctoral Science Foundation (2020M681115) and the National Natural Science Foundation of China (Grant No. 82103164).
Authors’ Contributions
R.Z. and D.S. did the conceptualization and designing, and had written—reviewed and edited the manuscript. Z.Q. and L.F. did the designing, data analysis, interpretation and had written—original draft of the manuscript. J.Z. and L.J. did methodology and project administration. X.C. did data acquisition and conducted research. All authors contributed to the article and approved the submitted version.
References
Author notes
Indicates Zhipeng Qian and Desi Shang authors contributed equally to this work.
