We downloaded gene expression profiles and clinical information of all 546 HNSCC patients from the TCGA database. Among them, 63(11.54%) patients were diagnosed with Grade 1, 310(56.78%) patients were Grade 2, 125(22.89%) patients were Grade 3, 7(1.28%) patients were Grade 4, and 41(7.51%) patients were of unknown grade. Based on ESTIMATE algorithm, stromal scores ranged from -1941.57 to 1947.28, and immune scores were distributed between -1057.57 to 2785.18, respectively (Figure 1a, b).
To find out the potential correlation of overall survival with immune scores and/or stromal scores, we divided the 500 HNSCC cases (only 500 of 502 cancer tissues samples have their grade clinical information) into top and bottom halves (high vs. low score groups) based on their scores. Kaplan-Meier survival curves (Figure 1c) showed that median overall survival of cases with the low score group of immune scores are shorter than the cases in the high score group (p = 0.107 in log-rank test). Consistently, cases with lower stromal scores also showed shorter median overall survival compared to patients with higher stromal scores (Figure 1d, p= 0.814 in log-rank test), although it was not statistically significant.
Comparison of gene expression profile with immune scores and stromal scores in HNSCC
To reveal the correlation of global gene expression profiles with immune scores and/or stromal scores, we compared Affymetrix microarray data of these 500 HNCSS cases obtained in TCGA database. Heatmaps in Figure 2 showed distinct gene expression profiles of cases belong to high vs. low immune scores/stromal scores groups. For comparison based on immune scores, 777 genes were upregulated and 161 genes downregulated in the high score than the low score group (fold change > 1.5, p < 0.05). Similarly, for the high and low groups based on stromal scores, 986 genes were upregulated and 63 genes were downregulated in the high score group (fold change > 1.5, p < 0.05). Moreover, Venn diagrams (Figure 2c, d) showed that 260 genes were commonly upregulated in the high scores groups, and 15 genes were commonly downregulated. It is worth mentioning that the DEGs extracted from the comparison of high vs. low immune scores groups covered the majority of genes extracted from the comparison based on stromal scores. Thus, we decided to focus on these DEGs for all subsequent analysis in this manuscript. To outline the potential function of the DEGs, we performed functional enrichment analysis of the 260 upregulated genes and 15 downregulated genes in high-immune scores group. Functional enrichment clustering of these genes showed strong association with immune response as well. Top GO terms identified neutrophil activation involved in and mediated immune response, extracellular matrix, and immunoglobulin binding (Figure 2e). In addition, all the pathways that were yielded from the KEGG identified Cytokine-cytokine receptor interaction, Complement and coagulation cascades and B cell receptor signaling pathway (Figure 2f)
Correlation of expression of individual DEGs in overall survival
To explore the potential roles of individual DEGs in overall survival, we generated Kaplan-Meier survival curves from TCGA database. Among the 275 DEGs in the high-immune scores group, a total of 59 DEGs were shown to significantly predict poor overall survival in log-rank test (p < 0.001, selected 6 genes are shown in Figure 3).