Differentially Infiltrated Identification of Novel Diagnostic Biomarkers Associated with Immune Infiltration in Nasopharyngeal Carcinoma

Abstract

Background: The prognostic value of tumor-infiltrating immune cells has been widely studied in nasopharyngeal carcinoma (NPC). However, the role of tumor-infiltrating immune cells in the diagnosis of NPC has not been fully elucidated. Thus, tumor-infiltrating immune cell-related biomarkers in the diagnosis of NPC patients were explored in the current study.

Methods: Gene expression profiles of NPC patients were downloaded from the Gene Expression Omnibus (GEO) database. Differentially infiltrating immune cells (DDICs) between NPC and control samples were analyzed by the CIBERSORT algorithm. Weighted gene coexpression network analysis (WGCNA) was performed to screen hub genes significantly correlated with DDIC. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of hub genes were performed with R package clusterProfiler. The diagnostic value of hub genes was evaluated by receiver operating characteristic (ROC) curves. RT-qPCR was conducted to validate the expression patterns of diagnostic markers in NPC and adjacent control tissues. The correlations between diagnostic markers and immunomodulators were analyzed using the TISIDB. The protein-protein interaction (PPI) network based on immunomodulators significantly associated with diagnostic biomarkers was constructed and visualized by STRING. The functional enrichment analysis of genes in the PPI network was analyzed by the WebGestalt online tool.

Results: The abundances of memory B cells, plasma cells, follicular helper T cells, activated NK cells, M0 macrophages, M1 macrophages, M2 macrophages, resting mast cells, and activated mast cells were significantly different between NPC and control samples. Dark orange was identified as the hub module, with a total of 371 genes associated with memory B cells, plasma cells, and M0 and M1 macrophages defined as hub genes, which were enriched into immune-related biological processes and pathways. FCER2, KHDRBS2, and IGSF9 were considered diagnostic biomarkers with areas under ROC curves as 0.985, 0.978, and 0.975, respectively. Moreover, real-time reverse transcriptase-polymerase chain reaction (RT-qPCR) suggested that the expression patterns of FCER2, KHDRBS2, and IGSF9 were consistent with the results in GEO datasets. TISIDB analysis revealed that FCER2, KHDRBS2, and IGSF9 had a strong association with 8 immunoinhibitors (BTLA, CD160, CD96, LAG3, PDCD1, TIGIT, CD244, and TGFB1) and 11 immunostimulators (CD27, CD28, CD40LG, CD48, ICOS, KLRC1, KLRK1, TMIGD2, TNFRSF13C, CXCR4, and C10 or f54). The PPI network implied that these 19 immunomodulators had interactions with other 50 genes. WebGestalt analysis demonstrated that 69 genes in the PPI network were enriched into cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, and pathways in cancer.

Conclusion: Our study identified novel diagnostic biomarkers and revealed potential immune-related mechanisms in NPC. These findings enlighten the investigation of the molecular mechanisms of tumor-infiltrating immune cells regulating NPC.

Figure 1

(a) The distribution of the immune cells in NCNT. (b) The distribution of the immune cells in NPC. (c) Comparison of 9 immune cells with significant difference between NPC and NCNT (p < 0.01).

Figure 2

(a) The sample clustering result. (b) The sample dendrogram and trait heatmap. (c) The pick soft threshold function of WGCNA. (d) Cluster dendrogram. (e) Screening of the hub module associated with immune infiltration in NPC. Each row represents a color-coded module eigengene; each column represents a type of infiltrating immune cells. The number in each cell means the correlation coefficient and p value.

Figure 3

(a–d) Genes related to memory B cells, plasma cells, M0 macrophages, and M1 macrophages. (e) 371 hub genes associated with immune cell infiltration were identified.

Figure 4

(a, b) B cell activation, lymphocyte differentiation, immune response activating cell surface receptor signaling pathway, immune response activating signal transduction, and B cell differentiation, and corresponding genes are visualized in these terms. (c, d) The first three pathways of KEGG, B cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, primary immunodeficiency, and hub genes involved in these three pathways are shown in the circular diagram.

Figure 5

(a) Compared with NCNT, 7 hub genes were upregulated and 43 hub genes were downregulated in NPC samples. (b) The diagnostic effectiveness of candidate biomarkers was evaluated by the ROC curve. (c) Differential expression of candidate diagnostic markers in NPC and NCNT in the GSE53819 dataset. (d) In the GSE12452 dataset, the expression of candidate diagnostic biomarkers. (e) The expression levels of candidate genes isolated from NPC and NCNT patients were analyzed by RT-qPCR.

Figure 6

(a) The heatmap of correlations between immunostimulators and target genes. (b) The heatmap of correlations between immunoinhibitors and target genes. ^∗ stands for p < 0.05, and ^∗∗ stands for p < 0.01. (c) PPI network composed of 19 immunomodulators and 50 genes interacting with them. Purple represents 11 immune enhancers, green represents 8 immunosuppressants, and pink represents 50 mRNA. (d) Bar charts of biological process categories, cellular component categories, and molecular function categories. (e)Bubble plot of 69 genes significantly enriched cancer and immune-related biological processes and pathways.