Immune and stromal scores were not correlate with clinical data and prognosis of OSCC
To explore the relationship between Immune and stromal scores and the clinical characteristics and prognosis of OSCC patients. In this study, a total of 19467 RNAs were extracted from the OSCC-TCGA RNA-seq database. The immune and stromal scores of 319 OSCC patients were determined using the ESTIMATE algorithm. ESTIMATE scores ranged between -2340.425 and 4674.735, and immune scores ranged between -392.987 and 2705.005, and stromal scores ranged between -1947.438 and 1969.731. To explore the relationship between the immune, stromal scores, and overall survival (OS) in OSCC patients, the patients were categorized into the high and a low score group based on the ESTIMATE, immune, and stromal scores. Kaplan–Meier survival curves were used to calculate the OS of OSCC patients, and we found that the ESTIMATE, immune, and stromal scores are positively correlated with the OS of OSCC patients (P = 0.635, P = 0.373, P = 0.668, respectively) (Fig. 1A, 1B and 1C). Furthermore, for the clinical data extracted from the TCGA-OSCC database, we evaluated the association of the ESTIMATE, immune and stromal scores, and T stage, and we found that ESTIMATE (P = 0.272), immune (P=0.208) and stromal (P=0.163) scores were not significantly related to T stage (Fig. 2A-C). Moreover, we also compared the immune, stromal, and ESTIMATE scores with the gender of OSCC patients, and we found that there was no significant difference between ESTIMATE (Fig. 2D) or stromal (Fig. 2E) scores with gender. However, male patients had lower immune scores compared with the females (Figure 2E). We also noticed that there were no significant differences in ESTIMATE, immune and stromal scores with tumor grade (Fig. 2G-I). However, ESTIMATE (Fig. 2G) and immune scores (Figure 2I) of high-grade tumors were higher than low-grade tumors. The above results indicate that immune and stromal scores are not significantly related to the clinical characteristics and prognosis of OSCC patients.
Identification of immune-related differentially expressed genes
To uncover TME-related differentially expressed genes (DEGs), we performed “limma” in the R package. As a result, a total of 1,589 DEGs were found to be differently expressed between low- and high-stromal score groups (Fig. 3A), while the heat map (Fig. 3B) showed that 1,625 genes were differentially expressed in OSCC patients based on the ESTIMATE, immune and stromal scores. Compared with the low-TME score group, there were 2,171 (82.83%) (Fig. 3C) upregulated genes and 450 (17.17%) (Fig. 3D) down-regulated genes in the high-immune score group. These data indicate that immune-related differentially expressed genes had been obtained and identified.
Functional annotation of immune-related DEGs
To explore the biological functions of the above DEGs, GO term and KEGG pathway enrichment were performed using the R package “clusterProfiler”. As shown in Fig. 4A, the top 10 GO terms of the up and down-regulated DEGs are related to different biological process, cellular component, and molecular functions. Interestingly, the top 5 terms related to the biological process included lymphocyte activation, regulation of lymphocyte-mediated immunity, phagocytosis, complement activation, and B cell-mediated immunity (Fig. 4A, top panel). In addition, the top 5 terms related to a cellular component included immunoglobulin complex, external side of plasma membrane, plasma membrane signaling receptor complex, T cell receptor, and immunoglobulin complex and circulating (Fig. 4A, middle panel). The top 5 terms related to molecular functions included antigen binding, immunoglobulin receptor binding, carbohydrate-binding, immune receptor activity, and cytokine receptor activity (Fig. 4A, bottom panel). On the other hand, the top 30 KEGG pathways for the DEGs related to the up and down-regulated DEGs included immune/inflammation related pathways, including B cell receptor signaling pathway, T cell receptor signaling pathway, intestinal immune network for IgA production, natural killer cell-mediated cytotoxicity, and C-type lectin receptor signaling pathway (Fig. 4B). The above data indicate that the signaling pathways and biological behaviors enriched by immune related DEGs are related to tumor microenvironment or immune functions.
Protein-Protein Interaction (PPI) analysis of the immune-related DEGs
To study the interaction between DEGs, and systematically conduct a comprehensive study on the crossover gene, all 167 intersection genes were uploaded to the STRING database to construct a PPI network. In our present study, the Search Tool Retrieval of Interacting Genes (STRING) database was used to draw the PPI network (Fig. 5). There was a total of 23 genes in the network with more than 20 interconnected nodes, including SUCNR1, P2RY13, GPR83, GALR2, FPR3, CCR8, CCL13, C5AR1, ADORA3, P2RY12, CXCR3, CX3CR1, CCR4, ADRA2A, CCR1, FPR1, CCR2, CCR5, ITGB2, FCER1G, C3AR1, C3 and FPR2 (Figure 6).
Construction of the risk score prognosis model
To analyze the survival-related immune-DGEs, Univariate Cox proportional hazards regression model survival analyses was performed. a total of 23 prognostic genes were identified as survival related DEGs: AC093278.2, AC103563.1, AL365361.1, BTLA, CCL22, CCR4, CD5, CELF2, F5, FCRL3, FLT3, GALR2, IGKV1D-8, IGLV1-36, IGLV4-60, IL10, LINC00861, LINC01508, MS4A2, P2RY14, P2RY8, PPP1R16B and RUBCNL(Fig. 7A). The least absolute shrinkage and selection operator (LASSO) was then used to exclude the over-fitting false positives data (Figure 7B-C). Finally, a prognostic model with 11 survival related DEGs was constructed for predicting the OS of OSCC patients. The risk scores were calculated for each patient as follows: risk score = -0.0418* AC103563.1 - 0.0289* CCL22 - 0.2829* FLT3 + 0.03681* GALR2+ 0.0048 * IGKV1D.8 + 0.0067* IGLV1.36 - 0.0139* IGLV4.60 - 0.0456* IL10 - 0.2040* LINC00861 + 0.1421* LINC01508 - 0.0245* MS4A2. The median value was chosen as the cut-off value and the entire cohort was classified into two risk groups accordingly.
Survival analysis of immune-related model
To test the effect of the above model on the evaluation of OSCC patients in the future, we analyzed the outcome survival (OS) of the high-risk group and the low-risk group OSCC patients. The risk distribution, survival status, and gene expression pattern are shown in Fig. 8 A-C. The scatter plot (Fig. 8A) shows that most of the patients in the high-risk score group died and most of the patients in the low risk group survived during 15 years follow-up. The heat map (Fig. 8C) shows that 7 ARGs were highly expressed in the low risk group while 4 ARGs were highly expressed in the high-risk group. To determine the relationship between immune-related DEGs and OS in OSCC patients, we used the Kaplan-Meier method to plot the survival curves using data obtained from the OSCC-TCGA database. Kaplan-Meier plots show that patients in the high-risk score group have a significantly poorer OS than those in the low risk score group (Fig. 8D). Moreover, the AUC for 1-year, 2-year, and 3- year PFS were 0.64, 0.63, and 0.65 respectively (Fig. 8E). Thus, these data suggest that the 11 DEGs prognostic model performs well in predicting OS in OSCC patients.