Data mining on identifying diagnosis and prognosis biomarkers in head and neck squamous carcinoma

<jats:title>Abstract</jats:title>         <jats:p><jats:bold>Background: </jats:bold>Head and neck squamous carcinoma (HNSC) is cancer with a poor prognosis that induces high cancer-related death worldwide. The biomarker screening on diagnosis and prognosis is of great importance. This research is aimed to explore the specific diagnostic and prognostic biomarkers for HNSC through bioinformatics analysis.<jats:bold>Results: </jats:bold>The mutation and dysregulation data were acquired from UCSC Xena and TCGA databases. The top ten genes with mutation frequency in HNSC were <jats:italic>TP53</jats:italic> (66%), <jats:italic>TTN</jats:italic> (35%), <jats:italic>FAT1</jats:italic> (21%), CDKN2A (20%), MUC16 (17%), CSMD3 (16%), PIK3CA (16%), NOTCH1 (16%), SYNE1 (15%), LRP1B (14%). A total of 1,060 DEGs were identified, with 396 up-regulated and 665 downregulated in HNSC patients. Patients with lower expression of <jats:italic>ACTN2</jats:italic> (<jats:italic>P</jats:italic>=0.039, HR=1.3), <jats:italic>MYH1</jats:italic> (<jats:italic>P</jats:italic>=0.005, HR=1.5), <jats:italic>MYH2</jats:italic> (<jats:italic>P</jats:italic>=0.035, HR=1.3), <jats:italic>MYH7</jats:italic> (<jats:italic>P</jats:italic>=0.053, HR=1.3), and <jats:italic>NEB</jats:italic> (<jats:italic>P</jats:italic>=0.0043, HR=1.5) exhibit longer overall survival time in HNSC patients. The main DEGs were further analyzed by pan-cancer expression and immune cell infiltration analyses. We can observe that MYH1, MYH2, and MYH7 were dysregulated in the cancers. Compared with HNSC, their expression levels are lower in the other types of cancers. All the five DEGs have a significant positive correlation with CD4+ T cells and macrophages.<jats:bold>Conclusion: </jats:bold><jats:italic>ACTN2</jats:italic>, <jats:italic>MYH1</jats:italic>, <jats:italic>MYH2</jats:italic>, <jats:italic>MYH7</jats:italic>, and<jats:italic> NEB</jats:italic> genes were significantly associated with poor prognosis in HNSC. Moreover, <jats:italic>MYH1</jats:italic>,<jats:italic> MYH2</jats:italic>, and <jats:italic>MYH7</jats:italic> were expected to be the specific diagnostic and prognostic molecular biomarkers of HNSC.</jats:p>

section_content: Mutations in HNSC patients. The mutation profile of HNSC patients from the TCGA database was obtained. Figure 1A showed that the frequency distribution of different variant types was divided into seven according to the impact on protein-coding. We found that missense mutation accounts for the majority of the variant types. In these variant types, the number of SNPs is significantly more than that of INS and DEL (Fig. 1B ). Besides, C > T transversion was the primary type of single nucleotide variants (SNV) in HNSC (Fig. 1C ). We show genes with a mutation frequency of top25 in Fig. 1D. The top ten genes with a mutation frequency were TP53 (66%), TTN (35%), FAT1 (21%), CDKN2A (20%), MUC16 (17%), CSMD3 (16%), PIK3CA (16%), NOTCH1 (16%), SYNE1 (15%), LRP1B (14%). 

Differentially expression genes calculation. In these patients, the expression levels in two groups were compared, and a total of 1061 DEGs were identified, with 396 up-regulated and 665 downregulated (Fig. 2A, B ). Then, we retained the DEGs with more than 1% mutation frequency. After the screening, 187 common DEGs were obtained (Fig. 2C ), including 66 up-regulated and 121 down-regulated. 

PPI and key module analyses on DEGs. PPI analysis was performed on the common DEGs with more mutation frequency. The PPI network showed 144 nodes and 489 interaction pairs (Fig. 3 ). The DEGs with more than one connection with other DEGs were retained (Supplementary File 1). In addition, two sub-network modules were screened, with 15 nodes and 102 interaction pairs in module A (score = 14. 571), and 14 nodes and 75 interaction pairs in module B (score = 11. 538). 

In addition, we listed the key DEGs in sub-network modules in Table 1.
section_name: Results

section_content: We further analyzed the function of DEGs by GO-BP and KEGG. We found that the DEGs in module A are mainly involved in GO-BP terms, including extracellular matrix organization, extracellular structure organization, collagen fibril organization, endodermal cell differentiation, endoderm formation, endoderm development, formation of the primary germ layer, cellular response to amino acid stimulus, collagen-activated tyrosine kinase receptor signaling pathway, and collagen-activated signaling pathway (Fig. 4A ). They are also involved in KEGG pathways of protein digestion and absorption, ECM-receptor interaction, AGE-RAGE signaling pathway in diabetic complications, amoebiasis, focal adhesion, human papillomavirus infection, PI3K-Akt signaling pathway, Relaxin signaling pathway, small cell lung cancer, and platelet activation (Fig. 4B ). The DEGs in module B mainly participated in muscle contraction, muscle system process, myofibril assembly, cellular component assembly involved in morphogenesis, striated muscle cell development, cardiac muscle cell development, cardiac cell development, muscle filament sliding, actin-myosin filament sliding, and cardiac muscle fiber development (Fig. 4C ). In addition, KEGG pathways of hypertrophic cardiomyopathy, dilated cardiomyopathy, cardiac muscle contraction, thyroid hormone signaling pathway, adrenergic signaling in cardiomyocytes, cGMP-PKG signaling pathway, viral myocarditis, and arrhythmogenic right ventricular cardiomyopathy were enriched (Fig. 4D ). The yellow circle is the down-regulated gene, the blue prism is the up-regulated gene, the node size is based on the degree value, the higher the degree value, the larger the node. lated in HNSC patients (P < 0. 05; Fig. 5 ). Patients with lower expression of ACTN2 (P = 0. 039, HR = 1. 3), MYH1 (P = 0. 005, HR = 1. 5), MYH2 (P = 0. 035, HR = 1. 3), MYH7 (P = 0. 053, HR = 1. 3), and NEB (P = 0. 0043, HR = 1. 5) exhibit longer overall survival time. These results demonstrated that these DEGs might be essential for predicting the overall survival of HNSC.
section_name: Function analysis on DEGs in sub-network modules.

section_content: Then, we evaluated their expression in other cancers (Fig. 6A-E ). ACTN2 was significantly downregulated in bladder urothelial carcinoma (BLCA), kidney chromophobe (KICH), kidney renal papillary cell carcinoma (KIRP), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC) (P < 0. 001); up-regulated in invasive breast carcinoma (BRCA), cholangiocarcinoma (CHOL), and colon adenocarcinoma (COAD). We can observe that although MYH1, MYH2, and MYH7 were dysregulated in some of the cancers, their expression levels are very low, except for HNSC. Therefore, MYH1, MYH2, and MYH7 can be considered specific HNSC biomarkers. As for NEB, it was significantly downregulated in BRCA, HNSC, KICH, KIRP, and THCA; up-regulated in BLCA, CHOL, COAD, esophageal carcinoma (ESCA), liver hepatocellular carcinoma (LIHC), LUAD, LUSC, rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), and UCEC (P < 0. 05; P < 0. 01; P < 0. 001). Moreover, we evaluated the expression levels of ACTN2, MYH1, MYH2, MYH7, and NEB in HPV-positive and -negative HNSC samples. Their expression levels showed no significant differences in HPV-positive and negative HNSC samples (Fig. 6F ).
section_name: Key DEG expression validation, survival curve, and UALCAN

section_content: The correlation between key DEGs and immune infiltrate cells in total HNSC, HNSC-positive, and HNSC-negative samples was evaluated. As for total HNSC samples, we found that ACTN2, MYH1, MYH2, and MYH7 were negatively correlated with B cell, CD8+T cell, and positively correlated with CD4+T cell, macrophage, neutrophil, and dendritic cell (Fig. 7A-E ). All five DEGs have a significant positive correlation with CD4+T cells and macrophages. ACTN2 had positive correlation with CD4+T cell (partial correlation = 0. 253; P = 2. 03e-08), macrophage (partial correlation = 0. 181; P = 6. 03e-05), neutrophil (partial correlation = 0. 157; P = 5. 82e-04), and dendritic cell (partial correlation = 0. 161; P = 4. 02e-04) in HNSC samples. It showed an opposed correlation in HNSC-HPV positive and HNSC-HPV negative in immune cell infiltration of B cell and CD8+T cell, but not significant (Fig. 7A ). Therefore, the HPV infection had no significant influence on immune cell infiltration of ACTN2. 

MYH1 showed positive correlation with CD4+T cell (partial correlation = 0. 201; P = 8. 71e-06), macrophage (partial correlation = 0. 146; P = 1. 31e-03), neutrophil (partial correlation = 0. 095; P = 3. 76e-02), and dendritic cell (partial correlation = 0. 104; P = 2. 30e-02; Fig. 7B ). The MYH1 showed a consistent correlation in the immune cells of HNSC-HPV positive and negative samples, except in neutrophils. Although MYH1 in HNSC-HPV positive sample is negatively correlated with neutrophils, but is not significant. 

As for MYH2, MYH7, and NEB were significantly positively correlated with CD4+T cell (P < 0. 001; Fig. 7C-E ). They showed consistent correlation in the immune cells among the total HNSC, HNSC-HPV positive and negative samples; Some slight correlation opposition was observed in HNSC-HPV positive or negative with total HNSC, but not significant. Therefore, we concluded that no significant differences were identified in HNSC-HPV positive or negative with total HNSC in terms of immune cell infiltration of ACTN2, MYH1, MYH2, MYH7, and NEB.
section_name: Immune cell infiltration correlation analysis.

section_content: HNSC, a dreadful opponent, has become one of the world's most deadly diseases. Treatment for HNSC is costly, has a long cycle, and is not affordable. Therefore, the precise and early prediction of HNSC outcomes is critical for therapy options and prognosis improvement. Our present study focused on gene mutation and dysregulation expression in HNSC, expecting to identify essential biomarkers for HNSC. 

In the present study, we explored the landscape of mutation and dysregulation in HNSC patients, showing that TP53, TTN, FAT1, CDKN2A, and MUC16 were the most predominant mutated genes. The mutation identification was consistent with Jiang et al. 14. Among human genes, TP53 is a critical tumor suppressor gene, with low expression in normal cells and high expression in malignant tumors, regulating cell proliferation, apoptosis, angiogenesis, and DNA repair 15. Mutated TP53 was found in various cancers and is correlated with reduced O. S. Furthermore; it showed that HNSC patients with TP53 mutations have a bleak prognosis than TP53-wildtype HNSC 16. Our present study also found that mutated TP53 is observed in HNSC. 

Then, the dysregulated genes were identified, with 187 DEGs having more mutation sites. We conducted the PPI analysis on these DEGs and identified two modules. We further analyzed these DEGs in these modules by GO-BP and KEGG analyses. The DEGs in module A are involved in the PI3K-Akt signaling pathway, a classic intracellular signaling pathway regulating cell apoptosis, proan essential pathway that regulates cell apoptosis, prognosis, and metastasis in HNSC [17] [18] [19]. The DEGs in module B participated in the cGMP-PKG signaling pathway. In addition, Tuttle et al. 20 demonstrated that the cGMP-PKG signaling pathway is one of the therapeutic targets of HNSC. 

The DEGs of ACTN2, MYH1, MYH2, MYH7, and NEB in module B were significantly associated with O. S. in HNSC. ACTN2 is a member of the spectrin gene superfamily that includes varying groups of cytoskeletal proteins. In breast cancer patients, mutated ACTN2 was related to invasive ductal carcinoma and suggested a worse O. S. than ductal carcinoma in situ 21. Sun et al. 22 demonstrated that ACTN2 is one of the hub genes selected by bioinformatics methods in PTEN mutation prostate cancer. Xu et al. 23 revealed that negative ACTN2 expression contributed to a better O. S. in HNSC patients, consistent with our present study. Therefore, ACTN2 might be an essential biomarker for predicting O. S. of HNSC, which might be employed in clinical. 

Through pan-cancer analysis, we found that MYH1, MYH2, and MYH7 are dysregulated in cancers, including HNSC. However, compared with HNSC, they were lower expressed in other cancers, revealing that these three genes can be considered single HNSC biomarkers. So far, different mutations in multiple members of the MYH gene family are associated with human hereditary diseases 24. Among them, MYH2 mutations can cause a class of skeletal muscle diseases characterized by ophthalmoplegia. MYH7 mutations can cause skeletal muscle diseases, including myosin deposition myopathy, and Distal Laing myopathy is also closely related to hypertrophic cardiomyopathy. However, MYH1, MYH2, and MYH7 functions in HNSC are rarely mentioned. Further study into the interaction of the MYH gene family and HHSC appears to be an intriguing research topic. 

As we all know, HPV infection is the most common causative factor for HNSC. Considering the significant influence of HPV on HNSC, some molecular research on HNSC have divided the samples into HNSC-HPV positive and negative groups. For example, Zhang et al. 25 found that Dedicator of cytokinesis 8 (DOCK8) is identified as one prognosis biomarker associated with immune infiltration HPV-positive HNSC. Our present study identified five prognosis DEGs that are significantly related to the O. S. in HNSC. To determine whether the HPV infection had significant influences on the results. We further analyzed their expression in HPV-positive and -negative HNSCs, and no significant differences were observed. Moreover, we also tested the immune cell infiltration of DEGs in both HPV-positive and -negative HNSCs. Although studies have demonstrated that HPV influences the immune infiltration situation 26, it is not observed in ACTN2, MYH1, MYH2, MYH7, and NEB. As genes highly expressed in the tumor cells are expected to affect tumor purity positively, the association between the expression of prognosis DEGs and six immune infiltrates was evaluated. We found that all five DEGs have a significant positive correlation with CD4+T cells and macrophages. Therefore, we speculate that genetic mutations and differential expression of ACTN2, MYH1, MYH2, MYH7, and NEB genes in HNSC cells may be important drivers of CD4+T cell and macrophage infiltration.
section_name: Discussion

section_content: In our present study, we identified the biomarker of HNSC in terms of diagnosis and prognosis. TP53 was found as the most mutated gene in HNSC. ACTN2, MYH1, MYH2, MYH7, and NEB genes were significantly associated with poor prognosis. Moreover, MYH1, MYH2, and MYH7 were expected to be the specific diagnostic and prognostic biomarkers molecular biomarkers of HNSC.
section_name: Conclusions

section_name: Funding

section_content: This work was supported by the General Program of National Natural Science Foundation of China ( 81972555 to Fang-Zhou Liu ).
section_name: Funding

section_name: Data availability

section_content: The RNA-Seq and mutation data re-analyzed during the current study are available in the TCGA database with the links of https:// xenab rowser. net/ datap ages/? datas et= TCGA-HNSC. htseq_ fpkm. tsv& host= https% 3A% 2F% 2Fgdc. xenah ubs. net and https:// xenab rowser. net/ datap ages/? datas et= TCGA-HNSC. varsc an2_ snv. tsv& host= https% 3A% 2F% 2Fgdc. xenah ubs. net.
section_name: Data availability

section_name: Methods

section_content: Data collection. The clinical data and gene expression profiles of HNSC patients were downloaded from the UCSC Xena database (https:// xenab rowser. net/ datap ages/). As calculated, a total of 506 tumors and 44 adjacent HNSC samples were included in our present study. The gene expression levels of examples were provided in Supplementary File 2. The somatic mutation data was from the TCGA database by Varscan 27. 

Screening of high-frequency mutation genes. The maftools 28 was used to summarize the HNSC mutation information downloaded from the TCGA official website and screen the high-frequency mutation genes. The top 25 high-frequency genes with the mutation frequency were exhibited with a waterfall chart. 

Differentially expressed gene (DEG) identification. The DEGs in tumor tissues were identified according to the FPKM value by limma package 29 (Version 3. 10. 3), using P-value < 0. 05 and log 2 fold change (F. C. ) > 2 as thresholds. The Pheatmap R package was employed for drawing the heatmap and volcano. 

Then, the common DEGs with high-frequency mutation genes that exhibit a mutation frequency of more than 1% were retained for further study.
section_name: Methods

section_content: The PPI analysis was performed on all the common DEGs using STRING 30 (Version: 10. 0, http:// www. string-db. org/) database, with a 0. 4 (medium confidence) parameter. The network was constructed by Cytoscape (version: 3. 2. 0) 31. The most significant clustered modules in the PPI network were analyzed using the Cytoscape plug-in MCODE 32 (Version 1. 4. 2, http:// apps. cytos cape. org/ apps/ MCODE) method, with a threshold of score ≥ 10.
section_name: Protein-protein interaction (PPI) and module analyses.

section_content: The biological process (B. P. ) of Gene Ontology (G. O. ) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses on DEGs were performed by clusterProfiler 33, with adjusted P-value < 0. 05 and gene number count ≥ 2. The top ten G. O. terms and KEGG pathways were visualized.
section_name: Functional analysis.

section_content: The key DEGs in the module were subjected to expression verification and survival analysis by GEPIA 34. Survival analysis was grouped by median data and performed using OS survival information. The pan-cancer expression levels of key DEGs significantly associated with prognosis were assessed using UALCAN (http:// ualcan. path. uab. edu/ index. html). The expression levels of key DEGs in HNSC-positive and HNSC-negative samples were evaluated.
section_name: Key DEG expression validation and survival curve analysis.

section_content: As for the immune cell infiltration analysis, the HNSC samples were divided into HNSC-positive and HNSC-negative. Then, the infiltration of the key DEGs in six immune cells (B cells, CD 4+ T cells, CD 8+ T cells, neutrophils, macrophages, and dendritic cells) were identified through the TIMER database. 

Ethics approval and consent to participate. UCSC Xena and TCGA belongs to public databases. The patients involved in the database have obtained ethical approval. All methods were carried out in accordance with relevant guidelines and regulations.
section_name: Immune cell infiltration analysis.

section_content: G. J. all data analysis, and writing; Z. Y. data analysis, and writing; Y. Z. performed part of the analysis; X. Z. conceptualization, interpretation and discussion of results; F. L. conceptualization, advised all research, design of experiments, evaluation and discussion of results, revision of manuscript, and provided research funding. All authors have read and agreed to the published version of the manuscript.
section_name: Author contributions

section_content: The authors declare no competing interests.
section_name: Competing interests