Identification of a novel ceRNA network related to prognosis and immunity in HNSCC based on integrated bioinformatic investigation

Abstract

Head and neck squamous cell carcinoma (HNSCC) is characterized by an immunosuppression environment and necessitates the development of new immunotherapy response predictors. The study aimed to build a prognosis-related competing endogenous RNA (ceRNA) network based on immune-related genes (IRGs) and analyze its immunological signatures. Differentially expressed IRGs were identified by bioinformatics analysis with Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA) and ImmPort databases. Finally, via upstream prognosis-related microRNAs (miRNAs) and long noncoding RNAs (lncRNAs) prediction and co-expression analysis, we built an immune-related ceRNA network (LINC00052/hsa-miR-148a-3p/PLAU) related to HNSCC patient prognosis. CIBERSORT analysis demonstrated that there were substantial differences in 11 infiltrating immune cells in HNSCC, and PLAU was closely correlated with 10 type cells, including T cells CD8+ (R = - 0.329), T cells follicular helper (R = - 0.342) and macrophage M0 (R = 0.278). Methylation and Tumor Immune Dysfunction and Exclusion (TIDE) analyses revealed that PLAU upregulation was most likely caused by hypomethylation and that high PLAU expression may be associated with tumor immune evasion in HNSCC, respectively.

Figure 1

Flow chart of establishing and analyzing the immune-related ceRNA network.

Figure 2

Assessment of differentially expressed genes (DEGs) in our analysis. (A, B) Volcano plots showing DEGs in TCGA-HNSC and GSE6631 (|log2 (Fold Change)| > 1, adjust p-value < 0.05); blue dots: significantly down-regulated; red dots: significantly up-regulated; grey dots: no significant differences; the top 10 significant genes are denoted. (C, D) Heatmaps showing DEGs of TCGA-HNSC and GSE6631. (E) Venn diagram showing detail information about overlaps of DEGs across four datasets. (F) The hub immune-related genes (IRGs) in our analysis.

Figure 3

GO and KEGG pathway enrichment of IRGs. (A) GO enrichment; (B) KEGG pathway enrichment.

Figure 4

Hub IRGs expression and prognosis of the HNSCC patients. (A–F) Validation of PLUA, SPP1, S100A8, S100A9, SPINK5 and ACKR1 expression roles and prognosis values using GEPIA and KM plotter (*p < 0.05).

Figure 5

(A) The volcano plot of the differentially expressed miRNAs (|log (Fold Change)| > 0.5, adjust p-value < 0.05); blue dots: significantly down-regulated; red dots: significantly up-regulated; grey dots: no significant differences. (B) The mRNA-miRNA networks were identified using Tarbase and starBase; blue ellipse: miRNAs; red triangle: mRNAs. (C) The volcano plot of the differentially expressed lncRNAs (|log (Fold Change)| > 0.5, adjust p-value < 0.05); blue dots: significantly down-regulated; red dots: significantly up-regulated; grey dots: no significant differences. (D) The miRNA-lncRNA networks were identified using miRNet and starBase; blue ellipse: lncRNAs; red triangle: miRNAs. (E) Sankey diagram for the ceRNA network in HCC. Sankey diagram showing the potential lncRNA–miRNA–mRNA regulatory axes.

Figure 6

(A–C) Co-expression correlation analysis of qualified ceRNA network. (D) Schematic representations of the immune-related ceRNA network. (E) Kaplan–Meier survival curve by the risk score of the TCGA-HNSC dataset. (F) Time-dependent ROC curve analysis for survival prediction based on the risk score.

Figure 7

(A) The expression of LINC00052 in pan-cancer. (B) Forest plot of LINC00052 in multiple tumors was analyzed by univariate cox regression. (C) The correlation between expression of LINC00052 in different tumors and immune cell infiltration scores. (*p < 0.05) (D) The cellular localization for LINC00052 was predicted using lncLocator.

Figure 8

(A) Differential expression of three DNA methyltransferases (DNMT1, DNMT3A, and DNMT3B). (B) Methylation was evaluated using UALCAN. (C) The methylation site of PLAU DNA sequence association with gene expression was visualized using MEXPRESS. The top 4 most significant methylation sites are marked and were negatively correlated with PLAU expression. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 9

Analysis of immune infiltration. (A) The composition of 22 immune cells was estimated by CIBERSORT in HNSCC. (B) Difference in the proportions of 22 immune cells between normal and tumor tissues (*p < 0.05, **p < 0.01, ***p < 0.001). (C) The correlation index among immune cells and ceRNA in HNSCC. Positive and negative correlations are represented by the red and blue colors, respectively. The degree of correlation index is represented by the color depth. The yellow color shows the genes with the highest positive and negative correlations. (D) Kaplan–Meier plots were used to analyze the immune infiltration and overall survival rate of HNSCC.

Figure 10

(A) Immune cell score heatmap, clustered by their relative expression of the PLAU. (B) Correlation between PLAU expression and immune infiltration levels. (C) The distribution of 8 immune checkpoint genes’ expression in high- and low-PLAU expression of HNSCC. (D) TIDE scores between high- and low-PLAU expression of HNSCC. *p < 0.05, **p < 0.01, ***p < 0.001.