Skip to main content

Immune-related lncRNA classification of head and neck squamous cell carcinoma



Long noncoding RNAs (lncRNAs) play a critical role in innate and adaptive immune responses. Thus, we aimed to identify ideal subtypes for head and neck squamous cell carcinoma (HNSCC) based on immune-related lncRNAs.


TCGA HNSCC cohort was divided into two datasets (training and validation dataset), and 960 previously characterized immune-related lncRNAs were extracted for non-negative matrix factorization analysis. We characterized our HNSCC subtypes based on biological behaviors, immune landscape and response to immunotherapy in both training and validation cohort. A lncRNA-signature was generated to predict our HNSCC subtypes, and essential lncRNAs involved in tumor microenvironment (TME) were identified.


We developed and validated two HNSCC subtypes (C1 and C2) based on the 70 lncRNAs in the training and validation cohort. C2 subtype displayed good prognosis, high immune cell infiltration, immune-related genes expression and sensitivity to PD-1 blockade. C1 subtype was associated with high activity of mTORC1 signaling and glycolysis as well as high fraction of inactive immune cells. Finally, we generated a 31-lncRNA signature that could predict our above subtypes with high accurate. Additionally, TRG-AS1 was identified as the essential lncRNA involving TME formation. Knockdown of TRG-AS1 inhibited the expression of HLA-A, HLA-B, HLA-C, CXCL9, CXCL10 and CXCL11. High expression of TRG-AS1 indicated a favorable prognosis in HNSCC and anti-PD-L1 cohort (IMvigor210).


Our study establishes a novel HNSCC classification on the basis of 31-lncRNA, helping to identify beneficiaries for anti-PD-1 treatment. In addition, a critical lncRNA TRG-AS1 is identified as a new potential prognosis biomarker as well as therapeutic target.


Head and neck squamous cell carcinoma (HNSCC) is the sixth most common malignant tumor, with approximately 640,000 new cases worldwide each year [1]. Despite significant advancements in treatment, mortality rates for HNSCC remain at around 50%. Thus, it is essential to explore novel and effective therapeutic strategies to improve the clinical outcomes of HNSCC.

Recently, immunotherapy has received more and more attention in the area of cancer treatment owing to its remarkable and stable overall survival advantages. HNSCC might also be effective to immunotherapies since its frequent mutations and resulted neoantigens [2]. Indeed, platinum-pretreated metastatic and recurrent HNSCC receiving anti-programmed cell death (PD)-1 therapy showed durable clinical and survival benefit [3, 4]. However, the overall response rates of immunotherapy are less than 20% in unselect HNSCC patients [3, 5]. A better understanding of the tumor microenvironment (TME) formation and selection of potential beneficiaries may help increase survival benefit of immunotherapy.

Increasing evidence have proved that long noncoding RNAs (lncRNAs) played a critical role in innate and adaptive immune responses via regulating the differentiation and function of immune cells. For example, knockdown of lncRNA Pvt1 could significantly suppress the immunosuppressive function of granulocytic myeloid-derived suppressor cells [6]. In terms of adaptive immune responses, lncRNA NKILA and EGFR could reduce cytotoxic T lymphocytes (CTLs) infiltration and CTLs activity [7, 8]. The function of few lncRNAs in HNSCC immunology has been demonstrated [9, 10], however, a large number of immune-related lncRNAs have not been investigated thoroughly.

In this study, we systematically analyzed the immune-related lncRNAs in HNSCC. Two HNSCC subtypes (C1 and C2) based on the prognostic value of immune-related lncRNAs were identified in both training and validation cohort. C2 subtype exhibited higher immune cell infiltration, fractions of active immune cells, expression of immune-associated molecular and the response rate of anti-PD-1 treatment than C1 subtype in both training and validation cohort. In addition, TRG-AS1 acted an important role in regulating TME of HNSCC and might be a potential therapeutic target.


Data source

Level 3 RNA-Seq data consisting of 502 HNSC tissues and 44 normal controls were downloaded using TCGAbiolinks R package [11] (up to April 21, 2020). We also obtained corresponding clinical information, including age, gender, tumor grade, TNM stage, survival time, and survival status. After filtering out non-primary tumors and following up for less than 30 days, 490 HNSC samples were finally included and randomly split into two cohorts: training cohort (n = 343) and validation cohort (n = 147). The baseline information was presented in Table 1.

Table 1 Clinicopathological characteristics of training and validation cohort

Identification of HNSCC subtypes

A total of 960 immune-related lncRNAs were achieved from a previous study [12]. First, univariate Cox analysis was utilized to filter out lncRNAs without prognosis value (P < 0.05). Subsequently, non-negative matrix factorization (NMF) clustering method was employed in the training dataset using CancerSutypes R package [13]. At last, the same candidate lncRNAs were also utilized in the process of NMF clustering in the validation dataset to further verify the robustness of the above classifier.

Function enrichment and gene set variation analysis

We applied limma R package [14] to identify differentially expressed mRNAs (DEmRNAs) based on the cut-off criteria: absolute log2fold change (FC) ≥ 1 and adjusted P value < 0.05. Subsequently, gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using clusterProfiler R package [15]. Adjusted P < 0.05 was considered statistically significant. In addition, ‘Hallmark’ gene sets were downloaded using msigdbr R package for running gene set variation analysis (GSVA) [16], which is a commonly performed method for assessing the variation in biological process activity and pathway in the expression datasets samples. Statistically significance was defined as |log2FC| > 0.1 and adjusted P < 0.05.

Estimation of immune infiltration

A total of 23 types of immune cells signatures were extracted from a published study [17] and the immune cell infiltration was quantified by single-sample gene set enrichment analysis (ssGSEA) based on GSVA R package. The abundance of immune cell was verified by the MCP counter and CIBERSORT, which were obtained from TIMER website ( [18]. In addition, the stromal scores and immune scores was also calculated by ESTIMATE algorithm [19].

Immunotherapy response and subtype prediction

To indirectly predict the immunotherapy response of our subtypes, we measured the similarity of gene expression profiles between our subclasses and immunotherapy-treated melanoma patients based on subclass mapping [20, 21]. In addition, in order to facilitate clinical application, logistic least absolute shrinkage and selector operation (LASSO) algorithm [22] was employed to predict our subtypes classification.

Identification of critical lncRNAs involved in TME formation

Random forest R package [23] was utilized to rank the importance of lncRNAs filtered in LASSO. The correlation between the important lncRNAs and immune infiltration was evaluated.

Cell culture

The human HNSCC cell line CAL27 was obtained from ATCC (Manassas, VA, USA). Cells were cultured in Dulbecco’s Modified Eagle’s Medium (Gibco, USA) supplemented with 10% fetal bovine serum (Gibco, USA) at 37 °C in a humidified 5% CO2 incubator.

Cell transfection

A small interfering RNA (siRNA) targeting TRG-AS1 (si-TRG-AS1) and a negative control siRNA (si-NC) were designed and purchased from GenePharma (Suzhou, China). The sequence of si-TRG-AS1#1 was 5′- GGAGCUGGACUACAGUGAUdTdT-3′, and the sequence of si-TRG-AS1#2 was 5′-GCAACUACCUCAUAGAUUUdTdT-3′. Lipofectamine 3000 (Invitrogen, Carlsbad, CA, United States) was used to transfect si- TRG-AS1 and si-NC into HNSCC cells following the instructions of the manufacturer.

Quantitative real-time PCR

Total RNA was isolated from CAL27 cells by the RNA-quick purification kit (ESscience Biotech, China) following the manufacturer’s instructions. A total of 2 µg RNA was synthesized from cDNA using HiScript III-RT SuperMix (Vazyme, China). The obtained cDNA product was used for real-time quantification PCR (RT-qPCR). The following primers were used: TRG-AS1, forward 5′-CTCCTGGCTGATCCCACT-3′, and reverse 5′-CACTATGCCATCCTGTACCAC-3′; HLA-A, forward 5′-AGATACACCTGCCATGTGCAGC-3′, and reverse 5′- GATCACAGCTCCAAGGAGAACC-3′; HLA-B, forward 5′-CTGCTGTGATGTGTAGGAGGAAG-3′, and reverse 5′ -GCTGTGAGAGACACATCAGAGC-3′; HLA-C, forward 5′-GGAGACACAGAAGTACAAGCGC-3′, and reverse 5′- ACATCCTCTGGAGGGTGTGAGA-3′; CXCL9, forward 5′-CAATTTGCCCCAAGCCCTTC-3′, and reverse 5′-TTTTCTTTTGGCTGACCTGT-3′; CXCL10, forward 5′-GGATGGCTGTCCTAGCTCTG-3′, and reverse 5′- TGAGCTAGGGAGGACAAGGA-3′; CXCL11, forward 5′- AGCCTTGGCTGTGATATTGT-3′, and reverse 5′- GGGTACATTATGGAGGCTTTCT-3′; GAPDH, forward 5′-CTCCTCCTGTTCGACAGTCAGC-3′, and reverse 5′-CCCAATACGACCAAATCCGTT-3′.

Statistical analysis

Log-rank test was used to compare the overall survival under different conditions assessed by Kaplan Meier survival curve. The cut-off of expression was identified by survminer package. The hazard ratios (HRs) and 95% confidence intervals (CIs) of OSCC mortality risk was estimated by univariate and multivariate Cox proportional hazards models. R was used to perform all statistical analyses.


NMF identifies two subclasses in HNSCC

A flowchart was developed to summarize our study (Fig. 1a). We obtained 70 prognosis-related candidate lncRNAs in training cohort by univariate Cox analysis, and according to their expression profile, we clustered the 343 HNSCC samples in the training cohort using NMF clustering. The results indicated that, when subtypes = 2, the maximum value of the average silhouette width (ASW) could be reached and the heatmap also maintained a highest level of similarity within a cluster, consequently, the training set should be divided into 2 subtypes (Fig. 1b). Subsequently, we applied NMF to validation cohort with 147 samples and the results also indicated that there were two distinct molecular subtypes of HNSCC (Additional file 1: Fig. S1). Based on above classification, C2 had better prognosis than C1 in training and validation cohort (Fig. 1c, d). After adjust potential confounds, C2 subtype still showed lower mortality risk compared with C1 subtype (0.52 [0.33, 0.80] in the training cohort; 0,49 [0.25, 0.94] in the validation cohort) (Table 2).

Fig. 1
figure 1

Identification of HNSCC subtype based on NMF analysis in the training cohort. a Flow diagram of the study. b NMF analysis based on 70 immune-related lncRNAs. c, d Survival analysis of the two HNSCC subtypes in training and validation cohort

Table 2 Relationship between HNSCC subtypes and overall survival in different models

The biological behaviors in distinct lncRNA-related patterns

To better understand the biological behaviors between two subtypes of HNSCC, we first performed differential analyses. A total of 309 DEmRNAs were obtained between C1 and C2 group and these genes were mainly enriched in immune related biological processes, such as T cell activation, regulation of lymphocyte/T cell activation and cell-cell adhesion (Fig. 2a). Besides, GSVA was applied to identify different pathways between C1 and C2 group. The results showed that immune related pathways (Interferon alpha/gamma response, IL6 JAK STAT3 signaling, Complement, inflammatory response and IL2 STAT5 signaling) were activate in C2 subtype, while MYC targets, mTORC1 signaling, glycolysis and cholesterol homeostasis etc. were higher in C1 group (Fig. 2b). Validation cohort also found similar results, and thus further verified the robustness of HNSCC classification (Fig. 2c, d).

Fig. 2
figure 2

Biological characteristics of HNSCC subtypes. a, c GO analysis in training and validation cohort. b, d GSVA analysis in training and validation cohort

TME cell infiltration characteristics in distinct lncRNA-related patterns

Since C2 subtype activated multiple immune-related pathways, we further explored the difference in TME between C1 and C2 subtype. C2 group gained more immune score and stromal score than C1 subtype in the training cohort, while a higher stromal score was not observed in the validation cohort (Fig. 3a, b). Next, immunologic landscape was characterized based on immune cell infiltration. Corresponding to the immune score, most of immune cells infiltration was higher in C2 subtype, including innate and adaptive immunity cells (Fig. 3c, d). Furthermore, MCP analysis was also performed to validate the above results. Fortunately, we obtained similar results, that is, the infiltration of B cell, T cell, CD8+ T cell, monocyte, macrophage/monocyte and myeloid dendritic cell was higher in C2 group. In addition, cytotoxicity score was also elevated in C2 group (Additional file 2: Fig. S2). Relative proportions of 22 immune cells subsets were estimated based on CIBERSORT method. In both cohorts, the infiltration fractions of native B cell, macrophage M1, activated mast cell, activated NK cell, activated memory CD4+ T cell, CD8+ T cell, T follicular helper (TFH) and regulatory T cells (Tregs) were significantly upregulated, while eosinophil, macrophage M0, resting master cells and resting NK cells were obviously downregulated in C2 subtype than C1 subtype (Fig. 3e, f). These results indicated that C2 group might be ‘hot tumor’.

Fig. 3
figure 3

Immune characteristics of HNSCC subtypes. a, b Boxplot of immune score in training cohort and validation cohort. c, d Boxplot of immune cell infiltration in training cohort and validation cohort. e, f Boxplot of immune cell fraction in training cohort and validation cohort. (*P < 0.05, **P < 0.01, ****P < 0.0001, ns represents no significance)

In addition to immune infiltration, we also assessed the expression of immune-related genes between the two groups, and the results were largely consistent with immune infiltration (Fig. 4). Most HLA molecules were generally upregulated in the C2 group. Immune checkpoints, such as LAG3, PDCD1, CD274, CTLA4, and TIGIT, were observed increased in the C2 group. Moreover, we also found the expression of most immunomodulators were higher in C2 group than C1 group, including chemokines and cytolytic activity related genes, such as CXCL9, CXCL10, IL2, GZMA and RPF1. Five of the above dysregulated immune-associated genes are the member of “IFNg signature” (IDO1, CXCL10, CXCL9, HLA-DRA, STAT1 and IFNG), suggesting their positive clinical response to anti-PD-1 treatment. In summary, patients in C2 group might be more sensitive to immunotherapy.

Fig. 4
figure 4

Heatmap of immune-related genes in training and validation cohort

Distinct response to immunotherapy for lncRNA-related patterns

To assess the response to immunotherapy for the two subtypes, we compared their gene expression profile with a published dataset which included 47 melanoma patients who received immunotherapies. As expected, significant association was observed between the expression profile of C2 group and PD-1-response group, indicating that C2 group was more promising to respond to anti-PD-1 (Fig. 5a). Therefore, our classification of HNSCC was robustness and had potential ability to seek general susceptibilities to anti-PD-1 therapy. To better advance their clinical application, we applied LASSO to construct a prediction model in the training cohort and obtained a 31-lncRNA signature. The formula of the lncRNA-signature was shown in (Additional file 3: Table S1). Based on the cut-off value (C1: < 3.68, C2: > 3.68) of training cohort, we further explored the subtypes of validation cohort. Fortunately, we observed a high concordance (92.4% in subtype C1, 98.2% in subtype C2) in the validation cohort (Fig. 5b).

Fig. 5
figure 5

Immunotherapeutic response and identification of predictive classifier. a C2 may be more response to the PD-1 inhibitor (nominal and Bonferroni corrected P < 0.05) by SubMap analysis in training and validation cohort. b Concordance of HNSCC subtypes prediction between original classification based on NMF and the 31-lncRNA classifier

TRG-AS1 as essential lncRNA in regulating TME

To further explore the essential lncRNAs that regulated TME, the 31 lncRNAs from the above obtained lncRNA-signature were subjected to random forest analysis. The results indicated that TRG-AS1 was the most important lncRNA in our clustering (Fig. 6a). TRG-AS1 was highly expression in C2 group and it was significantly positively correlated with most types of immune cell infiltration (Fig. 6b, c). The MHC-I molecules and chemokines (e.g. CXCL9/10/11) paly an important role in immune cell infiltration, and thus we further assessed the relationship between TRG-AS1 and above genes. The results indicated that knockdown of TRG-AS1 inhibited the expression of HLA-A, HLA-B, HLA-C, CXCL9, CXCL10 and CXCL11 (Fig. 6d). As expected, the high expression of TRG-AS1 indicated a favor prognosis (Fig. 6e). Thus, TRG-AS1 might play an essential role in the TME formation of HNSCC. In addition, we also found a strong positively relationship between TRG-AS1 and immune cells in other types of tumors, especially BRCA and CHOL (Additional file 4: Table S2), which was obtained from a previous study. Moreover, OS benefit was also observed in the high TRG-AS1 expression group in the anti PD-L1 treatment cohort (IMvigor210) (Additional file 5: Fig. S3). Collectively, TRG-AS1 might be of great importance in the formation of TME and tumor development.

Fig. 6
figure 6

Identification of TRG-AS1 as an essential lncRNA. a 31-lncRNA contribution to HNSCC subtype. b Boxplot of TRG-AS1 expression in training cohort and validation cohort. c The correlation between TRG-AS1 and immune cell infiltration in training cohort and validation cohort. d Expression levels of TRG-AS1, HLA-A and HLA-B, HLA-C, CXCL9, CXCL10 and CXCL11 in CAL27 cells after transfection with control siRNA or TRG-AS1 siRNA. e Survival analysis of TRG-AS1 in training and validation cohort


Immunotherapy is a promising treatment because it shows significant and durable clinical benefits in advanced HNSCC patients, however, less than 20% HNSCC could benefit from it, which highlights the demand to identify ideal subtypes for immunotherapy in HNSC. With the deepening in the molecular mechanisms of tumor immunity, lncRNAs attracted more and more attention due to their functions of regulating innate and adaptive immune cells response. Therefore, immune-related lncRNAs may help to explore the subtype which is more sensitive to immunotherapy.

In this study, we identified two robust subtypes (C1 and C2 subtype) based on the 70 immune-related lncRNAs with distinct TME characteristics. Compared with C1 subtype, C2 subtype presented increased immune cell infiltration and HLA molecules, suggesting initial recognition by host immune system and high level of cancer antigen presentation [24]. Higher fractions of activated immune cells (e.g. macrophage M1, activated mast/NK/memory CD4+ T cell, CD8+ T cell, TFH cell), upregulated expression of chemokines and cytolytic activity related genes as well as cytotoxicity score in C2 subtype indicated an active antitumor immune response [25,26,27]. Therefore, C2 subtype could be characterized by immunoinflammatory phenotype (also known as ‘hot tumor’), which maybe a possible explanation for the longer survival time of C2 subtype, and was consistent with previous studies [17, 25, 28]. Accordingly, it is reasonable to find that C2 subtype might more likely to benefit from immunotherapy owing to its TME characteristics.

C1 subtype characterized by low immune cell infiltration, decreased immune-related genes and high fraction of inactive immune cells, corresponding to immune-desert phenotype. GSVA analysis found that C1 subtype had increased activity of mTORC1 signaling, which is a metabolic master regulator. Active mTORC1 leads to the elevated glycolysis [29, 30], which is in line with our results, that is C1 subtype had higher activity of both mTORC1 signaling and glycolysis. Increased glycolysis in tumor cells are associated with immunosuppressive TME. On the one hand, as a result of glucose deprivation, tumor-infiltrating lymphocytes have decreases in antitumor effector molecules production and myeloid-derived suppressor cells were recruited to the TME [31,32,33]. On the other hand, an acidic TME formation hampers the expression of IFNγ in cytotoxic cells (T cells and NK cells) and promotes IL-17-mediated and IL-23-mediated inflammation [34, 35]. These might partly explain the TME characteristics as well as the poor prognosis of C1 subtype. Accordingly, it is reasonable to speculate the combination of glycolysis inhibitor and immune-checkpoint inhibitor might improve clinical outcome of C1 subtype.

In addition to identify and validate robust subtypes with distinct TME characteristics in HNSCC, we still found that lncRNA TRG-AS1, T cell receptor gamma locus antisense RNA 1, might play a critical role in regulating TME and might serve as potential prognosis biomarker as well as therapeutic target. Knockdown of TRG-AS1 suppressed the expression of HLA-A, HLA-B, HLA-C, CXCL9, CXCL10 and CXCL11. HLA-A, HLA-B and HLA-C play an important role in antigen presentation, thereby initiating an immune response. The expression of CXCL9, CXCL10 and CXCL11 is positively correlated with the density of tumor infiltrating NK and T cells [36]. Thus, it was not surprising to find that TRG-AS1 was significant positively correlated with multiple immune cell infiltration and long survival time. In addition, high expression of TRG-AS1 was also correlated with more immune cell infiltration in multiple tumors as well as OS benefit in the anti PD-L1 treatment cohort, mainly including bladder and kidney cancer. Thus, the above results highlighted the necessary to explore the function of TRG-AS1 in pan-cancer.


In summary, we generate a novel HNSCC classifier based on 31-lncRNA, which helps to identify ideal candidates for anti-PD-1 treatment. In addition, an important lncRNA TRG-AS1 is identified to a novel potential prognosis biomarker as well as therapeutic target.

Availability of data and materials

The datasets used during the current study are available from the Cancer Genome Atlas (TCGA) dataset (



Long noncoding RNAs


Head and neck squamous cell carcinoma


Tumor microenvironment


Programmed cell death


Cytotoxic T lymphocytes


Non-negative matrix factorization


Differentially expressed mRNAs


Kyoto Encyclopedia of Genes and Genomes


Gene set variation analysis


Single-sample gene set enrichment analysis


Logistic least absolute shrinkage and selector operation


Weighted gene co-expression network analysis


Hazard ratios


Average silhouette width


Follicular helper


Regulatory T cells


  1. Saenz-Ponce N, Pillay R, de Long LM, Kashyap T, Argueta C, Landesman Y, Hazar-Rethinam M, Boros S, Panizza B, Jacquemyn M, et al. Targeting the XPO1-dependent nuclear export of E2F7 reverses anthracycline resistance in head and neck squamous cell carcinomas. Sci Transl Med. 2018;10(447):eaar7223.

    Article  Google Scholar 

  2. Leemans CR, Snijders PJF, Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. 2018;18(5):269–82.

    Article  CAS  Google Scholar 

  3. Ferris RL, Blumenschein G Jr, Fayette J, Guigay J, Colevas AD, Licitra L, Harrington K, Kasper S, Vokes EE, Even C, et al. Nivolumab for recurrent squamous-cell carcinoma of the head and neck. N Engl J Med. 2016;375(19):1856–67.

    Article  Google Scholar 

  4. Oliva M, Spreafico A, Taberna M, Alemany L, Coburn B, Mesia R, Siu LL. Immune biomarkers of response to immune-checkpoint inhibitors in head and neck squamous cell carcinoma. Ann Oncol. 2019;30(1):57–67.

    Article  CAS  Google Scholar 

  5. Hanna GJ, Lizotte P, Cavanaugh M, Kuo FC, Shivdasani P, Frieden A, Chau NG, Schoenfeld JD, Lorch JH, Uppaluri R, et al. Frameshift events predict anti-PD-1/L1 response in head and neck cancer. JCI Insight. 2018;3(4):e98811.

    Article  Google Scholar 

  6. Zheng Y, Tian X, Wang T, Xia X, Cao F, Tian J, Xu P, Ma J, Xu H, Wang S. Long noncoding RNA Pvt1 regulates the immunosuppression activity of granulocytic myeloid-derived suppressor cells in tumor-bearing mice. Mol Cancer. 2019;18(1):61.

    Article  Google Scholar 

  7. Huang D, Chen J, Yang L, Ouyang Q, Li J, Lao L, Zhao J, Liu J, Lu Y, Xing Y, et al. NKILA lncRNA promotes tumor immune evasion by sensitizing T cells to activation-induced cell death. Nat Immunol. 2018;19(10):1112–25.

    Article  CAS  Google Scholar 

  8. Jiang R, Tang J, Chen Y, Deng L, Ji J, Xie Y, Wang K, Jia W, Chu WM, Sun B. The long noncoding RNA lnc-EGFR stimulates T-regulatory cells differentiation thus promoting hepatocellular carcinoma immune evasion. Nat Commun. 2017;8:15129.

    Article  CAS  Google Scholar 

  9. Ma H, Chang H, Yang W, Lu Y, Hu J, Jin S. A novel IFNα-induced long noncoding RNA negatively regulates immunosuppression by interrupting H3K27 acetylation in head and neck squamous cell carcinoma. Mol Cancer. 2020;19(1):4.

    Article  CAS  Google Scholar 

  10. Li H, Xiong HG, Xiao Y, Yang QC, Yang SC, Tang HC, Zhang WF, Sun ZJ. Long non-coding RNA LINC02195 as a regulator of MHC I molecules and favorable prognostic marker for head and neck squamous cell carcinoma. Front Oncol. 2020;10:615.

    Article  Google Scholar 

  11. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44(8):e71.

    Article  Google Scholar 

  12. Li Y, Jiang T, Zhou W, Li J, Li X, Wang Q, Jin X, Yin J, Chen L, Zhang Y, et al. Pan-cancer characterization of immune-related lncRNAs identifies potential oncogenic biomarkers. Nature communications. 2020;11(1):1000.

    Article  CAS  Google Scholar 

  13. Xu T, Le TD, Liu L, Su N, Wang R, Sun B, Colaprico A, Bontempi G, Li J. Cancer subtypes: an R/Bioconductor package for molecular cancer subtype identification, validation and visualization. Bioinformatics. 2017;33(19):3131–33.

    Article  CAS  Google Scholar 

  14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  Google Scholar 

  15. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.

    Article  CAS  Google Scholar 

  16. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7.

    Article  Google Scholar 

  17. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53.

    Article  CAS  Google Scholar 

  18. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108-10.

    Article  CAS  Google Scholar 

  19. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  Google Scholar 

  20. Roh W, Chen PL, Reuben A, Spencer CN, Prieto PA, Miller JP, Gopalakrishnan V, Wang F, Cooper ZA, Reddy SM, et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci Transl Med. 2017;9(379):eaah3560

    Article  Google Scholar 

  21. Hoshida Y, Brunet JP, Tamayo P, Golub TR, Mesirov JP. Subclass mapping: identifying common subtypes in independent disease data sets. PLoS ONE. 2007;2(11):e1195.

    Article  Google Scholar 

  22. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Statis Softw. 2010;33(1):1–22.

    Google Scholar 

  23. Liaw A, Wiener M. Classification and regression by random forest. R news. 2002;2(3):18–22.

    Google Scholar 

  24. Zeestraten EC, Reimers MS, Saadatmand S, Goossens-Beumer IJ, Dekker JW, Liefers GJ, van den Elsen PJ, van de Velde CJ, Kuppen PJ. Combined analysis of HLA class I, HLA-E and HLA-G predicts prognosis in colon cancer patients. Br J Cancer. 2014;110(2):459–68.

    Article  CAS  Google Scholar 

  25. Chen YP, Wang YQ, Lv JW, Li YQ, Chua MLK, Le QT, Lee N, Colevas AD, Seiwert T, Hayes DN, et al. Identification and validation of novel microenvironment-based immune molecular subgroups of head and neck squamous cell carcinoma: implications for immunotherapy. Ann Oncol. 2019;30(1):68–75.

    Article  Google Scholar 

  26. Narayanan S, Kawaguchi T, Yan L, Peng X, Qi Q, Takabe K. Cytolytic activity score to assess anticancer immunity in colorectal cancer. Ann Surg Oncol. 2018;25(8):2323–31.

    Article  Google Scholar 

  27. Pruneri G, Vingiani A, Denkert C. Tumor infiltrating lymphocytes in early breast cancer. Breast. 2018;37:207–14.

    Article  Google Scholar 

  28. Li ZX, Zheng ZQ, Wei ZH, Zhang LL, Li F, Lin L, Liu RQ, Huang XD, Lv JW, Chen FP, et al. Comprehensive characterization of the alternative splicing landscape in head and neck squamous cell carcinoma reveals novel events associated with tumorigenesis and the immune microenvironment. Theranostics. 2019;9(25):7648–65.

    Article  CAS  Google Scholar 

  29. Pusapati RV, Daemen A, Wilson C, Sandoval W, Gao M, Haley B, Baudy AR, Hatzivassiliou G, Evangelista M, Settleman J. mTORC1-dependent metabolic reprogramming underlies escape from glycolysis addiction in cancer cells. Cancer Cell. 2016;29(4):548–62.

    Article  CAS  Google Scholar 

  30. Malakar P, Stein I, Saragovi A, Winkler R, Stern-Ginossar N, Berger M, Pikarsky E, Karni R. Long noncoding RNA MALAT1 regulates cancer glucose metabolism by enhancing mTOR-mediated translation of TCF7L2. Cancer Res. 2019;79(10):2480–93.

    Article  CAS  Google Scholar 

  31. Ho PC, Bihuniak JD, Macintyre AN, Staron M, Liu X, Amezquita R, Tsui YC, Cui G, Micevic G, Perales JC, et al. Phosphoenolpyruvate is a metabolic checkpoint of anti-tumor t cell responses. Cell. 2015;162(6):1217–28.

    Article  CAS  Google Scholar 

  32. Li W, Tanikawa T, Kryczek I, Xia H, Li G, Wu K, Wei S, Zhao L, Vatan L, Wen B, et al. Aerobic glycolysis controls myeloid-derived suppressor cells and tumor immunity via a specific CEBPB isoform in triple-negative breast cancer. Cell Metab. 2018;28(1):87-103.e106.

    Article  CAS  Google Scholar 

  33. Li X, Wenes M, Romero P, Huang SC, Fendt SM, Ho PC. Navigating metabolic pathways to enhance antitumour immunity and immunotherapy. Nat Rev Clin Oncol. 2019;16(7):425–41.

    Article  CAS  Google Scholar 

  34. Brand A, Singer K, Koehl GE, Kolitzus M, Schoenhammer G, Thiel A, Matos C, Bruss C, Klobuch S, Peter K, et al. LDHA-associated lactic acid production blunts tumor immunosurveillance by T and NK cells. Cell Metab. 2016;24(5):657–71.

    Article  CAS  Google Scholar 

  35. Shime H, Yabu M, Akazawa T, Kodama K, Matsumoto M, Seya T, Inoue N. Tumor-secreted lactic acid promotes IL-23/IL-17 proinflammatory pathway. J Immunol. 2008;180(11):7175–83.

    Article  CAS  Google Scholar 

  36. Stoll G, Pol J, Soumelis V, Zitvogel L, Kroemer G. Impact of chemotactic factors and receptors on the cancer immune infiltrate: a bioinformatics study revealing homogeneity and heterogeneity among patient cohorts. Oncoimmunology. 2018;7(10):e1484980.

    Article  Google Scholar 

Download references


Not applicable.


This work was supported by the National Natural Science Foundation of China (Nos. 81671000, 81870769); Guangdong Financial fund for High-Caliber Hospital Construction (174-2018-XMZC-0001-03-0125/D-05).

Author information

Authors and Affiliations



XJ and CB designed the study. CRY performed all the bioinformatics analysis described here. CL performed in vitro experiments. CRY, CL and ZJY wrote and edited the manuscript. CL, ZJY and RXY collected and examined the data. XJ and CB supervised the project. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Bin Cheng or Juan Xia.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

Validating the HNSCC subtypes in validation cohort.

Additional file 2: Fig. S2.

Boxplot of Immune cell infiltration based on MCP method in training cohort and validation cohort.

Additional file 3: Table S1.

The formula of the 31-lncRNA signature.

Additional file 4: Table S2.

The correlation between TRG-AS1 and immune cell infiltration in multiple tumor.

Additional file 5: Fig. S3.

Survival analysis of TRG-AS1 in the anti PD-L1 treatment cohort (IMvigor210).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cao, R., Cui, L., Zhang, J. et al. Immune-related lncRNA classification of head and neck squamous cell carcinoma. Cancer Cell Int 22, 25 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: