Immune cell infiltration was abundant in HNSCC
We performed tumor immune analysis and evaluation on the selected samples through the CIBERSORTx deconvolution algorithm. CIBERSORTx calculated and displayed the infiltration abundance of 22 types of immune cells in the samples. The results indicated that compared with normal samples, the infiltration abundance of immune cells in tumor samples was more obvious (Fig. 1A). At the same time, the infiltration abundance of naive B cells, plasma cells, and monocytes in normal samples were significantly higher than that in tumor samples. The infiltration abundance of the memory CD4 T cells, resting NK cells, M0 macrophages, M1 macrophages, and dendritic cells activated in tumor samples was significantly higher than that of the normal samples. The correlation heatmap indicates a positive correlation between memory CD8 and CD4 T cells, and a positive correlation between CD4 naïve T cells and memory B cells. And a negative correlation between M0 macrophages and both CD8 and CD4 T cells (Fig. 1A-D).
B Lymphocyte Infiltration Of Tumors Affected The Patient Survival Rate
We used the consensus clustering algorithm in the ConsensusClusterPlus package to perform consistent clustering analysis. Through the ConsensusClusterPlus package, we set reps = 1000, pItem = 0.8, and pfeature = 1. Using the delta plots, consensus matrix (CM) plots, item-consensus (IC) plots, and cumulative distribution function (CDF) plots, CC algorithm was visualized to facilitate the finding of the optimal number of classifications. The approximate maximum k value of the distribution reached under the condition of maximum stability was found through the CDF plots, and the relative change in area under the CDF plot curve was represented by delta plots. The results indicated that when the k value was 4 (Fig. 2A), the change in the area under the CDF curve was no longer significant. At this point, k = 4 was concluded to be the optimal number of classifications (Fig. 2B). Then, all cases were divided into four categories using the optimal method and represented by CM plots (Fig. 2C). Each colored rectangular bar in the IC diagram represents each sample, and the height of the rectangle corresponds to the IC value. The consistent clusters were marked with colored asterisks on the rectangle, which is convenient for viewing which sample heights represent one group and which samples are mixed with other groups. The less mixed the color of the color rectangle, the more consistent the group is. The IC plots reveal that when k is greater than or equal to 5, the color of each color bar becomes increasingly mixed; similarly, k = 4 is the optimal classification number (Fig. 2D-E). According to the results, we divided the samples into four groups and performed cluster analysis as the optimal plan.
After analyzing the survival of these four groups using the Survminer package, we found that the survival prognoses of the first and third groups were significantly better than those of the second and fourth groups (Fig. 3A). To explore the reasons for the significant difference in survival prognosis, we classified the first and third groups as cluster A and the second and fourth groups as cluster B (Fig. 3B). Using the edge R package, we performed the differential analysis of the cluster A and cluster B gene expression levels, and 30 genes with significant differential expression between the two clusters were screened with |logFC| > 1 and p < 0.05.
Enrichment Analysis
After differential analysis, the 30 genes with significant differential expression were screened for online gene enrichment analysis using the Metascape website. Five genes, CD79A, IGHD, IGHG3, IGHM, and IGHV2-70, were enriched in the B-cell receptor signaling pathway (Fig. 4A-B, Table 1). We also found that the infiltration abundance of memory B cells and naive B cells in cluster A was significantly higher than that in cluster B. (Fig. 3C-D)
Table 1
Result of multivariate Cox regression
| coef | exp(coef) | se(coef) | z Pr(>|z|) |
CD79Alow | 1.6202 | 5.0541 | 0.8797 | 1.842 | 0.0655. |
IGHMlow | 0.9601 | 2.6119 | 0.5858 | 1.639 | 0.1012 |
Signif. codes: 0 ‘***’0.001 ‘**’0.01 ‘*’0.05 ‘.’0.1 ‘ ’1 |
| exp(coef) | exp(-coef) | lower .95 | upper .95 | |
CD79Alow | 5.054 | 0.1979 | 0.9012 | 28.345 | |
IGHMlow | 2.612 | 0.3829 | 0.8286 | 8.233 | |
Concordance0.615(se = 0.062 ) |
Likelihood ratio test = 7.9 on 2 dfp = 0.02 |
Wald test = 11.64 on 2 dfp = 0.003 |
Score (logrank) test = 17.64 on 2 dfp = 1e-04 |
Survival Prognosis Analysis
Through the "Surv_cutpoint" function in survminer package, an optimal cut point of the infiltration abundance of naive B cells and the memory B-cell infiltration in all samples were selected: 0.000621206 for naive B cells and 0.003937039 for memory B cells (Fig. 5A, C). According to the cut points, the infiltration abundance of cells was divided into high and low abundance, and K-M survival analysis was performed based on the results. The results indicated that patients with high naive B cell infiltration abundance had a better prognosis than those with low infiltration abundance, while patients with a high infiltration abundance of memory B cells had a poorer prognosis than those with low infiltration abundance (Fig. 5B, D).
Similarly, optimal cut points for the gene expression abundance of the five genes enriched in the B-cell receptor signaling pathway of each sample were selected using the “surv_cutpoint” function of the survminer package to create a division between high and low gene expression. Through the forward/backward multivariate Cox regression analysis, we found among the five genes enriched in the B-cell receptor signaling pathway, low expression of CD79A and IGHM significantly affected the survival prognosis of patients (Table 2). The K-M survival analysis was performed based on the high and low gene expression levels. The results indicated that the prognosis of patients with low expression of CD79A and IGHM was significantly lower than that of patients with high expression of IGHM (Fig. 6A-D).
Table 2
five genes with significant differential expression enriched in the B-cell receptor signaling pathway
Gene ID | Log2FoldChange | Log2CPM | P value | FDR |
CD79A | -1.420365627 | 4.929473725 | 0.00000209 | 0.0007769 |
IGHD | -1.214312563 | 5.298169419 | 0.001075435 | 0.042693127 |
IGHG3 | -1.038873442 | 8.403678554 | 0.000975217 | 0.04112342 |
IGHM | -1.464552247 | 8.375177082 | 1.23E-05 | 0.002856521 |
IGHV2-70 | -1.03617177 | 5.143032756 | 0.004168784 | 0.083538477 |
Table 3 Clinical correlation analysis of CD79A and oral squamous cell carcinoma |
| | + | - | p value |
Age | Mean ± SD | 57.68 ± 11.30 | 57.50 ± 11.13 | 0.4948 |
Sex | Male | 48 | 19 | |
| Female | 32 | 18 | 0.9837 |
T staging | 1 | 2 | 7 | |
| 2 | 7 | 19 | |
| 3 | 63 | 9 | |
| 4 | 8 | 2 | 0.0542 |
Degree of differentiation | High | 31 | 14 | |
| Medium | 27 | 11 | |
| Low | 22 | 12 | 0.5434 |
Lymph node metastasis | Yes | 68 | 6 | |
| None | 12 | 31 | 0.011* |
Primary lesion location | Lingual margin | 31 | 15 | |
| Dorsum | 5 | 7 | |
| Gingiva | 11 | 2 | |
| Mouth floor | 9 | 7 | |
| Cheek | 16 | 5 | |
| Palate | 3 | 0 | |
| Lip | 2 | 0 | |
| Oropharynx | 3 | 1 | 0.9373 |
Cd79a Was Closely Related To The Infiltration Abundance Of B Lymphocytes In Hnscc Tissues
We entered CD79A and IGHM into the TIMER2.0 database and used Spearman correlation to perform association analysis. The results revealed that CD79A had a significant correlation with the infiltration abundance of B lymphocytes of a variety of tumors in a variety of databases, such as TIMER, EPIC, QUANTISEQ, XCELL, MEPCOUNTER, and CIBERSORT (Fig. 7A). Among them, we noticed that in the CIBERSORT group, the naive B lymphocytes, memory B lymphocytes, and plasma cells in the HNSCC samples were positively correlated with CD79A (Fig. 7B-D). In the database analysis, no correlation was found between IGHM and immune cell infiltration in tumor tissues.
Clinical Samples And Follow-up Information Indicated That Cd79a Affected The Survival Prognosis Of Patients
From 2015 to 2018, a total of 117 patients with oral squamous cell carcinoma were admitted and treated in the Department of Oral and Maxillofacial Surgery, Second Affiliated Hospital of Shantou University Medical College, and underwent surgical treatment. Intraoperative frozen sections were used to confirm that the surgical boundary was expanded to the normal range. Immunohistochemistry of the pathological tissues after resection revealed that CD79A could be expressed to different degrees in oral squamous cell carcinoma tissues (Fig. 8A-B). According to the follow-up data analysis, the survival prognosis of patients with high CD79A expression was significantly lower than that of patients with low CD79A expression (p < 0.0001) (Fig. 8C). Survival was not significantly correlated with the patient’s sex, age, primary tumor location, clinical stage or degree of tumor differentiation but significantly correlated with lymph node metastasis (p = 0.011) (Table 3). Although the result showed that the clinical stage has no statistical differences with the expression level of CD79A, we still found out that more patients with high CD79A expression were in the stage III (n = 63) compared with patients with low CD79A expression.
In summary, we believe that CD79A expression not only affects the infiltration abundance of B lymphocytes in tumor tissues but can also be used as a potential target to analyze and predict the clinical prognosis of patients.