Log2(FPKM+1) gene expression RNA-seq data, somatic mutation (MuTect2) data, copy number (gene-level) data, and log2(RPM+1) miRNA mature strand expression RNA-seq Data, of Head and Neck Cancer samples of The Cancer Genome Atlas (TCGA), were downloaded from the UCSC Xena browser (https://gdc.xenahubs.net) as well as the corresponding clinical and survival information (8–10). The inclusion criteria for LSCC patients in the TCGA database were the patients whose sites of rection were larynx. Normal samples, recurrent tumor samples, or paraffin-embedded samples were excluded. We have added the inclusion and exclusion criteria in the method section as suggested. RNA-Seq data of GSE127165 and miRNA-Seq data of GSE133632 in the GEO database (https://www.ncbi.nlm.nih.gov/gds/) were used as the validation cohorts, and the data were also log2 conversed before being analyzed (11).

As previously reported (12–14), the expression profile of 585 immune cell infiltration related genes was used to calculate the infiltrations of the 24 immune cell populations based on the Gene Set Variation Analysis (GSVA) method using R software (version 4.0.5). A matrix containing the infiltration enrichment scores ranging from -1 to 1 for each immune cell type in every tumor sample was obtained. Then tumor samples with different immune cell infiltrating patterns were clustered and grouped by an unsupervised clustering method.

The clinical characteristics of the immune infiltration groups were compared using t-test, Wilcoxon, and chi-square tests. Kaplan-Meier method, log-rank test, univariate and multivariate cox analyses were used to investigate the survival. All the statistical analyses here were performed using R.

Somatic mutations and copy number variation were analyzed using the maftools of R. The threshold of p-value < 0.01 was considered significant. Differentially expressed genes (DEGs), miRNAs, and lncRNAs were identified in different immune infiltration phenotypes using the package limma of R, with a threshold of Log (fold change) > 0.5, p-value < 0.05, and adjusted p-value < 0.05.

Functional enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed with the clusterProfiler package of R, with a cutoff of adjusted P < 0.01 and false discovery rate (FDR) < 0.05. ceRNA network was built and visualized by Cytoscape (version 3.7.1) to further explore the relationship among DEGs, the differentially expressed miRNAs and lncRNAs. The regulations in the ceRNA network were predicted by miRwalk and miRcode. Only the targets which were down-regulated when the miRNAs were up-regulated were included and vice versa.

The LASSO regression was used in this study to select the optimal set of genes to calculate to which immune infiltration phenotype a patient belongs. Ten-fold cross-validation was adopted using the glmnet package in R to determine the optimal parameter λ and the corresponding set of genes. Then binary logistic regression was applied to construct the model using the expression of selected genes. At the same time, the ROC curve was used to determine the best cutoff value of the model and calculate its accuracy and specificity.

A total of 111 LSCC patients in the TCGA head and neck cancer cohort were enrolled in this study. Using the GSVA enrichment method, the infiltration of each immunity cells in each patient was calculated and then clustered into two different groups (high-infiltration group, n = 61, and low-infiltration group, n = 50) ( Figure 1A ). The infiltration scores of 21/24 immune cells were higher in the high-infiltration group than those in the low-infiltration group ( Figure 1B ).

The clinical characteristics were compared between the two immune infiltration phenotypes. As shown in Table 1 and Figure 1A , the N stage of the patients in the high-infiltration group was significantly lower (p = 0.048) as well as the TNM stage (p = 0.025), while there was no significant difference between the two infiltration groups regarding race, sex, age, grade, T stage, M stage, smoking or alcohol status. The 1-year, 3-year, and 5-year survival of patients with high-infiltration phenotype were 87.7 ± 4.4%, 68.0 ± 6.7%, and 61.5 ± 7.5%, respectively, significantly better than those of patients with low-infiltration (85.2 ± 5.2%, 46.1 ± 8.2%, 41.0 ± 8.7%, respectively) (p = 0.042) ( Figure 1C ).

One hundred and nine LSCC patients with gene mutation data were included in this analysis. As Figure 2 shows, the gene mutation status was quite different between patients with different immune infiltration phenotypes. There was a median of 120 and 161.5 mutations in each patient with a high-infiltration phenotype and a low-infiltration group, respectively. The most common DNA base mutation type was C>T in the high-infiltration group, but it was C>A in the low infiltration group ( Supplementary Figure 1 ).

The most commonly mutated genes in the high-infiltration group were TP53 (85%), TTN (46%), CSMD3 (36%), KMT2D (27%), and PCLO (27%), while in the low-infiltration group they were TP53 (84%), TTN (64%), CSMD3 (38%), and MUC16 (34%). In the top 20 mutated genes of either group, COL11A1 (7% vs. 28%, p = 0.004), MUC17 (8% vs. 26%, p = 0.019), FAM135B (12% vs. 30%, p = 0.030), and XIRP2 (10% vs.26%, p = 0.042) were significantly less mutated in the high-infiltration group ( Figure 2 and Supplementary Table 1 ). These differentially mutated genes might contribute to the distinct immune infiltration phenotypes.

DEGs between the high-infiltration and low-infiltration groups were firstly investigated using the TCGA data and further analyzed to identify the potential biological functions. As Figure 3A shows, a total of 558 up-regulated and 82 down-regulated DEGs were identified in the high-infiltration group compared with the low-infiltration group. The top-3 up-regulated genes with the lowest p-values were TNFRSF1B (logFC = 1.33), ABI3 (logFC = 1.04), and SELPLG (logFC = 1.36) (all p < 0.001), while the top-3 down-regulated genes were SNRPG (logFC = -0.55), FKBP1B (logFC = -0.56), and STK26 (logFC = -0.76) (all p < 0.001) ( Figure 3A ). Meanwhile, four targets of immunotherapies, PDCD1 (PD1), CD274 (PD-L1), LAG3, and CTLA4, were all significantly ug-regulated in patients with high immune infiltration phenotype, indiciting those patients might be more suitable for immunotherapy ( Figure 3B ).

The results of GO and KEGG analyses showed that the DEGs were significantly enriched in immune-related gene ontologies and pathways ( Figure 3C ), indicating the huge divergence between the two immune phenotypes.

Meanwhile, seven ug-regulated and nine down-regulated miRNAs were identified, including miR-150-5p (logFC = 1.21), miR-146a-5p (logFC = 0.80), miR-18a-5p (logFC = -0.77) and miR-96-5p (logFC = -0.83) (all p < 0.001) ( Figure 3D ). One hundred and eighty-three up-regulated and nine down-regulated lncRNAs were also identified, including TRBC2 (logFC = 1.57), TRAC (logFC = 1.56), RP11-174G6.1 (logFC = -0.50), and DANCR (logFC = -0.69) (all p < 0.001) ( Figure 3E ).

We searched the GEO database and chose GSE127165 and GSE133632, RNA-Seq and miRNA-Seq data of the same 57 LSCC patients, the largest cohort in the GEO database, as the validation cohort. There is a lack of clinical, survival, or gene mutation information in the GEO dataset. Therefore, unfortunately, we couldn’t validate the results of clinical, survival, and gene mutation analyses in the GEO data.

The 57 patients were also clustered into high-infiltration (n = 35) and low-infiltration (n = 22) groups ( Figure 4A ). The infiltration scores of most immune cells were higher in the high-infiltration group than those in the low-infiltration group, consistent with the results of the TCGA analyses ( Figure 4B ). As Figure 4B shows, 309 up-regulated DEGs and 43 down-regulated DEGs were indetified, including IL10RA (logFC = 1.06), NCKAP1L (logFC = 0.57), STAB1 (logFC = 1.00), NUSAP1 (logFC = -0.52), DEPDC1 (logFC = -0.55), and CENPF (logFC = -0.76) (all p < 0.001) ( Figure 4C ). The expression of PDCD1, LAG3, and CTLA4 were significantly raised in high-infiltration group, too ( Figure 4D ). The DEGs were also siginicantly enriched in lots of immue-related gene ontologies and pathways ( Supplementary Figure 2 ), quite similar as the GO and KEGG results of the TCGA data.

Six up-regulated miRNAs, including miR-150-5p (logFC = 1.68), miR-146a-5p (logFC = 0.94), and miR-148a-3p (logFC = 0.92), as well as 108 down-regulated miRNAs, including miR-2277-5p (logFC = -1.33), miR-301a-5p (logFC= -1.36), and miR-454-5p (logFC = -0.88) (all p < 0.001) were identified ( Figure 4E ). Meanwhile, 21 up-regulated lncRNAs and 12 down-regulated lncRNAs, including RN7SL116P (logFC = 1.11), RP11-284N8.3 (logFC = 1.00), GABPB1-AS1 (logFC = -0.65), and RP11-443B20.1 (logFC = -0.52) (all p < 0.001) were significantly differentially expressed between the two immune infiltration phenotypes ( Figure 4F ).

The intersections of differentially expressed genes, miRNAs, and lncRNAs obtained from the TCGA and GEO data were analyzed. As Supplementary Figures 3A , B shows, 202 genes (197 up-regulated and five down-regulated) and eight miRNAs (two up-regulated and six down-regulated) were differentially expressed both in the TCGA and GEO data, accounting for a large part of obtained DEGs and differentially expressed miRNAs, which demonstrates that the classification of immune infiltration phenotypes in our study is reasonable and universal. However, there was only one lncRNA, AL928768.3, differentially expressed in both two data cohorts ( Supplementary Figure 3C ), possibly due to the low and unstable expression of lncRNAs.

To investigate the regulatory network associated with immune infiltration phenotypes in LSCC patients, we constructed a ceRNA network using the genes, miRNAs, and lncRNAs those differentially expressed both in the TCGA and the GEO data ( Figure 5 and Supplementary Table 2 ). In this net, miR-149-5p regulated 28 DEGs while miR-17-5p, miR-18a-5p regulated 17 and 16 DEGs, respectively. Seven DEGs (BTN3A1, CCR1, CYBA, IRF1, LAIR1, MPEG1, and MS4A6A) were regulated by three miRNAs together.

In the above, we used GSVA and cluster analysis to determine the immune infiltration subtypes of the patients. However, this method is not suitable for analyzing the expression data of only one patient or quite a few patients because cluster analysis cannot be properly done at that time. To solve this problem, we constructed a predictive model using LASSO and binary logistic regression based on the expression of a set of genes to calculate which immune infiltration phenotype a patient belongs to.

Genes differentially expressed both in the TCGA and GEO cohorts were included in the LASSO analysis. The TCGA data were used as the training cohort while the GEO data were the validation cohort. Sixteen DEGs were selected in the logistic regression when lambda.min was chosen as the optimal parameter λ in the LASSO analysis ( Figures 6A , B ). Finally, a model was constructed as follows.

Here the expression value of each gene in this score is log2(FPKM+1). Most of these genes were significantly correlated with survival ( Supplementary Figure 4 ). Based on the ROC curve, the best cutoff value of the score was 0.541, with a 96.3% accuracy, 93.3% specificity, and 0.983 AUC ( Figure 6C ). When the gene expression data of an LSCC sample is substituted into the above formula, and the calculation result is less than 0.541, it indicates that the patient belongs to the high-infiltration phenotype. Otherwise, he belongs to the low-infiltration group.

We calculated the score’s performance in the GEO cohort, and found it also had an excellent performance with a 91.0% accuracy and 97.1% specificity, which demonstrated the universal usability of this score.