The preferred role of nurses in family physician team

1.   Introduction

Head and neck squamous cell cancer (HNSCC) ranks the 6^th place among all cancers globally, which occurs in the oral cavity, nasal cavity, sinus, throat and pharynx. The pathological type of HNSCC is squamous cell carcinoma (except thyroid tumor). In America, the incidence of HNSCC was 7.863% in 2006 ^[1], in China, the incidence of HNSCC was up to 3.268% ^[2]. At present, the standard treatment of HNSCC is still surgery, intensity-modulated radiotherapy and platinum-based chemotherapy. However, there is no better treatment plan when first-line therapy is ineffective. Therefore, new treatment is urgently needed for HNSCC to prolong the survival time of patient ^[3]. In recent years, the emergence of immunotherapy has brought new hope to cancer patients. Nivolumab and Pembrolizumab are approved by the FDA for advanced HNSCC after first-line treatment ^[4,5]. However, the effectiveness of immunotherapy of HNSCC is still low, in this regard, searching for novel therapeutic targets to treat HNSCC is of great necessity.

With the development of computer science and bioinformatics, network-based methods have become an effective tool for the research of pathogenic mechanism ^[6]. Robust Rank Aggregation (RRA) method can maximally reduce errors or biases between multiple data sets ^[7]. Nowadays, RRA method has been widely used in cancer research ^[7,8,9]. This method is better than RemoveBatchEffect method, which was widely used on analyzing GEO datasets now.

The cibersort algorithm is a general calculation method to quantify cell components from a large number of tissues expression profiles (GEPs). Combined with support vector regression and prior knowledge from the expression profile of purified leukocyte subsets, cibersort algorithm can accurately estimate the immune components of tumor tissues.

The protein-protein interaction (PPI) network exerts a vital part in cancer biology, and it is an effective method for screening cancer related Hub genes. At present, some studies have pointed out that the method based on PPI network can successfully predict the Hub gene of breast cancer ^[10], liver cancer ^[11] and gastric cancer ^[12].

In this study, we downloaded the mircoarray datasets GSE686 ^[13], GSE2379 ^[14], GSE6631 ^[15], GSE13399 ^[16], GSE33205 ^[17,18], GSE33493, GSE39376, GSE107591 ^[19,20] were downloaded from the GEO database (www.ncbi.nlm.nih.gov/geo/). In addition, the RRA approach was used to discover the potent DEGs in normal compared with HNSCC samples. This study discovered altogether 240 potent DEGs, among which, 105 showed up-regulation and 135 showed down-regulation. Further, GO functional annotation and KEGG analyses were conducted to explore the functions of those identified DEGs. Cibersort algorithm was used to analyze the infiltration of immune cells. Then, a PPI network was established, meanwhile, the key modules were also built. At last, we screened 5 hub genes from the whole network using cytoHubba. Hub gene survival was analyzed by the R packages. To sum up, the potent DEGs along with hub genes were discovered in this study by the combined bioinformatic approach, which might be a new and potential prognostic biomarker.

2.   Materials and methods

2.1.   Data collection and processing

We downloaded 8 head and neck squamous cell carcinoma chips from GEO databases, which contained GSE686 ^[13], GSE2379 ^[14], GSE6631 ^[15], GSE13399 ^[16], GSE33205 ^[17,18], GSE33493, GSE39376, GSE107591 ^[19,20]. These datasets covered both human normal and tumor samples, and each dataset contained at least 20 samples. We downloaded the matrix files of all the microarray along with the platform annotation profiles, so that the names of microarray probes were easily converted into genetic symbols by the use of Perl. We identified the normal tissue from the HNSCC tissues from the 8 datasets using the R package limma function by the thresholds of P-value < 0.05 and log2-fold change (FC) > 0.5.

2.2.   Robust rank aggregation analysis

The RRA method is a standard approach that can minimizes bias and errors between multiple datasets. It can detect genes whose ranking is consistently better than expected under null hypothesis of uncorrelated inputs, and assign a significance score to each gene, the specific algorithmic computations are shown by Eqs (1) and (2). The underlying probability model makes the algorithm parameter independent and robust to outliers, noise and errors. Significance sores also provide a rigorous way to keep only the statistically relevant genes in the final list.

In this study, in order to integrate eight microarray datasets, we used RRA method to determine the robust DEGs. Before RRA analysis, we sequenced the up-regulated and down-regulated were sequenced from each dataset according to FC, then the sequencing results of 8 datasets were combined, and the R package RobustRankAggreg function was used to select robust DEGs according to the above thresholds.

2.3.   Functional enrichment analysis: GO and KEGG

For exploring functions of those selected DEGs, the R package "Clusterprofiler" was applied to obtain GO enrichment results, which included biological processes (BPs), cell components (CCs), together with molecular functions (MFs). and we also used the R package to analyzed the KEGG pathway of the robust DEGs. A difference of P < 0.05 indicated statistical significance.

2.4.   Immune infiltration calculated by cibersort analysis

The cibersort algorithm ^[21] has been developed as the machine learning approach on the basis of linear support vector regression (SVR), as shown by Eq (3), which shows high robustness against noise. This algorithm outperforms others in terms of noise, tightly associated cell types, along with unclear mixture content. This algorithm was used in the present work for predicting the infiltrating degrees of 22 immunocyte types within HNSCC tissues. LM22 is called immune cell gene expression tag matrix. It contained 547 genes and can distinguish 22 kinds of human hematopoietic cell phenotypes. The HNSCC data together with LM22 matrix were used to be the input of the above algorithm to obtain the proportions of those 22 immunocyte types within HNSCC. As a result, the cell composition related to HNSCC response was quantified. Using P < 0.05 as the standard to screen immune cell matrix, the relative expression of immune cells in normal compared with HNSCC samples was detected by R software package. Differences in normal compared with cancer tissues were determined by principal component analysis (PCA).

2.5.   Establishment of a PPI network and module analysis

We uploaded the differentially expressed genes to string online database, and selected the confidence level > 0.7 as the screening criteria, and removed the free nodes to get the PPI (protein-protein interaction network), and downloaded the gene interaction files. The PPI network was visualized using Cytoscape (Version 3.6.1). In addition, the MCODE plug-in of Cytoscape was used to screen those significant modules from the as-constructed PPI network ^[22].

2.6.   Identification of hub gene

The Cytoscape plug-in CytoHubba can be used to sort those network-derived nodes based on network features. Cytohubba offers 11 topological analysis methods, such as edge penetration component, degree, maximum cluster centrality, maximum neighborhood component density, maximum neighborhood component, and six centralities (bottleneck, eccentricity, compactness, radius, intermediate degree and stress) based on the shortest path. Of these 11 approaches, the modified Maximal Clique Centrality (MCC) approach is more effective in predicting the essential proteins in PPI network.

2.7.   Survival analysis

We obtained RNA-seq and clinical data of HNSCC cases from The Cancer Genome Atlas (TCGA) database. Thereafter, hub gene survival was analyzed by survminer and survival functions of R package. A difference of P < 0.05 indicated statistical significance.

3.   Results

3.1.   Discovery of DEGs from every dataset

The bioinformatic approaches were utilized in the present work for identifying DEGs, and analyze the biological characteristics of these DEGs. We conducted this study according to the workflow shown in Figure 1.

The microarray data of HNSCC GSE686, GSE2379, GSE6631, GSE13399, GSE33205, GSE33493, GSE39376, GSE107591 were downloaded and analyzed by R package limma software. The distribution of DEGs is shown in the volcanic map (Figure 2). A total of 336 samples including 108 normal samples and 258 tumor samples were analyzed in our study. On the basis of the exclusion standard of log2FC > 1 and P < 0.05, there were 4 downregulated genes in GSE686. In GSE2379, there were 443 DEGs selected, among which, 238 showed up-regulation while 205 showed down-regulation. In the GSE6631 dataset, we selected 142 GSEs, among which, 53 showed up-regulation while 89 showed down-regulation. Altogether 322 DEGs were selected from the GSE13399 dataset, among which, 149 showed up-regulation while 173 showed down-regulation. There were 419 DEGs selected from the GSE33205 dataset, of them, 178 showed up-regulation while 241 showed down-regulation. Altogether 3330 DEGs were found from GSE33493, of them, 1773 showed up-regulation while 1557 showed down-regulation. Altogether 940 DEGs were selected from GSE39376, of them, 414 showed up-regulation while 526 showed down-regulation. In addition, altogether 466 DEGs were selected from GSE107591, of them, 200 showed up-regulation while 266 showed down-regulation. The volcano map for each GSE is shown in Figure 2, where the green and red dots indicate genes with down-regulation and up-regulation, separately.

3.2.   Identification of robust DEGs by RRA method

We used the RRA method to integrate eight datasets. We selected altogether 240 DEGs, among which, 105 showed up-regulation while 135 showed down-regulation. According to the P-value threshold for selecting potent DEGs, the top 20 potent DEGs with up-regulation and down-regulation were distributed within the heat map.

3.3.   Functional enrichment analysis

Functions of those identified potent DEGs were analyzed by GO as well as KEGG analysis. Three function types were included in the GO analysis results: including BP, CC and MF. For BP, the robust DEGs showed enrichment in extracellular matrix (ECM) organization, skin development and extracellular organization. In the CC term, these genes were enriched in cornified envelope, collagen trimer complex as well as collagen-containing ECM. And in MF sections, the most significantly enriched terms were activity of serine-type endopeptidase, activity of serine-type peptidase, along with activity of serine hydrolase. The Figure 3 also showed the KEGG pathway enrichment analysis. Among them, IL-17 signaling pathway, cell cycle and ECM-receptor interaction are highly related to tumor growth and progress.

3.4.   Immune cell infiltration analysis

Using cibersort algorithm, the barplot of immune cells between normal and tumor samples in HNSCC tissues was shown in Figure 4(A). The heapmap of immune cells between normal and tumor samples in HNSCC tissues was shown in Figure 4(B). From the visualized violin plot (Figure 4C), we could find that B cells memory, B cells naïve, Macrophages M0, M1, T cells CD4 memory activated, Dendritic cells activated and Dendritic cells resting showed significant differences in normal compared with cancer tissues. PCA revealed no difference in normal compared with cancer tissues (Figure 4D).

3.5.   Establishment of a PPI network and module analysis

For better exploring the associations among the potent DEGs, the String database was used to establish the PPI network. The Cytoscape was used to establish the visual PPI network when there were hidden disconnected nodes with the confidence value of > 0.7 (Figure 5A). The final network contained 125 edges along with 115 nodes, and there were 85 upregulated genes and 67 downregulated genes. We selected three key networks from the whole network through MCODE plugin (Figure 5B–D).

3.6.   The identification of hub genes

CytoHubba was applied for predicting and exploring the key nodes from the as-constructed PPI network. It can score each node in PPI network with topological algorithms. The genes with high scores are identified to be the Hub genes. The present work adopted ten topological algorithms ((Density of Maximum Neighhorhood Component (DMNC), Maximal Clique Centrality (MCC), Degrees, Maximum Neighborhood Component (MNC), BottleNeck, Edge Percolated Component (EPC), Closeness, EcCentricity, Betweenness and Radiality) to determine genes in the entire network. After calculation, there were 5 hub genes selected by the 10 algorithms, there were BIRC5, AURKA, UBE2C, CDC20, COL4A1.

3.7.   Survival analysis

We analyzed the relationship between 5 hub genes and the overall survival rate of patients with HNSCC using "survminer" and "survival" packages of R software. Based on the best cutoff value calculated by the "surv_cutpoint" function for all Hub genes, we classified HNSCC samples as 2 groups (namely, high or low expression group), then acquired the respective Kaplan Meier (K-M) survival curves. As a result, BIRC5, AURKA (P = 0.017), UBE2C (P = 0.015) and (P = 0.045) expression was markedly related to the HNSCC survival.

4.   Discussion

As a result of the development of bioinformatics, there are more and more studies on HNSCC biomarkers in public databases such as GEO and TCGA database. For example, Ding Y et al. ^[23] and Chamorro Petronacci CM et al. ^[24] screened potential prognosis biomarkers of HNSCC. Nevertheless, due to the fact that the DEGs in the above research were selected using one dataset only with a small sample size, the results are unstable. In our study, compared with other HNSCC researches, we conformitied eight datasets by RRA methods.

In this study, we selected altogether 240 DEGs, including 105 with up-regulation while 135 with down-regulation. As suggested by GO as well as KEGG enrichment analysis, those potent DEGs were mainly enriched to extracellular matrix organization, skin development, extracellular organization, cornified envelope, complex of collagen trimers, serine-type activity and IL-17 signaling pathway, cell cycle and ECM-receptor interaction, which were related to tumor growth and progress. Then, we performed analysis of the immune cell infiltration between the normal and tumor HNSCC samples by cibersort algorithm. Module analysis was adopted to construct a PPI network based on the STRING database. Finally, 5 hub genes were screened from the whole network by cytoHubba including BIRC5, AURKA, UBE2C, CDC20, COL4A1. And the sample survival was analyzed using the R packages according to hub gene expression.

We also examined the GO terms together with KEGG pathways within the HNSCC samples in enrichment analysis. It has been extensively suggested that, the epithelial-mesenchymal transition (EMT) indicates a metastatic process. In this process, the epithelial cells acquire migratory and invasive mesenchymal phenotypes ^[25]. Extracellular matrix also exerts a vital part during cancer development, which can regulate cell growth, metabolism, migration, proliferation and differentiation through integrin or other cell surface receptors. Therefore, the degradation of ECM is indispensable for the invasion and metastasis of malignant tumors ^[26,27]. The high IL-17B level is verified to be tightly associated with dismal prognostic outcomes for cancer patients, such as breast cancer ^[28,29], gastric cancer ^[30], colon cancer ^[31], lung cancer ^[32] and so on. The IL-17 pathway may also be involved in the development of HNSCC. DNA damage may occur during normal cell division or in the presence of external stimuli. If the repair mechanism of DNA damage is defective, genomic instability can be caused. Cell cycle disorder is an important part of genomic instability, which may promote the further malignant transformation of unstable cells, thus accelerating the occurrence and development of tumors ^[33,34]. In line with the analysis above, the potent DEGs were tightly related to the HNSCC pathogenic mechanism as well as progression.

An increasing number of recent articles suggest that, tumor microenvironment (TME) exerts a vital part in tumor genesis and progression ^[35,36]. Typically, TME mainly includes endothelial cells, mesenchymal cells, immune cells, and ECM ^[37]. From previous studies, we found that the B7 protein family contains seven members: CD80, CD86, ICOS-L, PD-L1, PD-L2, B7-H3 and B7-H4 ^[38]. All ligands of the B7 family can be detected in dendritic cells, B cells, together with macrophages. To be specific, the B7 family mainly plays a role in regulating immune response. If B7 family gene is knocked out, the mice will suffer from immune deficiency and autoimmune disease ^[39]. PD-L1 can lead to tumor escape in the immune system by weakening the specific response of T cells to tumor cells ^[40]. Based on two clinical trials of PD-1 and PD-L1 in HNSCC (KEYNOTE-012 and CheckMate-141), Nivolumab and Pembrolizumab are approved by the FDA to treat advanced HNSCC. In our study, we found that B cells, dendritic cells and macrophages were the principal infiltrating immune cells in HNSCC tissues. The immunotherapy for HNSCC needs to be further explored.

In our study, we discovered 5 hub genes according to the as-constructed PPI network. Among them, 3 key genes were screened for further exploration. AURKA is one of the mitotic serine/threonine kinase family members, which exerts a vital part in a variety of biological events, such as centrosome separation and maturation, chromosome alignment, spindle assembly, as well as G2-to-M transition ^[41,42]. Besides, AURKA is previously suggested to show over-expression in diverse cancers ^{[43,44,45,46,47]}, such as neuroblastoma ^[48], gastric cancer ^[49,50,51] and so on. AURKA can promote tumor development by inhibiting tumor suppressor genes such as p73 ^[52,53] and p53 ^[54], activating β-catenin ^[55], NF-κB ^[56] and cap-dependent translation of oncogenes ^[57]. Nonetheless, it remains unknown about the function of AURKA within HNSCC. BIRC5, which is also referred to as survivin, belongs to the IAP family ^[58]. Its expression can be detected in different cancers, including breast cancer (BC) ^[59], colorectal cancer (CRC) ^[60], liver cancer ^[61] and so on. It inhibited the Caspase activity and apoptosis by inhibiting the binding of Caspase-3 and Caspase-7, thus leading to the survival of cancer cells in the process of tumorigenesis ^[62]. UBE2C belongs to the E2 ubiquitin-conjugating enzyme family ^[63], which exerts an important part in regulating the cell cycle, and this is achieved through the catalysis of polyubiquitination-induced APC/C substrate degradation ^[64]. Many recent studies report the abnormal expression of UBE2C in different human cancers, like hepatocellular carcinoma (HCC) ^[65], CRC ^[66,67], lung cancer ^[68], BC ^[69] and so on.

These abundant evidences showed that UBE2C played an important part in tumor genesis and development. According to our results, the high expression of AURKA, BIRC5 and UBE2C genes in HNSCC tumor showed poor survival of patients in TCGA database, suggesting that these genes are related to the prognosis of HNSCC, and they may act as therapeutic targets of HNSCC.

5.   Conclusions

In conclusion, through Robust Rank Aggregation method and cibersort algorithm method, a series of robust DEFs and gene modules were identified in HNSCC. The identified genes were subjected to functional analyses, which revealed their close relationship with the occurrence and development of HNSCC. We not only screened five hub genes, but also analyzed the immune cell infiltration in HNSCC. From above discussion, we found that AURKA, BIRC5 and UBE2C may be considered as new biomarker and therapeutic targets of HNSCC. The Robust Rank Aggregation method and cibersort algorithm method can accurately predict the potential prognostic biomarker and therapeutic targets. Further investigation is warranted to explore their roles in HNSCC in further research.

Acknowledgement

This work was supported by grants from the National Natural Science Foundation of China (82073344, 81874217, 81703027, 81703028, 81672983).

Conflict of interests

The authors declare that there is no conflict of interests.